diff --git a/sbp/src/operators/algos.rs b/sbp/src/operators/algos.rs index 6c579e5..a504b4c 100644 --- a/sbp/src/operators/algos.rs +++ b/sbp/src/operators/algos.rs @@ -178,6 +178,71 @@ pub(crate) mod constmatrix { let m3 = m1 * m2; assert_eq!(m3, Matrix::new([[74, 80, 86, 92], [173, 188, 203, 218]])); } + #[test] + fn iter() { + let m = Matrix::new([[1_u8, 2, 3], [4, 5, 6]]); + let mut iter = m.iter(); + assert_eq!(iter.next(), Some(&1)); + assert_eq!(iter.next(), Some(&2)); + assert_eq!(iter.next(), Some(&3)); + assert_eq!(iter.next(), Some(&4)); + assert_eq!(iter.next(), Some(&5)); + assert_eq!(iter.next(), Some(&6)); + assert_eq!(iter.next(), None); + } + } + + mod approx { + use super::Matrix; + use ::approx::{AbsDiffEq, RelativeEq, UlpsEq}; + + impl AbsDiffEq for Matrix + where + T: AbsDiffEq, + { + type Epsilon = T::Epsilon; + fn default_epsilon() -> Self::Epsilon { + T::default_epsilon() + } + fn abs_diff_eq(&self, other: &Self, epsilon: Self::Epsilon) -> bool { + self.iter() + .zip(other.iter()) + .all(|(r, l)| r.abs_diff_eq(l, T::default_epsilon())) + } + } + impl RelativeEq for Matrix + where + T: RelativeEq, + Self::Epsilon: Copy, + { + fn default_max_relative() -> Self::Epsilon { + T::default_max_relative() + } + fn relative_eq( + &self, + other: &Self, + epsilon: Self::Epsilon, + max_relative: Self::Epsilon, + ) -> bool { + self.iter() + .zip(other.iter()) + .all(|(r, l)| r.relative_eq(l, epsilon, max_relative)) + } + } + impl UlpsEq for Matrix + where + T: UlpsEq, + Self::Epsilon: Copy, + { + fn default_max_ulps() -> u32 { + T::default_max_ulps() + } + fn ulps_eq(&self, other: &Self, epsilon: Self::Epsilon, max_ulps: u32) -> bool { + self.iter() + .zip(other.iter()) + .all(|(r, l)| r.ulps_eq(l, epsilon, max_ulps)) + } + } } } @@ -228,6 +293,7 @@ pub(crate) fn diff_op_1d_matrix( *f = diff * idx; } + let prev = prev.slice(ndarray::s![nx - N..]); for (bl, f) in blockend.iter_rows().zip(fut.iter_mut().rev().take(M).rev()) { let diff = bl .iter() diff --git a/sbp/src/operators/traditional4.rs b/sbp/src/operators/traditional4.rs index 57afd35..db5802f 100644 --- a/sbp/src/operators/traditional4.rs +++ b/sbp/src/operators/traditional4.rs @@ -240,8 +240,5 @@ fn block_equality() { let mut flipped_inverted = SBP4::BLOCK_MATRIX.flip(); flipped_inverted *= -1.0; - assert!(flipped_inverted - .iter() - .zip(SBP4::BLOCKEND_MATRIX.iter()) - .all(|(x, y)| (x - y).abs() < 1e-3)) + approx::assert_ulps_eq!(SBP4::BLOCKEND_MATRIX, flipped_inverted, max_ulps = 1); } diff --git a/sbp/src/operators/traditional8.rs b/sbp/src/operators/traditional8.rs index 88631b3..48891f8 100644 --- a/sbp/src/operators/traditional8.rs +++ b/sbp/src/operators/traditional8.rs @@ -221,8 +221,5 @@ fn block_equality() { let mut flipped_inverted = SBP8::BLOCK_MATRIX.flip(); flipped_inverted *= -1.0; - assert!(flipped_inverted - .iter() - .zip(SBP8::BLOCKEND_MATRIX.iter()) - .all(|(x, y)| (x - y).abs() < 1e-12)) + approx::assert_ulps_eq!(SBP8::BLOCKEND_MATRIX, flipped_inverted, max_ulps = 1); } diff --git a/sbp/src/operators/upwind4.rs b/sbp/src/operators/upwind4.rs index 8e577b6..66b9852 100644 --- a/sbp/src/operators/upwind4.rs +++ b/sbp/src/operators/upwind4.rs @@ -1,5 +1,6 @@ use super::{ - diff_op_col, diff_op_row, SbpOperator1d, SbpOperator2d, UpwindOperator1d, UpwindOperator2d, + diff_op_col, diff_op_row, Matrix, RowVector, SbpOperator1d, SbpOperator2d, UpwindOperator1d, + UpwindOperator2d, }; use crate::Float; use ndarray::{ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2, Axis}; @@ -174,12 +175,30 @@ impl Upwind4 { -1.0 / 24.0, 1.0 / 4.0, -7.0 / 8.0, 0.0, 7.0 / 8.0, -1.0 / 4.0, 1.0 / 24.0 ]; #[rustfmt::skip] + const DIAG_MATRIX : RowVector = RowVector::new([ + [-1.0 / 24.0, 1.0 / 4.0, -7.0 / 8.0, 0.0, 7.0 / 8.0, -1.0 / 4.0, 1.0 / 24.0] + ]); + #[rustfmt::skip] const BLOCK: &'static [&'static [Float]] = &[ &[ -72.0 / 49.0, 187.0 / 98.0, -20.0 / 49.0, -3.0 / 98.0, 0.0, 0.0, 0.0], &[-187.0 / 366.0, 0.0, 69.0 / 122.0, -16.0 / 183.0, 2.0 / 61.0, 0.0, 0.0], &[ 20.0 / 123.0, -69.0 / 82.0, 0.0, 227.0 / 246.0, -12.0 / 41.0, 2.0 / 41.0, 0.0], &[ 3.0 / 298.0, 16.0 / 149.0, -227.0 / 298.0, 0.0, 126.0 / 149.0, -36.0 / 149.0, 6.0 / 149.0], ]; + #[rustfmt::skip] + const BLOCK_MATRIX: Matrix = Matrix::new([ + [ -72.0 / 49.0, 187.0 / 98.0, -20.0 / 49.0, -3.0 / 98.0, 0.0, 0.0, 0.0], + [-187.0 / 366.0, 0.0, 69.0 / 122.0, -16.0 / 183.0, 2.0 / 61.0, 0.0, 0.0], + [ 20.0 / 123.0, -69.0 / 82.0, 0.0, 227.0 / 246.0, -12.0 / 41.0, 2.0 / 41.0, 0.0], + [ 3.0 / 298.0, 16.0 / 149.0, -227.0 / 298.0, 0.0, 126.0 / 149.0, -36.0 / 149.0, 6.0 / 149.0], + ]); + #[rustfmt::skip] + const BLOCKEND_MATRIX: Matrix = Matrix::new([ + [ -6.0 / 149.0, 36.0 / 149.0, -126.0 / 149.0, 0.0, 227.0 / 298.0, -16.0 / 149.0, -3.0 / 298.0], + [ 0.0, -2.0/41.0, 12.0/41.0, -227.0/246.0, 0.0, 69.0/82.0, -20.0/123.0], + [ 0.0, 0.0, -2.0/61.0, 16.0/183.0, -69.0/122.0, 0.0, 187.0/366.0], + [ 0.0, 0.0, 0.0, 3.0 / 98.0, 20.0 / 49.0, -187.0 / 98.0, 72.0/49.0], + ]); #[rustfmt::skip] const DISS_BLOCK: &'static [&'static [Float]] = &[ @@ -197,10 +216,10 @@ impl Upwind4 { impl SbpOperator1d for Upwind4 { fn diff(&self, prev: ArrayView1, fut: ArrayViewMut1) { - super::diff_op_1d( - Self::BLOCK, - Self::DIAG, - super::Symmetry::AntiSymmetric, + super::diff_op_1d_matrix( + &Self::BLOCK_MATRIX, + &Self::BLOCKEND_MATRIX, + &Self::DIAG_MATRIX, super::OperatorType::Normal, prev, fut, @@ -229,18 +248,29 @@ impl SbpOperator1d for Upwind4 { } } +fn diff_op_row_local(prev: ndarray::ArrayView2, mut fut: ndarray::ArrayViewMut2) { + for (p, mut f) in prev + .axis_iter(ndarray::Axis(0)) + .zip(fut.axis_iter_mut(ndarray::Axis(0))) + { + super::diff_op_1d_slice_matrix( + &Upwind4::BLOCK_MATRIX, + &Upwind4::BLOCKEND_MATRIX, + &Upwind4::DIAG_MATRIX, + super::OperatorType::Normal, + p.as_slice().unwrap(), + f.as_slice_mut().unwrap(), + ) + } +} + impl SbpOperator2d for Upwind4 { fn diffxi(&self, prev: ArrayView2, mut fut: ArrayViewMut2) { assert_eq!(prev.shape(), fut.shape()); assert!(prev.shape()[1] >= 2 * Upwind4::BLOCK.len()); match (prev.strides(), fut.strides()) { - ([_, 1], [_, 1]) => diff_op_row( - Upwind4::BLOCK, - Upwind4::DIAG, - super::Symmetry::AntiSymmetric, - super::OperatorType::Normal, - )(prev, fut), + ([_, 1], [_, 1]) => diff_op_row_local(prev, fut), ([1, _], [1, _]) if prev.len_of(Axis(0)) % SimdT::lanes() == 0 => { diff_simd_col(prev, fut) } @@ -465,3 +495,11 @@ fn upwind4_test2() { 1e-1, ); } + +#[test] +fn block_equality() { + let mut flipped_inverted = Upwind4::BLOCK_MATRIX.flip(); + flipped_inverted *= -1.0; + + approx::assert_ulps_eq!(Upwind4::BLOCKEND_MATRIX, flipped_inverted, max_ulps = 1); +}