mulimoen
/
SummationByParts
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

test diff operators more thoroughly

This commit is contained in:
Magnus Ulimoen 2020-02-29 10:53:20 +01:00
parent 64e1aec294
commit 2b7380ea3b
5 changed files with 221 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

44
sbp/src/operators.rs
View File

@ -78,3 +78,47 @@ mod traditional4;
pub use traditional4::SBP4;
mod traditional8;
pub use traditional8::SBP8;
#[cfg(test)]
pub(crate) mod testing {
use super::*;
use ndarray::prelude::*;
pub(crate) fn grid_eval<F: Fn(Float, Float) -> Float>(
n: (usize, usize),
f: F,
) -> Array2<Float> {
let nx = n.1;
let dx = 1.0 / (nx - 1) as Float;
let ny = n.0;
let dy = 1.0 / (ny - 1) as Float;
Array2::from_shape_fn(n, |(j, i)| {
let x = dx * i as Float;
let y = dy * j as Float;
f(x, y)
})
}
pub(crate) fn check_operator_on<SBP, F, FX, FY>(
n: (usize, usize),
f: F,
dfdx: FX,
dfdy: FY,
eps: Float,
) where
SBP: SbpOperator,
F: Fn(Float, Float) -> Float,
FX: Fn(Float, Float) -> Float,
FY: Fn(Float, Float) -> Float,
{
let mut y = Array2::zeros(n);
let x = grid_eval(n, f);
y.fill(0.0);
SBP::diffxi(x.view(), y.view_mut());
approx::assert_abs_diff_eq!(&y, &grid_eval(n, dfdx), epsilon = eps);
y.fill(0.0);
SBP::diffeta(x.view(), y.view_mut());
approx::assert_abs_diff_eq!(&y, &grid_eval(n, dfdy), epsilon = eps);
}
}

24
sbp/src/operators/traditional4.rs
View File

@ -45,3 +45,27 @@ impl SbpOperator for SBP4 {
Self::HBLOCK
}
}
#[test]
fn test_trad4() {
use super::testing::*;
use super::*;
let nx = 20;
let ny = 13;
check_operator_on::<SBP4, _, _, _>((ny, nx), |x, y| x + 2.0 * y, |_, _| 1.0, |_, _| 2.0, 1e-4);
check_operator_on::<SBP4, _, _, _>(
(ny, nx),
|x, y| x * x + 2.0 * x * y + 3.0 * y * y,
|x, y| 2.0 * x + 2.0 * y,
|x, y| 2.0 * x + 6.0 * y,
1e-3,
);
check_operator_on::<SBP4, _, _, _>(
(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-1,
);
}

58
sbp/src/operators/traditional8.rs
View File

@ -49,3 +49,61 @@ impl SbpOperator for SBP8 {
Self::HBLOCK
}
}
#[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,
);
}

37
sbp/src/operators/upwind4.rs
View File

@ -451,3 +451,40 @@ impl UpwindOperator for Upwind4 {
Self::dissxi(prev.reversed_axes(), fut.reversed_axes());
}
}
#[test]
fn upwind4_test2() {
use super::testing::*;
use super::*;
let nx = 32;
let ny = 16;
check_operator_on::<Upwind4, _, _, _>(
(ny, nx),
|x, y| x + 2.0 * y,
|_, _| 1.0,
|_, _| 2.0,
1e-4,
);
check_operator_on::<Upwind4, _, _, _>(
(ny, nx),
|x, y| x * x + 2.0 * x * y + 3.0 * y * y,
|x, y| 2.0 * x + 2.0 * y,
|x, y| 2.0 * x + 6.0 * y,
1e-3,
);
check_operator_on::<Upwind4, _, _, _>(
(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-1,
);
check_operator_on::<Upwind4, _, _, _>(
(32, 32),
|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-1,
);
}

58
sbp/src/operators/upwind9.rs
View File

@ -89,3 +89,61 @@ impl UpwindOperator for Upwind9 {
Self::dissxi(prev.reversed_axes(), fut.reversed_axes());
}
}
#[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,
);
}