diff --git a/sbp/Cargo.toml b/sbp/Cargo.toml index 2fe2450..a84bb47 100644 --- a/sbp/Cargo.toml +++ b/sbp/Cargo.toml @@ -10,6 +10,11 @@ approx = "0.3.2" packed_simd = "0.3.3" json = "0.12.1" +[features] +# Internal feature flag to gate the expensive tests +# which should be run only in release builds +expensive_tests = [] + [dev-dependencies] criterion = "0.3.0" diff --git a/sbp/src/euler.rs b/sbp/src/euler.rs index 4370898..a9d2512 100644 --- a/sbp/src/euler.rs +++ b/sbp/src/euler.rs @@ -270,7 +270,8 @@ impl Field { } impl Field { - pub fn err_diff(&self, other: &Self) -> f32 { + /// sqrt((self-other)^T*H*(self-other)) + pub fn h2_err(&self, other: &Self) -> f32 { assert_eq!(self.nx(), other.nx()); assert_eq!(self.ny(), other.ny()); @@ -312,6 +313,7 @@ impl Field { .zip(other.0.iter()) .map(|(((hx, hy), r0), r1)| (*r0 - *r1).powi(2) * hx * hy) .sum::() + .sqrt() } } @@ -323,10 +325,10 @@ fn h2_diff() { } let field1 = Field::new(20, 21); - assert!((field0.err_diff::(&field1) - 4.0).abs() < 1e-3); - assert!((field0.err_diff::(&field1) - 4.0).abs() < 1e-3); - assert!((field0.err_diff::(&field1) - 4.0).abs() < 1e-3); - assert!((field0.err_diff::(&field1) - 4.0).abs() < 1e-3); + assert!((field0.h2_err::(&field1).powi(2) - 4.0).abs() < 1e-3); + assert!((field0.h2_err::(&field1).powi(2) - 4.0).abs() < 1e-3); + assert!((field0.h2_err::(&field1).powi(2) - 4.0).abs() < 1e-3); + assert!((field0.h2_err::(&field1).powi(2) - 4.0).abs() < 1e-3); } #[derive(Copy, Clone)] diff --git a/sbp/tests/convergence.rs b/sbp/tests/convergence.rs index 082b246..aee9f46 100644 --- a/sbp/tests/convergence.rs +++ b/sbp/tests/convergence.rs @@ -1,9 +1,10 @@ +#![cfg(feature = "expensive_tests")] use ndarray::prelude::*; use sbp::euler::*; fn run_with_size(size: usize) -> f32 { - let nx = 4 * size; - let ny = 4 * size; + let nx = size; + let ny = size; let x = Array1::linspace(-5.0, 5.0, nx); let y = Array1::linspace(-5.0, 5.0, ny); @@ -26,7 +27,7 @@ fn run_with_size(size: usize) -> f32 { let mut sys = System::::new(x, y); sys.vortex(0.0, vortex_params); - let time = 0.1; + let time = 0.2; let dt = 0.2 * f32::min(1.0 / (nx - 1) as f32, 1.0 / (ny - 1) as f32); let nsteps = (time / dt) as usize; @@ -37,17 +38,30 @@ fn run_with_size(size: usize) -> f32 { let mut verifield = Field::new(ny, nx); verifield.vortex(sys.x(), sys.y(), nsteps as f32 * dt, vortex_params); - verifield.err_diff::(sys.field()) + verifield.h2_err::(sys.field()) } #[test] -fn test() { - let sizes = [25, 35, 50, 71, 100]; - let mut errs = Vec::with_capacity(sizes.len()); +fn convergence() { + let sizes = [25, 35, 50, 71, 100, 150]; + let mut prev: Option<(usize, f32)> = None; + println!("Size\tError(h2)\tq"); for size in &sizes { - println!("Size: {}", size); - let e = run_with_size::(*size); - errs.push(e); + print!("{:3}x{:3}", size, size); + let e = run_with_size::(*size); + print!("\t{:.10}", e); + if let Some(prev) = prev.take() { + let m0 = size * size; + let e0 = e; + + let (size1, e1) = prev; + let m1 = size1 * size1; + + let q = f32::log10(e0 / e1) / f32::log10((m0 as f32 / m1 as f32).powf(1.0 / 2.0)); + print!("\t{}", q); + } + println!(); + prev = Some((*size, e)); } - panic!("{:?}", errs); + panic!(); }