diff --git a/sbp/src/operators/algos/constmatrix.rs b/sbp/src/operators/algos/constmatrix.rs index 7c457cc..1f38d02 100644 --- a/sbp/src/operators/algos/constmatrix.rs +++ b/sbp/src/operators/algos/constmatrix.rs @@ -113,7 +113,6 @@ impl RowVector { } impl Matrix { - #[inline(always)] pub fn matmul_into(&mut self, lhs: &Matrix, rhs: &Matrix) where T: Default + Copy + core::ops::Mul + core::ops::Add, @@ -302,6 +301,28 @@ impl Matrix { } m } + + /// Zero extends if larger than self + pub(crate) const fn resize( + &self, + ) -> Matrix { + let mut m = Matrix::new([[0.0; N2]; M2]); + + let m_min = if M < M2 { M } else { M2 }; + let n_min = if N < N2 { N } else { N2 }; + + let mut i = 0; + while i < m_min { + let mut j = 0; + while j < n_min { + m.data[i][j] = self.data[i][j]; + j += 1; + } + i += 1; + } + + m + } } mod flipping { diff --git a/sbp/src/operators/traditional4.rs b/sbp/src/operators/traditional4.rs index 5cac183..2bb2ed6 100644 --- a/sbp/src/operators/traditional4.rs +++ b/sbp/src/operators/traditional4.rs @@ -18,16 +18,16 @@ impl SBP4 { 1.0 / 12.0, -2.0 / 3.0, 0.0, 2.0 / 3.0, -1.0 / 12.0 ]]); #[rustfmt::skip] - const DIFF_BLOCK: Matrix = Matrix::new([ + const MINIMAL_DIFF_BLOCK: Matrix = Matrix::new([ [-24.0/17.0, 59.0/34.0, -4.0/17.0, -3.0/34.0, 0.0, 0.0], [-1.0/2.0, 0.0, 1.0/2.0, 0.0, 0.0, 0.0], [4.0/43.0, -59.0/86.0, 0.0, 59.0/86.0, -4.0/43.0, 0.0], [3.0/98.0, 0.0, -59.0/98.0, 0.0, 32.0/49.0, -4.0/49.0] ]); - const DIFF_BLOCKEND: super::Matrix = - Self::DIFF_BLOCK.flip_lr().flip_ud().flip_sign(); + const DIFF_BLOCK: Matrix = Self::MINIMAL_DIFF_BLOCK.resize(); + const DIFF_BLOCKEND: Matrix = Self::DIFF_BLOCK.flip_lr().flip_ud().flip_sign(); - const DIFF: BlockMatrix = + const DIFF: BlockMatrix = BlockMatrix::new(Self::DIFF_BLOCK, Self::DIFF_DIAG, Self::DIFF_BLOCKEND); #[rustfmt::skip] diff --git a/sbp/src/operators/upwind4.rs b/sbp/src/operators/upwind4.rs index a24aec7..640d8ec 100644 --- a/sbp/src/operators/upwind4.rs +++ b/sbp/src/operators/upwind4.rs @@ -14,35 +14,38 @@ impl Upwind4 { 49.0 / 144.0, 61.0 / 48.0, 41.0 / 48.0, 149.0 / 144.0 ]); #[rustfmt::skip] - const DIFF_BLOCK: Matrix = Matrix::new([ + const MINIMAL_DIFF_BLOCK: 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], ]); + const DIFF_BLOCK: Matrix = Self::MINIMAL_DIFF_BLOCK.resize(); + #[rustfmt::skip] const DIFF_DIAG: 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 ]]); - const DIFF_BLOCKEND: Matrix = Self::DIFF_BLOCK.flip_lr().flip_ud().flip_sign(); + const DIFF_BLOCKEND: Matrix = Self::DIFF_BLOCK.flip_lr().flip_ud().flip_sign(); - const DIFF: BlockMatrix = + const DIFF: BlockMatrix = BlockMatrix::new(Self::DIFF_BLOCK, Self::DIFF_DIAG, Self::DIFF_BLOCKEND); #[rustfmt::skip] - const DISS_BLOCK: Matrix = Matrix::new([ + const MINIMAL_DISS_BLOCK: Matrix = Matrix::new([ [-3.0 / 49.0, 9.0 / 49.0, -9.0 / 49.0, 3.0 / 49.0, 0.0, 0.0, 0.0], [ 3.0 / 61.0, -11.0 / 61.0, 15.0 / 61.0, -9.0 / 61.0, 2.0 / 61.0, 0.0, 0.0], [-3.0 / 41.0, 15.0 / 41.0, -29.0 / 41.0, 27.0 / 41.0, -12.0 / 41.0, 2.0 / 41.0, 0.0], [3.0 / 149.0, -27.0 / 149.0, 81.0 / 149.0, -117.0 / 149.0, 90.0 / 149.0, -36.0 / 149.0, 6.0 / 149.0], ]); + const DISS_BLOCK: Matrix = Self::MINIMAL_DISS_BLOCK.resize(); #[rustfmt::skip] const DISS_DIAG: RowVector = Matrix::new([[ 1.0 / 24.0, -1.0 / 4.0, 5.0 / 8.0, -5.0 / 6.0, 5.0 / 8.0, -1.0 / 4.0, 1.0 / 24.0 ]]); - const DISS_BLOCKEND: Matrix = Self::DISS_BLOCK.flip_lr().flip_ud(); + const DISS_BLOCKEND: Matrix = Self::DISS_BLOCK.flip_lr().flip_ud(); - const DISS: BlockMatrix = + const DISS: BlockMatrix = BlockMatrix::new(Self::DISS_BLOCK, Self::DISS_DIAG, Self::DISS_BLOCKEND); } diff --git a/sbp/src/operators/upwind4h2.rs b/sbp/src/operators/upwind4h2.rs index 711c9b3..ab98006 100644 --- a/sbp/src/operators/upwind4h2.rs +++ b/sbp/src/operators/upwind4h2.rs @@ -14,35 +14,37 @@ impl Upwind4h2 { 0.91e2/0.720e3, 0.325e3/0.384e3, 0.595e3/0.576e3, 0.1909e4/0.1920e4, ]); #[rustfmt::skip] - const DIFF_BLOCK: Matrix = Matrix::new([ + const MINIMAL_DIFF_BLOCK: Matrix = Matrix::new([ [-3.95604395604396e+00, 5.41758241758242e+00, -1.94505494505495e+00, 4.83516483516484e-01, 0.0, 0.0, 0.0], [-8.09025641025641e-01, 0.00000000000000e+00, 1.03948717948718e+00, -2.47384615384615e-01, 1.69230769230769e-02, 0.0, 0.0], [2.37983193277311e-01, -8.51680672268908e-01, 0.00000000000000e+00, 7.35966386554622e-01, -1.36134453781513e-01, 1.38655462184874e-02, 0.0], [-6.14632442814737e-02, 2.10581456259822e-01, -7.64623712240265e-01, 0.00000000000000e+00, 7.42535358826611e-01, -1.41435306443164e-01, 1.44054478784704e-02], ]); + const DIFF_BLOCK: Matrix = Self::MINIMAL_DIFF_BLOCK.resize(); #[rustfmt::skip] const DIFF_DIAG: RowVector = RowVector::new([[ -1.43229166666667e-02, 1.40625000000000e-01, -7.38281250000000e-01, 0.00000000000000e+00, 7.38281250000000e-01, -1.40625000000000e-01, 1.43229166666667e-02 ]]); - const DIFF: BlockMatrix = BlockMatrix::new( + const DIFF: BlockMatrix = BlockMatrix::new( Self::DIFF_BLOCK, Self::DIFF_DIAG, Self::DIFF_BLOCK.flip_lr().flip_ud().flip_sign(), ); #[rustfmt::skip] - const DISS_BLOCK: Matrix = Matrix::new([ + const MINIMAL_DISS_BLOCK: Matrix = Matrix::new([ [-2.76989342315976e-01, 5.19355016842454e-01, -3.46236677894969e-01, 1.03871003368491e-01, 0.0, 0.0, 0.0], [7.75570158484731e-02, -1.62342481638964e-01, 1.47715500579822e-01, -7.98531117124082e-02, 1.69230769230769e-02, 0.0, 0.0], [-4.23630758836198e-02, 1.21027405937249e-01, -1.91609307039399e-01, 1.82272708078206e-01, -8.31932773109244e-02, 1.38655462184874e-02, 0.0], [1.32037874021759e-02, -6.79734450144910e-02, 1.89370108794365e-01, -2.78654929966754e-01, 2.16081718177056e-01, -8.64326872708224e-02, 1.44054478784704e-02], ]); + const DISS_BLOCK: Matrix = Self::MINIMAL_DISS_BLOCK.resize(); #[rustfmt::skip] const DISS_DIAG: RowVector = RowVector::new([[ 1.43229166666667e-02, -8.59375000000000e-02, 2.14843750000000e-01, -2.86458333333333e-01, 2.14843750000000e-01, -8.59375000000000e-02, 1.43229166666667e-02, ]]); - const DISS: BlockMatrix = BlockMatrix::new( + const DISS: BlockMatrix = BlockMatrix::new( Self::DISS_BLOCK, Self::DISS_DIAG, Self::DISS_BLOCK.flip_lr().flip_ud(),