From 3bffd184881697685b5e593c22890e3908f8ffd5 Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Mon, 30 Mar 2020 23:15:28 +0200 Subject: [PATCH] more tests --- sbp/src/euler.rs | 157 +++++++++++++++++++++++++++---------- sbp/tests/convergence.rs | 16 +++- sbp/tests/single_period.rs | 47 +++++++++++ 3 files changed, 174 insertions(+), 46 deletions(-) create mode 100644 sbp/tests/single_period.rs diff --git a/sbp/src/euler.rs b/sbp/src/euler.rs index a2aea48..1da5a4e 100644 --- a/sbp/src/euler.rs +++ b/sbp/src/euler.rs @@ -31,13 +31,14 @@ impl System { } pub fn advance(&mut self, dt: Float) { - let rhs_trad = |k: &mut Field, y: &Field, grid: &_, wb: &mut _| { - let boundaries = BoundaryTerms { - north: y.south(), - south: y.north(), - west: y.east(), - east: y.west(), - }; + let bc = BoundaryCharacteristics { + north: BoundaryCharacteristic::This, + south: BoundaryCharacteristic::This, + east: BoundaryCharacteristic::This, + west: BoundaryCharacteristic::This, + }; + let rhs_trad = |k: &mut Field, y: &Field, grid: &Grid<_>, wb: &mut _| { + let boundaries = boundary_extractor(y, grid, &bc); RHS_trad(k, y, grid, &boundaries, wb) }; integrate::rk4( @@ -91,13 +92,14 @@ impl System { impl System { pub fn advance_upwind(&mut self, dt: Float) { - let rhs_upwind = |k: &mut Field, y: &Field, grid: &_, wb: &mut _| { - let boundaries = BoundaryTerms { - north: y.south(), - south: y.north(), - west: y.east(), - east: y.west(), - }; + let bc = BoundaryCharacteristics { + north: BoundaryCharacteristic::This, + south: BoundaryCharacteristic::This, + east: BoundaryCharacteristic::This, + west: BoundaryCharacteristic::This, + }; + let rhs_upwind = |k: &mut Field, y: &Field, grid: &Grid<_>, wb: &mut _| { + let boundaries = boundary_extractor(y, grid, &bc); RHS_upwind(k, y, grid, &boundaries, wb) }; integrate::rk4( @@ -240,34 +242,18 @@ impl Field { assert_eq!(x.shape()[0], self.ny()); let (rho, rhou, rhov, e) = self.components_mut(); + let n = rho.len(); - let eps = vortex_param.eps; - let m = vortex_param.mach; - let rstar = vortex_param.rstar; - let p_inf = 1.0 / (GAMMA * m * m); - azip!((rho in rho, - rhou in rhou, - rhov in rhov, - e in e, - x in x, - y in y) - { - use crate::consts::PI; - - let dx = (x - vortex_param.x0) - t; - let dy = y - vortex_param.y0; - let f = (1.0 - (dx*dx + dy*dy))/(rstar*rstar); - - *rho = Float::powf(1.0 - eps*eps*(GAMMA - 1.0)*m*m/(8.0*PI*PI*p_inf*rstar*rstar)*f.exp(), 1.0/(GAMMA - 1.0)); - assert!(*rho > 0.0); - let p = Float::powf(*rho, GAMMA)*p_inf; - assert!(p > 0.0); - let u = 1.0 - eps*dy/(2.0*PI*p_inf.sqrt()*rstar*rstar)*(f/2.0).exp(); - let v = eps*dx/(2.0*PI*p_inf.sqrt()*rstar*rstar)*(f/2.0).exp(); - *rhou = *rho*u; - *rhov = *rho*v; - *e = p/(GAMMA - 1.0) + *rho*(u*u + v*v)/2.0; - }); + vortex( + rho.into_shape((n,)).unwrap(), + rhou.into_shape((n,)).unwrap(), + rhov.into_shape((n,)).unwrap(), + e.into_shape((n,)).unwrap(), + x.into_shape((n,)).unwrap(), + y.into_shape((n,)).unwrap(), + t, + vortex_param, + ) } } @@ -333,7 +319,7 @@ fn h2_diff() { assert!((field0.h2_err::(&field1).powi(2) - 4.0).abs() < 1e-3); } -#[derive(Copy, Clone)] +#[derive(Copy, Clone, Debug)] pub struct VortexParameters { pub x0: Float, pub y0: Float, @@ -342,6 +328,52 @@ pub struct VortexParameters { pub mach: Float, } +pub fn vortex( + rho: ArrayViewMut1, + rhou: ArrayViewMut1, + rhov: ArrayViewMut1, + e: ArrayViewMut1, + x: ArrayView1, + y: ArrayView1, + t: Float, + vortex_param: VortexParameters, +) { + assert_eq!(rho.len(), rhou.len()); + assert_eq!(rho.len(), rhov.len()); + assert_eq!(rho.len(), e.len()); + assert_eq!(rho.len(), x.len()); + assert_eq!(rho.len(), y.len()); + assert_eq!(x.shape(), y.shape()); + + let eps = vortex_param.eps; + let m = vortex_param.mach; + let rstar = vortex_param.rstar; + let p_inf = 1.0 / (GAMMA * m * m); + azip!((rho in rho, + rhou in rhou, + rhov in rhov, + e in e, + x in x, + y in y) + { + use crate::consts::PI; + + let dx = (x - vortex_param.x0) - t; + let dy = y - vortex_param.y0; + let f = (1.0 - (dx*dx + dy*dy))/(rstar*rstar); + + *rho = Float::powf(1.0 - eps*eps*(GAMMA - 1.0)*m*m/(8.0*PI*PI*p_inf*rstar*rstar)*f.exp(), 1.0/(GAMMA - 1.0)); + assert!(*rho > 0.0); + let p = Float::powf(*rho, GAMMA)*p_inf; + assert!(p > 0.0); + let u = 1.0 - eps*dy/(2.0*PI*p_inf.sqrt()*rstar*rstar)*(f/2.0).exp(); + let v = eps*dx/(2.0*PI*p_inf.sqrt()*rstar*rstar)*(f/2.0).exp(); + *rhou = *rho*u; + *rhov = *rho*v; + *e = p/(GAMMA - 1.0) + *rho*(u*u + v*v)/2.0; + }); +} + fn pressure(gamma: Float, rho: Float, rhou: Float, rhov: Float, e: Float) -> Float { (gamma - 1.0) * (e - (rhou * rhou + rhov * rhov) / (2.0 * rho)) } @@ -542,6 +574,47 @@ pub struct BoundaryTerms<'a> { pub west: ArrayView2<'a, Float>, } +#[derive(Clone, Debug)] +pub enum BoundaryCharacteristic { + This, + // Grid(usize), + Vortex(VortexParameters), + // Vortices(Vec), +} + +#[derive(Clone, Debug)] +pub struct BoundaryCharacteristics { + north: BoundaryCharacteristic, + south: BoundaryCharacteristic, + east: BoundaryCharacteristic, + west: BoundaryCharacteristic, +} + +fn boundary_extractor<'a, SBP: SbpOperator>( + field: &'a Field, + _grid: &Grid, + bc: &BoundaryCharacteristics, +) -> BoundaryTerms<'a> { + BoundaryTerms { + north: match bc.north { + BoundaryCharacteristic::This => field.south(), + BoundaryCharacteristic::Vortex(_params) => todo!(), + }, + south: match bc.south { + BoundaryCharacteristic::This => field.north(), + BoundaryCharacteristic::Vortex(_params) => todo!(), + }, + west: match bc.west { + BoundaryCharacteristic::This => field.east(), + BoundaryCharacteristic::Vortex(_params) => todo!(), + }, + east: match bc.east { + BoundaryCharacteristic::This => field.west(), + BoundaryCharacteristic::Vortex(_params) => todo!(), + }, + } +} + #[allow(non_snake_case)] /// Boundary conditions (SAT) fn SAT_characteristics( diff --git a/sbp/tests/convergence.rs b/sbp/tests/convergence.rs index 1d5759b..83f8866 100644 --- a/sbp/tests/convergence.rs +++ b/sbp/tests/convergence.rs @@ -42,14 +42,13 @@ fn run_with_size(size: usize) -> Float { verifield.h2_err::(sys.field()) } -#[test] -fn convergence() { +fn convergence() { let sizes = [25, 35, 50, 71, 100, 150, 200]; let mut prev: Option<(usize, Float)> = None; println!("Size\tError(h2)\tq"); for size in &sizes { print!("{:3}x{:3}", size, size); - let e = run_with_size::(*size); + let e = run_with_size::(*size); print!("\t{:.10}", e); if let Some(prev) = prev.take() { let m0 = size * size; @@ -65,5 +64,14 @@ fn convergence() { println!(); prev = Some((*size, e)); } - panic!(); +} + +#[test] +fn convergence_upwind4() { + convergence::(); +} + +#[test] +fn convergence_upwind9() { + convergence::(); } diff --git a/sbp/tests/single_period.rs b/sbp/tests/single_period.rs new file mode 100644 index 0000000..b5cce08 --- /dev/null +++ b/sbp/tests/single_period.rs @@ -0,0 +1,47 @@ +#![cfg(feature = "expensive_tests")] +use ndarray::prelude::*; +use sbp::euler::*; +use sbp::Float; + +#[test] +#[ignore] +fn single_period_upwind4() { + let nx = 100; + let ny = 100; + + let x = Array1::linspace(-5.0, 5.0, nx); + let y = Array1::linspace(-5.0, 5.0, ny); + + let x = x.broadcast((ny, nx)).unwrap().to_owned(); + let y = y + .reversed_axes() + .broadcast((nx, ny)) + .unwrap() + .reversed_axes() + .to_owned(); + + let vortex_params = VortexParameters { + x0: -1.0, + y0: 0.0, + mach: 0.5, + rstar: 0.5, + eps: 1.0, + }; + + let mut sys = System::::new(x, y); + sys.vortex(0.0, vortex_params); + + let time = 10.0; + let dt = 0.2 * Float::min(1.0 / (nx - 1) as Float, 1.0 / (ny - 1) as Float); + + let nsteps = (time / dt) as usize; + for _ in 0..nsteps { + sys.advance_upwind(dt); + } + + let mut verifield = Field::new(ny, nx); + verifield.vortex(sys.x(), sys.y(), nsteps as Float * dt - 10.0, vortex_params); + + let err = verifield.h2_err::(sys.field()); + panic!("{}", err); +}