use super::{ BlockMatrix, DiagonalMatrix, Matrix, OperatorType, RowVector, SbpOperator1d, SbpOperator2d, UpwindOperator1d, UpwindOperator2d, }; use crate::Float; use ndarray::{ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2}; #[derive(Debug, Copy, Clone)] pub struct Upwind9; impl Upwind9 { #[rustfmt::skip] const H: DiagonalMatrix<8> = DiagonalMatrix::new([ 1070017.0/3628800.0, 5537111.0/3628800.0, 103613.0/403200.0, 261115.0/145152.0, 298951.0/725760.0, 515677.0/403200.0, 3349879.0/3628800.0, 3662753.0/3628800.0 ]); #[rustfmt::skip] const DIFF_DIAG: RowVector = RowVector::new([[ -1.0/1260.0, 5.0/504.0, -5.0/84.0, 5.0/21.0, -5.0/6.0, 0.0, 5.0/6.0, -5.0/21.0, 5.0/84.0, -5.0/504.0, 1.0/1260.0, ]]); #[rustfmt::skip] const DIFF_BLOCK: Matrix = Matrix::new([ [-1.69567399396458e+00, 2.29023358159400e+00, -2.16473500425698e-01, -5.05879766354449e-01, -1.01161106778154e-01, 2.59147072064383e-01, 1.93922119400659e-02, -4.95844980755642e-02, 0.0, 0.0, 0.0, 0.0, 0.0], [-4.42575354959737e-01, 0.00000000000000e+00, 1.91582959381899e-01, 2.82222626681305e-01, 1.12083989713257e-01, -1.51334868892111e-01, -2.23600502721044e-02, 3.03806983474913e-02, 0.0, 0.0, 0.0, 0.0,0.0], [2.48392603571843e-01, -1.13758367065272e+00, 0.00000000000000e+00, 1.95334726810969e+00, -1.58879011773212e+00, 3.93797129320378e-01, 2.52140821030291e-01, -1.21304033647356e-01, 0.0, 0.0, 0.0, 0.0, 0.0], [8.29213105268236e-02, -2.39388470313226e-01, -2.79038666398460e-01, 0.00000000000000e+00, 3.43018053395471e-01, 1.10370852514749e-01, 1.72029988649808e-03, -2.00445645303789e-02, 4.41184918522490e-04, 0.0, 0.0, 0.0, 0.0], [7.24159504343116e-02, -4.15199475743626e-01, 9.91181694804303e-01, -1.49802407438608e+00, 0.00000000000000e+00, 1.30188867830442e+00, -6.03535071819214e-01, 1.73429775718218e-01, -2.40842144699299e-02, 1.92673715759439e-03, 0.0, 0.0, 0.0], [-5.97470838462221e-02, 1.80551858630298e-01, -7.91241454636765e-02, -1.55240829877729e-01, -4.19298775383066e-01, 0.00000000000000e+00, 6.42287612546289e-01, -1.48833147569152e-01, 4.65407609802260e-02, -7.75679349670433e-03, 6.20543479736347e-04, 0.0, 0.0], [-6.19425252179959e-03, 3.69595678895333e-02, -7.01892820620398e-02, -3.35233082197107e-03, 2.69304373763091e-01, -8.89857974743355e-01, 0.00000000000000e+00, 8.66656645522330e-01, -2.57919763669076e-01, 6.44799409172690e-02, -1.07466568195448e-02, 8.59732545563586e-04, 0.0], [1.44853491014330e-02, -4.59275574977554e-02, 3.08833474560615e-02, 3.57240610228828e-02, -7.07760049349999e-02, 1.88587240076292e-01, -7.92626447113877e-01, 0.00000000000000e+00, 8.25608497215073e-01, -2.35888142061449e-01, 5.89720355153623e-02, -9.82867258589373e-03, 7.86293806871498e-04], ]); 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([ [-3.99020778658945e-04, 2.05394169917502e-03, -4.24493243399805e-03, 4.38126393542801e-03, -2.18883813216888e-03, 2.98565988131608e-04, 1.38484104084115e-04, -3.94643819928825e-05, 0.0, 0.0, 0.0, 0.0, 0.0], [3.96913216138553e-04, -2.28230530115522e-03, 5.43069719436758e-03, -6.81086901935894e-03, 4.69064759201504e-03, -1.61429862514855e-03, 1.62083873811316e-04, 2.71310693302277e-05, 0.0, 0.0, 0.0, 0.0, 0.0 ], [-4.87084939816571e-03, 3.22464611075207e-02, -9.06094757860846e-02, 1.39830191253413e-01, -1.27675500367419e-01, 6.87310321912961e-02, -2.00917702215270e-02, 2.43991122096699e-03, 0.0, 0.0, 0.0, 0.0, 0.0], [7.18155125886276e-04, -5.77715378536685e-03, 1.99749582302141e-02, -3.87940986951101e-02, 4.62756436981388e-02, -3.46770570075288e-02, 1.59058082995305e-02, -4.06744078428648e-03, 4.41184918522490e-04, 0.0, 0.0, 0.0, 0.0], [-1.56687484682703e-03, 1.73758484693946e-02, -7.96515646886111e-02, 2.02094401829054e-01, -3.16098733124618e-01, 3.17999240131250e-01, -2.06522928911140e-01, 8.37112455598470e-02, -1.92673715759439e-02, 1.92673715759439e-03, 0.0, 0.0, 0.0], [6.88352254356072e-05, -1.92595810396278e-03, 1.38098624496279e-02, -4.87746083763075e-02, 1.02417890394006e-01, -1.38292226669620e-01, 1.23829022892659e-01, -7.34723830823462e-02, 2.79244565881356e-02, -6.20543479736347e-03, 6.20543479736347e-04, 0.0, 0.0], [4.42345367100640e-05, 2.67913080025652e-04, -5.59301314813691e-03, 3.09954862110834e-02, -9.21529346596015e-02, 1.71559035817103e-01, -2.12738289547735e-01, 1.79835101537893e-01, -1.03167905467630e-01, 3.86879645503614e-02, -8.59732545563586e-03, 8.59732545563586e-04, 0.0], [-1.15289127131636e-05, 4.10149803795578e-05, 6.21188131452618e-04, -7.24912245235322e-03, 3.41622279353287e-02, -9.30972311856124e-02, 1.64473506705108e-01, -1.98013074867399e-01, 1.65121699443015e-01, -9.43552568245798e-02, 3.53832213092174e-02, -7.86293806871498e-03, 7.86293806871498e-04] ]); #[rustfmt::skip] const DISS_DIAG: RowVector = RowVector::new([[ 1.0/1260.0, -1.0/126.0, 1.0/28.0, -2.0/21.0, 1.0/6.0, -1.0/5.0, 1.0/6.0, -2.0/21.0, 1.0/28.0, -1.0/126.0, 1.0/1260.0, ]]); const DISS: BlockMatrix = BlockMatrix::new( Self::DISS_BLOCK, Self::DISS_DIAG, Self::DISS_BLOCK.flip_lr().flip_ud(), ); } impl SbpOperator1d for Upwind9 { fn diff(&self, prev: ArrayView1, fut: ArrayViewMut1) { super::diff_op_1d(&Self::DIFF, super::OperatorType::Normal, prev, fut); } fn h(&self) -> &'static [Float] { &Self::H.start } #[cfg(feature = "sparse")] fn diff_matrix(&self, n: usize) -> sprs::CsMat { super::sparse_from_block(&Self::DIFF, super::OperatorType::Normal, n) } #[cfg(feature = "sparse")] fn h_matrix(&self, n: usize) -> sprs::CsMat { super::h_matrix(&Self::H, n, self.is_h2()) } fn upwind(&self) -> Option<&dyn UpwindOperator1d> { Some(&Self) } } impl SbpOperator2d for Upwind9 { fn diffxi(&self, prev: ArrayView2, fut: ArrayViewMut2) { assert_eq!(prev.shape(), fut.shape()); super::diff_op_2d(&Self::DIFF, OperatorType::Normal, prev, fut); } fn op_xi(&self) -> &dyn SbpOperator1d { &Self } fn upwind(&self) -> Option> { Some(Box::new(Self)) } } impl UpwindOperator1d for Upwind9 { fn diss(&self, prev: ArrayView1, fut: ArrayViewMut1) { super::diff_op_1d(&Self::DISS, super::OperatorType::Normal, prev, fut); } fn as_sbp(&self) -> &dyn SbpOperator1d { self } #[cfg(feature = "sparse")] fn diss_matrix(&self, n: usize) -> sprs::CsMat { super::sparse_from_block(&Self::DISS, super::OperatorType::Normal, n) } } impl UpwindOperator2d for Upwind9 { fn dissxi(&self, prev: ArrayView2, fut: ArrayViewMut2) { assert_eq!(prev.shape(), fut.shape()); super::diff_op_2d(&Self::DISS, OperatorType::Normal, prev, fut); } fn op_xi(&self) -> &dyn UpwindOperator1d { &Self } } #[test] fn test_upwind9() { use super::testing::*; let nx = 32; let ny = 16; // Order one polynomial check_operator_on( Upwind9, (ny, nx), |x, y| x + 2.0 * y, |_x, _y| 1.0, |_x, _y| 2.0, 1e-4, ); // Order two polynomial check_operator_on( Upwind9, (ny, nx), |x, y| x * x + 0.5 * y * y, |x, _y| 2.0 * x, |_x, y| y, 1e-4, ); check_operator_on(Upwind9, (ny, nx), |x, y| x * y, |_x, y| y, |x, _y| x, 1e-4); // Order three polynomials check_operator_on( Upwind9, (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( Upwind9, (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( Upwind9, (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( Upwind9, (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, ); }