add the trad 8/upwind 9 operator

This commit is contained in:
Magnus Ulimoen 2019-12-15 17:49:41 +01:00
parent 03a4de5ab0
commit 231d218965
3 changed files with 142 additions and 0 deletions

View File

@ -67,5 +67,9 @@ macro_rules! diff_op_1d {
mod upwind4; mod upwind4;
pub use upwind4::Upwind4; pub use upwind4::Upwind4;
mod upwind9;
pub use upwind9::Upwind9;
mod traditional4; mod traditional4;
pub use traditional4::SBP4; pub use traditional4::SBP4;
mod traditional8;
pub use traditional8::SBP8;

View File

@ -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<f32>, mut fut: ArrayViewMut2<f32>) {
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<f32>, fut: ArrayViewMut2<f32>) {
// transpose then use diffxi
Self::diffxi(prev.reversed_axes(), fut.reversed_axes());
}
fn h() -> &'static [f32] {
Self::HBLOCK
}
}

89
src/operators/upwind9.rs Normal file
View File

@ -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<f32>, mut fut: ArrayViewMut2<f32>) {
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<f32>, fut: ArrayViewMut2<f32>) {
// 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<f32>, mut fut: ArrayViewMut2<f32>) {
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<f32>, fut: ArrayViewMut2<f32>) {
// diffeta = transpose then use dissxi
Self::dissxi(prev.reversed_axes(), fut.reversed_axes());
}
}