mulimoen
/
SummationByParts
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

separate metrics and grid

This commit is contained in:
Magnus Ulimoen 2020-04-04 00:29:02 +02:00
parent 6b94990015
commit 358f831513
6 changed files with 215 additions and 173 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

27
sbp/examples/multigrid.rs
View File

@ -14,12 +14,13 @@ struct System<T: operators::UpwindOperator> {
euler::Field, euler::Field,
)>, )>,
k: [Vec<euler::Field>; 4], k: [Vec<euler::Field>; 4],
grids: Vec<grid::Grid<T>>, grids: Vec<grid::Grid>,
metrics: Vec<grid::Metrics<T>>,
bt: Vec<euler::BoundaryCharacteristics>, bt: Vec<euler::BoundaryCharacteristics>,
} }
impl<T: operators::UpwindOperator> System<T> { impl<T: operators::UpwindOperator> System<T> {
fn new(grids: Vec<grid::Grid<T>>, bt: Vec<euler::BoundaryCharacteristics>) -> Self { fn new(grids: Vec<grid::Grid>, bt: Vec<euler::BoundaryCharacteristics>) -> Self {
let fnow = grids let fnow = grids
.iter() .iter()
.map(|g| euler::Field::new(g.ny(), g.nx())) .map(|g| euler::Field::new(g.ny(), g.nx()))
@ -33,6 +34,7 @@ impl<T: operators::UpwindOperator> System<T> {
}) })
.collect(); .collect();
let k = [fnow.clone(), fnow.clone(), fnow.clone(), fnow.clone()]; let k = [fnow.clone(), fnow.clone(), fnow.clone(), fnow.clone()];
let metrics = grids.iter().map(|g| g.metrics().unwrap()).collect();
Self { Self {
fnow, fnow,
@ -40,6 +42,7 @@ impl<T: operators::UpwindOperator> System<T> {
k, k,
wb, wb,
grids, grids,
metrics,
bt, bt,
} }
} }
@ -142,14 +145,14 @@ impl<T: operators::UpwindOperator> System<T> {
}) })
.collect::<Vec<_>>(); .collect::<Vec<_>>();
for ((((prev, fut), grid), wb), bt) in fields for ((((prev, fut), metrics), wb), bt) in fields
.iter() .iter()
.zip(fnext) .zip(fnext)
.zip(&self.grids) .zip(&self.metrics)
.zip(&mut self.wb) .zip(&mut self.wb)
.zip(bt) .zip(bt)
{ {
euler::RHS_upwind(fut, prev, grid, &bt, wb) euler::RHS_upwind(fut, prev, metrics, &bt, wb)
} }
} }
} }
@ -257,14 +260,14 @@ impl<T: operators::UpwindOperator> System<T> {
}, },
}) })
.collect::<Vec<_>>(); .collect::<Vec<_>>();
for ((((prev, fut), grid), wb), bt) in fields for ((((prev, fut), metrics), wb), bt) in fields
.iter() .iter()
.zip(&mut self.k[i]) .zip(&mut self.k[i])
.zip(&self.grids) .zip(&self.metrics)
.zip(&mut self.wb) .zip(&mut self.wb)
.zip(bt) .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()), west: determine_bc(grid.dirw.as_ref()),
}); });
} }
let mut grids: Vec<grid::Grid<operators::Upwind4>> = Vec::with_capacity(jgrids.len()); let grids = jgrids.into_iter().map(|egrid| egrid.grid).collect();
for grid in jgrids {
grids.push(grid::Grid::new(grid.x, grid.y).unwrap());
}
let integration_time: f64 = json["integration_time"].as_number().unwrap().into(); let integration_time: f64 = json["integration_time"].as_number().unwrap().into();
let vortexparams = utils::json_to_vortex(json["vortex"].clone()); let vortexparams = utils::json_to_vortex(json["vortex"].clone());
let mut sys = System::new(grids, bt); let mut sys = System::<sbp::operators::Upwind4>::new(grids, bt);
sys.vortex(0.0, vortexparams); sys.vortex(0.0, vortexparams);
let max_n = { let max_n = {

111
sbp/src/euler.rs
View File

@ -1,4 +1,4 @@
use super::grid::Grid; use super::grid::{Grid, Metrics};
use super::integrate; use super::integrate;
use super::operators::{SbpOperator, UpwindOperator}; use super::operators::{SbpOperator, UpwindOperator};
use super::Float; use super::Float;
@ -13,7 +13,7 @@ pub const GAMMA: Float = 1.4;
pub struct System<SBP: SbpOperator> { pub struct System<SBP: SbpOperator> {
sys: (Field, Field), sys: (Field, Field),
wb: WorkBuffers, wb: WorkBuffers,
grid: Grid<SBP>, grid: (Grid, Metrics<SBP>),
} }
impl<SBP: SbpOperator> System<SBP> { impl<SBP: SbpOperator> System<SBP> {
@ -21,11 +21,12 @@ impl<SBP: SbpOperator> System<SBP> {
let grid = Grid::new(x, y).expect( let grid = Grid::new(x, y).expect(
"Could not create grid. Different number of elements compared to width*height?", "Could not create grid. Different number of elements compared to width*height?",
); );
let metrics = grid.metrics().unwrap();
let nx = grid.nx(); let nx = grid.nx();
let ny = grid.ny(); let ny = grid.ny();
Self { Self {
sys: (Field::new(ny, nx), Field::new(ny, nx)), sys: (Field::new(ny, nx), Field::new(ny, nx)),
grid, grid: (grid, metrics),
wb: WorkBuffers::new(ny, nx), wb: WorkBuffers::new(ny, nx),
} }
} }
@ -37,16 +38,17 @@ impl<SBP: SbpOperator> System<SBP> {
east: BoundaryCharacteristic::This, east: BoundaryCharacteristic::This,
west: 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); let boundaries = boundary_extractor(y, grid, &bc);
RHS_trad(k, y, grid, &boundaries, wb) RHS_trad(k, y, metrics, &boundaries, wb)
}; };
integrate::rk4( integrate::rk4(
rhs_trad, rhs_trad,
&self.sys.0, &self.sys.0,
&mut self.sys.1, &mut self.sys.1,
dt, dt,
&self.grid, &self.grid.0,
&self.grid.1,
&mut self.wb.k, &mut self.wb.k,
&mut self.wb.tmp, &mut self.wb.tmp,
); );
@ -56,7 +58,7 @@ impl<SBP: SbpOperator> System<SBP> {
pub fn vortex(&mut self, t: Float, vortex_parameters: VortexParameters) { pub fn vortex(&mut self, t: Float, vortex_parameters: VortexParameters) {
self.sys self.sys
.0 .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)] #[allow(clippy::many_single_char_names)]
@ -70,12 +72,9 @@ impl<SBP: SbpOperator> System<SBP> {
mach: 0.5, mach: 0.5,
}; };
self.sys.0.vortex( self.sys
self.grid.x.view(), .0
self.grid.y.view(), .vortex(self.grid.0.x(), self.grid.0.y(), 0.0, vortex_parameters)
0.0,
vortex_parameters,
)
} }
pub fn field(&self) -> &Field { pub fn field(&self) -> &Field {
@ -83,10 +82,10 @@ impl<SBP: SbpOperator> System<SBP> {
} }
pub fn x(&self) -> ArrayView2<Float> { pub fn x(&self) -> ArrayView2<Float> {
self.grid.x.view() self.grid.0.x.view()
} }
pub fn y(&self) -> ArrayView2<Float> { pub fn y(&self) -> ArrayView2<Float> {
self.grid.y.view() self.grid.0.y.view()
} }
} }
@ -98,16 +97,18 @@ impl<UO: UpwindOperator> System<UO> {
east: BoundaryCharacteristic::This, east: BoundaryCharacteristic::This,
west: BoundaryCharacteristic::This, west: BoundaryCharacteristic::This,
}; };
let rhs_upwind = |k: &mut Field, y: &Field, grid: &Grid<_>, wb: &mut _| { let rhs_upwind =
let boundaries = boundary_extractor(y, grid, &bc); |k: &mut Field, y: &Field, grid: &Grid, metrics: &Metrics<_>, wb: &mut _| {
RHS_upwind(k, y, grid, &boundaries, wb) let boundaries = boundary_extractor(y, grid, &bc);
}; RHS_upwind(k, y, metrics, &boundaries, wb)
};
integrate::rk4( integrate::rk4(
rhs_upwind, rhs_upwind,
&self.sys.0, &self.sys.0,
&mut self.sys.1, &mut self.sys.1,
dt, dt,
&self.grid, &self.grid.0,
&self.grid.1,
&mut self.wb.k, &mut self.wb.k,
&mut self.wb.tmp, &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<SBP: SbpOperator>( pub(crate) fn RHS_trad<SBP: SbpOperator>(
k: &mut Field, k: &mut Field,
y: &Field, y: &Field,
grid: &Grid<SBP>, metrics: &Metrics<SBP>,
boundaries: &BoundaryTerms, boundaries: &BoundaryTerms,
tmp: &mut (Field, Field, Field, Field, Field, Field), tmp: &mut (Field, Field, Field, Field, Field, Field),
) { ) {
let ehat = &mut tmp.0; let ehat = &mut tmp.0;
let fhat = &mut tmp.1; let fhat = &mut tmp.1;
fluxes((ehat, fhat), y, grid); fluxes((ehat, fhat), y, metrics);
let dE = &mut tmp.2; let dE = &mut tmp.2;
let dF = &mut tmp.3; let dF = &mut tmp.3;
@ -405,24 +406,24 @@ pub(crate) fn RHS_trad<SBP: SbpOperator>(
azip!((out in &mut k.0, azip!((out in &mut k.0,
eflux in &dE.0, eflux in &dE.0,
fflux in &dF.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 *out = (-eflux - fflux)/detj
}); });
SAT_characteristics(k, y, grid, boundaries); SAT_characteristics(k, y, metrics, boundaries);
} }
#[allow(non_snake_case)] #[allow(non_snake_case)]
pub fn RHS_upwind<UO: UpwindOperator>( pub fn RHS_upwind<UO: UpwindOperator>(
k: &mut Field, k: &mut Field,
y: &Field, y: &Field,
grid: &Grid<UO>, metrics: &Metrics<UO>,
boundaries: &BoundaryTerms, boundaries: &BoundaryTerms,
tmp: &mut (Field, Field, Field, Field, Field, Field), tmp: &mut (Field, Field, Field, Field, Field, Field),
) { ) {
let ehat = &mut tmp.0; let ehat = &mut tmp.0;
let fhat = &mut tmp.1; let fhat = &mut tmp.1;
fluxes((ehat, fhat), y, grid); fluxes((ehat, fhat), y, metrics);
let dE = &mut tmp.2; let dE = &mut tmp.2;
let dF = &mut tmp.3; let dF = &mut tmp.3;
@ -438,25 +439,25 @@ pub fn RHS_upwind<UO: UpwindOperator>(
let ad_xi = &mut tmp.4; let ad_xi = &mut tmp.4;
let ad_eta = &mut tmp.5; 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, azip!((out in &mut k.0,
eflux in &dE.0, eflux in &dE.0,
fflux in &dF.0, fflux in &dF.0,
ad_xi in &ad_xi.0, ad_xi in &ad_xi.0,
ad_eta in &ad_eta.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 *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)] #[allow(clippy::many_single_char_names)]
fn upwind_dissipation<UO: UpwindOperator>( fn upwind_dissipation<UO: UpwindOperator>(
k: (&mut Field, &mut Field), k: (&mut Field, &mut Field),
y: &Field, y: &Field,
grid: &Grid<UO>, metrics: &Metrics<UO>,
tmp: (&mut Field, &mut Field), tmp: (&mut Field, &mut Field),
) { ) {
let n = y.nx() * y.ny(); let n = y.nx() * y.ny();
@ -471,11 +472,11 @@ fn upwind_dissipation<UO: UpwindOperator>(
.axis_iter(ndarray::Axis(1)) .axis_iter(ndarray::Axis(1))
.zip(tmp0.axis_iter_mut(ndarray::Axis(1))) .zip(tmp0.axis_iter_mut(ndarray::Axis(1)))
.zip(tmp1.axis_iter_mut(ndarray::Axis(1))) .zip(tmp1.axis_iter_mut(ndarray::Axis(1)))
.zip(grid.detj.iter()) .zip(metrics.detj.iter())
.zip(grid.detj_dxi_dx.iter()) .zip(metrics.detj_dxi_dx.iter())
.zip(grid.detj_dxi_dy.iter()) .zip(metrics.detj_dxi_dy.iter())
.zip(grid.detj_deta_dx.iter()) .zip(metrics.detj_deta_dx.iter())
.zip(grid.detj_deta_dy.iter()) .zip(metrics.detj_deta_dy.iter())
{ {
let rho = y[0]; let rho = y[0];
assert!(rho > 0.0); assert!(rho > 0.0);
@ -520,11 +521,11 @@ fn upwind_dissipation<UO: UpwindOperator>(
UO::disseta(tmp.1.e(), k.1.e_mut()); UO::disseta(tmp.1.e(), k.1.e_mut());
} }
fn fluxes<SBP: SbpOperator>(k: (&mut Field, &mut Field), y: &Field, grid: &Grid<SBP>) { fn fluxes<SBP: SbpOperator>(k: (&mut Field, &mut Field), y: &Field, metrics: &Metrics<SBP>) {
let j_dxi_dx = grid.detj_dxi_dx.view(); let j_dxi_dx = metrics.detj_dxi_dx.view();
let j_dxi_dy = grid.detj_dxi_dy.view(); let j_dxi_dy = metrics.detj_dxi_dy.view();
let j_deta_dx = grid.detj_deta_dx.view(); let j_deta_dx = metrics.detj_deta_dx.view();
let j_deta_dy = grid.detj_deta_dy.view(); let j_deta_dy = metrics.detj_deta_dy.view();
let rho = y.rho(); let rho = y.rho();
let rhou = y.rhou(); let rhou = y.rhou();
@ -590,9 +591,9 @@ pub struct BoundaryCharacteristics {
pub west: BoundaryCharacteristic, pub west: BoundaryCharacteristic,
} }
fn boundary_extractor<'a, SBP: SbpOperator>( fn boundary_extractor<'a>(
field: &'a Field, field: &'a Field,
_grid: &Grid<SBP>, _grid: &Grid,
bc: &BoundaryCharacteristics, bc: &BoundaryCharacteristics,
) -> BoundaryTerms<'a> { ) -> BoundaryTerms<'a> {
BoundaryTerms { BoundaryTerms {
@ -624,7 +625,7 @@ fn boundary_extractor<'a, SBP: SbpOperator>(
fn SAT_characteristics<SBP: SbpOperator>( fn SAT_characteristics<SBP: SbpOperator>(
k: &mut Field, k: &mut Field,
y: &Field, y: &Field,
grid: &Grid<SBP>, metrics: &Metrics<SBP>,
boundaries: &BoundaryTerms, boundaries: &BoundaryTerms,
) { ) {
// North boundary // North boundary
@ -640,9 +641,9 @@ fn SAT_characteristics<SBP: SbpOperator>(
hi, hi,
sign, sign,
tau, tau,
grid.detj.slice(slice), metrics.detj.slice(slice),
grid.detj_deta_dx.slice(slice), metrics.detj_deta_dx.slice(slice),
grid.detj_deta_dy.slice(slice), metrics.detj_deta_dy.slice(slice),
); );
} }
// South boundary // South boundary
@ -658,9 +659,9 @@ fn SAT_characteristics<SBP: SbpOperator>(
hi, hi,
sign, sign,
tau, tau,
grid.detj.slice(slice), metrics.detj.slice(slice),
grid.detj_deta_dx.slice(slice), metrics.detj_deta_dx.slice(slice),
grid.detj_deta_dy.slice(slice), metrics.detj_deta_dy.slice(slice),
); );
} }
// West Boundary // West Boundary
@ -676,9 +677,9 @@ fn SAT_characteristics<SBP: SbpOperator>(
hi, hi,
sign, sign,
tau, tau,
grid.detj.slice(slice), metrics.detj.slice(slice),
grid.detj_dxi_dx.slice(slice), metrics.detj_dxi_dx.slice(slice),
grid.detj_dxi_dy.slice(slice), metrics.detj_dxi_dy.slice(slice),
); );
} }
// East Boundary // East Boundary
@ -694,9 +695,9 @@ fn SAT_characteristics<SBP: SbpOperator>(
hi, hi,
sign, sign,
tau, tau,
grid.detj.slice(slice), metrics.detj.slice(slice),
grid.detj_dxi_dx.slice(slice), metrics.detj_dxi_dx.slice(slice),
grid.detj_dxi_dy.slice(slice), metrics.detj_dxi_dy.slice(slice),
); );
} }
} }

60
sbp/src/grid.rs
View File

@ -2,13 +2,16 @@ use crate::Float;
use ndarray::Array2; use ndarray::Array2;
#[derive(Debug, Clone)] #[derive(Debug, Clone)]
pub struct Grid<SBP> pub struct Grid {
pub(crate) x: Array2<Float>,
pub(crate) y: Array2<Float>,
}
#[derive(Debug, Clone)]
pub struct Metrics<SBP>
where where
SBP: super::operators::SbpOperator, SBP: super::operators::SbpOperator,
{ {
pub(crate) x: Array2<Float>,
pub(crate) y: Array2<Float>,
pub(crate) detj: Array2<Float>, pub(crate) detj: Array2<Float>,
pub(crate) detj_dxi_dx: Array2<Float>, pub(crate) detj_dxi_dx: Array2<Float>,
pub(crate) detj_dxi_dy: Array2<Float>, pub(crate) detj_dxi_dy: Array2<Float>,
@ -18,11 +21,39 @@ where
operator: std::marker::PhantomData<SBP>, operator: std::marker::PhantomData<SBP>,
} }
impl<SBP: super::operators::SbpOperator> Grid<SBP> { impl Grid {
pub fn new(x: Array2<Float>, y: Array2<Float>) -> Result<Self, ndarray::ShapeError> { pub fn new(x: Array2<Float>, y: Array2<Float>) -> Result<Self, ndarray::ShapeError> {
assert_eq!(x.shape(), y.shape()); 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<Float> {
self.x.view()
}
pub fn y(&self) -> ndarray::ArrayView2<Float> {
self.y.view()
}
pub fn metrics<SBP: super::operators::SbpOperator>(
&self,
) -> Result<Metrics<SBP>, ndarray::ShapeError> {
Metrics::new(self)
}
}
impl<SBP: super::operators::SbpOperator> Metrics<SBP> {
fn new(grid: &Grid) -> Result<Self, ndarray::ShapeError> {
let ny = grid.ny();
let nx = grid.nx();
let x = &grid.x;
let y = &grid.y;
let mut dx_dxi = Array2::zeros((ny, nx)); let mut dx_dxi = Array2::zeros((ny, nx));
SBP::diffxi(x.view(), dx_dxi.view_mut()); SBP::diffxi(x.view(), dx_dxi.view_mut());
@ -57,8 +88,6 @@ impl<SBP: super::operators::SbpOperator> Grid<SBP> {
let detj_deta_dy = dx_dxi; let detj_deta_dy = dx_dxi;
Ok(Self { Ok(Self {
x,
y,
detj, detj,
detj_dxi_dx, detj_dxi_dx,
detj_dxi_dy, detj_dxi_dy,
@ -67,17 +96,4 @@ impl<SBP: super::operators::SbpOperator> Grid<SBP> {
operator: std::marker::PhantomData, 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<Float> {
self.x.view()
}
pub fn y(&self) -> ndarray::ArrayView2<Float> {
self.y.view()
}
} }

9
sbp/src/integrate.rs
View File

@ -1,4 +1,4 @@
use super::grid::Grid; use super::grid::{Grid, Metrics};
use super::operators::SbpOperator; use super::operators::SbpOperator;
use super::Float; use super::Float;
use ndarray::{Array3, Zip}; use ndarray::{Array3, Zip};
@ -8,13 +8,14 @@ pub(crate) fn rk4<'a, F: 'a, SBP, RHS, WB>(
prev: &F, prev: &F,
fut: &mut F, fut: &mut F,
dt: Float, dt: Float,
grid: &Grid<SBP>, grid: &Grid,
metrics: &Metrics<SBP>,
k: &mut [F; 4], k: &mut [F; 4],
mut wb: &mut WB, mut wb: &mut WB,
) where ) where
F: std::ops::Deref<Target = Array3<Float>> + std::ops::DerefMut<Target = Array3<Float>>, F: std::ops::Deref<Target = Array3<Float>> + std::ops::DerefMut<Target = Array3<Float>>,
SBP: SbpOperator, SBP: SbpOperator,
RHS: Fn(&mut F, &F, &Grid<SBP>, &mut WB), RHS: Fn(&mut F, &F, &Grid, &Metrics<SBP>, &mut WB),
{ {
assert_eq!(prev.shape(), fut.shape()); 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);
} }
} }

91
sbp/src/maxwell.rs
View File

@ -1,4 +1,4 @@
use super::grid::Grid; use super::grid::{Grid, Metrics};
use super::integrate; use super::integrate;
use super::operators::{SbpOperator, UpwindOperator}; use super::operators::{SbpOperator, UpwindOperator};
use crate::Float; use crate::Float;
@ -77,7 +77,8 @@ impl Field {
pub struct System<SBP: SbpOperator> { pub struct System<SBP: SbpOperator> {
sys: (Field, Field), sys: (Field, Field),
wb: WorkBuffers, wb: WorkBuffers,
grid: Grid<SBP>, grid: Grid,
metrics: Metrics<SBP>,
} }
impl<SBP: SbpOperator> System<SBP> { impl<SBP: SbpOperator> System<SBP> {
@ -87,10 +88,12 @@ impl<SBP: SbpOperator> System<SBP> {
let nx = x.shape()[1]; let nx = x.shape()[1];
let grid = Grid::new(x, y).unwrap(); let grid = Grid::new(x, y).unwrap();
let metrics = grid.metrics().unwrap();
Self { Self {
sys: (Field::new(ny, nx), Field::new(ny, nx)), sys: (Field::new(ny, nx), Field::new(ny, nx)),
grid, grid,
metrics,
wb: WorkBuffers::new(ny, nx), wb: WorkBuffers::new(ny, nx),
} }
} }
@ -118,6 +121,7 @@ impl<SBP: SbpOperator> System<SBP> {
&mut self.sys.1, &mut self.sys.1,
dt, dt,
&self.grid, &self.grid,
&self.metrics,
&mut self.wb.k, &mut self.wb.k,
&mut self.wb.tmp, &mut self.wb.tmp,
); );
@ -134,6 +138,7 @@ impl<UO: UpwindOperator> System<UO> {
&mut self.sys.1, &mut self.sys.1,
dt, dt,
&self.grid, &self.grid,
&self.metrics,
&mut self.wb.k, &mut self.wb.k,
&mut self.wb.tmp, &mut self.wb.tmp,
); );
@ -171,10 +176,11 @@ fn gaussian(x: Float, x0: Float, y: Float, y0: Float) -> Float {
fn RHS<SBP: SbpOperator>( fn RHS<SBP: SbpOperator>(
k: &mut Field, k: &mut Field,
y: &Field, y: &Field,
grid: &Grid<SBP>, _grid: &Grid,
metrics: &Metrics<SBP>,
tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>), tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
) { ) {
fluxes(k, y, grid, tmp); fluxes(k, y, metrics, tmp);
let boundaries = BoundaryTerms { let boundaries = BoundaryTerms {
north: Boundary::This, north: Boundary::This,
@ -182,10 +188,10 @@ fn RHS<SBP: SbpOperator>(
west: Boundary::This, west: Boundary::This,
east: Boundary::This, east: Boundary::This,
}; };
SAT_characteristics(k, y, grid, &boundaries); SAT_characteristics(k, y, metrics, &boundaries);
azip!((k in &mut k.0, 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; *k /= detj;
}); });
} }
@ -194,11 +200,12 @@ fn RHS<SBP: SbpOperator>(
fn RHS_upwind<UO: UpwindOperator>( fn RHS_upwind<UO: UpwindOperator>(
k: &mut Field, k: &mut Field,
y: &Field, y: &Field,
grid: &Grid<UO>, _grid: &Grid,
metrics: &Metrics<UO>,
tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>), tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
) { ) {
fluxes(k, y, grid, tmp); fluxes(k, y, metrics, tmp);
dissipation(k, y, grid, tmp); dissipation(k, y, metrics, tmp);
let boundaries = BoundaryTerms { let boundaries = BoundaryTerms {
north: Boundary::This, north: Boundary::This,
@ -206,10 +213,10 @@ fn RHS_upwind<UO: UpwindOperator>(
west: Boundary::This, west: Boundary::This,
east: Boundary::This, east: Boundary::This,
}; };
SAT_characteristics(k, y, grid, &boundaries); SAT_characteristics(k, y, metrics, &boundaries);
azip!((k in &mut k.0, 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; *k /= detj;
}); });
} }
@ -217,20 +224,20 @@ fn RHS_upwind<UO: UpwindOperator>(
fn fluxes<SBP: SbpOperator>( fn fluxes<SBP: SbpOperator>(
k: &mut Field, k: &mut Field,
y: &Field, y: &Field,
grid: &Grid<SBP>, metrics: &Metrics<SBP>,
tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>), tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
) { ) {
// ex = hz_y // ex = hz_y
{ {
ndarray::azip!((a in &mut tmp.0, 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()) &hz in &y.hz())
*a = dxi_dy * hz *a = dxi_dy * hz
); );
SBP::diffxi(tmp.0.view(), tmp.1.view_mut()); SBP::diffxi(tmp.0.view(), tmp.1.view_mut());
ndarray::azip!((b in &mut tmp.2, 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()) &hz in &y.hz())
*b = deta_dy * hz *b = deta_dy * hz
); );
@ -244,8 +251,8 @@ fn fluxes<SBP: SbpOperator>(
{ {
// hz = -ey_x + ex_y // hz = -ey_x + ex_y
ndarray::azip!((a in &mut tmp.0, ndarray::azip!((a in &mut tmp.0,
&dxi_dx in &grid.detj_dxi_dx, &dxi_dx in &metrics.detj_dxi_dx,
&dxi_dy in &grid.detj_dxi_dy, &dxi_dy in &metrics.detj_dxi_dy,
&ex in &y.ex(), &ex in &y.ex(),
&ey in &y.ey()) &ey in &y.ey())
*a = dxi_dx * -ey + dxi_dy * ex *a = dxi_dx * -ey + dxi_dy * ex
@ -253,8 +260,8 @@ fn fluxes<SBP: SbpOperator>(
SBP::diffxi(tmp.0.view(), tmp.1.view_mut()); SBP::diffxi(tmp.0.view(), tmp.1.view_mut());
ndarray::azip!((b in &mut tmp.2, ndarray::azip!((b in &mut tmp.2,
&deta_dx in &grid.detj_deta_dx, &deta_dx in &metrics.detj_deta_dx,
&deta_dy in &grid.detj_deta_dy, &deta_dy in &metrics.detj_deta_dy,
&ex in &y.ex(), &ex in &y.ex(),
&ey in &y.ey()) &ey in &y.ey())
*b = deta_dx * -ey + deta_dy * ex *b = deta_dx * -ey + deta_dy * ex
@ -269,14 +276,14 @@ fn fluxes<SBP: SbpOperator>(
// ey = -hz_x // ey = -hz_x
{ {
ndarray::azip!((a in &mut tmp.0, 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()) &hz in &y.hz())
*a = dxi_dx * -hz *a = dxi_dx * -hz
); );
SBP::diffxi(tmp.0.view(), tmp.1.view_mut()); SBP::diffxi(tmp.0.view(), tmp.1.view_mut());
azip!((b in &mut tmp.2, azip!((b in &mut tmp.2,
&deta_dx in &grid.detj_deta_dx, &deta_dx in &metrics.detj_deta_dx,
&hz in &y.hz()) &hz in &y.hz())
*b = deta_dx * -hz *b = deta_dx * -hz
); );
@ -291,14 +298,14 @@ fn fluxes<SBP: SbpOperator>(
fn dissipation<UO: UpwindOperator>( fn dissipation<UO: UpwindOperator>(
k: &mut Field, k: &mut Field,
y: &Field, y: &Field,
grid: &Grid<UO>, metrics: &Metrics<UO>,
tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>), tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
) { ) {
// ex component // ex component
{ {
ndarray::azip!((a in &mut tmp.0, ndarray::azip!((a in &mut tmp.0,
&kx in &grid.detj_dxi_dx, &kx in &metrics.detj_dxi_dx,
&ky in &grid.detj_dxi_dy, &ky in &metrics.detj_dxi_dy,
&ex in &y.ex(), &ex in &y.ex(),
&ey in &y.ey()) { &ey in &y.ey()) {
let r = Float::hypot(kx, ky); let r = Float::hypot(kx, ky);
@ -307,8 +314,8 @@ fn dissipation<UO: UpwindOperator>(
UO::dissxi(tmp.0.view(), tmp.1.view_mut()); UO::dissxi(tmp.0.view(), tmp.1.view_mut());
ndarray::azip!((b in &mut tmp.2, ndarray::azip!((b in &mut tmp.2,
&kx in &grid.detj_deta_dx, &kx in &metrics.detj_deta_dx,
&ky in &grid.detj_deta_dy, &ky in &metrics.detj_deta_dy,
&ex in &y.ex(), &ex in &y.ex(),
&ey in &y.ey()) { &ey in &y.ey()) {
let r = Float::hypot(kx, ky); let r = Float::hypot(kx, ky);
@ -324,8 +331,8 @@ fn dissipation<UO: UpwindOperator>(
// hz component // hz component
{ {
ndarray::azip!((a in &mut tmp.0, ndarray::azip!((a in &mut tmp.0,
&kx in &grid.detj_dxi_dx, &kx in &metrics.detj_dxi_dx,
&ky in &grid.detj_dxi_dy, &ky in &metrics.detj_dxi_dy,
&hz in &y.hz()) { &hz in &y.hz()) {
let r = Float::hypot(kx, ky); let r = Float::hypot(kx, ky);
*a = r * hz; *a = r * hz;
@ -333,8 +340,8 @@ fn dissipation<UO: UpwindOperator>(
UO::dissxi(tmp.0.view(), tmp.1.view_mut()); UO::dissxi(tmp.0.view(), tmp.1.view_mut());
ndarray::azip!((b in &mut tmp.2, ndarray::azip!((b in &mut tmp.2,
&kx in &grid.detj_deta_dx, &kx in &metrics.detj_deta_dx,
&ky in &grid.detj_deta_dy, &ky in &metrics.detj_deta_dy,
&hz in &y.hz()) { &hz in &y.hz()) {
let r = Float::hypot(kx, ky); let r = Float::hypot(kx, ky);
*b = r * hz; *b = r * hz;
@ -349,8 +356,8 @@ fn dissipation<UO: UpwindOperator>(
// ey // ey
{ {
ndarray::azip!((a in &mut tmp.0, ndarray::azip!((a in &mut tmp.0,
&kx in &grid.detj_dxi_dx, &kx in &metrics.detj_dxi_dx,
&ky in &grid.detj_dxi_dy, &ky in &metrics.detj_dxi_dy,
&ex in &y.ex(), &ex in &y.ex(),
&ey in &y.ey()) { &ey in &y.ey()) {
let r = Float::hypot(kx, ky); let r = Float::hypot(kx, ky);
@ -359,8 +366,8 @@ fn dissipation<UO: UpwindOperator>(
UO::dissxi(tmp.0.view(), tmp.1.view_mut()); UO::dissxi(tmp.0.view(), tmp.1.view_mut());
ndarray::azip!((b in &mut tmp.2, ndarray::azip!((b in &mut tmp.2,
&kx in &grid.detj_deta_dx, &kx in &metrics.detj_deta_dx,
&ky in &grid.detj_deta_dy, &ky in &metrics.detj_deta_dy,
&ex in &y.ex(), &ex in &y.ex(),
&ey in &y.ey()) { &ey in &y.ey()) {
let r = Float::hypot(kx, ky); let r = Float::hypot(kx, ky);
@ -392,7 +399,7 @@ pub struct BoundaryTerms {
fn SAT_characteristics<SBP: SbpOperator>( fn SAT_characteristics<SBP: SbpOperator>(
k: &mut Field, k: &mut Field,
y: &Field, y: &Field,
grid: &Grid<SBP>, metrics: &Metrics<SBP>,
boundaries: &BoundaryTerms, boundaries: &BoundaryTerms,
) { ) {
let ny = y.ny(); let ny = y.ny();
@ -427,8 +434,8 @@ fn SAT_characteristics<SBP: SbpOperator>(
.into_iter() .into_iter()
.zip(y.slice(s![.., .., nx - 1]).gencolumns()) .zip(y.slice(s![.., .., nx - 1]).gencolumns())
.zip(g.gencolumns()) .zip(g.gencolumns())
.zip(grid.detj_dxi_dx.slice(s![.., nx - 1])) .zip(metrics.detj_dxi_dx.slice(s![.., nx - 1]))
.zip(grid.detj_dxi_dy.slice(s![.., nx - 1])) .zip(metrics.detj_dxi_dy.slice(s![.., nx - 1]))
{ {
// East boundary, positive flux // East boundary, positive flux
let tau = -1.0; let tau = -1.0;
@ -461,8 +468,8 @@ fn SAT_characteristics<SBP: SbpOperator>(
.into_iter() .into_iter()
.zip(y.slice(s![.., .., 0]).gencolumns()) .zip(y.slice(s![.., .., 0]).gencolumns())
.zip(g.gencolumns()) .zip(g.gencolumns())
.zip(grid.detj_dxi_dx.slice(s![.., 0])) .zip(metrics.detj_dxi_dx.slice(s![.., 0]))
.zip(grid.detj_dxi_dy.slice(s![.., 0])) .zip(metrics.detj_dxi_dy.slice(s![.., 0]))
{ {
let tau = 1.0; let tau = 1.0;
@ -500,8 +507,8 @@ fn SAT_characteristics<SBP: SbpOperator>(
.into_iter() .into_iter()
.zip(y.slice(s![.., ny - 1, ..]).gencolumns()) .zip(y.slice(s![.., ny - 1, ..]).gencolumns())
.zip(g.gencolumns()) .zip(g.gencolumns())
.zip(grid.detj_deta_dx.slice(s![ny - 1, ..])) .zip(metrics.detj_deta_dx.slice(s![ny - 1, ..]))
.zip(grid.detj_deta_dy.slice(s![ny - 1, ..])) .zip(metrics.detj_deta_dy.slice(s![ny - 1, ..]))
{ {
// North boundary, positive flux // North boundary, positive flux
let tau = -1.0; let tau = -1.0;
@ -533,8 +540,8 @@ fn SAT_characteristics<SBP: SbpOperator>(
.into_iter() .into_iter()
.zip(y.slice(s![.., 0, ..]).gencolumns()) .zip(y.slice(s![.., 0, ..]).gencolumns())
.zip(g.gencolumns()) .zip(g.gencolumns())
.zip(grid.detj_deta_dx.slice(s![0, ..])) .zip(metrics.detj_deta_dx.slice(s![0, ..]))
.zip(grid.detj_deta_dy.slice(s![0, ..])) .zip(metrics.detj_deta_dy.slice(s![0, ..]))
{ {
// South boundary, negative flux // South boundary, negative flux

90
sbp/src/utils.rs
View File

@ -1,10 +1,10 @@
use crate::grid::Grid;
use crate::Float; use crate::Float;
use json::JsonValue; use json::JsonValue;
#[derive(Debug, Clone)] #[derive(Debug, Clone)]
pub struct SimpleGrid { pub struct ExtendedGrid {
pub x: ndarray::Array2<Float>, pub grid: Grid,
pub y: ndarray::Array2<Float>,
pub name: Option<String>, pub name: Option<String>,
pub dire: Option<String>, pub dire: Option<String>,
pub dirw: Option<String>, pub dirw: Option<String>,
@ -28,8 +28,8 @@ pub struct SimpleGrid {
/// Optional parameters: /// Optional parameters:
/// * name (for relating boundaries) /// * name (for relating boundaries)
/// * dir{e,w,n,s} (for boundary terms) /// * dir{e,w,n,s} (for boundary terms)
pub fn json_to_grids(json: JsonValue) -> Result<Vec<SimpleGrid>, String> { pub fn json_to_grids(json: JsonValue) -> Result<Vec<ExtendedGrid>, String> {
fn json_to_grid(mut grid: JsonValue) -> Result<SimpleGrid, String> { fn json_to_grid(mut grid: JsonValue) -> Result<ExtendedGrid, String> {
#[derive(Debug)] #[derive(Debug)]
enum ArrayForm { enum ArrayForm {
/// Only know the one dimension, will broadcast to /// Only know the one dimension, will broadcast to
@ -173,9 +173,8 @@ pub fn json_to_grids(json: JsonValue) -> Result<Vec<SimpleGrid>, String> {
} }
} }
Ok(SimpleGrid { Ok(ExtendedGrid {
x, grid: Grid::new(x, y).unwrap(),
y,
name, name,
dire, dire,
dirw, dirw,
@ -195,56 +194,73 @@ pub fn json_to_grids(json: JsonValue) -> Result<Vec<SimpleGrid>, String> {
#[test] #[test]
fn parse_linspace() { fn parse_linspace() {
let grids = let grids = json_to_grids(
json_to_grids(r#"[{"name": "main", "x": "linspace:0:10:20", "y": "linspace:0:10:21"}]"#) json::parse(r#"[{"name": "main", "x": "linspace:0:10:20", "y": "linspace:0:10:21"}]"#)
.unwrap(); .unwrap(),
)
.unwrap();
assert_eq!(grids.len(), 1); assert_eq!(grids.len(), 1);
assert_eq!(grids[0].x.shape(), [21, 20]); assert_eq!(grids[0].grid.x.shape(), [21, 20]);
assert_eq!(grids[0].y.shape(), [21, 20]); assert_eq!(grids[0].grid.y.shape(), [21, 20]);
assert_eq!(grids[0].name.as_ref().unwrap(), "main"); assert_eq!(grids[0].name.as_ref().unwrap(), "main");
let grids = let grids = json_to_grids(
json_to_grids(r#"{"name": "main", "x": "linspace:0:10:20", "y": "linspace:0:10:21"}"#) json::parse(r#"{"name": "main", "x": "linspace:0:10:20", "y": "linspace:0:10:21"}"#)
.unwrap(); .unwrap(),
)
.unwrap();
assert_eq!(grids.len(), 1); assert_eq!(grids.len(), 1);
assert_eq!(grids[0].x.shape(), [21, 20]); assert_eq!(grids[0].grid.x.shape(), [21, 20]);
assert_eq!(grids[0].y.shape(), [21, 20]); assert_eq!(grids[0].grid.y.shape(), [21, 20]);
assert_eq!(grids[0].name.as_ref().unwrap(), "main"); assert_eq!(grids[0].name.as_ref().unwrap(), "main");
} }
#[test] #[test]
fn parse_1d() { 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); assert_eq!(grids.len(), 1);
let grid = &grids[0]; let grid = &grids[0];
assert_eq!(grid.x.shape(), &[2, 6]); assert_eq!(grid.grid.x.shape(), &[2, 6]);
assert_eq!(grid.x.shape(), grid.y.shape()); assert_eq!(grid.grid.x.shape(), grid.grid.y.shape());
} }
#[test] #[test]
fn parse_2d() { 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 = 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(); .unwrap();
assert_eq!(grids.len(), 1); 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] #[test]
fn parse_err() { fn parse_err() {
json_to_grids(r#"[{"#).unwrap_err(); json_to_grids(json::parse(r#"{}"#).unwrap()).unwrap_err();
json_to_grids(r#"{}"#).unwrap_err(); json_to_grids(json::parse(r#"0.45"#).unwrap()).unwrap_err();
json_to_grids(r#"0.45"#).unwrap_err(); json_to_grids(json::parse(r#"{"x": "linspace", "y": [0.1, 0.2]}"#).unwrap()).unwrap_err();
json_to_grids(r#"{"x": "linspace", "y": [0.1, 0.2]}"#).unwrap_err(); json_to_grids(json::parse(r#"{"x": "linspace:::", "y": [0.1, 0.2]}"#).unwrap()).unwrap_err();
json_to_grids(r#"{"x": "linspace:::", "y": [0.1, 0.2]}"#).unwrap_err(); json_to_grids(json::parse(r#"{"x": "linspace:1.2:3.1:412.2", "y": [0.1, 0.2]}"#).unwrap())
json_to_grids(r#"{"x": "linspace:1.2:3.1:412.2", "y": [0.1, 0.2]}"#).unwrap_err(); .unwrap_err();
json_to_grids(r#"{"x": [-2, -3, "dfd"], "y": [0.1, 0.2]}"#).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 { pub fn json_to_vortex(mut json: JsonValue) -> super::euler::VortexParameters {