use super::{diff_op_col, diff_op_row, Matrix, RowVector, SbpOperator1d, SbpOperator2d}; use crate::Float; use ndarray::{ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2}; #[derive(Debug, Copy, Clone)] pub struct SBP8; impl SBP8 { #[rustfmt::skip] const HBLOCK: &'static [Float] = &[ 2.94890676177879e-01, 1.52572062389771e+00, 2.57452876984127e-01, 1.79811370149912e+00, 4.12708057760141e-01, 1.27848462301587e+00, 9.23295579805997e-01, 1.00933386085916e+00 ]; #[rustfmt::skip] const DIAG: &'static [Float] = &[ 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], &[3.88218088704253e-02, -5.44329744454984e-01, 0.00000000000000e+00, 3.89516189693211e-01, 1.36433486528546e-01, 1.03290582800845e-01, -1.79720579323281e-01, 5.59882558852296e-02], &[1.28088632226564e-01, -4.19172130036008e-01, -5.57707021445779e-02, 0.00000000000000e+00, 1.24714160903055e-01, 2.81285212519100e-01, -3.94470423942641e-02, -1.96981310738430e-02], &[-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], &[-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], &[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], ]); #[rustfmt::skip] const BLOCKEND_MATRIX: Matrix = Matrix::new([ [0.0035384016230155199, -0.037742950645498902, 0.19815049088886899, -0.79260196355547696, -0.0, 0.75166824372712304, -0.11941174319355199, 0.012224423573111201, -0.035091936128702303, 0.0142810403117803, 0.010833666120933699, -0.0058476727216045096], [-0.0, 0.0038681313433548601, -0.041260067662451799, 0.21661535522787201, -0.821713248844682, -0.0, 0.75587358166357999, -0.046018121340651898, -0.076822925360096897, -0.050113507756358401, 0.087129564256063705, -0.027558761526617698], [-0.0, -0.0, 0.0027934857464329702, -0.029797181295285, 0.094272792663829805, -0.54587651996612696, -0.0, 0.205756876210586, 0.39561116774840099, 0.020799982439143602, -0.22287532317150199, 0.079314719624520302], [-0.0, -0.0, -0.0, 0.0086536439119011006, -0.029896495621603899, 0.102950081118829, -0.63739245543855805, -0.0, 0.54336288636530095, 0.085109057027750595, -0.090998664615454999, 0.01821194725190089], [-0.0, -0.0, -0.0, -0.0, 0.019698131073842998, 0.039447042394264102, -0.28128521251909999, -0.124714160903055, -0.0, 0.055770702144577897, 0.41917213003600801, -0.128088632226564], [-0.0, -0.0, -0.0, -0.0, -0.055988255885229599, 0.179720579323281, -0.10329058280084499, -0.13643348652854601, -0.389516189693211, -0.0, 0.54432974445498405, -0.038821808870425301], [-0.0, -0.0, -0.0, -0.0, -0.0071669648308011497, -0.052726783847581303, 0.18675940343293501, 0.024615176293723499, -0.49400862680798402, -0.091851192507295606, -0.0, 0.43437898826698501], [-0.0, -0.0, -0.0, -0.0, 0.020015058331576099, 0.086285816263333495, -0.34386522738887298, -0.0254881486107905, 0.78102816812674902, 0.033893192260150001, -2.2474124634140402, 1.6955436044319001], ]); } impl SbpOperator1d for SBP8 { fn diff(&self, prev: ArrayView1, fut: ArrayViewMut1) { super::diff_op_1d( Self::BLOCK, Self::DIAG, super::Symmetry::AntiSymmetric, super::OperatorType::Normal, prev, fut, ) } fn h(&self) -> &'static [Float] { Self::HBLOCK } #[cfg(feature = "sparse")] fn diff_matrix(&self, n: usize) -> sprs::CsMat { super::sparse_from_block( Self::BLOCK, Self::DIAG, super::Symmetry::AntiSymmetric, super::OperatorType::Normal, n, ) } #[cfg(feature = "sparse")] fn h_matrix(&self, n: usize) -> sprs::CsMat { super::h_matrix(Self::HBLOCK, n, self.is_h2()) } } 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, &SBP8::BLOCKEND_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()); assert!(prev.shape()[1] >= 2 * SBP8::BLOCK.len()); let symmetry = super::Symmetry::AntiSymmetric; let optype = super::OperatorType::Normal; match (prev.strides(), fut.strides()) { ([_, 1], [_, 1]) => { //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); } ([_, _], [_, _]) => { // Fallback, work row by row for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) { SBP8.diff(r0, r1); } } _ => unreachable!("Should only be two elements in the strides vectors"), } } fn op_xi(&self) -> &dyn SbpOperator1d { &Self } fn op_eta(&self) -> &dyn SbpOperator1d { &Self } } #[test] fn test_trad8() { use super::testing::*; let nx = 32; let ny = 16; // Order one polynomial check_operator_on( SBP8, (ny, nx), |x, y| x + 2.0 * y, |_x, _y| 1.0, |_x, _y| 2.0, 1e-4, ); // Order two polynomial 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(SBP8, (ny, nx), |x, y| x * y, |_x, y| y, |x, _y| x, 1e-4); // Order three polynomials 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( 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( 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), |x, y| 2.0 * x.powi(2) + 6.0 * x * y + 12.0 * y.powi(2), 1e-4, ); // Order four polynomials 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), |x, y| x.powi(3) + 2.0 * x.powi(2) * y + 3.0 * x * y.powi(2) + 4.0 * y.powi(3), 1e-4, ); } #[test] 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)) }