convergence test

sbp/src/euler.rs

@@ -51,6 +51,12 @@ impl<SBP: SbpOperator> System<SBP> {
         std::mem::swap(&mut self.sys.0, &mut self.sys.1);
     }
 
+    pub fn vortex(&mut self, t: f32, vortex_parameters: VortexParameters) {
+        self.sys
+            .0
+            .vortex(self.grid.x.view(), self.grid.y.view(), t, vortex_parameters);
+    }
+
     #[allow(clippy::many_single_char_names)]
     pub fn init_with_vortex(&mut self, x0: f32, y0: f32) {
         // Should parametrise such that we have radius, drop in pressure at center, etc
@@ -73,6 +79,13 @@ impl<SBP: SbpOperator> System<SBP> {
     pub fn field(&self) -> &Field {
         &self.sys.0
     }
+
+    pub fn x(&self) -> ArrayView2<f32> {
+        self.grid.x.view()
+    }
+    pub fn y(&self) -> ArrayView2<f32> {
+        self.grid.y.view()
+    }
 }
 
 impl<UO: UpwindOperator> System<UO> {
@@ -214,7 +227,7 @@ impl Field {
         self.slice_mut(s![.., .., 0])
     }
 
-    fn vortex(
+    pub fn vortex(
         &mut self,
         x: ArrayView2<f32>,
         y: ArrayView2<f32>,
@@ -257,7 +270,7 @@ impl Field {
 }
 
 impl Field {
-    fn err_diff<SBP: SbpOperator>(&self, other: &Self) -> f32 {
+    pub fn err_diff<SBP: SbpOperator>(&self, other: &Self) -> f32 {
         assert_eq!(self.nx(), other.nx());
         assert_eq!(self.ny(), other.ny());
 
@@ -318,11 +331,11 @@ fn h2_diff() {
 
 #[derive(Copy, Clone)]
 pub struct VortexParameters {
-    x0: f32,
+    pub x0: f32,
-    y0: f32,
+    pub y0: f32,
-    rstar: f32,
+    pub rstar: f32,
-    eps: f32,
+    pub eps: f32,
-    mach: f32,
+    pub mach: f32,
 }
 
 fn pressure(gamma: f32, rho: f32, rhou: f32, rhov: f32, e: f32) -> f32 {

sbp/tests/convergence.rs

@@ -0,0 +1,53 @@
+use ndarray::prelude::*;
+use sbp::euler::*;
+
+fn run_with_size<SBP: sbp::operators::UpwindOperator>(size: usize) -> f32 {
+    let nx = 4 * size;
+    let ny = 4 * size;
+    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>::new(x, y);
+    sys.vortex(0.0, vortex_params);
+
+    let time = 0.1;
+    let dt = 0.2 * f32::min(1.0 / (nx - 1) as f32, 1.0 / (ny - 1) as f32);
+
+    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 f32 * dt, vortex_params);
+
+    verifield.err_diff::<SBP>(sys.field())
+}
+
+#[test]
+fn test() {
+    let sizes = [25, 35, 50, 71, 100];
+    let mut errs = Vec::with_capacity(sizes.len());
+    for size in &sizes {
+        println!("Size: {}", size);
+        let e = run_with_size::<sbp::operators::Upwind9>(*size);
+        errs.push(e);
+    }
+    panic!("{:?}", errs);
+}