more tests

sbp/src/euler.rs

@@ -31,13 +31,14 @@ impl<SBP: SbpOperator> System<SBP> {
     }
 
     pub fn advance(&mut self, dt: Float) {
-        let rhs_trad = |k: &mut Field, y: &Field, grid: &_, wb: &mut _| {
+        let bc = BoundaryCharacteristics {
-            let boundaries = BoundaryTerms {
+            north: BoundaryCharacteristic::This,
-                north: y.south(),
+            south: BoundaryCharacteristic::This,
-                south: y.north(),
+            east: BoundaryCharacteristic::This,
-                west: y.east(),
+            west: BoundaryCharacteristic::This,
-                east: y.west(),
         };
+        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<SBP: SbpOperator> System<SBP> {
 
 impl<UO: UpwindOperator> System<UO> {
     pub fn advance_upwind(&mut self, dt: Float) {
-        let rhs_upwind = |k: &mut Field, y: &Field, grid: &_, wb: &mut _| {
+        let bc = BoundaryCharacteristics {
-            let boundaries = BoundaryTerms {
+            north: BoundaryCharacteristic::This,
-                north: y.south(),
+            south: BoundaryCharacteristic::This,
-                south: y.north(),
+            east: BoundaryCharacteristic::This,
-                west: y.east(),
+            west: BoundaryCharacteristic::This,
-                east: y.west(),
         };
+        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;
+        vortex(
-        let m = vortex_param.mach;
+            rho.into_shape((n,)).unwrap(),
-        let rstar = vortex_param.rstar;
+            rhou.into_shape((n,)).unwrap(),
-        let p_inf = 1.0 / (GAMMA * m * m);
+            rhov.into_shape((n,)).unwrap(),
-        azip!((rho in rho,
+            e.into_shape((n,)).unwrap(),
-               rhou in rhou,
+            x.into_shape((n,)).unwrap(),
-               rhov in rhov,
+            y.into_shape((n,)).unwrap(),
-               e in e,
+            t,
-               x in x,
+            vortex_param,
-               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;
-        });
     }
 }
 
@@ -333,7 +319,7 @@ fn h2_diff() {
     assert!((field0.h2_err::<super::operators::SBP8>(&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<Float>,
+    rhou: ArrayViewMut1<Float>,
+    rhov: ArrayViewMut1<Float>,
+    e: ArrayViewMut1<Float>,
+    x: ArrayView1<Float>,
+    y: ArrayView1<Float>,
+    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<VortexParameters>),
+}
+
+#[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<SBP>,
+    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<SBP: SbpOperator>(

sbp/tests/convergence.rs

@@ -42,14 +42,13 @@ fn run_with_size<SBP: sbp::operators::UpwindOperator>(size: usize) -> Float {
     verifield.h2_err::<SBP>(sys.field())
 }
 
-#[test]
+fn convergence<SBP: sbp::operators::UpwindOperator>() {
-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::<sbp::operators::Upwind4>(*size);
+        let e = run_with_size::<SBP>(*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::<sbp::operators::Upwind4>();
+}
+
+#[test]
+fn convergence_upwind9() {
+    convergence::<sbp::operators::Upwind9>();
 }

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::<sbp::operators::Upwind4>::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::<sbp::operators::Upwind4>(sys.field());
+    panic!("{}", err);
+}