From 85f3c46430a9ac48eab541d4d0c4c2d0d0945e4c Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Tue, 14 Apr 2020 21:59:02 +0200 Subject: [PATCH] SbpOperator takes &self --- sbp/benches/euler.rs | 6 +- sbp/benches/maxwell.rs | 6 +- sbp/src/euler.rs | 142 +++++++++++++++++------------- sbp/src/grid.rs | 37 ++++---- sbp/src/maxwell.rs | 125 ++++++++++++++------------ sbp/src/operators.rs | 31 +++---- sbp/src/operators/traditional4.rs | 23 +++-- sbp/src/operators/traditional8.rs | 28 +++--- sbp/src/operators/upwind4.rs | 48 +++++----- sbp/src/operators/upwind4h2.rs | 18 ++-- sbp/src/operators/upwind9.rs | 30 ++++--- sbp/src/operators/upwind9h2.rs | 18 ++-- webfront/lib.rs | 6 +- 13 files changed, 285 insertions(+), 233 deletions(-) diff --git a/sbp/benches/euler.rs b/sbp/benches/euler.rs index 07fa2a4..2ed59b8 100644 --- a/sbp/benches/euler.rs +++ b/sbp/benches/euler.rs @@ -25,7 +25,7 @@ fn performance_benchmark(c: &mut Criterion) { let y = ndarray::Array1::linspace(-10.0, 10.0, h); let y = y.broadcast((w, h)).unwrap().reversed_axes(); - let mut universe = System::::new(x.into_owned(), y.into_owned()); + let mut universe = System::new(x.into_owned(), y.into_owned(), Upwind4); group.bench_function("advance", |b| { b.iter(|| { universe.init_with_vortex(0.0, 0.0); @@ -33,7 +33,7 @@ fn performance_benchmark(c: &mut Criterion) { }) }); - let mut universe = System::::new(x.into_owned(), y.into_owned()); + let mut universe = System::new(x.into_owned(), y.into_owned(), Upwind4); group.bench_function("advance_upwind", |b| { b.iter(|| { universe.init_with_vortex(0.0, 0.0); @@ -41,7 +41,7 @@ fn performance_benchmark(c: &mut Criterion) { }) }); - let mut universe = System::::new(x.into_owned(), y.into_owned()); + let mut universe = System::new(x.into_owned(), y.into_owned(), SBP4); group.bench_function("advance_trad4", |b| { b.iter(|| { universe.init_with_vortex(0.0, 0.0); diff --git a/sbp/benches/maxwell.rs b/sbp/benches/maxwell.rs index f532e5e..70005ca 100644 --- a/sbp/benches/maxwell.rs +++ b/sbp/benches/maxwell.rs @@ -24,7 +24,7 @@ fn performance_benchmark(c: &mut Criterion) { let x = ndarray::Array2::from_shape_fn((h, w), |(_, i)| i as Float / (w - 1) as Float); let y = ndarray::Array2::from_shape_fn((h, w), |(j, _)| j as Float / (h - 1) as Float); - let mut universe = System::::new(x.clone(), y.clone()); + let mut universe = System::new(x.clone(), y.clone(), Upwind4); group.bench_function("advance", |b| { b.iter(|| { universe.set_gaussian(0.5, 0.5); @@ -32,7 +32,7 @@ fn performance_benchmark(c: &mut Criterion) { }) }); - let mut universe = System::::new(x.clone(), y.clone()); + let mut universe = System::new(x.clone(), y.clone(), Upwind4); group.bench_function("advance_upwind", |b| { b.iter(|| { universe.set_gaussian(0.5, 0.5); @@ -40,7 +40,7 @@ fn performance_benchmark(c: &mut Criterion) { }) }); - let mut universe = System::::new(x, y); + let mut universe = System::new(x, y, SBP4); group.bench_function("advance_trad4", |b| { b.iter(|| { universe.set_gaussian(0.5, 0.5); diff --git a/sbp/src/euler.rs b/sbp/src/euler.rs index 08b824a..a988159 100644 --- a/sbp/src/euler.rs +++ b/sbp/src/euler.rs @@ -14,15 +14,16 @@ pub struct System { sys: (Field, Field), k: [Field; 4], wb: WorkBuffers, - grid: (Grid, Metrics), + grid: (Grid, Metrics), + op: SBP, } impl System { - pub fn new(x: ndarray::Array2, y: ndarray::Array2) -> Self { + pub fn new(x: ndarray::Array2, y: ndarray::Array2, op: SBP) -> Self { 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 metrics = grid.metrics(op, op).unwrap(); let nx = grid.nx(); let ny = grid.ny(); Self { @@ -35,6 +36,7 @@ impl System { Field::new(ny, nx), ], wb: WorkBuffers::new(ny, nx), + op, } } @@ -45,10 +47,11 @@ impl System { east: BoundaryCharacteristic::This, west: BoundaryCharacteristic::This, }; + let op = self.op; let rhs_trad = |k: &mut Field, y: &Field, _time: Float, gm: &(_, _), wb: &mut _| { let (grid, metrics) = gm; let boundaries = boundary_extractor(y, grid, &bc); - RHS_trad(k, y, metrics, &boundaries, wb) + RHS_trad((op, op), k, y, metrics, &boundaries, wb) }; integrate::integrate::( rhs_trad, @@ -105,10 +108,11 @@ impl System { east: BoundaryCharacteristic::This, west: BoundaryCharacteristic::This, }; + let op = self.op; let rhs_upwind = |k: &mut Field, y: &Field, _time: Float, gm: &(_, _), wb: &mut _| { let (grid, metrics) = gm; let boundaries = boundary_extractor(y, grid, &bc); - RHS_upwind(k, y, metrics, &boundaries, wb) + RHS_upwind((op, op), k, y, metrics, &boundaries, wb) }; integrate::integrate::( rhs_upwind, @@ -268,12 +272,14 @@ impl Field { impl Field { /// sqrt((self-other)^T*H*(self-other)) - pub fn h2_err(&self, other: &Self) -> Float { + pub fn h2_err( + &self, + (opeta, opxi): (SBPeta, SBPxi), + other: &Self, + ) -> Float { assert_eq!(self.nx(), other.nx()); assert_eq!(self.ny(), other.ny()); - let h = SBP::h(); - // Resulting structure should be // serialized(F0 - F1)^T (Hx kron Hy) serialized(F0 - F1) // @@ -282,7 +288,7 @@ impl Field { // This chains the h block into the form [h, 1, 1, 1, rev(h)], // and multiplies with a factor - let itermaker = move |n: usize, factor: Float| { + let itermaker = move |h: &'static [Float], n: usize, factor: Float| { h.iter() .copied() .chain(std::iter::repeat(1.0).take(n - 2 * h.len())) @@ -291,8 +297,9 @@ impl Field { }; let hxiterator = itermaker( + opxi.h(), self.nx(), - if SBP::is_h2() { + if opxi.is_h2() { 1.0 / (self.nx() - 2) as Float } else { 1.0 / (self.nx() - 1) as Float @@ -303,8 +310,9 @@ impl Field { let hxiterator = hxiterator.cycle().take(self.nx() * self.ny()); let hyiterator = itermaker( + opeta.h(), self.ny(), - 1.0 / if SBP::is_h2() { + 1.0 / if opeta.is_h2() { (self.ny() - 2) as Float } else { (self.ny() - 1) as Float @@ -333,10 +341,12 @@ fn h2_diff() { } let field1 = Field::new(20, 21); - assert!((field0.h2_err::(&field1).powi(2) - 4.0).abs() < 1e-3); - assert!((field0.h2_err::(&field1).powi(2) - 4.0).abs() < 1e-3); - assert!((field0.h2_err::(&field1).powi(2) - 4.0).abs() < 1e-3); - assert!((field0.h2_err::(&field1).powi(2) - 4.0).abs() < 1e-3); + use super::operators::{Upwind4, Upwind9, SBP4, SBP8}; + + assert!((field0.h2_err((Upwind4, Upwind4), &field1).powi(2) - 4.0).abs() < 1e-3); + assert!((field0.h2_err((Upwind9, Upwind9), &field1).powi(2) - 4.0).abs() < 1e-3); + assert!((field0.h2_err((SBP4, SBP4), &field1).powi(2) - 4.0).abs() < 1e-3); + assert!((field0.h2_err((SBP8, SBP8), &field1).powi(2) - 4.0).abs() < 1e-3); } #[derive(Copy, Clone, Debug)] @@ -401,9 +411,10 @@ fn pressure(gamma: Float, rho: Float, rhou: Float, rhov: Float, e: Float) -> Flo #[allow(non_snake_case)] pub fn RHS_trad( + (opeta, opxi): (SBPeta, SBPxi), k: &mut Field, y: &Field, - metrics: &Metrics, + metrics: &Metrics, boundaries: &BoundaryTerms, tmp: &mut (Field, Field, Field, Field, Field, Field), ) { @@ -413,15 +424,15 @@ pub fn RHS_trad( let dE = &mut tmp.2; let dF = &mut tmp.3; - SBPxi::diffxi(ehat.rho(), dE.rho_mut()); - SBPxi::diffxi(ehat.rhou(), dE.rhou_mut()); - SBPxi::diffxi(ehat.rhov(), dE.rhov_mut()); - SBPxi::diffxi(ehat.e(), dE.e_mut()); + opxi.diffxi(ehat.rho(), dE.rho_mut()); + opxi.diffxi(ehat.rhou(), dE.rhou_mut()); + opxi.diffxi(ehat.rhov(), dE.rhov_mut()); + opxi.diffxi(ehat.e(), dE.e_mut()); - SBPeta::diffeta(fhat.rho(), dF.rho_mut()); - SBPeta::diffeta(fhat.rhou(), dF.rhou_mut()); - SBPeta::diffeta(fhat.rhov(), dF.rhov_mut()); - SBPeta::diffeta(fhat.e(), dF.e_mut()); + opeta.diffeta(fhat.rho(), dF.rho_mut()); + opeta.diffeta(fhat.rhou(), dF.rhou_mut()); + opeta.diffeta(fhat.rhov(), dF.rhov_mut()); + opeta.diffeta(fhat.e(), dF.e_mut()); azip!((out in &mut k.0, eflux in &dE.0, @@ -430,14 +441,15 @@ pub fn RHS_trad( *out = (-eflux - fflux)/detj }); - SAT_characteristics::(k, y, metrics, boundaries); + SAT_characteristics((opeta, opxi), k, y, metrics, boundaries); } #[allow(non_snake_case)] pub fn RHS_upwind( + (opeta, opxi): (UOeta, UOxi), k: &mut Field, y: &Field, - metrics: &Metrics, + metrics: &Metrics, boundaries: &BoundaryTerms, tmp: &mut (Field, Field, Field, Field, Field, Field), ) { @@ -447,19 +459,25 @@ pub fn RHS_upwind( let dE = &mut tmp.2; let dF = &mut tmp.3; - UOxi::diffxi(ehat.rho(), dE.rho_mut()); - UOxi::diffxi(ehat.rhou(), dE.rhou_mut()); - UOxi::diffxi(ehat.rhov(), dE.rhov_mut()); - UOxi::diffxi(ehat.e(), dE.e_mut()); + opxi.diffxi(ehat.rho(), dE.rho_mut()); + opxi.diffxi(ehat.rhou(), dE.rhou_mut()); + opxi.diffxi(ehat.rhov(), dE.rhov_mut()); + opxi.diffxi(ehat.e(), dE.e_mut()); - UOeta::diffeta(fhat.rho(), dF.rho_mut()); - UOeta::diffeta(fhat.rhou(), dF.rhou_mut()); - UOeta::diffeta(fhat.rhov(), dF.rhov_mut()); - UOeta::diffeta(fhat.e(), dF.e_mut()); + opeta.diffeta(fhat.rho(), dF.rho_mut()); + opeta.diffeta(fhat.rhou(), dF.rhou_mut()); + opeta.diffeta(fhat.rhov(), dF.rhov_mut()); + opeta.diffeta(fhat.e(), dF.e_mut()); let ad_xi = &mut tmp.4; let ad_eta = &mut tmp.5; - upwind_dissipation::((ad_xi, ad_eta), y, metrics, (&mut tmp.0, &mut tmp.1)); + upwind_dissipation( + (opeta, opxi), + (ad_xi, ad_eta), + y, + metrics, + (&mut tmp.0, &mut tmp.1), + ); azip!((out in &mut k.0, eflux in &dE.0, @@ -470,14 +488,15 @@ pub fn RHS_upwind( *out = (-eflux - fflux + ad_xi + ad_eta)/detj }); - SAT_characteristics::(k, y, metrics, boundaries); + SAT_characteristics((opeta, opxi), k, y, metrics, boundaries); } #[allow(clippy::many_single_char_names)] fn upwind_dissipation( + (opeta, opxi): (UOeta, UOxi), k: (&mut Field, &mut Field), y: &Field, - metrics: &Metrics, + metrics: &Metrics, tmp: (&mut Field, &mut Field), ) { let n = y.nx() * y.ny(); @@ -530,22 +549,18 @@ fn upwind_dissipation( tmp1[3] = alpha_v * e * detj; } - UOxi::dissxi(tmp.0.rho(), k.0.rho_mut()); - UOxi::dissxi(tmp.0.rhou(), k.0.rhou_mut()); - UOxi::dissxi(tmp.0.rhov(), k.0.rhov_mut()); - UOxi::dissxi(tmp.0.e(), k.0.e_mut()); + opxi.dissxi(tmp.0.rho(), k.0.rho_mut()); + opxi.dissxi(tmp.0.rhou(), k.0.rhou_mut()); + opxi.dissxi(tmp.0.rhov(), k.0.rhov_mut()); + opxi.dissxi(tmp.0.e(), k.0.e_mut()); - UOeta::disseta(tmp.1.rho(), k.1.rho_mut()); - UOeta::disseta(tmp.1.rhou(), k.1.rhou_mut()); - UOeta::disseta(tmp.1.rhov(), k.1.rhov_mut()); - UOeta::disseta(tmp.1.e(), k.1.e_mut()); + opeta.disseta(tmp.1.rho(), k.1.rho_mut()); + opeta.disseta(tmp.1.rhou(), k.1.rhou_mut()); + opeta.disseta(tmp.1.rhov(), k.1.rhov_mut()); + opeta.disseta(tmp.1.e(), k.1.e_mut()); } -fn fluxes( - k: (&mut Field, &mut Field), - y: &Field, - metrics: &Metrics, -) { +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(); @@ -839,17 +854,18 @@ fn vortexify( #[allow(non_snake_case)] /// Boundary conditions (SAT) fn SAT_characteristics( + (opeta, opxi): (SBPeta, SBPxi), k: &mut Field, y: &Field, - metrics: &Metrics, + metrics: &Metrics, boundaries: &BoundaryTerms, ) { // North boundary { - let hi = if SBPeta::is_h2() { - (k.ny() - 2) as Float / SBPeta::h()[0] + let hi = if opeta.is_h2() { + (k.ny() - 2) as Float / opeta.h()[0] } else { - (k.ny() - 1) as Float / SBPeta::h()[0] + (k.ny() - 1) as Float / opeta.h()[0] }; let sign = -1.0; let tau = 1.0; @@ -868,10 +884,10 @@ fn SAT_characteristics( } // South boundary { - let hi = if SBPeta::is_h2() { - (k.ny() - 2) as Float / SBPeta::h()[0] + let hi = if opeta.is_h2() { + (k.ny() - 2) as Float / opeta.h()[0] } else { - (k.ny() - 1) as Float / SBPeta::h()[0] + (k.ny() - 1) as Float / opeta.h()[0] }; let sign = 1.0; let tau = -1.0; @@ -890,10 +906,10 @@ fn SAT_characteristics( } // West Boundary { - let hi = if SBPxi::is_h2() { - (k.nx() - 2) as Float / SBPxi::h()[0] + let hi = if opxi.is_h2() { + (k.nx() - 2) as Float / opxi.h()[0] } else { - (k.nx() - 1) as Float / SBPxi::h()[0] + (k.nx() - 1) as Float / opxi.h()[0] }; let sign = 1.0; let tau = -1.0; @@ -912,10 +928,10 @@ fn SAT_characteristics( } // East Boundary { - let hi = if SBPxi::is_h2() { - (k.nx() - 2) as Float / SBPxi::h()[0] + let hi = if opxi.is_h2() { + (k.nx() - 2) as Float / opxi.h()[0] } else { - (k.nx() - 1) as Float / SBPxi::h()[0] + (k.nx() - 1) as Float / opxi.h()[0] }; let sign = -1.0; let tau = 1.0; diff --git a/sbp/src/grid.rs b/sbp/src/grid.rs index 9db994e..c99e6a3 100644 --- a/sbp/src/grid.rs +++ b/sbp/src/grid.rs @@ -8,19 +8,12 @@ pub struct Grid { } #[derive(Debug, Clone)] -pub struct Metrics -where - SBPeta: super::operators::SbpOperator, - SBPxi: super::operators::SbpOperator, -{ +pub struct Metrics { pub(crate) detj: Array2, pub(crate) detj_dxi_dx: Array2, pub(crate) detj_dxi_dy: Array2, pub(crate) detj_deta_dx: Array2, pub(crate) detj_deta_dy: Array2, - - operatoreta: std::marker::PhantomData, - operatorxi: std::marker::PhantomData, } impl Grid { @@ -43,10 +36,12 @@ impl Grid { self.y.view() } - pub fn metrics( + pub fn metrics( &self, - ) -> Result, ndarray::ShapeError> { - Metrics::new(self) + opeta: SBPeta, + opxi: SBPxi, + ) -> Result { + Metrics::new(self, opeta, opxi) } pub fn north(&self) -> (ndarray::ArrayView1, ndarray::ArrayView1) { @@ -75,23 +70,25 @@ impl Grid { } } -impl - Metrics -{ - fn new(grid: &Grid) -> Result { +impl Metrics { + fn new( + grid: &Grid, + opeta: SBPeta, + opxi: SBPxi, + ) -> 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)); - SBPxi::diffxi(x.view(), dx_dxi.view_mut()); + opxi.diffxi(x.view(), dx_dxi.view_mut()); let mut dx_deta = Array2::zeros((ny, nx)); - SBPeta::diffeta(x.view(), dx_deta.view_mut()); + opeta.diffeta(x.view(), dx_deta.view_mut()); let mut dy_dxi = Array2::zeros((ny, nx)); - SBPxi::diffxi(y.view(), dy_dxi.view_mut()); + opxi.diffxi(y.view(), dy_dxi.view_mut()); let mut dy_deta = Array2::zeros((ny, nx)); - SBPeta::diffeta(y.view(), dy_deta.view_mut()); + opeta.diffeta(y.view(), dy_deta.view_mut()); let mut detj = Array2::zeros((ny, nx)); ndarray::azip!((detj in &mut detj, @@ -122,8 +119,6 @@ impl { sys: (Field, Field), wb: WorkBuffers, grid: Grid, - metrics: Metrics, + metrics: Metrics, + op: SBP, } impl System { - pub fn new(x: Array2, y: Array2) -> Self { + pub fn new(x: Array2, y: Array2, op: SBP) -> Self { assert_eq!(x.shape(), y.shape()); let ny = x.shape()[0]; let nx = x.shape()[1]; let grid = Grid::new(x, y).unwrap(); - let metrics = grid.metrics().unwrap(); + let metrics = grid.metrics(op, op).unwrap(); Self { + op, sys: (Field::new(ny, nx), Field::new(ny, nx)), grid, metrics, @@ -115,16 +117,20 @@ impl System { } pub fn advance(&mut self, dt: Float) { - fn rhs_adaptor( - fut: &mut Field, - prev: &Field, - _time: Float, - c: &(&Grid, &Metrics), - m: &mut (Array2, Array2, Array2, Array2), - ) { + let op = self.op; + let rhs_adaptor = move |fut: &mut Field, + prev: &Field, + _time: Float, + c: &(&Grid, &Metrics), + m: &mut ( + Array2, + Array2, + Array2, + Array2, + )| { let (grid, metrics) = c; - RHS(fut, prev, grid, metrics, m); - } + RHS(op, fut, prev, grid, metrics, m); + }; let mut _time = 0.0; integrate::integrate::( rhs_adaptor, @@ -143,16 +149,20 @@ impl System { impl System { /// Using artificial dissipation with the upwind operator pub fn advance_upwind(&mut self, dt: Float) { - fn rhs_adaptor( - fut: &mut Field, - prev: &Field, - _time: Float, - c: &(&Grid, &Metrics), - m: &mut (Array2, Array2, Array2, Array2), - ) { + let op = self.op; + let rhs_adaptor = move |fut: &mut Field, + prev: &Field, + _time: Float, + c: &(&Grid, &Metrics), + m: &mut ( + Array2, + Array2, + Array2, + Array2, + )| { let (grid, metrics) = c; - RHS_upwind(fut, prev, grid, metrics, m); - } + RHS_upwind(op, fut, prev, grid, metrics, m); + }; let mut _time = 0.0; integrate::integrate::( rhs_adaptor, @@ -196,13 +206,14 @@ fn gaussian(x: Float, x0: Float, y: Float, y0: Float) -> Float { /// /// This is used both in fluxes and SAT terms fn RHS( + op: SBP, k: &mut Field, y: &Field, _grid: &Grid, - metrics: &Metrics, + metrics: &Metrics, tmp: &mut (Array2, Array2, Array2, Array2), ) { - fluxes(k, y, metrics, tmp); + fluxes(op, k, y, metrics, tmp); let boundaries = BoundaryTerms { north: Boundary::This, @@ -210,7 +221,7 @@ fn RHS( west: Boundary::This, east: Boundary::This, }; - SAT_characteristics(k, y, metrics, &boundaries); + SAT_characteristics(op, k, y, metrics, &boundaries); azip!((k in &mut k.0, &detj in &metrics.detj.broadcast((3, y.ny(), y.nx())).unwrap()) { @@ -220,14 +231,15 @@ fn RHS( #[allow(non_snake_case)] fn RHS_upwind( + op: UO, k: &mut Field, y: &Field, _grid: &Grid, - metrics: &Metrics, + metrics: &Metrics, tmp: &mut (Array2, Array2, Array2, Array2), ) { - fluxes(k, y, metrics, tmp); - dissipation(k, y, metrics, tmp); + fluxes(op, k, y, metrics, tmp); + dissipation(op, k, y, metrics, tmp); let boundaries = BoundaryTerms { north: Boundary::This, @@ -235,7 +247,7 @@ fn RHS_upwind( west: Boundary::This, east: Boundary::This, }; - SAT_characteristics(k, y, metrics, &boundaries); + SAT_characteristics(op, k, y, metrics, &boundaries); azip!((k in &mut k.0, &detj in &metrics.detj.broadcast((3, y.ny(), y.nx())).unwrap()) { @@ -244,9 +256,10 @@ fn RHS_upwind( } fn fluxes( + op: SBP, k: &mut Field, y: &Field, - metrics: &Metrics, + metrics: &Metrics, tmp: &mut (Array2, Array2, Array2, Array2), ) { // ex = hz_y @@ -256,14 +269,14 @@ fn fluxes( &hz in &y.hz()) *a = dxi_dy * hz ); - SBP::diffxi(tmp.0.view(), tmp.1.view_mut()); + op.diffxi(tmp.0.view(), tmp.1.view_mut()); ndarray::azip!((b in &mut tmp.2, &deta_dy in &metrics.detj_deta_dy, &hz in &y.hz()) *b = deta_dy * hz ); - SBP::diffeta(tmp.2.view(), tmp.3.view_mut()); + op.diffeta(tmp.2.view(), tmp.3.view_mut()); ndarray::azip!((flux in &mut k.ex_mut(), &ax in &tmp.1, &by in &tmp.3) *flux = ax + by @@ -279,7 +292,7 @@ fn fluxes( &ey in &y.ey()) *a = dxi_dx * -ey + dxi_dy * ex ); - SBP::diffxi(tmp.0.view(), tmp.1.view_mut()); + op.diffxi(tmp.0.view(), tmp.1.view_mut()); ndarray::azip!((b in &mut tmp.2, &deta_dx in &metrics.detj_deta_dx, @@ -288,7 +301,7 @@ fn fluxes( &ey in &y.ey()) *b = deta_dx * -ey + deta_dy * ex ); - SBP::diffeta(tmp.2.view(), tmp.3.view_mut()); + op.diffeta(tmp.2.view(), tmp.3.view_mut()); ndarray::azip!((flux in &mut k.hz_mut(), &ax in &tmp.1, &by in &tmp.3) *flux = ax + by @@ -302,14 +315,14 @@ fn fluxes( &hz in &y.hz()) *a = dxi_dx * -hz ); - SBP::diffxi(tmp.0.view(), tmp.1.view_mut()); + op.diffxi(tmp.0.view(), tmp.1.view_mut()); azip!((b in &mut tmp.2, &deta_dx in &metrics.detj_deta_dx, &hz in &y.hz()) *b = deta_dx * -hz ); - SBP::diffeta(tmp.2.view(), tmp.3.view_mut()); + op.diffeta(tmp.2.view(), tmp.3.view_mut()); azip!((flux in &mut k.ey_mut(), &ax in &tmp.1, &by in &tmp.3) *flux = ax + by @@ -318,9 +331,10 @@ fn fluxes( } fn dissipation( + op: UO, k: &mut Field, y: &Field, - metrics: &Metrics, + metrics: &Metrics, tmp: &mut (Array2, Array2, Array2, Array2), ) { // ex component @@ -333,7 +347,7 @@ fn dissipation( let r = Float::hypot(kx, ky); *a = ky*ky/r * ex + -kx*ky/r*ey; }); - UO::dissxi(tmp.0.view(), tmp.1.view_mut()); + op.dissxi(tmp.0.view(), tmp.1.view_mut()); ndarray::azip!((b in &mut tmp.2, &kx in &metrics.detj_deta_dx, @@ -343,7 +357,7 @@ fn dissipation( let r = Float::hypot(kx, ky); *b = ky*ky/r * ex + -kx*ky/r*ey; }); - UO::disseta(tmp.2.view(), tmp.3.view_mut()); + op.disseta(tmp.2.view(), tmp.3.view_mut()); ndarray::azip!((flux in &mut k.ex_mut(), &ax in &tmp.1, &by in &tmp.3) *flux += ax + by @@ -359,7 +373,7 @@ fn dissipation( let r = Float::hypot(kx, ky); *a = r * hz; }); - UO::dissxi(tmp.0.view(), tmp.1.view_mut()); + op.dissxi(tmp.0.view(), tmp.1.view_mut()); ndarray::azip!((b in &mut tmp.2, &kx in &metrics.detj_deta_dx, @@ -368,7 +382,7 @@ fn dissipation( let r = Float::hypot(kx, ky); *b = r * hz; }); - UO::disseta(tmp.2.view(), tmp.3.view_mut()); + op.disseta(tmp.2.view(), tmp.3.view_mut()); ndarray::azip!((flux in &mut k.hz_mut(), &ax in &tmp.1, &by in &tmp.3) *flux += ax + by @@ -385,7 +399,7 @@ fn dissipation( let r = Float::hypot(kx, ky); *a = -kx*ky/r * ex + kx*kx/r*ey; }); - UO::dissxi(tmp.0.view(), tmp.1.view_mut()); + op.dissxi(tmp.0.view(), tmp.1.view_mut()); ndarray::azip!((b in &mut tmp.2, &kx in &metrics.detj_deta_dx, @@ -395,7 +409,7 @@ fn dissipation( let r = Float::hypot(kx, ky); *b = -kx*ky/r * ex + kx*kx/r*ey; }); - UO::disseta(tmp.2.view(), tmp.3.view_mut()); + op.disseta(tmp.2.view(), tmp.3.view_mut()); ndarray::azip!((flux in &mut k.hz_mut(), &ax in &tmp.1, &by in &tmp.3) *flux += ax + by @@ -419,9 +433,10 @@ pub struct BoundaryTerms { #[allow(non_snake_case)] /// Boundary conditions (SAT) fn SAT_characteristics( + op: SBP, k: &mut Field, y: &Field, - metrics: &Metrics, + metrics: &Metrics, boundaries: &BoundaryTerms, ) { let ny = y.ny(); @@ -449,10 +464,10 @@ fn SAT_characteristics( Boundary::This => y.slice(s![.., .., 0]), }; // East boundary - let hinv = if SBP::is_h2() { - (nx - 2) as Float / SBP::h()[0] + let hinv = if op.is_h2() { + (nx - 2) as Float / op.h()[0] } else { - (nx - 1) as Float / SBP::h()[0] + (nx - 1) as Float / op.h()[0] }; for ((((mut k, v), g), &kx), &ky) in k .slice_mut(s![.., .., nx - 1]) @@ -487,10 +502,10 @@ fn SAT_characteristics( let g = match boundaries.east { Boundary::This => y.slice(s![.., .., nx - 1]), }; - let hinv = if SBP::is_h2() { - (nx - 2) as Float / SBP::h()[0] + let hinv = if op.is_h2() { + (nx - 2) as Float / op.h()[0] } else { - (nx - 1) as Float / SBP::h()[0] + (nx - 1) as Float / op.h()[0] }; for ((((mut k, v), g), &kx), &ky) in k .slice_mut(s![.., .., 0]) @@ -530,10 +545,10 @@ fn SAT_characteristics( let g = match boundaries.north { Boundary::This => y.slice(s![.., 0, ..]), }; - let hinv = if SBP::is_h2() { - (ny - 2) as Float / SBP::h()[0] + let hinv = if op.is_h2() { + (ny - 2) as Float / op.h()[0] } else { - (ny - 1) as Float / SBP::h()[0] + (ny - 1) as Float / op.h()[0] }; for ((((mut k, v), g), &kx), &ky) in k .slice_mut(s![.., ny - 1, ..]) @@ -567,10 +582,10 @@ fn SAT_characteristics( let g = match boundaries.south { Boundary::This => y.slice(s![.., ny - 1, ..]), }; - let hinv = if SBP::is_h2() { - (ny - 2) as Float / SBP::h()[0] + let hinv = if op.is_h2() { + (ny - 2) as Float / op.h()[0] } else { - (ny - 1) as Float / SBP::h()[0] + (ny - 1) as Float / op.h()[0] }; for ((((mut k, v), g), &kx), &ky) in k .slice_mut(s![.., 0, ..]) diff --git a/sbp/src/operators.rs b/sbp/src/operators.rs index 9f5d040..5ad2742 100644 --- a/sbp/src/operators.rs +++ b/sbp/src/operators.rs @@ -5,33 +5,33 @@ use crate::Float; use ndarray::{ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2}; -pub trait SbpOperator: Send + Sync { - fn diff1d(prev: ArrayView1, fut: ArrayViewMut1); - fn diffxi(prev: ArrayView2, mut fut: ArrayViewMut2) { +pub trait SbpOperator: Send + Sync + Copy { + fn diff1d(&self, prev: ArrayView1, fut: ArrayViewMut1); + fn diffxi(&self, prev: ArrayView2, mut fut: ArrayViewMut2) { assert_eq!(prev.shape(), fut.shape()); for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) { - Self::diff1d(r0, r1); + self.diff1d(r0, r1); } } - fn diffeta(prev: ArrayView2, fut: ArrayViewMut2) { - Self::diffxi(prev.reversed_axes(), fut.reversed_axes()) + fn diffeta(&self, prev: ArrayView2, fut: ArrayViewMut2) { + self.diffxi(prev.reversed_axes(), fut.reversed_axes()) } - fn h() -> &'static [Float]; - fn is_h2() -> bool { + fn h(&self) -> &'static [Float]; + fn is_h2(&self) -> bool { false } } pub trait UpwindOperator: SbpOperator { - fn diss1d(prev: ArrayView1, fut: ArrayViewMut1); - fn dissxi(prev: ArrayView2, mut fut: ArrayViewMut2) { + fn diss1d(&self, prev: ArrayView1, fut: ArrayViewMut1); + fn dissxi(&self, prev: ArrayView2, mut fut: ArrayViewMut2) { assert_eq!(prev.shape(), fut.shape()); for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) { - Self::diss1d(r0, r1); + self.diss1d(r0, r1); } } - fn disseta(prev: ArrayView2, fut: ArrayViewMut2) { - Self::dissxi(prev.reversed_axes(), fut.reversed_axes()) + fn disseta(&self, prev: ArrayView2, fut: ArrayViewMut2) { + self.dissxi(prev.reversed_axes(), fut.reversed_axes()) } } @@ -133,6 +133,7 @@ pub(crate) mod testing { } pub(crate) fn check_operator_on( + op: SBP, n: (usize, usize), f: F, dfdx: FX, @@ -148,11 +149,11 @@ pub(crate) mod testing { let x = grid_eval(n, f); y.fill(0.0); - SBP::diffxi(x.view(), y.view_mut()); + op.diffxi(x.view(), y.view_mut()); approx::assert_abs_diff_eq!(&y, &grid_eval(n, dfdx), epsilon = eps); y.fill(0.0); - SBP::diffeta(x.view(), y.view_mut()); + op.diffeta(x.view(), y.view_mut()); approx::assert_abs_diff_eq!(&y, &grid_eval(n, dfdy), epsilon = eps); } } diff --git a/sbp/src/operators/traditional4.rs b/sbp/src/operators/traditional4.rs index 4d4eebb..4dd9917 100644 --- a/sbp/src/operators/traditional4.rs +++ b/sbp/src/operators/traditional4.rs @@ -2,8 +2,8 @@ use super::SbpOperator; use crate::Float; use ndarray::{ArrayView1, ArrayViewMut1}; -#[derive(Debug)] -pub struct SBP4 {} +#[derive(Debug, Copy, Clone)] +pub struct SBP4; impl SBP4 { #[rustfmt::skip] @@ -24,7 +24,7 @@ impl SBP4 { } impl SbpOperator for SBP4 { - fn diff1d(prev: ArrayView1, fut: ArrayViewMut1) { + fn diff1d(&self, prev: ArrayView1, fut: ArrayViewMut1) { super::diff_op_1d( ndarray::arr2(Self::BLOCK).view(), ndarray::arr1(Self::DIAG).view(), @@ -35,7 +35,7 @@ impl SbpOperator for SBP4 { ) } - fn h() -> &'static [Float] { + fn h(&self) -> &'static [Float] { Self::HBLOCK } } @@ -47,15 +47,24 @@ fn test_trad4() { let nx = 20; let ny = 13; - check_operator_on::((ny, nx), |x, y| x + 2.0 * y, |_, _| 1.0, |_, _| 2.0, 1e-4); - check_operator_on::( + check_operator_on( + SBP4, + (ny, nx), + |x, y| x + 2.0 * y, + |_, _| 1.0, + |_, _| 2.0, + 1e-4, + ); + check_operator_on( + SBP4, (ny, nx), |x, y| x * x + 2.0 * x * y + 3.0 * y * y, |x, y| 2.0 * x + 2.0 * y, |x, y| 2.0 * x + 6.0 * y, 1e-3, ); - check_operator_on::( + check_operator_on( + SBP4, (ny, nx), |x, y| x.powi(3) + 2.0 * x.powi(2) * y + 3.0 * x * y.powi(2) + 4.0 * y.powi(3), |x, y| 3.0 * x.powi(2) + 4.0 * x * y + 3.0 * y.powi(2), diff --git a/sbp/src/operators/traditional8.rs b/sbp/src/operators/traditional8.rs index d5f0c6f..8fe19b6 100644 --- a/sbp/src/operators/traditional8.rs +++ b/sbp/src/operators/traditional8.rs @@ -2,8 +2,8 @@ use super::SbpOperator; use crate::Float; use ndarray::{ArrayView1, ArrayViewMut1}; -#[derive(Debug)] -pub struct SBP8 {} +#[derive(Debug, Copy, Clone)] +pub struct SBP8; impl SBP8 { #[rustfmt::skip] @@ -28,7 +28,7 @@ impl SBP8 { } impl SbpOperator for SBP8 { - fn diff1d(prev: ArrayView1, fut: ArrayViewMut1) { + fn diff1d(&self, prev: ArrayView1, fut: ArrayViewMut1) { super::diff_op_1d( ndarray::arr2(Self::BLOCK).view(), ndarray::arr1(Self::DIAG).view(), @@ -39,7 +39,7 @@ impl SbpOperator for SBP8 { ) } - fn h() -> &'static [Float] { + fn h(&self) -> &'static [Float] { Self::HBLOCK } } @@ -51,7 +51,8 @@ fn test_trad8() { let ny = 16; // Order one polynomial - check_operator_on::( + check_operator_on( + SBP8, (ny, nx), |x, y| x + 2.0 * y, |_x, _y| 1.0, @@ -60,31 +61,35 @@ fn test_trad8() { ); // Order two polynomial - check_operator_on::( + check_operator_on( + SBP8, (ny, nx), |x, y| x * x + 0.5 * y * y, |x, _y| 2.0 * x, |_x, y| y, 1e-4, ); - check_operator_on::((ny, nx), |x, y| x * y, |_x, y| y, |x, _y| x, 1e-4); + check_operator_on(SBP8, (ny, nx), |x, y| x * y, |_x, y| y, |x, _y| x, 1e-4); // Order three polynomials - check_operator_on::( + check_operator_on( + SBP8, (ny, nx), |x, y| x * x * x + y * y * y / 6.0, |x, _y| 3.0 * x * x, |_x, y| y * y / 2.0, 1e-4, ); - check_operator_on::( + check_operator_on( + SBP8, (ny, nx), |x, y| x * x * y + x * y * y / 2.0, |x, y| 2.0 * x * y + y * y / 2.0, |x, y| x * x + x * y, 1e-4, ); - check_operator_on::( + check_operator_on( + SBP8, (ny, nx), |x, y| x.powi(3) + 2.0 * x.powi(2) * y + 3.0 * x * y.powi(2) + 4.0 * y.powi(3), |x, y| 3.0 * x.powi(2) + 4.0 * x * y + 3.0 * y.powi(2), @@ -93,7 +98,8 @@ fn test_trad8() { ); // Order four polynomials - check_operator_on::( + check_operator_on( + SBP8, (ny, nx), |x, y| x.powi(4) + x.powi(3) * y + x.powi(2) * y.powi(2) + x * y.powi(3) + y.powi(4), |x, y| 4.0 * x.powi(3) + 3.0 * x.powi(2) * y + 2.0 * x * y.powi(2) + y.powi(3), diff --git a/sbp/src/operators/upwind4.rs b/sbp/src/operators/upwind4.rs index 8efa7c6..834fbeb 100644 --- a/sbp/src/operators/upwind4.rs +++ b/sbp/src/operators/upwind4.rs @@ -2,8 +2,8 @@ use super::{SbpOperator, UpwindOperator}; use crate::Float; use ndarray::{ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2, Axis}; -#[derive(Debug)] -pub struct Upwind4 {} +#[derive(Debug, Copy, Clone)] +pub struct Upwind4; /// Simdtype used in diff_simd_col and diff_simd_row #[cfg(feature = "f32")] @@ -276,7 +276,7 @@ impl Upwind4 { } impl SbpOperator for Upwind4 { - fn diff1d(prev: ArrayView1, fut: ArrayViewMut1) { + fn diff1d(&self, prev: ArrayView1, fut: ArrayViewMut1) { super::diff_op_1d( ndarray::arr2(Self::BLOCK).view(), ndarray::arr1(Self::DIAG).view(), @@ -286,7 +286,7 @@ impl SbpOperator for Upwind4 { fut, ) } - fn diffxi(prev: ArrayView2, mut fut: ArrayViewMut2) { + fn diffxi(&self, prev: ArrayView2, mut fut: ArrayViewMut2) { assert_eq!(prev.shape(), fut.shape()); assert!(prev.shape()[1] >= 2 * Self::BLOCK.len()); @@ -300,14 +300,14 @@ impl SbpOperator for Upwind4 { ([_, _], [_, _]) => { // Fallback, work row by row for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) { - Self::diff1d(r0, r1); + Self.diff1d(r0, r1); } } _ => unreachable!("Should only be two elements in the strides vectors"), } } - fn h() -> &'static [Float] { + fn h(&self) -> &'static [Float] { Self::HBLOCK } } @@ -326,14 +326,14 @@ fn upwind4_test() { target[i] = 1.0; } res.fill(0.0); - Upwind4::diff1d(source.view(), res.view_mut()); + Upwind4.diff1d(source.view(), res.view_mut()); approx::assert_abs_diff_eq!(&res, &target, epsilon = 1e-4); { let source = source.to_owned().insert_axis(ndarray::Axis(0)); let mut res = res.to_owned().insert_axis(ndarray::Axis(0)); let target = target.to_owned().insert_axis(ndarray::Axis(0)); res.fill(0.0); - Upwind4::diffxi(source.view(), res.view_mut()); + Upwind4.diffxi(source.view(), res.view_mut()); approx::assert_abs_diff_eq!(&res, &target, epsilon = 1e-2); } @@ -342,7 +342,7 @@ fn upwind4_test() { let target = Array2::from_shape_fn((nx, 8), |(i, _)| target[i]); let mut res = Array2::zeros((nx, 8)); res.fill(0.0); - Upwind4::diffeta(source.view(), res.view_mut()); + Upwind4.diffeta(source.view(), res.view_mut()); approx::assert_abs_diff_eq!(&res.to_owned(), &target.to_owned(), epsilon = 1e-2); } @@ -352,14 +352,14 @@ fn upwind4_test() { target[i] = 2.0 * x; } res.fill(0.0); - Upwind4::diff1d(source.view(), res.view_mut()); + Upwind4.diff1d(source.view(), res.view_mut()); approx::assert_abs_diff_eq!(&res, &target, epsilon = 1e-4); { let source = source.to_owned().insert_axis(ndarray::Axis(0)); let mut res = res.to_owned().insert_axis(ndarray::Axis(0)); let target = target.to_owned().insert_axis(ndarray::Axis(0)); res.fill(0.0); - Upwind4::diffxi(source.view(), res.view_mut()); + Upwind4.diffxi(source.view(), res.view_mut()); approx::assert_abs_diff_eq!(&res, &target, epsilon = 1e-2); } @@ -368,7 +368,7 @@ fn upwind4_test() { let target = Array2::from_shape_fn((nx, 8), |(i, _)| target[i]); let mut res = Array2::zeros((nx, 8)); res.fill(0.0); - Upwind4::diffeta(source.view(), res.view_mut()); + Upwind4.diffeta(source.view(), res.view_mut()); approx::assert_abs_diff_eq!(&res.to_owned(), &target.to_owned(), epsilon = 1e-2); } @@ -378,7 +378,7 @@ fn upwind4_test() { target[i] = 3.0 * x * x; } res.fill(0.0); - Upwind4::diff1d(source.view(), res.view_mut()); + Upwind4.diff1d(source.view(), res.view_mut()); approx::assert_abs_diff_eq!(&res, &target, epsilon = 1e-2); { @@ -386,7 +386,7 @@ fn upwind4_test() { let mut res = res.to_owned().insert_axis(ndarray::Axis(0)); let target = target.to_owned().insert_axis(ndarray::Axis(0)); res.fill(0.0); - Upwind4::diffxi(source.view(), res.view_mut()); + Upwind4.diffxi(source.view(), res.view_mut()); approx::assert_abs_diff_eq!(&res, &target, epsilon = 1e-2); } @@ -395,13 +395,13 @@ fn upwind4_test() { let target = Array2::from_shape_fn((nx, 8), |(i, _)| target[i]); let mut res = Array2::zeros((nx, 8)); res.fill(0.0); - Upwind4::diffeta(source.view(), res.view_mut()); + Upwind4.diffeta(source.view(), res.view_mut()); approx::assert_abs_diff_eq!(&res.to_owned(), &target.to_owned(), epsilon = 1e-2); } } impl UpwindOperator for Upwind4 { - fn diss1d(prev: ArrayView1, fut: ArrayViewMut1) { + fn diss1d(&self, prev: ArrayView1, fut: ArrayViewMut1) { super::diff_op_1d( ndarray::arr2(Self::DISS_BLOCK).view(), ndarray::arr1(Self::DISS_DIAG).view(), @@ -411,7 +411,7 @@ impl UpwindOperator for Upwind4 { fut, ) } - fn dissxi(prev: ArrayView2, mut fut: ArrayViewMut2) { + fn dissxi(&self, prev: ArrayView2, mut fut: ArrayViewMut2) { assert_eq!(prev.shape(), fut.shape()); assert!(prev.shape()[1] >= 2 * Self::BLOCK.len()); @@ -425,7 +425,7 @@ impl UpwindOperator for Upwind4 { ([_, _], [_, _]) => { // Fallback, work row by row for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) { - Self::diss1d(r0, r1); + Self.diss1d(r0, r1); } } _ => unreachable!("Should only be two elements in the strides vectors"), @@ -440,28 +440,32 @@ fn upwind4_test2() { let nx = 32; let ny = 16; - check_operator_on::( + check_operator_on( + Upwind4, (ny, nx), |x, y| x + 2.0 * y, |_, _| 1.0, |_, _| 2.0, 1e-4, ); - check_operator_on::( + check_operator_on( + Upwind4, (ny, nx), |x, y| x * x + 2.0 * x * y + 3.0 * y * y, |x, y| 2.0 * x + 2.0 * y, |x, y| 2.0 * x + 6.0 * y, 1e-3, ); - check_operator_on::( + check_operator_on( + Upwind4, (ny, nx), |x, y| x.powi(3) + 2.0 * x.powi(2) * y + 3.0 * x * y.powi(2) + 4.0 * y.powi(3), |x, y| 3.0 * x.powi(2) + 4.0 * x * y + 3.0 * y.powi(2), |x, y| 2.0 * x.powi(2) + 6.0 * x * y + 12.0 * y.powi(2), 1e-1, ); - check_operator_on::( + check_operator_on( + Upwind4, (32, 32), |x, y| x.powi(3) + 2.0 * x.powi(2) * y + 3.0 * x * y.powi(2) + 4.0 * y.powi(3), |x, y| 3.0 * x.powi(2) + 4.0 * x * y + 3.0 * y.powi(2), diff --git a/sbp/src/operators/upwind4h2.rs b/sbp/src/operators/upwind4h2.rs index b6857a3..d9304ba 100644 --- a/sbp/src/operators/upwind4h2.rs +++ b/sbp/src/operators/upwind4h2.rs @@ -2,8 +2,8 @@ use super::{SbpOperator, UpwindOperator}; use crate::Float; use ndarray::{ArrayView1, ArrayViewMut1}; -#[derive(Debug)] -pub struct Upwind4h2 {} +#[derive(Debug, Copy, Clone)] +pub struct Upwind4h2; impl Upwind4h2 { #[rustfmt::skip] @@ -37,7 +37,7 @@ impl Upwind4h2 { } impl SbpOperator for Upwind4h2 { - fn diff1d(prev: ArrayView1, fut: ArrayViewMut1) { + fn diff1d(&self, prev: ArrayView1, fut: ArrayViewMut1) { super::diff_op_1d( ndarray::arr2(Self::BLOCK).view(), ndarray::arr1(Self::DIAG).view(), @@ -48,10 +48,10 @@ impl SbpOperator for Upwind4h2 { ) } - fn h() -> &'static [Float] { + fn h(&self) -> &'static [Float] { Self::HBLOCK } - fn is_h2() -> bool { + fn is_h2(&self) -> bool { true } } @@ -64,25 +64,25 @@ fn upwind4h2_test() { let mut res = ndarray::Array1::zeros(nx); - Upwind4h2::diff1d(x.view(), res.view_mut()); + Upwind4h2.diff1d(x.view(), res.view_mut()); let ans = &x * 0.0 + 1.0; approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-4); res.fill(0.0); let y = &x * &x / 2.0; - Upwind4h2::diff1d(y.view(), res.view_mut()); + Upwind4h2.diff1d(y.view(), res.view_mut()); let ans = &x; approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-4); res.fill(0.0); let y = &x * &x * &x / 3.0; - Upwind4h2::diff1d(y.view(), res.view_mut()); + Upwind4h2.diff1d(y.view(), res.view_mut()); let ans = &x * &x; approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-2); } impl UpwindOperator for Upwind4h2 { - fn diss1d(prev: ArrayView1, fut: ArrayViewMut1) { + fn diss1d(&self, prev: ArrayView1, fut: ArrayViewMut1) { super::diff_op_1d( ndarray::arr2(Self::DISS_BLOCK).view(), ndarray::arr1(Self::DISS_DIAG).view(), diff --git a/sbp/src/operators/upwind9.rs b/sbp/src/operators/upwind9.rs index 18009f9..9c11143 100644 --- a/sbp/src/operators/upwind9.rs +++ b/sbp/src/operators/upwind9.rs @@ -2,8 +2,8 @@ use super::{SbpOperator, UpwindOperator}; use crate::Float; use ndarray::{ArrayView1, ArrayViewMut1}; -#[derive(Debug)] -pub struct Upwind9 {} +#[derive(Debug, Copy, Clone)] +pub struct Upwind9; impl Upwind9 { #[rustfmt::skip] @@ -45,7 +45,7 @@ impl Upwind9 { } impl SbpOperator for Upwind9 { - fn diff1d(prev: ArrayView1, fut: ArrayViewMut1) { + fn diff1d(&self, prev: ArrayView1, fut: ArrayViewMut1) { super::diff_op_1d( ndarray::arr2(Self::BLOCK).view(), ndarray::arr1(Self::DIAG).view(), @@ -56,13 +56,13 @@ impl SbpOperator for Upwind9 { ) } - fn h() -> &'static [Float] { + fn h(&self) -> &'static [Float] { Self::HBLOCK } } impl UpwindOperator for Upwind9 { - fn diss1d(prev: ArrayView1, fut: ArrayViewMut1) { + fn diss1d(&self, prev: ArrayView1, fut: ArrayViewMut1) { super::diff_op_1d( ndarray::arr2(Self::DISS_BLOCK).view(), ndarray::arr1(Self::DISS_DIAG).view(), @@ -81,7 +81,8 @@ fn test_upwind9() { let ny = 16; // Order one polynomial - check_operator_on::( + check_operator_on( + Upwind9, (ny, nx), |x, y| x + 2.0 * y, |_x, _y| 1.0, @@ -90,31 +91,35 @@ fn test_upwind9() { ); // Order two polynomial - check_operator_on::( + check_operator_on( + Upwind9, (ny, nx), |x, y| x * x + 0.5 * y * y, |x, _y| 2.0 * x, |_x, y| y, 1e-4, ); - check_operator_on::((ny, nx), |x, y| x * y, |_x, y| y, |x, _y| x, 1e-4); + check_operator_on(Upwind9, (ny, nx), |x, y| x * y, |_x, y| y, |x, _y| x, 1e-4); // Order three polynomials - check_operator_on::( + check_operator_on( + Upwind9, (ny, nx), |x, y| x * x * x + y * y * y / 6.0, |x, _y| 3.0 * x * x, |_x, y| y * y / 2.0, 1e-4, ); - check_operator_on::( + check_operator_on( + Upwind9, (ny, nx), |x, y| x * x * y + x * y * y / 2.0, |x, y| 2.0 * x * y + y * y / 2.0, |x, y| x * x + x * y, 1e-4, ); - check_operator_on::( + check_operator_on( + Upwind9, (ny, nx), |x, y| x.powi(3) + 2.0 * x.powi(2) * y + 3.0 * x * y.powi(2) + 4.0 * y.powi(3), |x, y| 3.0 * x.powi(2) + 4.0 * x * y + 3.0 * y.powi(2), @@ -123,7 +128,8 @@ fn test_upwind9() { ); // Order four polynomials - check_operator_on::( + check_operator_on( + Upwind9, (ny, nx), |x, y| x.powi(4) + x.powi(3) * y + x.powi(2) * y.powi(2) + x * y.powi(3) + y.powi(4), |x, y| 4.0 * x.powi(3) + 3.0 * x.powi(2) * y + 2.0 * x * y.powi(2) + y.powi(3), diff --git a/sbp/src/operators/upwind9h2.rs b/sbp/src/operators/upwind9h2.rs index 06defed..b65296d 100644 --- a/sbp/src/operators/upwind9h2.rs +++ b/sbp/src/operators/upwind9h2.rs @@ -2,8 +2,8 @@ use super::{SbpOperator, UpwindOperator}; use crate::Float; use ndarray::{ArrayView1, ArrayViewMut1}; -#[derive(Debug)] -pub struct Upwind9h2 {} +#[derive(Debug, Copy, Clone)] +pub struct Upwind9h2; impl Upwind9h2 { #[rustfmt::skip] @@ -45,7 +45,7 @@ impl Upwind9h2 { } impl SbpOperator for Upwind9h2 { - fn diff1d(prev: ArrayView1, fut: ArrayViewMut1) { + fn diff1d(&self, prev: ArrayView1, fut: ArrayViewMut1) { super::diff_op_1d( ndarray::arr2(Self::BLOCK).view(), ndarray::arr1(Self::DIAG).view(), @@ -56,10 +56,10 @@ impl SbpOperator for Upwind9h2 { ) } - fn h() -> &'static [Float] { + fn h(&self) -> &'static [Float] { Self::HBLOCK } - fn is_h2() -> bool { + fn is_h2(&self) -> bool { true } } @@ -72,25 +72,25 @@ fn upwind9h2_test() { let mut res = ndarray::Array1::zeros(nx); - Upwind9h2::diff1d(x.view(), res.view_mut()); + Upwind9h2.diff1d(x.view(), res.view_mut()); let ans = &x * 0.0 + 1.0; approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-4); res.fill(0.0); let y = &x * &x / 2.0; - Upwind9h2::diff1d(y.view(), res.view_mut()); + Upwind9h2.diff1d(y.view(), res.view_mut()); let ans = &x; approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-4); res.fill(0.0); let y = &x * &x * &x / 3.0; - Upwind9h2::diff1d(y.view(), res.view_mut()); + Upwind9h2.diff1d(y.view(), res.view_mut()); let ans = &x * &x; approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-2); } impl UpwindOperator for Upwind9h2 { - fn diss1d(prev: ArrayView1, fut: ArrayViewMut1) { + fn diss1d(&self, prev: ArrayView1, fut: ArrayViewMut1) { super::diff_op_1d( ndarray::arr2(Self::DISS_BLOCK).view(), ndarray::arr1(Self::DISS_DIAG).view(), diff --git a/webfront/lib.rs b/webfront/lib.rs index 13f3a9a..aa44711 100644 --- a/webfront/lib.rs +++ b/webfront/lib.rs @@ -20,7 +20,7 @@ impl MaxwellUniverse { pub fn new(height: usize, width: usize, x: &[f32], y: &[f32]) -> Self { let x = ndarray::Array2::from_shape_vec((height, width), x.to_vec()).unwrap(); let y = ndarray::Array2::from_shape_vec((height, width), y.to_vec()).unwrap(); - Self(maxwell::System::new(x, y)) + Self(maxwell::System::new(x, y, operators::Upwind4)) } pub fn init(&mut self, x0: f32, y0: f32) { @@ -53,7 +53,7 @@ pub struct EulerUniverse(euler::System); impl EulerUniverse { pub fn new(x: ndarray::Array2, y: ndarray::Array2) -> Self { - Self(euler::System::new(x, y)) + Self(euler::System::new(x, y, operators::Upwind4)) } } @@ -63,7 +63,7 @@ impl EulerUniverse { pub fn new_with_slice(height: usize, width: usize, x: &[f32], y: &[f32]) -> Self { let x = ndarray::Array2::from_shape_vec((height, width), x.to_vec()).unwrap(); let y = ndarray::Array2::from_shape_vec((height, width), y.to_vec()).unwrap(); - Self(euler::System::new(x, y)) + Self(euler::System::new(x, y, operators::Upwind4)) } pub fn init(&mut self, x0: f32, y0: f32) {