From 90f72847fbbb8688334d9bc7d2840fafc306e211 Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Sun, 12 Apr 2020 19:27:18 +0200 Subject: [PATCH] add h2 operator --- sbp/src/euler.rs | 42 ++++++++++++-- sbp/src/lib.rs | 2 + sbp/src/maxwell.rs | 24 ++++++-- sbp/src/operators.rs | 14 ++++- sbp/src/operators/traditional4.rs | 1 + sbp/src/operators/traditional8.rs | 1 + sbp/src/operators/upwind4.rs | 2 + sbp/src/operators/upwind4h2.rs | 95 +++++++++++++++++++++++++++++++ sbp/src/operators/upwind9.rs | 2 + sbp/src/utils.rs | 42 ++++++++++++-- 10 files changed, 208 insertions(+), 17 deletions(-) create mode 100644 sbp/src/operators/upwind4h2.rs diff --git a/sbp/src/euler.rs b/sbp/src/euler.rs index 90b9ea3..d054f91 100644 --- a/sbp/src/euler.rs +++ b/sbp/src/euler.rs @@ -283,12 +283,26 @@ impl Field { .map(move |x| x * factor) }; - let hxiterator = itermaker(self.nx(), 1.0 / (self.nx() - 1) as Float); + let hxiterator = itermaker( + self.nx(), + if SBP::is_h2() { + 1.0 / (self.nx() - 2) as Float + } else { + 1.0 / (self.nx() - 1) as Float + }, + ); // Repeating to get the form // [[hx0, hx1, ..., hxn], [hx0, hx1, ..., hxn], ..., [hx0, hx1, ..., hxn]] let hxiterator = hxiterator.into_iter().cycle().take(self.nx() * self.ny()); - let hyiterator = itermaker(self.ny(), 1.0 / (self.ny() - 1) as Float); + let hyiterator = itermaker( + self.ny(), + 1.0 / if SBP::is_h2() { + (self.ny() - 2) as Float + } else { + (self.ny() - 1) as Float + }, + ); // Repeating to get the form // [[hy0, hy0, ..., hy0], [hy1, hy1, ..., hy1], ..., [hym, hym, ..., hym]] let hyiterator = hyiterator @@ -801,7 +815,11 @@ fn SAT_characteristics( ) { // North boundary { - let hi = (k.ny() - 1) as Float / SBP::h()[0]; + let hi = if SBP::is_h2() { + (k.ny() - 2) as Float / SBP::h()[0] + } else { + (k.ny() - 1) as Float / SBP::h()[0] + }; let sign = -1.0; let tau = 1.0; let slice = s![y.ny() - 1, ..]; @@ -819,7 +837,11 @@ fn SAT_characteristics( } // South boundary { - let hi = (k.ny() - 1) as Float / SBP::h()[0]; + let hi = if SBP::is_h2() { + (k.ny() - 2) as Float / SBP::h()[0] + } else { + (k.ny() - 1) as Float / SBP::h()[0] + }; let sign = 1.0; let tau = -1.0; let slice = s![0, ..]; @@ -837,7 +859,11 @@ fn SAT_characteristics( } // West Boundary { - let hi = (k.nx() - 1) as Float / SBP::h()[0]; + let hi = if SBP::is_h2() { + (k.nx() - 2) as Float / SBP::h()[0] + } else { + (k.nx() - 1) as Float / SBP::h()[0] + }; let sign = 1.0; let tau = -1.0; let slice = s![.., 0]; @@ -855,7 +881,11 @@ fn SAT_characteristics( } // East Boundary { - let hi = (k.nx() - 1) as Float / SBP::h()[0]; + let hi = if SBP::is_h2() { + (k.nx() - 2) as Float / SBP::h()[0] + } else { + (k.nx() - 1) as Float / SBP::h()[0] + }; let sign = -1.0; let tau = 1.0; let slice = s![.., y.nx() - 1]; diff --git a/sbp/src/lib.rs b/sbp/src/lib.rs index 9b361fd..2f37194 100644 --- a/sbp/src/lib.rs +++ b/sbp/src/lib.rs @@ -1,3 +1,5 @@ +#![feature(str_strip)] + #[cfg(feature = "f32")] pub type Float = f32; #[cfg(not(feature = "f32"))] diff --git a/sbp/src/maxwell.rs b/sbp/src/maxwell.rs index e253944..0dacf35 100644 --- a/sbp/src/maxwell.rs +++ b/sbp/src/maxwell.rs @@ -449,7 +449,11 @@ fn SAT_characteristics( Boundary::This => y.slice(s![.., .., 0]), }; // East boundary - let hinv = 1.0 / (SBP::h()[0] / (nx - 1) as Float); + let hinv = if SBP::is_h2() { + (nx - 2) as Float / SBP::h()[0] + } else { + (nx - 1) as Float / SBP::h()[0] + }; for ((((mut k, v), g), &kx), &ky) in k .slice_mut(s![.., .., nx - 1]) .gencolumns_mut() @@ -483,7 +487,11 @@ fn SAT_characteristics( let g = match boundaries.east { Boundary::This => y.slice(s![.., .., nx - 1]), }; - let hinv = 1.0 / (SBP::h()[0] / (nx - 1) as Float); + let hinv = if SBP::is_h2() { + (nx - 2) as Float / SBP::h()[0] + } else { + (nx - 1) as Float / SBP::h()[0] + }; for ((((mut k, v), g), &kx), &ky) in k .slice_mut(s![.., .., 0]) .gencolumns_mut() @@ -522,7 +530,11 @@ fn SAT_characteristics( let g = match boundaries.north { Boundary::This => y.slice(s![.., 0, ..]), }; - let hinv = 1.0 / (SBP::h()[0] / (ny - 1) as Float); + let hinv = if SBP::is_h2() { + (ny - 2) as Float / SBP::h()[0] + } else { + (ny - 1) as Float / SBP::h()[0] + }; for ((((mut k, v), g), &kx), &ky) in k .slice_mut(s![.., ny - 1, ..]) .gencolumns_mut() @@ -555,7 +567,11 @@ fn SAT_characteristics( let g = match boundaries.south { Boundary::This => y.slice(s![.., ny - 1, ..]), }; - let hinv = 1.0 / (SBP::h()[0] / (ny - 1) as Float); + let hinv = if SBP::is_h2() { + (ny - 2) as Float / SBP::h()[0] + } else { + (ny - 1) as Float / SBP::h()[0] + }; for ((((mut k, v), g), &kx), &ky) in k .slice_mut(s![.., 0, ..]) .gencolumns_mut() diff --git a/sbp/src/operators.rs b/sbp/src/operators.rs index f1b1b48..8264a42 100644 --- a/sbp/src/operators.rs +++ b/sbp/src/operators.rs @@ -17,6 +17,9 @@ pub trait SbpOperator: Send + Sync { Self::diffxi(prev.reversed_axes(), fut.reversed_axes()) } fn h() -> &'static [Float]; + fn is_h2() -> bool { + false + } } pub trait UpwindOperator: SbpOperator { @@ -42,6 +45,7 @@ pub(crate) fn diff_op_1d( block: ndarray::ArrayView2, diag: ndarray::ArrayView1, symmetric: bool, + is_h2: bool, prev: ArrayView1, mut fut: ArrayViewMut1, ) { @@ -49,7 +53,11 @@ pub(crate) fn diff_op_1d( let nx = prev.shape()[0]; assert!(nx >= 2 * block.len_of(ndarray::Axis(0))); - let dx = 1.0 / (nx - 1) as Float; + let dx = if is_h2 { + 1.0 / (nx - 2) as Float + } else { + 1.0 / (nx - 1) as Float + }; let idx = 1.0 / dx; let first_elems = prev.slice(::ndarray::s!(..block.len_of(::ndarray::Axis(1)))); @@ -91,6 +99,10 @@ mod upwind4; pub use upwind4::Upwind4; mod upwind9; pub use upwind9::Upwind9; + +mod upwind4h2; +pub use upwind4h2::Upwind4h2; + mod traditional4; pub use traditional4::SBP4; mod traditional8; diff --git a/sbp/src/operators/traditional4.rs b/sbp/src/operators/traditional4.rs index 1ea42c3..4d4eebb 100644 --- a/sbp/src/operators/traditional4.rs +++ b/sbp/src/operators/traditional4.rs @@ -29,6 +29,7 @@ impl SbpOperator for SBP4 { ndarray::arr2(Self::BLOCK).view(), ndarray::arr1(Self::DIAG).view(), false, + false, prev, fut, ) diff --git a/sbp/src/operators/traditional8.rs b/sbp/src/operators/traditional8.rs index 63fd017..d5f0c6f 100644 --- a/sbp/src/operators/traditional8.rs +++ b/sbp/src/operators/traditional8.rs @@ -33,6 +33,7 @@ impl SbpOperator for SBP8 { ndarray::arr2(Self::BLOCK).view(), ndarray::arr1(Self::DIAG).view(), false, + false, prev, fut, ) diff --git a/sbp/src/operators/upwind4.rs b/sbp/src/operators/upwind4.rs index e0a7bb8..8efa7c6 100644 --- a/sbp/src/operators/upwind4.rs +++ b/sbp/src/operators/upwind4.rs @@ -281,6 +281,7 @@ impl SbpOperator for Upwind4 { ndarray::arr2(Self::BLOCK).view(), ndarray::arr1(Self::DIAG).view(), false, + false, prev, fut, ) @@ -405,6 +406,7 @@ impl UpwindOperator for Upwind4 { ndarray::arr2(Self::DISS_BLOCK).view(), ndarray::arr1(Self::DISS_DIAG).view(), true, + false, prev, fut, ) diff --git a/sbp/src/operators/upwind4h2.rs b/sbp/src/operators/upwind4h2.rs new file mode 100644 index 0000000..b6857a3 --- /dev/null +++ b/sbp/src/operators/upwind4h2.rs @@ -0,0 +1,95 @@ +use super::{SbpOperator, UpwindOperator}; +use crate::Float; +use ndarray::{ArrayView1, ArrayViewMut1}; + +#[derive(Debug)] +pub struct Upwind4h2 {} + +impl Upwind4h2 { + #[rustfmt::skip] + const HBLOCK: &'static [Float] = &[ + 0.91e2/0.720e3, 0.325e3/0.384e3, 0.595e3/0.576e3, 0.1909e4/0.1920e4, + ]; + #[rustfmt::skip] + const DIAG: &'static [Float] = &[ + -1.43229166666667e-02, 1.40625000000000e-01, -7.38281250000000e-01, 0.00000000000000e+00, 7.38281250000000e-01, -1.40625000000000e-01, 1.43229166666667e-02 + ]; + #[rustfmt::skip] + const BLOCK: &'static [[Float; 7]] = &[ + [-3.95604395604396e+00, 5.41758241758242e+00, -1.94505494505495e+00, 4.83516483516484e-01, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00], + [-8.09025641025641e-01, 0.00000000000000e+00, 1.03948717948718e+00, -2.47384615384615e-01, 1.69230769230769e-02, 0.00000000000000e+00, 0.00000000000000e+00], + [2.37983193277311e-01, -8.51680672268908e-01, 0.00000000000000e+00, 7.35966386554622e-01, -1.36134453781513e-01, 1.38655462184874e-02, 0.00000000000000e+00], + [-6.14632442814737e-02, 2.10581456259822e-01, -7.64623712240265e-01, 0.00000000000000e+00, 7.42535358826611e-01, -1.41435306443164e-01, 1.44054478784704e-02], + ]; + + #[rustfmt::skip] + const DISS_BLOCK: &'static [[Float; 7]; 4] = &[ + [-2.76989342315976e-01, 5.19355016842454e-01, -3.46236677894969e-01, 1.03871003368491e-01, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00], + [7.75570158484731e-02, -1.62342481638964e-01, 1.47715500579822e-01, -7.98531117124082e-02, 1.69230769230769e-02, 0.00000000000000e+00, 0.00000000000000e+00], + [-4.23630758836198e-02, 1.21027405937249e-01, -1.91609307039399e-01, 1.82272708078206e-01, -8.31932773109244e-02, 1.38655462184874e-02, 0.00000000000000e+00], + [1.32037874021759e-02, -6.79734450144910e-02, 1.89370108794365e-01, -2.78654929966754e-01, 2.16081718177056e-01, -8.64326872708224e-02, 1.44054478784704e-02], + ]; + + #[rustfmt::skip] + const DISS_DIAG: &'static [Float; 7] = &[ + 1.43229166666667e-02, -8.59375000000000e-02, 2.14843750000000e-01, -2.86458333333333e-01, 2.14843750000000e-01, -8.59375000000000e-02, 1.43229166666667e-02, + ]; +} + +impl SbpOperator for Upwind4h2 { + fn diff1d(prev: ArrayView1, fut: ArrayViewMut1) { + super::diff_op_1d( + ndarray::arr2(Self::BLOCK).view(), + ndarray::arr1(Self::DIAG).view(), + false, + true, + prev, + fut, + ) + } + + fn h() -> &'static [Float] { + Self::HBLOCK + } + fn is_h2() -> bool { + true + } +} + +#[test] +fn upwind4h2_test() { + let nx = 20; + + let x = crate::utils::h2linspace(0.0, 1.0, nx); + + let mut res = ndarray::Array1::zeros(nx); + + 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()); + 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()); + let ans = &x * &x; + approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-2); +} + +impl UpwindOperator for Upwind4h2 { + fn diss1d(prev: ArrayView1, fut: ArrayViewMut1) { + super::diff_op_1d( + ndarray::arr2(Self::DISS_BLOCK).view(), + ndarray::arr1(Self::DISS_DIAG).view(), + true, + true, + prev, + fut, + ) + } +} diff --git a/sbp/src/operators/upwind9.rs b/sbp/src/operators/upwind9.rs index 522083b..18009f9 100644 --- a/sbp/src/operators/upwind9.rs +++ b/sbp/src/operators/upwind9.rs @@ -50,6 +50,7 @@ impl SbpOperator for Upwind9 { ndarray::arr2(Self::BLOCK).view(), ndarray::arr1(Self::DIAG).view(), false, + false, prev, fut, ) @@ -66,6 +67,7 @@ impl UpwindOperator for Upwind9 { ndarray::arr2(Self::DISS_BLOCK).view(), ndarray::arr1(Self::DISS_DIAG).view(), true, + false, prev, fut, ) diff --git a/sbp/src/utils.rs b/sbp/src/utils.rs index dd84fa4..df3c921 100644 --- a/sbp/src/utils.rs +++ b/sbp/src/utils.rs @@ -19,6 +19,7 @@ pub struct ExtendedGrid { /// x: [x0, x1, ..., xn] /// which results in x being broadcasted to nx/ny size /// x: linspace:start:end:num +/// x: linspace:h2:start:end:num /// where x will be from start to end inclusive, with num steps /// x: [[x00, x01, .., x0n], [x10, x11, .., x1n], ... [xm0, xm1, ..., xmn]] /// which is the full grid x @@ -49,11 +50,15 @@ pub fn json_to_grids(json: JsonValue) -> Result, String> { let to_array_form = |mut x: JsonValue| { if let Some(s) = x.take_string() { - if s.starts_with("linspace") { + if let Some(s) = s.strip_prefix("linspace:") { + let (s, h2) = if let Some(s) = s.strip_prefix("h2:") { + (s, true) + } else { + (s, false) + }; + // linspace:start:stop:steps let mut iter = s.split(':'); - let name = iter.next().unwrap(); - assert_eq!(name, "linspace"); let start = iter.next(); let start: Float = match start { @@ -73,9 +78,11 @@ pub fn json_to_grids(json: JsonValue) -> Result, String> { if iter.next().is_some() { return Err(format!("linspace: contained more than expected")); } - Ok(ArrayForm::Array1(ndarray::Array::linspace( - start, end, steps, - ))) + Ok(ArrayForm::Array1(if h2 { + h2linspace(start, end, steps) + } else { + ndarray::Array::linspace(start, end, steps) + })) } else { Err(format!("Could not parse gridline")) } @@ -285,3 +292,26 @@ pub fn json_to_vortex(mut json: JsonValue) -> super::euler::VortexParameters { eps, } } + +pub fn h2linspace(start: Float, end: Float, n: usize) -> ndarray::Array1 { + let h = 1.0 / (n - 2) as Float; + ndarray::Array1::from_shape_fn(n, |i| match i { + 0 => start, + i if i == n - 1 => end, + i => h * (i as Float - 0.5), + }) +} + +#[test] +fn test_h2linspace() { + let x = h2linspace(0.0, 1.0, 50); + approx::assert_abs_diff_eq!(x[0], 0.0, epsilon = 1e-6); + approx::assert_abs_diff_eq!(x[49], 1.0, epsilon = 1e-6); + let hend = x[1] - x[0]; + let h = x[2] - x[1]; + approx::assert_abs_diff_eq!(x[49] - x[48], hend, epsilon = 1e-6); + approx::assert_abs_diff_eq!(2.0 * hend, h, epsilon = 1e-6); + for i in 1..48 { + approx::assert_abs_diff_eq!(x[i + 1] - x[i], h, epsilon = 1e-6); + } +}