From 14fefe97abeaedc7dbb18d50d3a7ab20e50bc2df Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Fri, 29 Jan 2021 08:28:55 +0100 Subject: [PATCH] add Matrix approach for SBP8 --- sbp/src/operators/traditional8.rs | 39 +++++++++++++++++++++++++++++-- 1 file changed, 37 insertions(+), 2 deletions(-) diff --git a/sbp/src/operators/traditional8.rs b/sbp/src/operators/traditional8.rs index 6947022..14d0d18 100644 --- a/sbp/src/operators/traditional8.rs +++ b/sbp/src/operators/traditional8.rs @@ -1,4 +1,4 @@ -use super::{diff_op_col, diff_op_row, SbpOperator1d, SbpOperator2d}; +use super::{diff_op_col, diff_op_row, Matrix, RowVector, SbpOperator1d, SbpOperator2d}; use crate::Float; use ndarray::{ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2}; @@ -15,6 +15,10 @@ impl SBP8 { 3.57142857142857e-03, -3.80952380952381e-02, 2.00000000000000e-01, -8.00000000000000e-01, -0.00000000000000e+00, 8.00000000000000e-01, -2.00000000000000e-01, 3.80952380952381e-02, -3.57142857142857e-03 ]; #[rustfmt::skip] + const DIAG_MATRIX: RowVector = Matrix::new([[ + 3.57142857142857e-03, -3.80952380952381e-02, 2.00000000000000e-01, -8.00000000000000e-01, -0.00000000000000e+00, 8.00000000000000e-01, -2.00000000000000e-01, 3.80952380952381e-02, -3.57142857142857e-03 + ]]); + #[rustfmt::skip] const BLOCK: &'static [&'static [Float]] = &[ &[-1.69554360443190e+00, 2.24741246341404e+00, -3.38931922601500e-02, -7.81028168126749e-01, 2.54881486107905e-02, 3.43865227388873e-01, -8.62858162633335e-02, -2.00150583315761e-02], &[-4.34378988266985e-01, 0.00000000000000e+00, 9.18511925072956e-02, 4.94008626807984e-01, -2.46151762937235e-02, -1.86759403432935e-01, 5.27267838475813e-02, 7.16696483080115e-03], @@ -25,6 +29,17 @@ impl SBP8 { &[2.75587615266177e-02, -8.71295642560637e-02, 5.01135077563584e-02, 7.68229253600969e-02, 4.60181213406519e-02, -7.55873581663580e-01, 0.00000000000000e+00, 8.21713248844682e-01, -2.16615355227872e-01, 4.12600676624518e-02, -3.86813134335486e-03], &[5.84767272160451e-03, -1.08336661209337e-02, -1.42810403117803e-02, 3.50919361287023e-02, -1.22244235731112e-02, 1.19411743193552e-01, -7.51668243727123e-01, 0.00000000000000e+00, 7.92601963555477e-01, -1.98150490888869e-01, 3.77429506454989e-02, -3.53840162301552e-03], ]; + #[rustfmt::skip] + const BLOCK_MATRIX: Matrix = Matrix::new([ + [-1.69554360443190e+00, 2.24741246341404e+00, -3.38931922601500e-02, -7.81028168126749e-01, 2.54881486107905e-02, 3.43865227388873e-01, -8.62858162633335e-02, -2.00150583315761e-02, 0.0, 0.0, 0.0, 0.0], + [-4.34378988266985e-01, 0.00000000000000e+00, 9.18511925072956e-02, 4.94008626807984e-01, -2.46151762937235e-02, -1.86759403432935e-01, 5.27267838475813e-02, 7.16696483080115e-03, 0.0, 0.0, 0.0, 0.0], + [3.88218088704253e-02, -5.44329744454984e-01, 0.00000000000000e+00, 3.89516189693211e-01, 1.36433486528546e-01, 1.03290582800845e-01, -1.79720579323281e-01, 5.59882558852296e-02, 0.0, 0.0, 0.0, 0.0], + [1.28088632226564e-01, -4.19172130036008e-01, -5.57707021445779e-02, 0.00000000000000e+00, 1.24714160903055e-01, 2.81285212519100e-01, -3.94470423942641e-02, -1.96981310738430e-02, 0.0, 0.0, 0.0, 0.0], + [-1.82119472519009e-02, 9.09986646154550e-02, -8.51090570277506e-02, -5.43362886365301e-01, 0.00000000000000e+00, 6.37392455438558e-01, -1.02950081118829e-01, 2.98964956216039e-02, -8.65364391190110e-03, 0.0, 0.0, 0.0], + [-7.93147196245203e-02, 2.22875323171502e-01, -2.07999824391436e-02, -3.95611167748401e-01, -2.05756876210586e-01, 0.00000000000000e+00, 5.45876519966127e-01, -9.42727926638298e-02, 2.97971812952850e-02, -2.79348574643297e-03, 0.0, 0.0], + [2.75587615266177e-02, -8.71295642560637e-02, 5.01135077563584e-02, 7.68229253600969e-02, 4.60181213406519e-02, -7.55873581663580e-01, 0.00000000000000e+00, 8.21713248844682e-01, -2.16615355227872e-01, 4.12600676624518e-02, -3.86813134335486e-03, 0.0], + [5.84767272160451e-03, -1.08336661209337e-02, -1.42810403117803e-02, 3.50919361287023e-02, -1.22244235731112e-02, 1.19411743193552e-01, -7.51668243727123e-01, 0.00000000000000e+00, 7.92601963555477e-01, -1.98150490888869e-01, 3.77429506454989e-02, -3.53840162301552e-03], + ]); } impl SbpOperator1d for SBP8 { @@ -59,6 +74,25 @@ impl SbpOperator1d for SBP8 { } } +fn diff_op_row_local(prev: ndarray::ArrayView2, mut fut: ndarray::ArrayViewMut2) { + let mut flipmatrix = SBP8::BLOCK_MATRIX; + flipmatrix *= &-1.0; + + 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( + &SBP8::BLOCK_MATRIX, + &flipmatrix, + &SBP8::DIAG_MATRIX, + super::OperatorType::Normal, + p.as_slice().unwrap(), + f.as_slice_mut().unwrap(), + ) + } +} + impl SbpOperator2d for SBP8 { fn diffxi(&self, prev: ArrayView2, mut fut: ArrayViewMut2) { assert_eq!(prev.shape(), fut.shape()); @@ -69,7 +103,8 @@ impl SbpOperator2d for SBP8 { match (prev.strides(), fut.strides()) { ([_, 1], [_, 1]) => { - diff_op_row(SBP8::BLOCK, SBP8::DIAG, symmetry, optype)(prev, fut); + //diff_op_row(SBP8::BLOCK, SBP8::DIAG, symmetry, optype)(prev, fut); + diff_op_row_local(prev, fut); } ([1, _], [1, _]) => { diff_op_col(SBP8::BLOCK, SBP8::DIAG, symmetry, optype)(prev, fut);