From 41935728e1ab0ad5ec371cc1f623f65affb8a214 Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Wed, 15 Apr 2020 23:58:39 +0200 Subject: [PATCH] change json format --- multigrid/quad.json | 62 +++-- multigrid/sedecim.json | 212 ++++++++------- multigrid/src/main.rs | 68 +---- multigrid/src/parsing.rs | 544 +++++++++++++++++++++------------------ sbp/src/euler.rs | 90 +++---- sbp/src/operators.rs | 2 +- sbp/src/utils.rs | 26 ++ 7 files changed, 506 insertions(+), 498 deletions(-) diff --git a/multigrid/quad.json b/multigrid/quad.json index da8c17c..5f75041 100644 --- a/multigrid/quad.json +++ b/multigrid/quad.json @@ -1,42 +1,52 @@ { - "grids": [ - { - "name": "grid0", + "grids": { + "default": { + "operators": { + "xi": "upwind9", + "eta": "upwind9" + } + }, + "grid0": { "x": "linspace:-5:0:50", "y": "linspace:0:5:50", - "dirS": "grid1", - "dirN": "grid1", - "dirE": "grid3", - "dirW": "grid3" + "boundary_conditions": { + "south": "grid1", + "north": "grid1", + "east": "grid3", + "west": "grid3" + } }, - { - "name": "grid1", + "grid1": { "x": "linspace:-5:0:50", "y": "linspace:-5:0:50", - "dirS": "grid0", - "dirN": "grid0", - "dirE": "grid2", - "dirW": "grid2" + "boundary_conditions": { + "south": "grid0", + "north": "grid0", + "east": "grid2", + "west": "grid2" + } }, - { - "name": "grid2", + "grid2": { "x": "linspace:0:5:50", "y": "linspace:-5:0:50", - "dirS": "grid3", - "dirN": "grid3", - "dirE": "grid1", - "dirW": "grid1" + "boundary_conditions": { + "south": "grid3", + "north": "grid3", + "east": "grid1", + "west": "grid1" + } }, - { - "name": "grid3", + "grid3": { "x": "linspace:0:5:50", "y": "linspace:0:5:50", - "dirS": "grid2", - "dirN": "grid2", - "dirE": "grid0", - "dirW": "grid0" + "boundary_conditions": { + "south": "grid2", + "north": "grid2", + "east": "grid0", + "west": "grid0" + } } - ], + }, "integration_time": 2.0, "vortex": { "x0": -1.0, diff --git a/multigrid/sedecim.json b/multigrid/sedecim.json index 311ea8f..7dc0142 100644 --- a/multigrid/sedecim.json +++ b/multigrid/sedecim.json @@ -1,150 +1,166 @@ { - "grids": [ - { - "name": "grid00", + "grids": { + "grid00": { "x": "linspace:-5:-2.5:128", "y": "linspace:2.5:5:128", - "dirS": "grid01", - "dirN": "grid11", - "dirE": "grid03", - "dirW": "grid33" + "boundary_conditions": { + "south": "grid01", + "north": "grid11", + "east": "grid03", + "west": "grid33" + } }, - { - "name": "grid01", + "grid01": { "x": "linspace:-5:-2.5:128", "y": "linspace:0:2.5:128", - "dirS": "grid10", - "dirN": "grid00", - "dirE": "grid02", - "dirW": "grid32" + "boundary_conditions": { + "south": "grid10", + "north": "grid00", + "east": "grid02", + "west": "grid32" + } }, - { - "name": "grid02", + "grid02": { "x": "linspace:-2.5:0:128", "y": "linspace:0:2.5:128", - "dirS": "grid13", - "dirN": "grid03", - "dirE": "grid31", - "dirW": "grid01" + "boundary_conditions": { + "south": "grid13", + "north": "grid03", + "east": "grid31", + "west": "grid01" + } }, - { - "name": "grid03", + "grid03": { "x": "linspace:-2.5:0:128", "y": "linspace:2.5:5:128", - "dirS": "grid02", - "dirN": "grid12", - "dirE": "grid30", - "dirW": "grid00" + "boundary_conditions": { + "south": "grid02", + "north": "grid12", + "east": "grid30", + "west": "grid00" + } }, - { - "name": "grid10", + "grid10": { "x": "linspace:-5:-2.5:128", "y": "linspace:-2.5:0:128", - "dirS": "grid11", - "dirN": "grid01", - "dirE": "grid13", - "dirW": "grid23" + "boundary_conditions": { + "south": "grid11", + "north": "grid01", + "east": "grid13", + "west": "grid23" + } }, - { - "name": "grid11", + "grid11": { "x": "linspace:-5:-2.5:128", "y": "linspace:-5:-2.5:128", - "dirS": "grid00", - "dirN": "grid10", - "dirE": "grid12", - "dirW": "grid22" + "boundary_conditions": { + "south": "grid00", + "north": "grid10", + "east": "grid12", + "west": "grid22" + } }, - { - "name": "grid12", + "grid12": { "x": "linspace:-2.5:0:128", "y": "linspace:-5:-2.5:128", - "dirS": "grid03", - "dirN": "grid13", - "dirE": "grid21", - "dirW": "grid11" + "boundary_conditions": { + "south": "grid03", + "north": "grid13", + "east": "grid21", + "west": "grid11" + } }, - { - "name": "grid13", + "grid13": { "x": "linspace:-2.5:0:128", "y": "linspace:-2.5:0:128", - "dirS": "grid12", - "dirN": "grid02", - "dirE": "grid20", - "dirW": "grid10" + "boundary_conditions": { + "south": "grid12", + "north": "grid02", + "east": "grid20", + "west": "grid10" + } }, - { - "name": "grid20", + "grid20": { "x": "linspace:0:2.5:128", "y": "linspace:-2.5:0:128", - "dirS": "grid21", - "dirN": "grid31", - "dirE": "grid23", - "dirW": "grid13" + "boundary_conditions": { + "south": "grid21", + "north": "grid31", + "east": "grid23", + "west": "grid13" + } }, - { - "name": "grid21", + "grid21": { "x": "linspace:0:2.5:128", "y": "linspace:-5:-2.5:128", - "dirS": "grid30", - "dirN": "grid20", - "dirE": "grid22", - "dirW": "grid12" + "boundary_conditions": { + "south": "grid30", + "north": "grid20", + "east": "grid22", + "west": "grid12" + } }, - { - "name": "grid22", + "grid22": { "x": "linspace:2.5:5:128", "y": "linspace:-5:-2.5:128", - "dirS": "grid33", - "dirN": "grid23", - "dirE": "grid11", - "dirW": "grid21" + "boundary_conditions": { + "south": "grid33", + "north": "grid23", + "east": "grid11", + "west": "grid21" + } }, - { - "name": "grid23", + "grid23": { "x": "linspace:2.5:5:128", "y": "linspace:-2.5:0:128", - "dirS": "grid22", - "dirN": "grid32", - "dirE": "grid10", - "dirW": "grid20" + "boundary_conditions": { + "south": "grid22", + "north": "grid32", + "east": "grid10", + "west": "grid20" + } }, - { - "name": "grid30", + "grid30": { "x": "linspace:0:2.5:128", "y": "linspace:2.5:5:128", - "dirS": "grid31", - "dirN": "grid21", - "dirE": "grid33", - "dirW": "grid03" + "boundary_conditions": { + "south": "grid31", + "north": "grid21", + "east": "grid33", + "west": "grid03" + } }, - { - "name": "grid31", + "grid31": { "x": "linspace:0:2.5:128", "y": "linspace:0:2.5:128", - "dirS": "grid20", - "dirN": "grid30", - "dirE": "grid32", - "dirW": "grid02" + "boundary_conditions": { + "south": "grid20", + "north": "grid30", + "east": "grid32", + "west": "grid02" + } }, - { - "name": "grid32", + "grid32": { "x": "linspace:2.5:5:128", "y": "linspace:0:2.5:128", - "dirS": "grid23", - "dirN": "grid33", - "dirE": "grid01", - "dirW": "grid31" + "boundary_conditions": { + "south": "grid23", + "north": "grid33", + "east": "grid01", + "west": "grid31" + } }, - { - "name": "grid33", + "grid33": { "x": "linspace:2.5:5:128", "y": "linspace:2.5:5:128", - "dirS": "grid32", - "dirN": "grid22", - "dirE": "grid00", - "dirW": "grid30" + "boundary_conditions": { + "south": "grid32", + "north": "grid22", + "east": "grid00", + "west": "grid30" + } } - ], + }, "integration_time": 2.0, "vortex": { "x0": -1.0, diff --git a/multigrid/src/main.rs b/multigrid/src/main.rs index a924aed..65ed1cc 100644 --- a/multigrid/src/main.rs +++ b/multigrid/src/main.rs @@ -10,6 +10,8 @@ use parsing::{json_to_grids, json_to_vortex}; mod file; use file::*; +pub(crate) type DiffOp = Either, Box>; + struct System { fnow: Vec, fnext: Vec, @@ -20,16 +22,14 @@ struct System { bt: Vec, eb: Vec, time: Float, - operators: Vec, Box>>, - interpolation_operators: Vec, + operators: Vec, } impl System { fn new( grids: Vec, bt: Vec, - interpolation_operators: Vec, - operators: Vec, Box>>, + operators: Vec, ) -> Self { let fnow = grids .iter() @@ -66,7 +66,6 @@ impl System { bt, eb, time: 0.0, - interpolation_operators, operators, } } @@ -83,7 +82,6 @@ impl System { let bt = &self.bt; let wb = &mut self.wb; let mut eb = &mut self.eb; - let intops = &self.interpolation_operators; let operators = &self.operators; let rhs = move |fut: &mut [euler::Field], @@ -91,7 +89,7 @@ impl System { time: Float, _c: (), _mt: &mut ()| { - let bc = euler::extract_boundaries(prev, &bt, &mut eb, &grids, time, Some(intops)); + let bc = euler::extract_boundaries(prev, &bt, &mut eb, &grids, time); pool.scope(|s| { for (((((fut, prev), bc), wb), metrics), op) in fut .iter_mut() @@ -162,63 +160,13 @@ fn main() { let filecontents = std::fs::read_to_string(&opt.json).unwrap(); let json = json::parse(&filecontents).unwrap(); - let jgrids = json_to_grids(json["grids"].clone()).unwrap(); + let vortexparams = json_to_vortex(json["vortex"].clone()); - - let mut bt = Vec::with_capacity(jgrids.len()); - let determine_bc = |dir: Option<&String>| match dir { - Some(dir) => { - if dir == "vortex" { - euler::BoundaryCharacteristic::Vortex(vortexparams) - } else if let Some(grid) = dir.strip_prefix("interpolate:") { - euler::BoundaryCharacteristic::Interpolate( - jgrids - .iter() - .position(|other| other.name.as_ref().map_or(false, |name| name == grid)) - .unwrap(), - ) - } else { - euler::BoundaryCharacteristic::Grid( - jgrids - .iter() - .position(|other| other.name.as_ref().map_or(false, |name| name == dir)) - .unwrap(), - ) - } - } - None => euler::BoundaryCharacteristic::This, - }; - for grid in &jgrids { - bt.push(euler::BoundaryCharacteristics { - north: determine_bc(grid.dirn.as_ref()), - south: determine_bc(grid.dirs.as_ref()), - east: determine_bc(grid.dire.as_ref()), - west: determine_bc(grid.dirw.as_ref()), - }); - } - let interpolation_operators = jgrids - .iter() - .map(|_g| euler::InterpolationOperators { - north: Some(Box::new(operators::Interpolation4)), - south: Some(Box::new(operators::Interpolation4)), - east: Some(Box::new(operators::Interpolation4)), - west: Some(Box::new(operators::Interpolation4)), - }) - .collect::>(); - - let grids = jgrids - .into_iter() - .map(|egrid| egrid.grid) - .collect::>(); + let (_names, grids, bt, operators) = json_to_grids(json["grids"].clone(), vortexparams); let integration_time: Float = json["integration_time"].as_number().unwrap().into(); - let operators = grids - .iter() - .map(|_| Right(Box::new(operators::Upwind4) as Box)) - .collect::>(); - - let mut sys = System::new(grids, bt, interpolation_operators, operators); + let mut sys = System::new(grids, bt, operators); sys.vortex(0.0, vortexparams); let max_n = { diff --git a/multigrid/src/parsing.rs b/multigrid/src/parsing.rs index 3452e1d..d0170f6 100644 --- a/multigrid/src/parsing.rs +++ b/multigrid/src/parsing.rs @@ -1,17 +1,180 @@ +use super::DiffOp; use crate::grid::Grid; use crate::Float; +use either::*; +use json::JsonValue; +use sbp::utils::h2linspace; -#[derive(Debug, Clone)] -pub struct ExtendedGrid { - pub grid: Grid, - pub name: Option, - pub dire: Option, - pub dirw: Option, - pub dirn: Option, - pub dirs: Option, +pub fn json_to_grids( + mut jsongrids: JsonValue, + vortexparams: sbp::euler::VortexParameters, +) -> ( + Vec, + Vec, + Vec, + Vec, +) { + let default = jsongrids.remove("default"); + let default_operator = { + let operators = &default["operators"]; + let defaultxi = operators["xi"].as_str().unwrap_or("upwind4"); + let defaulteta = operators["eta"].as_str().unwrap_or("upwind4"); + + (defaulteta.to_string(), defaultxi.to_string()) + }; + /* + let default_bc: sbp::utils::Direction> = { + let bc = &default["boundary_conditions"]; + Direction { + south: bc["south"].as_str().map(|x| x.to_string()), + north: bc["north"].as_str().map(|x| x.to_string()), + west: bc["west"].as_str().map(|x| x.to_string()), + east: bc["east"].as_str().map(|x| x.to_string()), + } + }; + */ + let mut names = Vec::new(); + let mut grids = Vec::new(); + + let mut operators: Vec = Vec::new(); + for (name, grid) in jsongrids.entries() { + names.push(name.to_string()); + grids.push(json2grid(grid["x"].clone(), grid["y"].clone()).unwrap()); + + operators.push({ + use sbp::operators::*; + let opxi = grid["operators"]["xi"] + .as_str() + .unwrap_or(&default_operator.1); + + let opeta = grid["operators"]["eta"] + .as_str() + .unwrap_or(&default_operator.1); + + match (opeta, opxi) { + ("upwind4", "upwind4") => Left(Box::new(Upwind4) as Box), + ("upwind9", "upwind9") => Left(Box::new(Upwind9) as Box), + ("upwind4h2", "upwind4h2") => Left(Box::new(Upwind4h2) as Box), + ("upwind9h2", "upwind9h2") => Left(Box::new(Upwind9h2) as Box), + + ("upwind4", "upwind9") => { + Left(Box::new((&Upwind4, &Upwind9)) as Box) + } + ("upwind4", "upwind4h2") => { + Left(Box::new((&Upwind4, &Upwind4h2)) as Box) + } + ("upwind4", "upwind9h2") => { + Left(Box::new((&Upwind4, &Upwind9h2)) as Box) + } + + ("upwind9", "upwind4") => { + Left(Box::new((&Upwind9, &Upwind4)) as Box) + } + ("upwind9", "upwind4h2") => { + Left(Box::new((&Upwind9, &Upwind4h2)) as Box) + } + ("upwind9", "upwind9h2") => { + Left(Box::new((&Upwind9, &Upwind9h2)) as Box) + } + + ("upwind4h2", "upwind4") => { + Left(Box::new((&Upwind4h2, &Upwind4)) as Box) + } + ("upwind4h2", "upwind9") => { + Left(Box::new((&Upwind4h2, &Upwind9)) as Box) + } + ("upwind4h2", "upwind9h2") => { + Left(Box::new((&Upwind4h2, &Upwind9h2)) as Box) + } + + ("upwind9h2", "upwind4") => { + Left(Box::new((&Upwind9h2, &Upwind4)) as Box) + } + ("upwind9h2", "upwind9") => { + Left(Box::new((&Upwind9h2, &Upwind9)) as Box) + } + ("upwind9h2", "upwind4h2") => { + Left(Box::new((&Upwind9h2, &Upwind4h2)) as Box) + } + + (opeta, opxi) => panic!("combination {} {} not yet implemented", opeta, opxi), + } + }); + } + + let mut bcs = Vec::new(); + let determine_bc = |dir: Option<&str>| match dir { + Some(dir) => { + if dir == "vortex" { + sbp::euler::BoundaryCharacteristic::Vortex(vortexparams) + } else if let Some(grid) = dir.strip_prefix("interpolate:") { + use sbp::operators::*; + let (grid, int_op) = if let Some(rest) = grid.strip_prefix("4:") { + ( + rest, + Box::new(Interpolation4) as Box, + ) + } else if let Some(rest) = grid.strip_prefix("9:") { + ( + rest, + Box::new(Interpolation9) as Box, + ) + } else if let Some(rest) = grid.strip_prefix("8:") { + ( + rest, + Box::new(Interpolation8) as Box, + ) + } else if let Some(rest) = grid.strip_prefix("9h2:") { + ( + rest, + Box::new(Interpolation9h2) as Box, + ) + } else { + ( + grid, + Box::new(Interpolation4) as Box, + ) + }; + sbp::euler::BoundaryCharacteristic::Interpolate( + names.iter().position(|other| other == grid).unwrap(), + int_op, + ) + } else { + sbp::euler::BoundaryCharacteristic::Grid( + names.iter().position(|other| other == dir).unwrap(), + ) + } + } + None => sbp::euler::BoundaryCharacteristic::This, + }; + for name in &names { + let bc = &jsongrids[name]["boundary_conditions"]; + let bc_n = determine_bc(bc["north"].as_str()); + let bc_s = determine_bc(bc["south"].as_str()); + let bc_e = determine_bc(bc["east"].as_str()); + let bc_w = determine_bc(bc["west"].as_str()); + + let bc = sbp::euler::BoundaryCharacteristics { + north: bc_n, + south: bc_s, + east: bc_e, + west: bc_w, + }; + + bcs.push(bc); + } + + (names, grids, bcs, operators) +} +#[derive(Debug)] +enum ArrayForm { + /// Only know the one dimension, will broadcast to + /// two dimensions once we know about both dims + Array1(ndarray::Array1), + /// The usize is the inner dimension (nx) + Array2(ndarray::Array2), } -use json::JsonValue; /// Parsing json strings to some gridlike form /// /// Each grid should be an object with the descriptors on the form @@ -29,245 +192,136 @@ use json::JsonValue; /// Optional parameters: /// * name (for relating boundaries) /// * dir{e,w,n,s} (for boundary terms) -pub fn json_to_grids(json: JsonValue) -> Result, String> { - fn json_to_grid(mut grid: JsonValue) -> Result { - #[derive(Debug)] - enum ArrayForm { - /// Only know the one dimension, will broadcast to - /// two dimensions once we know about both dims - Array1(ndarray::Array1), - /// The usize is the inner dimension (nx) - Array2(ndarray::Array2), - } - if grid.is_empty() { - return Err("empty object".to_string()); - } - let name = grid.remove("name").take_string(); - let dire = grid.remove("dirE").take_string(); - let dirw = grid.remove("dirW").take_string(); - let dirn = grid.remove("dirN").take_string(); - let dirs = grid.remove("dirS").take_string(); - - let to_array_form = |mut x: JsonValue| { - if let Some(s) = x.take_string() { - if let Some(s) = s.strip_prefix("linspace:") { - let (s, h2) = if let Some(s) = s.strip_prefix("h2:") { - (s, true) - } else { - (s, false) - }; - - // linspace:start:stop:steps - let mut iter = s.split(':'); - - let start = iter.next(); - let start: Float = match start { - Some(x) => x.parse().map_err(|e| format!("linspace: {}", e))?, - None => return Err(format!("")), - }; - let end = iter.next(); - let end: Float = match end { - Some(x) => x.parse().map_err(|e| format!("linspace: {}", e))?, - None => return Err(format!("")), - }; - let steps = iter.next(); - let steps: usize = match steps { - Some(x) => x.parse().map_err(|e| format!("linspace: {}", e))?, - None => return Err(format!("")), - }; - if iter.next().is_some() { - return Err("linspace: contained more than expected".to_string()); - } - Ok(ArrayForm::Array1(if h2 { - h2linspace(start, end, steps) - } else { - ndarray::Array::linspace(start, end, steps) - })) +fn json2grid(x: JsonValue, y: JsonValue) -> Result { + let to_array_form = |mut x: JsonValue| { + if let Some(s) = x.take_string() { + if let Some(s) = s.strip_prefix("linspace:") { + let (s, h2) = if let Some(s) = s.strip_prefix("h2:") { + (s, true) } else { - Err("Could not parse gridline".to_string()) + (s, false) + }; + + // linspace:start:stop:steps + let mut iter = s.split(':'); + + let start = iter.next(); + let start: Float = match start { + Some(x) => x.parse().map_err(|e| format!("linspace: {}", e))?, + None => return Err(format!("")), + }; + let end = iter.next(); + let end: Float = match end { + Some(x) => x.parse().map_err(|e| format!("linspace: {}", e))?, + None => return Err(format!("")), + }; + let steps = iter.next(); + let steps: usize = match steps { + Some(x) => x.parse().map_err(|e| format!("linspace: {}", e))?, + None => return Err(format!("")), + }; + if iter.next().is_some() { + return Err("linspace: contained more than expected".to_string()); } - } else if x.is_array() { - let arrlen = x.len(); - if arrlen == 0 { + Ok(ArrayForm::Array1(if h2 { + h2linspace(start, end, steps) + } else { + ndarray::Array::linspace(start, end, steps) + })) + } else { + Err("Could not parse gridline".to_string()) + } + } else if x.is_array() { + let arrlen = x.len(); + if arrlen == 0 { + return Err("gridline does not have any members".to_string()); + } + if !x[0].is_array() { + let v = x + .members() + .map(|x: &JsonValue| -> Result { + Ok(x.as_number() + .ok_or_else(|| { + "Array contained something that could not be converted to an array" + .to_string() + })? + .into()) + }) + .collect::, _>>()?; + Ok(ArrayForm::Array1(ndarray::Array::from(v))) + } else { + let arrlen2 = x[0].len(); + if arrlen2 == 0 { return Err("gridline does not have any members".to_string()); } - if !x[0].is_array() { - let v = x - .members() - .map(|x: &JsonValue| -> Result { - Ok(x.as_number().ok_or_else(|| "Array contained something that could not be converted to an array".to_string())?.into()) - }) - .collect::, _>>()?; - Ok(ArrayForm::Array1(ndarray::Array::from(v))) - } else { - let arrlen2 = x[0].len(); - if arrlen2 == 0 { - return Err("gridline does not have any members".to_string()); + for member in x.members() { + if arrlen2 != member.len() { + return Err("some arrays seems to have differing lengths".to_string()); } - for member in x.members() { - if arrlen2 != member.len() { - return Err("some arrays seems to have differing lengths".to_string()); - } - } - let mut arr = ndarray::Array::zeros((arrlen, arrlen2)); - for (mut arr, member) in arr.outer_iter_mut().zip(x.members()) { - for (a, m) in arr.iter_mut().zip(member.members()) { - *a = m - .as_number() - .ok_or_else(|| { - "array contained something which was not a number".to_string() - })? - .into() - } - } - Ok(ArrayForm::Array2(arr)) } - } else { - Err("Inner object was not a string value, or an array".to_string()) + let mut arr = ndarray::Array::zeros((arrlen, arrlen2)); + for (mut arr, member) in arr.outer_iter_mut().zip(x.members()) { + for (a, m) in arr.iter_mut().zip(member.members()) { + *a = m + .as_number() + .ok_or_else(|| { + "array contained something which was not a number".to_string() + })? + .into() + } + } + Ok(ArrayForm::Array2(arr)) } - }; - - let x = grid.remove("x"); - if x.is_empty() { - return Err("x was empty".to_string()); + } else { + Err("Inner object was not a string value, or an array".to_string()) } - let x = to_array_form(x)?; + }; - let y = grid.remove("y"); - if y.is_empty() { - return Err("y was empty".to_string()); - } - let y = to_array_form(y)?; - - let (x, y) = match (x, y) { - (ArrayForm::Array1(x), ArrayForm::Array1(y)) => { - let xlen = x.len(); - let ylen = y.len(); - let x = x.broadcast((ylen, xlen)).unwrap().to_owned(); - let y = y - .broadcast((xlen, ylen)) - .unwrap() - .reversed_axes() - .to_owned(); - - (x, y) - } - (ArrayForm::Array2(x), ArrayForm::Array2(y)) => { - assert_eq!(x.shape(), y.shape()); - (x, y) - } - (ArrayForm::Array1(x), ArrayForm::Array2(y)) => { - assert_eq!(x.len(), y.shape()[1]); - let x = x.broadcast((y.shape()[1], x.len())).unwrap().to_owned(); - (x, y) - } - (ArrayForm::Array2(x), ArrayForm::Array1(y)) => { - assert_eq!(x.shape()[0], y.len()); - let y = y - .broadcast((x.shape()[1], y.len())) - .unwrap() - .reversed_axes() - .to_owned(); - (x, y) - } - }; - assert_eq!(x.shape(), y.shape()); - - if !grid.is_empty() { - eprintln!("Grid contains some unused entries"); - for i in grid.entries() { - eprintln!("{:#?}", i); - } - } - - Ok(ExtendedGrid { - grid: Grid::new(x, y).unwrap(), - name, - dire, - dirw, - dirn, - dirs, - }) + if x.is_empty() { + return Err("x was empty".to_string()); } + let x = to_array_form(x)?; - match json { - JsonValue::Array(a) => a - .into_iter() - .map(json_to_grid) - .collect::, _>>(), - grid => Ok(vec![json_to_grid(grid)?]), + if y.is_empty() { + return Err("y was empty".to_string()); } -} + let y = to_array_form(y)?; -#[test] -fn parse_linspace() { - let grids = json_to_grids( - json::parse(r#"[{"name": "main", "x": "linspace:0:10:20", "y": "linspace:0:10:21"}]"#) - .unwrap(), - ) - .unwrap(); - assert_eq!(grids.len(), 1); - assert_eq!(grids[0].grid.x.shape(), [21, 20]); - assert_eq!(grids[0].grid.y.shape(), [21, 20]); - assert_eq!(grids[0].name.as_ref().unwrap(), "main"); - let grids = json_to_grids( - json::parse(r#"{"name": "main", "x": "linspace:0:10:20", "y": "linspace:0:10:21"}"#) - .unwrap(), - ) - .unwrap(); - assert_eq!(grids.len(), 1); - assert_eq!(grids[0].grid.x.shape(), [21, 20]); - assert_eq!(grids[0].grid.y.shape(), [21, 20]); - assert_eq!(grids[0].name.as_ref().unwrap(), "main"); -} + let (x, y) = match (x, y) { + (ArrayForm::Array1(x), ArrayForm::Array1(y)) => { + let xlen = x.len(); + let ylen = y.len(); + let x = x.broadcast((ylen, xlen)).unwrap().to_owned(); + let y = y + .broadcast((xlen, ylen)) + .unwrap() + .reversed_axes() + .to_owned(); -#[test] -fn parse_1d() { - let grids = - json_to_grids(json::parse(r#"{"x": [1, 2, 3, 4, 5.1, 3], "y": [1, 2]}"#).unwrap()).unwrap(); - assert_eq!(grids.len(), 1); - let grid = &grids[0]; - assert_eq!(grid.grid.x.shape(), &[2, 6]); - assert_eq!(grid.grid.x.shape(), grid.grid.y.shape()); -} + (x, y) + } + (ArrayForm::Array2(x), ArrayForm::Array2(y)) => { + assert_eq!(x.shape(), y.shape()); + (x, y) + } + (ArrayForm::Array1(x), ArrayForm::Array2(y)) => { + assert_eq!(x.len(), y.shape()[1]); + let x = x.broadcast((y.shape()[1], x.len())).unwrap().to_owned(); + (x, y) + } + (ArrayForm::Array2(x), ArrayForm::Array1(y)) => { + assert_eq!(x.shape()[0], y.len()); + let y = y + .broadcast((x.shape()[1], y.len())) + .unwrap() + .reversed_axes() + .to_owned(); + (x, y) + } + }; + assert_eq!(x.shape(), y.shape()); -#[test] -fn parse_2d() { - let grids = - json_to_grids(json::parse(r#"{"x": [[1, 2], [3, 4], [5.1, 3]], "y": [1, 2, 3]}"#).unwrap()) - .unwrap(); - assert_eq!(grids.len(), 1); - let grid = &grids[0]; - assert_eq!(grid.grid.x.shape(), &[3, 2]); - assert_eq!(grid.grid.x.shape(), grid.grid.y.shape()); - json_to_grids( - json::parse(r#"{"x": [[1, 2], [3, 4], [5.1, 3], [1]], "y": [1, 2, 3]}"#).unwrap(), - ) - .unwrap_err(); - json_to_grids( - json::parse(r#"{"y": [[1, 2], [3, 4], [5.1, 3], [1]], "x": [1, 2, 3]}"#).unwrap(), - ) - .unwrap_err(); - let grids = json_to_grids( - json::parse(r#"{"x": [[1, 2], [3, 4], [5.1, 3]], "y": [[1, 2], [3, 4], [5, 6]]}"#).unwrap(), - ) - .unwrap(); - assert_eq!(grids.len(), 1); - json_to_grids( - json::parse(r#"{"x": [[1, 2], [3, 4], [5.1, 3]], "y": [[1, 2], [3, 4], [5]]}"#).unwrap(), - ) - .unwrap_err(); -} - -#[test] -fn parse_err() { - json_to_grids(json::parse(r#"{}"#).unwrap()).unwrap_err(); - json_to_grids(json::parse(r#"0.45"#).unwrap()).unwrap_err(); - json_to_grids(json::parse(r#"{"x": "linspace", "y": [0.1, 0.2]}"#).unwrap()).unwrap_err(); - json_to_grids(json::parse(r#"{"x": "linspace:::", "y": [0.1, 0.2]}"#).unwrap()).unwrap_err(); - json_to_grids(json::parse(r#"{"x": "linspace:1.2:3.1:412.2", "y": [0.1, 0.2]}"#).unwrap()) - .unwrap_err(); - json_to_grids(json::parse(r#"{"x": [-2, -3, "dfd"], "y": [0.1, 0.2]}"#).unwrap()).unwrap_err(); + Ok(Grid::new(x, y).unwrap()) } pub fn json_to_vortex(mut json: JsonValue) -> super::euler::VortexParameters { @@ -292,27 +346,3 @@ pub fn json_to_vortex(mut json: JsonValue) -> super::euler::VortexParameters { eps, } } - -pub fn h2linspace(start: Float, end: Float, n: usize) -> ndarray::Array1 { - let h = (end - start) / (n - 2) as Float; - ndarray::Array1::from_shape_fn(n, |i| match i { - 0 => start, - i if i == n - 1 => end, - i => start + h * (i as Float - 0.5), - }) -} - -#[test] -fn test_h2linspace() { - let x = h2linspace(-1.0, 1.0, 50); - println!("{}", x); - approx::assert_abs_diff_eq!(x[0], -1.0, epsilon = 1e-6); - approx::assert_abs_diff_eq!(x[49], 1.0, epsilon = 1e-6); - let hend = x[1] - x[0]; - let h = x[2] - x[1]; - approx::assert_abs_diff_eq!(x[49] - x[48], hend, epsilon = 1e-6); - approx::assert_abs_diff_eq!(2.0 * hend, h, epsilon = 1e-6); - for i in 1..48 { - approx::assert_abs_diff_eq!(x[i + 1] - x[i], h, epsilon = 1e-6); - } -} diff --git a/sbp/src/euler.rs b/sbp/src/euler.rs index 6ed6dce..bab75d9 100644 --- a/sbp/src/euler.rs +++ b/sbp/src/euler.rs @@ -597,13 +597,12 @@ fn fluxes(k: (&mut Field, &mut Field), y: &Field, metrics: &Metrics) { } } -#[derive(Clone, Debug)] pub enum BoundaryCharacteristic { This, Grid(usize), Vortex(VortexParameters), // Vortices(Vec), - Interpolate(usize), + Interpolate(usize, Box), } pub type BoundaryTerms<'a> = Direction>; @@ -618,68 +617,59 @@ fn boundary_extractor<'a>( north: match bc.north { BoundaryCharacteristic::This => field.south(), BoundaryCharacteristic::Vortex(_params) => todo!(), - BoundaryCharacteristic::Grid(_) | BoundaryCharacteristic::Interpolate(_) => { + BoundaryCharacteristic::Grid(_) | BoundaryCharacteristic::Interpolate(_, _) => { panic!("Only working on self grid") } }, south: match bc.south { BoundaryCharacteristic::This => field.north(), BoundaryCharacteristic::Vortex(_params) => todo!(), - BoundaryCharacteristic::Grid(_) | BoundaryCharacteristic::Interpolate(_) => { + BoundaryCharacteristic::Grid(_) | BoundaryCharacteristic::Interpolate(_, _) => { panic!("Only working on self grid") } }, west: match bc.west { BoundaryCharacteristic::This => field.east(), BoundaryCharacteristic::Vortex(_params) => todo!(), - BoundaryCharacteristic::Grid(_) | BoundaryCharacteristic::Interpolate(_) => { + BoundaryCharacteristic::Grid(_) | BoundaryCharacteristic::Interpolate(_, _) => { panic!("Only working on self grid") } }, east: match bc.east { BoundaryCharacteristic::This => field.west(), BoundaryCharacteristic::Vortex(_params) => todo!(), - BoundaryCharacteristic::Grid(_) | BoundaryCharacteristic::Interpolate(_) => { + BoundaryCharacteristic::Grid(_) | BoundaryCharacteristic::Interpolate(_, _) => { panic!("Only working on self grid") } }, } } -pub type InterpolationOperators = Direction>>; - pub fn extract_boundaries<'a>( fields: &'a [Field], bt: &[BoundaryCharacteristics], eb: &'a mut [BoundaryStorage], grids: &[Grid], time: Float, - interpolation_operators: Option<&[InterpolationOperators]>, ) -> Vec> { bt.iter() .zip(eb) .zip(grids) .zip(fields) - .enumerate() - .map(|(ig, (((bt, eb), grid), field))| BoundaryTerms { - north: match bt.north { + .map(|(((bt, eb), grid), field)| BoundaryTerms { + north: match &bt.north { BoundaryCharacteristic::This => field.south(), - BoundaryCharacteristic::Grid(g) => fields[g].south(), + BoundaryCharacteristic::Grid(g) => fields[*g].south(), BoundaryCharacteristic::Vortex(v) => { let field = eb.north.as_mut().unwrap(); - vortexify(field.view_mut(), grid.north(), v, time); + vortexify(field.view_mut(), grid.north(), *v, time); field.view() } - BoundaryCharacteristic::Interpolate(g) => { + BoundaryCharacteristic::Interpolate(g, operator) => { let to = eb.north.as_mut().unwrap(); - let fine2coarse = field.nx() < fields[g].nx(); + let fine2coarse = field.nx() < fields[*g].nx(); - let operator = interpolation_operators.as_ref().unwrap()[ig] - .north - .as_ref() - .unwrap(); - - for (mut to, from) in to.outer_iter_mut().zip(fields[g].south().outer_iter()) { + for (mut to, from) in to.outer_iter_mut().zip(fields[*g].south().outer_iter()) { if fine2coarse { operator.fine2coarse(from.view(), to.view_mut()); } else { @@ -689,23 +679,19 @@ pub fn extract_boundaries<'a>( to.view() } }, - south: match bt.south { + south: match &bt.south { BoundaryCharacteristic::This => field.north(), - BoundaryCharacteristic::Grid(g) => fields[g].north(), + BoundaryCharacteristic::Grid(g) => fields[*g].north(), BoundaryCharacteristic::Vortex(v) => { let field = eb.south.as_mut().unwrap(); - vortexify(field.view_mut(), grid.south(), v, time); + vortexify(field.view_mut(), grid.south(), *v, time); field.view() } - BoundaryCharacteristic::Interpolate(g) => { + BoundaryCharacteristic::Interpolate(g, operator) => { let to = eb.south.as_mut().unwrap(); - let fine2coarse = field.nx() < fields[g].nx(); - let operator = interpolation_operators.as_ref().unwrap()[ig] - .south - .as_ref() - .unwrap(); + let fine2coarse = field.nx() < fields[*g].nx(); - for (mut to, from) in to.outer_iter_mut().zip(fields[g].north().outer_iter()) { + for (mut to, from) in to.outer_iter_mut().zip(fields[*g].north().outer_iter()) { if fine2coarse { operator.fine2coarse(from.view(), to.view_mut()); } else { @@ -715,23 +701,19 @@ pub fn extract_boundaries<'a>( to.view() } }, - west: match bt.west { + west: match &bt.west { BoundaryCharacteristic::This => field.east(), - BoundaryCharacteristic::Grid(g) => fields[g].east(), + BoundaryCharacteristic::Grid(g) => fields[*g].east(), BoundaryCharacteristic::Vortex(v) => { let field = eb.west.as_mut().unwrap(); - vortexify(field.view_mut(), grid.west(), v, time); + vortexify(field.view_mut(), grid.west(), *v, time); field.view() } - BoundaryCharacteristic::Interpolate(g) => { + BoundaryCharacteristic::Interpolate(g, operator) => { let to = eb.west.as_mut().unwrap(); - let fine2coarse = field.ny() < fields[g].ny(); - let operator = interpolation_operators.as_ref().unwrap()[ig] - .west - .as_ref() - .unwrap(); + let fine2coarse = field.ny() < fields[*g].ny(); - for (mut to, from) in to.outer_iter_mut().zip(fields[g].east().outer_iter()) { + for (mut to, from) in to.outer_iter_mut().zip(fields[*g].east().outer_iter()) { if fine2coarse { operator.fine2coarse(from.view(), to.view_mut()); } else { @@ -741,23 +723,19 @@ pub fn extract_boundaries<'a>( to.view() } }, - east: match bt.east { + east: match &bt.east { BoundaryCharacteristic::This => field.west(), - BoundaryCharacteristic::Grid(g) => fields[g].west(), + BoundaryCharacteristic::Grid(g) => fields[*g].west(), BoundaryCharacteristic::Vortex(v) => { let field = eb.east.as_mut().unwrap(); - vortexify(field.view_mut(), grid.east(), v, time); + vortexify(field.view_mut(), grid.east(), *v, time); field.view() } - BoundaryCharacteristic::Interpolate(g) => { + BoundaryCharacteristic::Interpolate(g, operator) => { let to = eb.east.as_mut().unwrap(); - let fine2coarse = field.ny() < fields[g].ny(); - let operator = interpolation_operators.as_ref().unwrap()[ig] - .east - .as_ref() - .unwrap(); + let fine2coarse = field.ny() < fields[*g].ny(); - for (mut to, from) in to.outer_iter_mut().zip(fields[g].west().outer_iter()) { + for (mut to, from) in to.outer_iter_mut().zip(fields[*g].west().outer_iter()) { if fine2coarse { operator.fine2coarse(from.view(), to.view_mut()); } else { @@ -778,25 +756,25 @@ impl BoundaryStorage { pub fn new(bt: &BoundaryCharacteristics, grid: &Grid) -> Self { Self { north: match bt.north { - BoundaryCharacteristic::Vortex(_) | BoundaryCharacteristic::Interpolate(_) => { + BoundaryCharacteristic::Vortex(_) | BoundaryCharacteristic::Interpolate(_, _) => { Some(ndarray::Array2::zeros((4, grid.nx()))) } _ => None, }, south: match bt.south { - BoundaryCharacteristic::Vortex(_) | BoundaryCharacteristic::Interpolate(_) => { + BoundaryCharacteristic::Vortex(_) | BoundaryCharacteristic::Interpolate(_, _) => { Some(ndarray::Array2::zeros((4, grid.nx()))) } _ => None, }, east: match bt.east { - BoundaryCharacteristic::Vortex(_) | BoundaryCharacteristic::Interpolate(_) => { + BoundaryCharacteristic::Vortex(_) | BoundaryCharacteristic::Interpolate(_, _) => { Some(ndarray::Array2::zeros((4, grid.ny()))) } _ => None, }, west: match bt.west { - BoundaryCharacteristic::Vortex(_) | BoundaryCharacteristic::Interpolate(_) => { + BoundaryCharacteristic::Vortex(_) | BoundaryCharacteristic::Interpolate(_, _) => { Some(ndarray::Array2::zeros((4, grid.ny()))) } _ => None, diff --git a/sbp/src/operators.rs b/sbp/src/operators.rs index 502357d..5e07370 100644 --- a/sbp/src/operators.rs +++ b/sbp/src/operators.rs @@ -185,7 +185,7 @@ mod traditional8; pub use traditional8::SBP8; mod interpolation; -pub use interpolation::Interpolation4; +pub use interpolation::{Interpolation4, Interpolation8, Interpolation9, Interpolation9h2}; #[cfg(test)] pub(crate) mod testing { diff --git a/sbp/src/utils.rs b/sbp/src/utils.rs index d96ccdd..ab9c227 100644 --- a/sbp/src/utils.rs +++ b/sbp/src/utils.rs @@ -1,6 +1,32 @@ +use crate::Float; + pub struct Direction { pub north: T, pub south: T, pub west: T, pub east: T, } + +pub fn h2linspace(start: Float, end: Float, n: usize) -> ndarray::Array1 { + let h = (end - start) / (n - 2) as Float; + ndarray::Array1::from_shape_fn(n, |i| match i { + 0 => start, + i if i == n - 1 => end, + i => start + h * (i as Float - 0.5), + }) +} + +#[test] +fn test_h2linspace() { + let x = h2linspace(-1.0, 1.0, 50); + println!("{}", x); + approx::assert_abs_diff_eq!(x[0], -1.0, epsilon = 1e-6); + approx::assert_abs_diff_eq!(x[49], 1.0, epsilon = 1e-6); + let hend = x[1] - x[0]; + let h = x[2] - x[1]; + approx::assert_abs_diff_eq!(x[49] - x[48], hend, epsilon = 1e-6); + approx::assert_abs_diff_eq!(2.0 * hend, h, epsilon = 1e-6); + for i in 1..48 { + approx::assert_abs_diff_eq!(x[i + 1] - x[i], h, epsilon = 1e-6); + } +}