From 231d218965499634644ae53c0b2a860410e5111d Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Sun, 15 Dec 2019 17:49:41 +0100 Subject: [PATCH] add the trad 8/upwind 9 operator --- src/operators.rs | 4 ++ src/operators/traditional8.rs | 49 +++++++++++++++++++ src/operators/upwind9.rs | 89 +++++++++++++++++++++++++++++++++++ 3 files changed, 142 insertions(+) create mode 100644 src/operators/traditional8.rs create mode 100644 src/operators/upwind9.rs diff --git a/src/operators.rs b/src/operators.rs index 79b4f28..12e3bbd 100644 --- a/src/operators.rs +++ b/src/operators.rs @@ -67,5 +67,9 @@ macro_rules! diff_op_1d { mod upwind4; pub use upwind4::Upwind4; +mod upwind9; +pub use upwind9::Upwind9; mod traditional4; pub use traditional4::SBP4; +mod traditional8; +pub use traditional8::SBP8; diff --git a/src/operators/traditional8.rs b/src/operators/traditional8.rs new file mode 100644 index 0000000..27c30d9 --- /dev/null +++ b/src/operators/traditional8.rs @@ -0,0 +1,49 @@ +use super::SbpOperator; +use crate::diff_op_1d; +use ndarray::{s, ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2}; + +pub struct SBP8 {} + +diff_op_1d!(SBP8, diff_1d, SBP8::BLOCK, SBP8::DIAG, false); + +impl SBP8 { + #[rustfmt::skip] + const HBLOCK: &'static [f32] = &[ + 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 [f32] = &[ + 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 [[f32; 12]] = &[ + [-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.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00], + [-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.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00], + [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.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00], + [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.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00], + [-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.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00], + [-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.00000000000000e+00, 0.00000000000000e+00], + [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.00000000000000e+00], + [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 SbpOperator for SBP8 { + fn diffxi(prev: ArrayView2, mut fut: ArrayViewMut2) { + assert_eq!(prev.shape(), fut.shape()); + assert!(prev.shape()[1] >= 2 * Self::BLOCK.len()); + + for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) { + Self::diff_1d(r0, r1); + } + } + + fn diffeta(prev: ArrayView2, fut: ArrayViewMut2) { + // transpose then use diffxi + Self::diffxi(prev.reversed_axes(), fut.reversed_axes()); + } + + fn h() -> &'static [f32] { + Self::HBLOCK + } +} diff --git a/src/operators/upwind9.rs b/src/operators/upwind9.rs new file mode 100644 index 0000000..8053e40 --- /dev/null +++ b/src/operators/upwind9.rs @@ -0,0 +1,89 @@ +use super::{SbpOperator, UpwindOperator}; +use crate::diff_op_1d; +use ndarray::{s, ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2}; + +pub struct Upwind9 {} + +diff_op_1d!(Upwind9, diff_1d, Upwind9::BLOCK, Upwind9::DIAG, false); +diff_op_1d!( + Upwind9, + diss_1d, + Upwind9::DISS_BLOCK, + Upwind9::DISS_DIAG, + true +); + +impl Upwind9 { + #[rustfmt::skip] + const HBLOCK: &'static [f32] = &[ + 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 DIAG: &'static [f32] = &[ + -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 BLOCK: &'static [[f32; 13]] = &[ + [-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.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00], + [-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.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00], + [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.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00], + [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.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00], + [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.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00], + [-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.00000000000000e+00, 0.00000000000000e+00], + [-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.00000000000000e+00], + [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], + ]; + + #[rustfmt::skip] + const DISS_BLOCK: &'static [[f32; 13]] = &[ + [-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.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00], + [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.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00], + [-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.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00], + [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.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00], + [-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.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00], + [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.00000000000000e+00, 0.00000000000000e+00], + [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.00000000000000e+00], + [-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: &'static [f32] = &[ + 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, + ]; +} + +impl SbpOperator for Upwind9 { + fn diffxi(prev: ArrayView2, mut fut: ArrayViewMut2) { + assert_eq!(prev.shape(), fut.shape()); + assert!(prev.shape()[1] >= 2 * Self::BLOCK.len()); + + for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) { + Self::diff_1d(r0, r1); + } + } + + fn diffeta(prev: ArrayView2, fut: ArrayViewMut2) { + // transpose then use diffxi + Self::diffxi(prev.reversed_axes(), fut.reversed_axes()); + } + + fn h() -> &'static [f32] { + Self::HBLOCK + } +} + +impl UpwindOperator for Upwind9 { + fn dissxi(prev: ArrayView2, mut fut: ArrayViewMut2) { + assert_eq!(prev.shape(), fut.shape()); + assert!(prev.shape()[1] >= 2 * Self::BLOCK.len()); + + for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) { + Self::diss_1d(r0, r1); + } + } + + fn disseta(prev: ArrayView2, fut: ArrayViewMut2) { + // diffeta = transpose then use dissxi + Self::dissxi(prev.reversed_axes(), fut.reversed_axes()); + } +}