From 358f83151343895e5e5f70b23ba868da8f80d0c5 Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Sat, 4 Apr 2020 00:29:02 +0200 Subject: [PATCH] separate metrics and grid --- sbp/examples/multigrid.rs | 27 +++++----- sbp/src/euler.rs | 111 +++++++++++++++++++------------------- sbp/src/grid.rs | 60 +++++++++++++-------- sbp/src/integrate.rs | 9 ++-- sbp/src/maxwell.rs | 91 ++++++++++++++++--------------- sbp/src/utils.rs | 90 ++++++++++++++++++------------- 6 files changed, 215 insertions(+), 173 deletions(-) diff --git a/sbp/examples/multigrid.rs b/sbp/examples/multigrid.rs index 42f508d..71c54ce 100644 --- a/sbp/examples/multigrid.rs +++ b/sbp/examples/multigrid.rs @@ -14,12 +14,13 @@ struct System { euler::Field, )>, k: [Vec; 4], - grids: Vec>, + grids: Vec, + metrics: Vec>, bt: Vec, } impl System { - fn new(grids: Vec>, bt: Vec) -> Self { + fn new(grids: Vec, bt: Vec) -> Self { let fnow = grids .iter() .map(|g| euler::Field::new(g.ny(), g.nx())) @@ -33,6 +34,7 @@ impl System { }) .collect(); let k = [fnow.clone(), fnow.clone(), fnow.clone(), fnow.clone()]; + let metrics = grids.iter().map(|g| g.metrics().unwrap()).collect(); Self { fnow, @@ -40,6 +42,7 @@ impl System { k, wb, grids, + metrics, bt, } } @@ -142,14 +145,14 @@ impl System { }) .collect::>(); - for ((((prev, fut), grid), wb), bt) in fields + for ((((prev, fut), metrics), wb), bt) in fields .iter() .zip(fnext) - .zip(&self.grids) + .zip(&self.metrics) .zip(&mut self.wb) .zip(bt) { - euler::RHS_upwind(fut, prev, grid, &bt, wb) + euler::RHS_upwind(fut, prev, metrics, &bt, wb) } } } @@ -257,14 +260,14 @@ impl System { }, }) .collect::>(); - for ((((prev, fut), grid), wb), bt) in fields + for ((((prev, fut), metrics), wb), bt) in fields .iter() .zip(&mut self.k[i]) - .zip(&self.grids) + .zip(&self.metrics) .zip(&mut self.wb) .zip(bt) { - s.spawn(move |_| euler::RHS_upwind(fut, prev, grid, &bt, wb)); + s.spawn(move |_| euler::RHS_upwind(fut, prev, metrics, &bt, wb)); } }); } @@ -307,15 +310,13 @@ fn main() { west: determine_bc(grid.dirw.as_ref()), }); } - let mut grids: Vec> = Vec::with_capacity(jgrids.len()); - for grid in jgrids { - grids.push(grid::Grid::new(grid.x, grid.y).unwrap()); - } + let grids = jgrids.into_iter().map(|egrid| egrid.grid).collect(); + let integration_time: f64 = json["integration_time"].as_number().unwrap().into(); let vortexparams = utils::json_to_vortex(json["vortex"].clone()); - let mut sys = System::new(grids, bt); + let mut sys = System::::new(grids, bt); sys.vortex(0.0, vortexparams); let max_n = { diff --git a/sbp/src/euler.rs b/sbp/src/euler.rs index 6c03857..3132f71 100644 --- a/sbp/src/euler.rs +++ b/sbp/src/euler.rs @@ -1,4 +1,4 @@ -use super::grid::Grid; +use super::grid::{Grid, Metrics}; use super::integrate; use super::operators::{SbpOperator, UpwindOperator}; use super::Float; @@ -13,7 +13,7 @@ pub const GAMMA: Float = 1.4; pub struct System { sys: (Field, Field), wb: WorkBuffers, - grid: Grid, + grid: (Grid, Metrics), } impl System { @@ -21,11 +21,12 @@ impl System { let grid = Grid::new(x, y).expect( "Could not create grid. Different number of elements compared to width*height?", ); + let metrics = grid.metrics().unwrap(); let nx = grid.nx(); let ny = grid.ny(); Self { sys: (Field::new(ny, nx), Field::new(ny, nx)), - grid, + grid: (grid, metrics), wb: WorkBuffers::new(ny, nx), } } @@ -37,16 +38,17 @@ impl System { east: BoundaryCharacteristic::This, west: BoundaryCharacteristic::This, }; - let rhs_trad = |k: &mut Field, y: &Field, grid: &Grid<_>, wb: &mut _| { + let rhs_trad = |k: &mut Field, y: &Field, grid: &Grid, metrics: &Metrics<_>, wb: &mut _| { let boundaries = boundary_extractor(y, grid, &bc); - RHS_trad(k, y, grid, &boundaries, wb) + RHS_trad(k, y, metrics, &boundaries, wb) }; integrate::rk4( rhs_trad, &self.sys.0, &mut self.sys.1, dt, - &self.grid, + &self.grid.0, + &self.grid.1, &mut self.wb.k, &mut self.wb.tmp, ); @@ -56,7 +58,7 @@ impl System { pub fn vortex(&mut self, t: Float, vortex_parameters: VortexParameters) { self.sys .0 - .vortex(self.grid.x.view(), self.grid.y.view(), t, vortex_parameters); + .vortex(self.grid.0.x(), self.grid.0.y(), t, vortex_parameters); } #[allow(clippy::many_single_char_names)] @@ -70,12 +72,9 @@ impl System { mach: 0.5, }; - self.sys.0.vortex( - self.grid.x.view(), - self.grid.y.view(), - 0.0, - vortex_parameters, - ) + self.sys + .0 + .vortex(self.grid.0.x(), self.grid.0.y(), 0.0, vortex_parameters) } pub fn field(&self) -> &Field { @@ -83,10 +82,10 @@ impl System { } pub fn x(&self) -> ArrayView2 { - self.grid.x.view() + self.grid.0.x.view() } pub fn y(&self) -> ArrayView2 { - self.grid.y.view() + self.grid.0.y.view() } } @@ -98,16 +97,18 @@ impl System { east: BoundaryCharacteristic::This, west: BoundaryCharacteristic::This, }; - let rhs_upwind = |k: &mut Field, y: &Field, grid: &Grid<_>, wb: &mut _| { - let boundaries = boundary_extractor(y, grid, &bc); - RHS_upwind(k, y, grid, &boundaries, wb) - }; + let rhs_upwind = + |k: &mut Field, y: &Field, grid: &Grid, metrics: &Metrics<_>, wb: &mut _| { + let boundaries = boundary_extractor(y, grid, &bc); + RHS_upwind(k, y, metrics, &boundaries, wb) + }; integrate::rk4( rhs_upwind, &self.sys.0, &mut self.sys.1, dt, - &self.grid, + &self.grid.0, + &self.grid.1, &mut self.wb.k, &mut self.wb.tmp, ); @@ -382,13 +383,13 @@ fn pressure(gamma: Float, rho: Float, rhou: Float, rhov: Float, e: Float) -> Flo pub(crate) fn RHS_trad( k: &mut Field, y: &Field, - grid: &Grid, + metrics: &Metrics, boundaries: &BoundaryTerms, tmp: &mut (Field, Field, Field, Field, Field, Field), ) { let ehat = &mut tmp.0; let fhat = &mut tmp.1; - fluxes((ehat, fhat), y, grid); + fluxes((ehat, fhat), y, metrics); let dE = &mut tmp.2; let dF = &mut tmp.3; @@ -405,24 +406,24 @@ pub(crate) fn RHS_trad( azip!((out in &mut k.0, eflux in &dE.0, fflux in &dF.0, - detj in &grid.detj.broadcast((4, y.ny(), y.nx())).unwrap()) { + detj in &metrics.detj.broadcast((4, y.ny(), y.nx())).unwrap()) { *out = (-eflux - fflux)/detj }); - SAT_characteristics(k, y, grid, boundaries); + SAT_characteristics(k, y, metrics, boundaries); } #[allow(non_snake_case)] pub fn RHS_upwind( k: &mut Field, y: &Field, - grid: &Grid, + metrics: &Metrics, boundaries: &BoundaryTerms, tmp: &mut (Field, Field, Field, Field, Field, Field), ) { let ehat = &mut tmp.0; let fhat = &mut tmp.1; - fluxes((ehat, fhat), y, grid); + fluxes((ehat, fhat), y, metrics); let dE = &mut tmp.2; let dF = &mut tmp.3; @@ -438,25 +439,25 @@ pub fn RHS_upwind( let ad_xi = &mut tmp.4; let ad_eta = &mut tmp.5; - upwind_dissipation((ad_xi, ad_eta), y, grid, (&mut tmp.0, &mut tmp.1)); + upwind_dissipation((ad_xi, ad_eta), y, metrics, (&mut tmp.0, &mut tmp.1)); azip!((out in &mut k.0, eflux in &dE.0, fflux in &dF.0, ad_xi in &ad_xi.0, ad_eta in &ad_eta.0, - detj in &grid.detj.broadcast((4, y.ny(), y.nx())).unwrap()) { + detj in &metrics.detj.broadcast((4, y.ny(), y.nx())).unwrap()) { *out = (-eflux - fflux + ad_xi + ad_eta)/detj }); - SAT_characteristics(k, y, grid, boundaries); + SAT_characteristics(k, y, metrics, boundaries); } #[allow(clippy::many_single_char_names)] fn upwind_dissipation( k: (&mut Field, &mut Field), y: &Field, - grid: &Grid, + metrics: &Metrics, tmp: (&mut Field, &mut Field), ) { let n = y.nx() * y.ny(); @@ -471,11 +472,11 @@ fn upwind_dissipation( .axis_iter(ndarray::Axis(1)) .zip(tmp0.axis_iter_mut(ndarray::Axis(1))) .zip(tmp1.axis_iter_mut(ndarray::Axis(1))) - .zip(grid.detj.iter()) - .zip(grid.detj_dxi_dx.iter()) - .zip(grid.detj_dxi_dy.iter()) - .zip(grid.detj_deta_dx.iter()) - .zip(grid.detj_deta_dy.iter()) + .zip(metrics.detj.iter()) + .zip(metrics.detj_dxi_dx.iter()) + .zip(metrics.detj_dxi_dy.iter()) + .zip(metrics.detj_deta_dx.iter()) + .zip(metrics.detj_deta_dy.iter()) { let rho = y[0]; assert!(rho > 0.0); @@ -520,11 +521,11 @@ fn upwind_dissipation( UO::disseta(tmp.1.e(), k.1.e_mut()); } -fn fluxes(k: (&mut Field, &mut Field), y: &Field, grid: &Grid) { - let j_dxi_dx = grid.detj_dxi_dx.view(); - let j_dxi_dy = grid.detj_dxi_dy.view(); - let j_deta_dx = grid.detj_deta_dx.view(); - let j_deta_dy = grid.detj_deta_dy.view(); +fn fluxes(k: (&mut Field, &mut Field), y: &Field, metrics: &Metrics) { + let j_dxi_dx = metrics.detj_dxi_dx.view(); + let j_dxi_dy = metrics.detj_dxi_dy.view(); + let j_deta_dx = metrics.detj_deta_dx.view(); + let j_deta_dy = metrics.detj_deta_dy.view(); let rho = y.rho(); let rhou = y.rhou(); @@ -590,9 +591,9 @@ pub struct BoundaryCharacteristics { pub west: BoundaryCharacteristic, } -fn boundary_extractor<'a, SBP: SbpOperator>( +fn boundary_extractor<'a>( field: &'a Field, - _grid: &Grid, + _grid: &Grid, bc: &BoundaryCharacteristics, ) -> BoundaryTerms<'a> { BoundaryTerms { @@ -624,7 +625,7 @@ fn boundary_extractor<'a, SBP: SbpOperator>( fn SAT_characteristics( k: &mut Field, y: &Field, - grid: &Grid, + metrics: &Metrics, boundaries: &BoundaryTerms, ) { // North boundary @@ -640,9 +641,9 @@ fn SAT_characteristics( hi, sign, tau, - grid.detj.slice(slice), - grid.detj_deta_dx.slice(slice), - grid.detj_deta_dy.slice(slice), + metrics.detj.slice(slice), + metrics.detj_deta_dx.slice(slice), + metrics.detj_deta_dy.slice(slice), ); } // South boundary @@ -658,9 +659,9 @@ fn SAT_characteristics( hi, sign, tau, - grid.detj.slice(slice), - grid.detj_deta_dx.slice(slice), - grid.detj_deta_dy.slice(slice), + metrics.detj.slice(slice), + metrics.detj_deta_dx.slice(slice), + metrics.detj_deta_dy.slice(slice), ); } // West Boundary @@ -676,9 +677,9 @@ fn SAT_characteristics( hi, sign, tau, - grid.detj.slice(slice), - grid.detj_dxi_dx.slice(slice), - grid.detj_dxi_dy.slice(slice), + metrics.detj.slice(slice), + metrics.detj_dxi_dx.slice(slice), + metrics.detj_dxi_dy.slice(slice), ); } // East Boundary @@ -694,9 +695,9 @@ fn SAT_characteristics( hi, sign, tau, - grid.detj.slice(slice), - grid.detj_dxi_dx.slice(slice), - grid.detj_dxi_dy.slice(slice), + metrics.detj.slice(slice), + metrics.detj_dxi_dx.slice(slice), + metrics.detj_dxi_dy.slice(slice), ); } } diff --git a/sbp/src/grid.rs b/sbp/src/grid.rs index ed50d43..0cbce71 100644 --- a/sbp/src/grid.rs +++ b/sbp/src/grid.rs @@ -2,13 +2,16 @@ use crate::Float; use ndarray::Array2; #[derive(Debug, Clone)] -pub struct Grid +pub struct Grid { + pub(crate) x: Array2, + pub(crate) y: Array2, +} + +#[derive(Debug, Clone)] +pub struct Metrics where SBP: super::operators::SbpOperator, { - pub(crate) x: Array2, - pub(crate) y: Array2, - pub(crate) detj: Array2, pub(crate) detj_dxi_dx: Array2, pub(crate) detj_dxi_dy: Array2, @@ -18,11 +21,39 @@ where operator: std::marker::PhantomData, } -impl Grid { +impl Grid { pub fn new(x: Array2, y: Array2) -> Result { assert_eq!(x.shape(), y.shape()); - let ny = x.shape()[0]; - let nx = x.shape()[1]; + + Ok(Self { x, y }) + } + pub fn nx(&self) -> usize { + self.x.shape()[1] + } + pub fn ny(&self) -> usize { + self.x.shape()[0] + } + + pub fn x(&self) -> ndarray::ArrayView2 { + self.x.view() + } + pub fn y(&self) -> ndarray::ArrayView2 { + self.y.view() + } + + pub fn metrics( + &self, + ) -> Result, ndarray::ShapeError> { + Metrics::new(self) + } +} + +impl Metrics { + fn new(grid: &Grid) -> Result { + let ny = grid.ny(); + let nx = grid.nx(); + let x = &grid.x; + let y = &grid.y; let mut dx_dxi = Array2::zeros((ny, nx)); SBP::diffxi(x.view(), dx_dxi.view_mut()); @@ -57,8 +88,6 @@ impl Grid { let detj_deta_dy = dx_dxi; Ok(Self { - x, - y, detj, detj_dxi_dx, detj_dxi_dy, @@ -67,17 +96,4 @@ impl Grid { operator: std::marker::PhantomData, }) } - pub fn nx(&self) -> usize { - self.x.shape()[1] - } - pub fn ny(&self) -> usize { - self.x.shape()[0] - } - - pub fn x(&self) -> ndarray::ArrayView2 { - self.x.view() - } - pub fn y(&self) -> ndarray::ArrayView2 { - self.y.view() - } } diff --git a/sbp/src/integrate.rs b/sbp/src/integrate.rs index 7969561..7cce642 100644 --- a/sbp/src/integrate.rs +++ b/sbp/src/integrate.rs @@ -1,4 +1,4 @@ -use super::grid::Grid; +use super::grid::{Grid, Metrics}; use super::operators::SbpOperator; use super::Float; use ndarray::{Array3, Zip}; @@ -8,13 +8,14 @@ pub(crate) fn rk4<'a, F: 'a, SBP, RHS, WB>( prev: &F, fut: &mut F, dt: Float, - grid: &Grid, + grid: &Grid, + metrics: &Metrics, k: &mut [F; 4], mut wb: &mut WB, ) where F: std::ops::Deref> + std::ops::DerefMut>, SBP: SbpOperator, - RHS: Fn(&mut F, &F, &Grid, &mut WB), + RHS: Fn(&mut F, &F, &Grid, &Metrics, &mut WB), { assert_eq!(prev.shape(), fut.shape()); @@ -48,6 +49,6 @@ pub(crate) fn rk4<'a, F: 'a, SBP, RHS, WB>( } }; - rhs(&mut k[i], &fut, grid, &mut wb); + rhs(&mut k[i], &fut, grid, metrics, &mut wb); } } diff --git a/sbp/src/maxwell.rs b/sbp/src/maxwell.rs index 21e3a06..1e7c287 100644 --- a/sbp/src/maxwell.rs +++ b/sbp/src/maxwell.rs @@ -1,4 +1,4 @@ -use super::grid::Grid; +use super::grid::{Grid, Metrics}; use super::integrate; use super::operators::{SbpOperator, UpwindOperator}; use crate::Float; @@ -77,7 +77,8 @@ impl Field { pub struct System { sys: (Field, Field), wb: WorkBuffers, - grid: Grid, + grid: Grid, + metrics: Metrics, } impl System { @@ -87,10 +88,12 @@ impl System { let nx = x.shape()[1]; let grid = Grid::new(x, y).unwrap(); + let metrics = grid.metrics().unwrap(); Self { sys: (Field::new(ny, nx), Field::new(ny, nx)), grid, + metrics, wb: WorkBuffers::new(ny, nx), } } @@ -118,6 +121,7 @@ impl System { &mut self.sys.1, dt, &self.grid, + &self.metrics, &mut self.wb.k, &mut self.wb.tmp, ); @@ -134,6 +138,7 @@ impl System { &mut self.sys.1, dt, &self.grid, + &self.metrics, &mut self.wb.k, &mut self.wb.tmp, ); @@ -171,10 +176,11 @@ fn gaussian(x: Float, x0: Float, y: Float, y0: Float) -> Float { fn RHS( k: &mut Field, y: &Field, - grid: &Grid, + _grid: &Grid, + metrics: &Metrics, tmp: &mut (Array2, Array2, Array2, Array2), ) { - fluxes(k, y, grid, tmp); + fluxes(k, y, metrics, tmp); let boundaries = BoundaryTerms { north: Boundary::This, @@ -182,10 +188,10 @@ fn RHS( west: Boundary::This, east: Boundary::This, }; - SAT_characteristics(k, y, grid, &boundaries); + SAT_characteristics(k, y, metrics, &boundaries); azip!((k in &mut k.0, - &detj in &grid.detj.broadcast((3, y.ny(), y.nx())).unwrap()) { + &detj in &metrics.detj.broadcast((3, y.ny(), y.nx())).unwrap()) { *k /= detj; }); } @@ -194,11 +200,12 @@ fn RHS( fn RHS_upwind( k: &mut Field, y: &Field, - grid: &Grid, + _grid: &Grid, + metrics: &Metrics, tmp: &mut (Array2, Array2, Array2, Array2), ) { - fluxes(k, y, grid, tmp); - dissipation(k, y, grid, tmp); + fluxes(k, y, metrics, tmp); + dissipation(k, y, metrics, tmp); let boundaries = BoundaryTerms { north: Boundary::This, @@ -206,10 +213,10 @@ fn RHS_upwind( west: Boundary::This, east: Boundary::This, }; - SAT_characteristics(k, y, grid, &boundaries); + SAT_characteristics(k, y, metrics, &boundaries); azip!((k in &mut k.0, - &detj in &grid.detj.broadcast((3, y.ny(), y.nx())).unwrap()) { + &detj in &metrics.detj.broadcast((3, y.ny(), y.nx())).unwrap()) { *k /= detj; }); } @@ -217,20 +224,20 @@ fn RHS_upwind( fn fluxes( k: &mut Field, y: &Field, - grid: &Grid, + metrics: &Metrics, tmp: &mut (Array2, Array2, Array2, Array2), ) { // ex = hz_y { ndarray::azip!((a in &mut tmp.0, - &dxi_dy in &grid.detj_dxi_dy, + &dxi_dy in &metrics.detj_dxi_dy, &hz in &y.hz()) *a = dxi_dy * hz ); SBP::diffxi(tmp.0.view(), tmp.1.view_mut()); ndarray::azip!((b in &mut tmp.2, - &deta_dy in &grid.detj_deta_dy, + &deta_dy in &metrics.detj_deta_dy, &hz in &y.hz()) *b = deta_dy * hz ); @@ -244,8 +251,8 @@ fn fluxes( { // hz = -ey_x + ex_y ndarray::azip!((a in &mut tmp.0, - &dxi_dx in &grid.detj_dxi_dx, - &dxi_dy in &grid.detj_dxi_dy, + &dxi_dx in &metrics.detj_dxi_dx, + &dxi_dy in &metrics.detj_dxi_dy, &ex in &y.ex(), &ey in &y.ey()) *a = dxi_dx * -ey + dxi_dy * ex @@ -253,8 +260,8 @@ fn fluxes( SBP::diffxi(tmp.0.view(), tmp.1.view_mut()); ndarray::azip!((b in &mut tmp.2, - &deta_dx in &grid.detj_deta_dx, - &deta_dy in &grid.detj_deta_dy, + &deta_dx in &metrics.detj_deta_dx, + &deta_dy in &metrics.detj_deta_dy, &ex in &y.ex(), &ey in &y.ey()) *b = deta_dx * -ey + deta_dy * ex @@ -269,14 +276,14 @@ fn fluxes( // ey = -hz_x { ndarray::azip!((a in &mut tmp.0, - &dxi_dx in &grid.detj_dxi_dx, + &dxi_dx in &metrics.detj_dxi_dx, &hz in &y.hz()) *a = dxi_dx * -hz ); SBP::diffxi(tmp.0.view(), tmp.1.view_mut()); azip!((b in &mut tmp.2, - &deta_dx in &grid.detj_deta_dx, + &deta_dx in &metrics.detj_deta_dx, &hz in &y.hz()) *b = deta_dx * -hz ); @@ -291,14 +298,14 @@ fn fluxes( fn dissipation( k: &mut Field, y: &Field, - grid: &Grid, + metrics: &Metrics, tmp: &mut (Array2, Array2, Array2, Array2), ) { // ex component { ndarray::azip!((a in &mut tmp.0, - &kx in &grid.detj_dxi_dx, - &ky in &grid.detj_dxi_dy, + &kx in &metrics.detj_dxi_dx, + &ky in &metrics.detj_dxi_dy, &ex in &y.ex(), &ey in &y.ey()) { let r = Float::hypot(kx, ky); @@ -307,8 +314,8 @@ fn dissipation( UO::dissxi(tmp.0.view(), tmp.1.view_mut()); ndarray::azip!((b in &mut tmp.2, - &kx in &grid.detj_deta_dx, - &ky in &grid.detj_deta_dy, + &kx in &metrics.detj_deta_dx, + &ky in &metrics.detj_deta_dy, &ex in &y.ex(), &ey in &y.ey()) { let r = Float::hypot(kx, ky); @@ -324,8 +331,8 @@ fn dissipation( // hz component { ndarray::azip!((a in &mut tmp.0, - &kx in &grid.detj_dxi_dx, - &ky in &grid.detj_dxi_dy, + &kx in &metrics.detj_dxi_dx, + &ky in &metrics.detj_dxi_dy, &hz in &y.hz()) { let r = Float::hypot(kx, ky); *a = r * hz; @@ -333,8 +340,8 @@ fn dissipation( UO::dissxi(tmp.0.view(), tmp.1.view_mut()); ndarray::azip!((b in &mut tmp.2, - &kx in &grid.detj_deta_dx, - &ky in &grid.detj_deta_dy, + &kx in &metrics.detj_deta_dx, + &ky in &metrics.detj_deta_dy, &hz in &y.hz()) { let r = Float::hypot(kx, ky); *b = r * hz; @@ -349,8 +356,8 @@ fn dissipation( // ey { ndarray::azip!((a in &mut tmp.0, - &kx in &grid.detj_dxi_dx, - &ky in &grid.detj_dxi_dy, + &kx in &metrics.detj_dxi_dx, + &ky in &metrics.detj_dxi_dy, &ex in &y.ex(), &ey in &y.ey()) { let r = Float::hypot(kx, ky); @@ -359,8 +366,8 @@ fn dissipation( UO::dissxi(tmp.0.view(), tmp.1.view_mut()); ndarray::azip!((b in &mut tmp.2, - &kx in &grid.detj_deta_dx, - &ky in &grid.detj_deta_dy, + &kx in &metrics.detj_deta_dx, + &ky in &metrics.detj_deta_dy, &ex in &y.ex(), &ey in &y.ey()) { let r = Float::hypot(kx, ky); @@ -392,7 +399,7 @@ pub struct BoundaryTerms { fn SAT_characteristics( k: &mut Field, y: &Field, - grid: &Grid, + metrics: &Metrics, boundaries: &BoundaryTerms, ) { let ny = y.ny(); @@ -427,8 +434,8 @@ fn SAT_characteristics( .into_iter() .zip(y.slice(s![.., .., nx - 1]).gencolumns()) .zip(g.gencolumns()) - .zip(grid.detj_dxi_dx.slice(s![.., nx - 1])) - .zip(grid.detj_dxi_dy.slice(s![.., nx - 1])) + .zip(metrics.detj_dxi_dx.slice(s![.., nx - 1])) + .zip(metrics.detj_dxi_dy.slice(s![.., nx - 1])) { // East boundary, positive flux let tau = -1.0; @@ -461,8 +468,8 @@ fn SAT_characteristics( .into_iter() .zip(y.slice(s![.., .., 0]).gencolumns()) .zip(g.gencolumns()) - .zip(grid.detj_dxi_dx.slice(s![.., 0])) - .zip(grid.detj_dxi_dy.slice(s![.., 0])) + .zip(metrics.detj_dxi_dx.slice(s![.., 0])) + .zip(metrics.detj_dxi_dy.slice(s![.., 0])) { let tau = 1.0; @@ -500,8 +507,8 @@ fn SAT_characteristics( .into_iter() .zip(y.slice(s![.., ny - 1, ..]).gencolumns()) .zip(g.gencolumns()) - .zip(grid.detj_deta_dx.slice(s![ny - 1, ..])) - .zip(grid.detj_deta_dy.slice(s![ny - 1, ..])) + .zip(metrics.detj_deta_dx.slice(s![ny - 1, ..])) + .zip(metrics.detj_deta_dy.slice(s![ny - 1, ..])) { // North boundary, positive flux let tau = -1.0; @@ -533,8 +540,8 @@ fn SAT_characteristics( .into_iter() .zip(y.slice(s![.., 0, ..]).gencolumns()) .zip(g.gencolumns()) - .zip(grid.detj_deta_dx.slice(s![0, ..])) - .zip(grid.detj_deta_dy.slice(s![0, ..])) + .zip(metrics.detj_deta_dx.slice(s![0, ..])) + .zip(metrics.detj_deta_dy.slice(s![0, ..])) { // South boundary, negative flux diff --git a/sbp/src/utils.rs b/sbp/src/utils.rs index aa2d6c7..dd84fa4 100644 --- a/sbp/src/utils.rs +++ b/sbp/src/utils.rs @@ -1,10 +1,10 @@ +use crate::grid::Grid; use crate::Float; use json::JsonValue; #[derive(Debug, Clone)] -pub struct SimpleGrid { - pub x: ndarray::Array2, - pub y: ndarray::Array2, +pub struct ExtendedGrid { + pub grid: Grid, pub name: Option, pub dire: Option, pub dirw: Option, @@ -28,8 +28,8 @@ pub struct SimpleGrid { /// 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 { +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 @@ -173,9 +173,8 @@ pub fn json_to_grids(json: JsonValue) -> Result, String> { } } - Ok(SimpleGrid { - x, - y, + Ok(ExtendedGrid { + grid: Grid::new(x, y).unwrap(), name, dire, dirw, @@ -195,56 +194,73 @@ pub fn json_to_grids(json: JsonValue) -> Result, String> { #[test] fn parse_linspace() { - let grids = - json_to_grids(r#"[{"name": "main", "x": "linspace:0:10:20", "y": "linspace:0:10:21"}]"#) - .unwrap(); + 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].x.shape(), [21, 20]); - assert_eq!(grids[0].y.shape(), [21, 20]); + 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(r#"{"name": "main", "x": "linspace:0:10:20", "y": "linspace:0:10:21"}"#) - .unwrap(); + 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].x.shape(), [21, 20]); - assert_eq!(grids[0].y.shape(), [21, 20]); + 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"); } #[test] fn parse_1d() { - let grids = json_to_grids(r#"{"x": [1, 2, 3, 4, 5.1, 3], "y": [1, 2]}"#).unwrap(); + 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.x.shape(), &[2, 6]); - assert_eq!(grid.x.shape(), grid.y.shape()); + assert_eq!(grid.grid.x.shape(), &[2, 6]); + assert_eq!(grid.grid.x.shape(), grid.grid.y.shape()); } #[test] fn parse_2d() { - let grids = json_to_grids(r#"{"x": [[1, 2], [3, 4], [5.1, 3]], "y": [1, 2, 3]}"#).unwrap(); - assert_eq!(grids.len(), 1); - let grid = &grids[0]; - assert_eq!(grid.x.shape(), &[3, 2]); - assert_eq!(grid.x.shape(), grid.y.shape()); - json_to_grids(r#"{"x": [[1, 2], [3, 4], [5.1, 3], [1]], "y": [1, 2, 3]}"#).unwrap_err(); - json_to_grids(r#"{"y": [[1, 2], [3, 4], [5.1, 3], [1]], "x": [1, 2, 3]}"#).unwrap_err(); let grids = - json_to_grids(r#"{"x": [[1, 2], [3, 4], [5.1, 3]], "y": [[1, 2], [3, 4], [5, 6]]}"#) + 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); - json_to_grids(r#"{"x": [[1, 2], [3, 4], [5.1, 3]], "y": [[1, 2], [3, 4], [5]]}"#).unwrap_err(); + 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(r#"[{"#).unwrap_err(); - json_to_grids(r#"{}"#).unwrap_err(); - json_to_grids(r#"0.45"#).unwrap_err(); - json_to_grids(r#"{"x": "linspace", "y": [0.1, 0.2]}"#).unwrap_err(); - json_to_grids(r#"{"x": "linspace:::", "y": [0.1, 0.2]}"#).unwrap_err(); - json_to_grids(r#"{"x": "linspace:1.2:3.1:412.2", "y": [0.1, 0.2]}"#).unwrap_err(); - json_to_grids(r#"{"x": [-2, -3, "dfd"], "y": [0.1, 0.2]}"#).unwrap_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(); } pub fn json_to_vortex(mut json: JsonValue) -> super::euler::VortexParameters {