allow several vortices

multigrid/src/main.rs

@@ -70,9 +70,9 @@ impl System {
         }
     }
 
-    fn vortex(&mut self, t: Float, vortex_params: euler::VortexParameters) {
+    fn vortex(&mut self, t: Float, vortex_params: &euler::VortexParameters) {
         for (f, g) in self.fnow.iter_mut().zip(&self.grids) {
-            f.vortex(g.x(), g.y(), t, vortex_params);
+            f.vortex(g.x(), g.y(), t, &vortex_params);
         }
     }
 
@@ -161,12 +161,12 @@ fn main() {
     let json = json::parse(&filecontents).unwrap();
 
     let vortexparams = json_to_vortex(json["vortex"].clone());
-    let (names, grids, bt, operators) = json_to_grids(json["grids"].clone(), vortexparams);
+    let (names, grids, bt, operators) = json_to_grids(json["grids"].clone(), vortexparams.clone());
 
     let integration_time: Float = json["integration_time"].as_number().unwrap().into();
 
     let mut sys = System::new(grids, bt, operators);
-    sys.vortex(0.0, vortexparams);
+    sys.vortex(0.0, &vortexparams);
 
     let max_n = {
         let max_nx = sys.grids.iter().map(|g| g.nx()).max().unwrap();
@@ -233,7 +233,7 @@ fn main() {
         let mut e = 0.0;
         for ((fmod, grid), op) in sys.fnow.iter().zip(&sys.grids).zip(&sys.operators) {
             let mut fvort = fmod.clone();
-            fvort.vortex(grid.x(), grid.y(), time, vortexparams);
+            fvort.vortex(grid.x(), grid.y(), time, &vortexparams);
             let sbpop: &dyn SbpOperator2d = op.as_ref().either(|op| &**op, |uo| uo.as_sbp());
             e += fmod.h2_err(&fvort, sbpop);
         }

multigrid/src/parsing.rs

@@ -110,7 +110,7 @@ pub fn json_to_grids(
     let determine_bc = |dir: Option<&str>| match dir {
         Some(dir) => {
             if dir == "vortex" {
-                sbp::euler::BoundaryCharacteristic::Vortex(vortexparams)
+                sbp::euler::BoundaryCharacteristic::Vortex(vortexparams.clone())
             } else if let Some(grid) = dir.strip_prefix("interpolate:") {
                 use sbp::operators::*;
                 let (grid, int_op) = if let Some(rest) = grid.strip_prefix("4:") {
@@ -363,11 +363,13 @@ pub fn json_to_vortex(mut json: JsonValue) -> super::euler::VortexParameters {
         }
     }
 
-    super::euler::VortexParameters {
+    let vortice = super::euler::Vortice { x0, y0, eps, rstar };
-        x0,
+
-        y0,
+    let vortices = {
-        mach,
+        let mut a = super::euler::ArrayVec::new();
-        rstar,
+        a.push(vortice);
-        eps,
+        a
-    }
+    };
+
+    super::euler::VortexParameters { vortices, mach }
 }

sbp/Cargo.toml

@@ -9,6 +9,7 @@ ndarray = { version = "0.13.1", features = ["approx"] }
 approx = "0.3.2"
 packed_simd = "0.3.3"
 rayon = { version = "1.3.0", optional = true }
+arrayvec = "0.5.1"
 
 [features]
 # Internal feature flag to gate the expensive tests

sbp/src/euler.rs

@@ -3,6 +3,7 @@ use super::integrate;
 use super::operators::{InterpolationOperator, SbpOperator2d, UpwindOperator2d};
 use super::utils::Direction;
 use super::Float;
+pub use arrayvec::ArrayVec;
 use ndarray::azip;
 use ndarray::prelude::*;
 
@@ -70,23 +71,29 @@ impl<SBP: SbpOperator2d> System<SBP> {
     pub fn vortex(&mut self, t: Float, vortex_parameters: VortexParameters) {
         self.sys
             .0
-            .vortex(self.grid.0.x(), self.grid.0.y(), t, vortex_parameters);
+            .vortex(self.grid.0.x(), self.grid.0.y(), t, &vortex_parameters);
     }
 
     #[allow(clippy::many_single_char_names)]
     pub fn init_with_vortex(&mut self, x0: Float, y0: Float) {
         // Should parametrise such that we have radius, drop in pressure at center, etc
         let vortex_parameters = VortexParameters {
-            x0,
+            vortices: {
-            y0,
+                let mut v = ArrayVec::new();
-            rstar: 1.0,
+                v.push(Vortice {
-            eps: 3.0,
+                    x0,
+                    y0,
+                    rstar: 1.0,
+                    eps: 3.0,
+                });
+                v
+            },
             mach: 0.5,
         };
 
         self.sys
             .0
-            .vortex(self.grid.0.x(), self.grid.0.y(), 0.0, vortex_parameters)
+            .vortex(self.grid.0.x(), self.grid.0.y(), 0.0, &vortex_parameters)
     }
 
     pub fn field(&self) -> &Field {
@@ -297,7 +304,7 @@ impl Field {
         x: ArrayView2<Float>,
         y: ArrayView2<Float>,
         time: Float,
-        vortex_param: VortexParameters,
+        vortex_param: &VortexParameters,
     ) {
         assert_eq!(x.shape(), y.shape());
         assert_eq!(x.shape()[1], self.nx());
@@ -314,7 +321,7 @@ impl Field {
             x.into_shape((n,)).unwrap(),
             y.into_shape((n,)).unwrap(),
             time,
-            vortex_param,
+            &vortex_param,
         )
     }
 }
@@ -395,11 +402,16 @@ fn h2_diff() {
 }
 
 #[derive(Copy, Clone, Debug)]
-pub struct VortexParameters {
+pub struct Vortice {
     pub x0: Float,
     pub y0: Float,
     pub rstar: Float,
     pub eps: Float,
+}
+
+#[derive(Clone, Debug)]
+pub struct VortexParameters {
+    pub vortices: ArrayVec<[Vortice; 5]>,
     pub mach: Float,
 }
 
@@ -412,7 +424,7 @@ pub fn vortex(
     x: ArrayView1<Float>,
     y: ArrayView1<Float>,
     time: Float,
-    vortex_param: VortexParameters,
+    vortex_param: &VortexParameters,
 ) {
     assert_eq!(rho.len(), rhou.len());
     assert_eq!(rho.len(), rhov.len());
@@ -421,10 +433,14 @@ pub fn vortex(
     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);
+
+    let rho_inf: Float = 1.0;
+    let u_inf: Float = 1.0;
+    let v_inf: Float = 0.0;
+    let e_inf = p_inf / (GAMMA - 1.0) + rho_inf * (u_inf.powi(2) + v_inf.powi(2)) / 2.0;
+
     azip!((rho in rho,
            rhou in rhou,
            rhov in rhov,
@@ -432,21 +448,36 @@ pub fn vortex(
            x in x,
            y in y)
     {
-        use crate::consts::PI;
 
-        let dx = (x - vortex_param.x0) - time;
+        *rho = rho_inf;
-        let dy = y - vortex_param.y0;
+        *rhou = rho_inf*u_inf;
-        let f = (1.0 - (dx*dx + dy*dy))/(rstar*rstar);
+        *rhou = rho_inf*v_inf;
+        *e = e_inf;
 
-        *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));
+        for vortice in vortex_param.vortices.iter() {
-        assert!(*rho > 0.0);
+            use crate::consts::PI;
-        let p = Float::powf(*rho, GAMMA)*p_inf;
+
-        assert!(p > 0.0);
+            let rstar = vortice.rstar;
-        let u = 1.0 - eps*dy/(2.0*PI*p_inf.sqrt()*rstar*rstar)*(f/2.0).exp();
+            let eps = vortice.eps;
-        let v =       eps*dx/(2.0*PI*p_inf.sqrt()*rstar*rstar)*(f/2.0).exp();
+
-        *rhou = *rho*u;
+            let dx = (x - vortice.x0) - time;
-        *rhov = *rho*v;
+            let dy = y - vortice.y0;
-        *e = p/(GAMMA - 1.0) + *rho*(u*u + v*v)/2.0;
+            let f = (1.0 - (dx*dx + dy*dy))/(rstar*rstar);
+
+            let rho_vortice = 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));
+            let p = Float::powf(rho_vortice, GAMMA)*p_inf;
+            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();
+
+            assert!(rho_vortice > 0.0);
+            assert!(p > 0.0);
+
+            *rho += rho_vortice - rho_inf;
+            assert!(*rho > 0.0);
+            *rhou += rho_vortice*u - rho_inf*u_inf;
+            *rhov += rho_vortice*v - rho_inf*v_inf;
+            *e += (p/(GAMMA - 1.0) + rho_vortice*(u*u + v*v)/2.0) - e_inf;
+        }
     });
 }
 
@@ -663,28 +694,28 @@ fn boundary_extractor<'a>(
     bc: &BoundaryCharacteristics,
 ) -> BoundaryTerms<'a> {
     BoundaryTerms {
-        north: match bc.north {
+        north: match &bc.north {
             BoundaryCharacteristic::This => field.south(),
             BoundaryCharacteristic::Vortex(_params) => todo!(),
             BoundaryCharacteristic::Grid(_)
             | BoundaryCharacteristic::Interpolate(_, _)
             | BoundaryCharacteristic::MultiGrid(_) => panic!("Only working on self grid"),
         },
-        south: match bc.south {
+        south: match &bc.south {
             BoundaryCharacteristic::This => field.north(),
             BoundaryCharacteristic::Vortex(_params) => todo!(),
             BoundaryCharacteristic::Grid(_)
             | BoundaryCharacteristic::Interpolate(_, _)
             | BoundaryCharacteristic::MultiGrid(_) => panic!("Only working on self grid"),
         },
-        west: match bc.west {
+        west: match &bc.west {
             BoundaryCharacteristic::This => field.east(),
             BoundaryCharacteristic::Vortex(_params) => todo!(),
             BoundaryCharacteristic::Grid(_)
             | BoundaryCharacteristic::Interpolate(_, _)
             | BoundaryCharacteristic::MultiGrid(_) => panic!("Only working on self grid"),
         },
-        east: match bc.east {
+        east: match &bc.east {
             BoundaryCharacteristic::This => field.west(),
             BoundaryCharacteristic::Vortex(_params) => todo!(),
             BoundaryCharacteristic::Grid(_)
@@ -708,7 +739,7 @@ fn boundary_extract<'a>(
         BoundaryCharacteristic::Grid(g) => seldir(&fields[*g]),
         BoundaryCharacteristic::Vortex(v) => {
             let field = eb.unwrap();
-            vortexify(field.view_mut(), grid, *v, time);
+            vortexify(field.view_mut(), grid, v, time);
             field.view()
         }
         BoundaryCharacteristic::Interpolate(g, operator) => {
@@ -854,7 +885,7 @@ impl BoundaryStorage {
 fn vortexify(
     mut field: ndarray::ArrayViewMut2<Float>,
     yx: (ndarray::ArrayView1<Float>, ndarray::ArrayView1<Float>),
-    vparams: VortexParameters,
+    vparams: &VortexParameters,
     time: Float,
 ) {
     let mut fiter = field.outer_iter_mut();
@@ -865,7 +896,7 @@ fn vortexify(
         fiter.next().unwrap(),
     );
     let (y, x) = yx;
-    vortex(rho, rhou, rhov, e, x, y, time, vparams);
+    vortex(rho, rhou, rhov, e, x, y, time, &vparams);
 }
 
 #[allow(non_snake_case)]