move from lib to euler

benches/euler.rs

@@ -1,14 +1,14 @@
 use criterion::{black_box, criterion_group, criterion_main, Criterion};
+use sbp::euler::System;
 use sbp::operators::{SbpOperator, Upwind4, UpwindOperator, SBP4};
-use sbp::EulerSystem;
 
-fn advance_system<SBP: SbpOperator>(universe: &mut EulerSystem<SBP>, n: usize) {
+fn advance_system<SBP: SbpOperator>(universe: &mut System<SBP>, n: usize) {
     for _ in 0..n {
         universe.advance(1.0 / 40.0 * 0.2);
     }
 }
 
-fn advance_system_upwind<UO: UpwindOperator>(universe: &mut EulerSystem<UO>, n: usize) {
+fn advance_system_upwind<UO: UpwindOperator>(universe: &mut System<UO>, n: usize) {
     for _ in 0..n {
         universe.advance_upwind(1.0 / 40.0 * 0.2);
     }
@@ -25,26 +25,26 @@ fn performance_benchmark(c: &mut Criterion) {
     let y = ndarray::Array1::linspace(-10.0, 10.0, h);
     let y = y.broadcast((w, h)).unwrap().reversed_axes();
 
-    let mut universe = EulerSystem::<Upwind4>::new(x.into_owned(), y.into_owned());
+    let mut universe = System::<Upwind4>::new(x.into_owned(), y.into_owned());
     group.bench_function("advance", |b| {
         b.iter(|| {
-            universe.vortex(0.0, 0.0);
+            universe.init_with_vortex(0.0, 0.0);
             advance_system(&mut universe, black_box(20))
         })
     });
 
-    let mut universe = EulerSystem::<Upwind4>::new(x.into_owned(), y.into_owned());
+    let mut universe = System::<Upwind4>::new(x.into_owned(), y.into_owned());
     group.bench_function("advance_upwind", |b| {
         b.iter(|| {
-            universe.vortex(0.0, 0.0);
+            universe.init_with_vortex(0.0, 0.0);
             advance_system_upwind(&mut universe, black_box(20))
         })
     });
 
-    let mut universe = EulerSystem::<SBP4>::new(x.into_owned(), y.into_owned());
+    let mut universe = System::<SBP4>::new(x.into_owned(), y.into_owned());
     group.bench_function("advance_trad4", |b| {
         b.iter(|| {
-            universe.vortex(0.0, 0.0);
+            universe.init_with_vortex(0.0, 0.0);
             advance_system(&mut universe, black_box(20))
         })
     });

src/euler.rs

@@ -5,6 +5,106 @@ use ndarray::{azip, Zip};
 
 pub const GAMMA: f32 = 1.4;
 
+// A collection of buffers that allows one to efficiently
+// move to the next state
+pub struct System<SBP: SbpOperator> {
+    sys: (Field, Field),
+    wb: WorkBuffers,
+    grid: Grid<SBP>,
+}
+
+impl<SBP: SbpOperator> System<SBP> {
+    pub fn new(x: ndarray::Array2<f32>, y: ndarray::Array2<f32>) -> Self {
+        let grid = Grid::new(x, y).expect(
+            "Could not create grid. Different number of elements compared to width*height?",
+        );
+        let nx = grid.nx();
+        let ny = grid.ny();
+        Self {
+            sys: (Field::new(ny, nx), Field::new(ny, nx)),
+            grid,
+            wb: WorkBuffers::new(ny, nx),
+        }
+    }
+
+    pub fn advance(&mut self, dt: f32) {
+        advance(
+            RHS_trad,
+            &self.sys.0,
+            &mut self.sys.1,
+            dt,
+            &self.grid,
+            Some(&mut self.wb),
+        );
+        std::mem::swap(&mut self.sys.0, &mut self.sys.1);
+    }
+
+    pub fn init_with_vortex(&mut self, x0: f32, y0: f32) {
+        // Should parametrise such that we have radius, drop in pressure at center, etc
+        let rstar = 1.0;
+        let eps = 3.0;
+        #[allow(non_snake_case)]
+        let M = 0.5;
+
+        let p_inf = 1.0 / (GAMMA * M * M);
+        let t = 0.0;
+
+        let nx = self.grid.nx();
+        let ny = self.grid.ny();
+
+        for j in 0..ny {
+            for i in 0..nx {
+                let x = self.grid.x[(j, i)];
+                let y = self.grid.y[(j, i)];
+
+                let dx = (x - x0) - t;
+                let dy = y - y0;
+                let f = (1.0 - (dx * dx + dy * dy)) / (rstar * rstar);
+
+                use std::f32::consts::PI;
+                let u =
+                    1.0 - eps * dy / (2.0 * PI * p_inf.sqrt() * rstar * rstar) * (f / 2.0).exp();
+                let v =
+                    0.0 + eps * dx / (2.0 * PI * p_inf.sqrt() * rstar * rstar) * (f / 2.0).exp();
+                let rho = f32::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 = p_inf * rho.powf(GAMMA);
+                assert!(p > 0.0);
+                let e = p / (GAMMA - 1.0) + rho * (u * u + v * v) / 2.0;
+                assert!(e > 0.0);
+
+                self.sys.0[(0, j, i)] = rho;
+                self.sys.0[(1, j, i)] = rho * u;
+                self.sys.0[(2, j, i)] = rho * v;
+                self.sys.0[(3, j, i)] = e;
+            }
+        }
+    }
+
+    pub fn field(&self) -> &Field {
+        &self.sys.0
+    }
+}
+
+impl<UO: UpwindOperator> System<UO> {
+    pub fn advance_upwind(&mut self, dt: f32) {
+        advance(
+            RHS_upwind,
+            &self.sys.0,
+            &mut self.sys.1,
+            dt,
+            &self.grid,
+            Some(&mut self.wb),
+        );
+        std::mem::swap(&mut self.sys.0, &mut self.sys.1);
+    }
+}
+
 #[derive(Clone, Debug)]
 /// A 4 x ny x nx array
 pub struct Field(pub(crate) Array3<f32>);

src/grid.rs

@@ -1,4 +1,6 @@
 use ndarray::Array2;
+
+#[derive(Debug, Clone)]
 pub struct Grid<SBP>
 where
     SBP: super::operators::SbpOperator,

src/lib.rs

@@ -1,6 +1,6 @@
 use wasm_bindgen::prelude::*;
 
-mod euler;
+pub mod euler;
 mod grid;
 mod maxwell;
 pub mod operators;
@@ -121,18 +121,12 @@ fn gaussian(x: f32, x0: f32, y: f32, y0: f32) -> f32 {
     1.0 / (2.0 * f32::consts::PI * sigma * sigma) * (-(x * x + y * y) / (2.0 * sigma * sigma)).exp()
 }
 
-pub struct EulerSystem<SBP: operators::SbpOperator> {
-    sys: (euler::Field, euler::Field),
-    wb: euler::WorkBuffers,
-    grid: Grid<SBP>,
-}
-
 #[wasm_bindgen]
-pub struct EulerUniverse(EulerSystem<operators::Upwind4>);
+pub struct EulerUniverse(euler::System<operators::Upwind4>);
 
 impl EulerUniverse {
     pub fn new(x: ndarray::Array2<f32>, y: ndarray::Array2<f32>) -> Self {
-        Self(EulerSystem::new(x, y))
+        Self(euler::System::new(x, y))
     }
 }
 
@@ -142,11 +136,11 @@ impl EulerUniverse {
     pub fn new_with_slice(height: usize, width: usize, x: &[f32], y: &[f32]) -> Self {
         let x = ndarray::Array2::from_shape_vec((height, width), x.to_vec()).unwrap();
         let y = ndarray::Array2::from_shape_vec((height, width), y.to_vec()).unwrap();
-        Self(EulerSystem::new(x, y))
+        Self(euler::System::new(x, y))
     }
 
     pub fn init(&mut self, x0: f32, y0: f32) {
-        self.0.vortex(x0, y0)
+        self.0.init_with_vortex(x0, y0)
     }
 
     pub fn advance(&mut self, dt: f32) {
@@ -158,105 +152,16 @@ impl EulerUniverse {
     }
 
     pub fn get_rho_ptr(&self) -> *const u8 {
-        self.0.sys.0.rho().as_ptr() as *const u8
+        self.0.field().rho().as_ptr() as *const u8
     }
     pub fn get_rhou_ptr(&self) -> *const u8 {
-        self.0.sys.0.rhou().as_ptr() as *const u8
+        self.0.field().rhou().as_ptr() as *const u8
     }
     pub fn get_rhov_ptr(&self) -> *const u8 {
-        self.0.sys.0.rhov().as_ptr() as *const u8
+        self.0.field().rhov().as_ptr() as *const u8
     }
     pub fn get_e_ptr(&self) -> *const u8 {
-        self.0.sys.0.e().as_ptr() as *const u8
+        self.0.field().e().as_ptr() as *const u8
-    }
-}
-
-impl<SBP: operators::SbpOperator> EulerSystem<SBP> {
-    pub fn new(x: ndarray::Array2<f32>, y: ndarray::Array2<f32>) -> Self {
-        let grid = Grid::new(x, y).expect(
-            "Could not create grid. Different number of elements compared to width*height?",
-        );
-        let nx = grid.nx();
-        let ny = grid.ny();
-        Self {
-            sys: (euler::Field::new(ny, nx), euler::Field::new(ny, nx)),
-            grid,
-            wb: euler::WorkBuffers::new(ny, nx),
-        }
-    }
-
-    pub fn advance(&mut self, dt: f32) {
-        euler::advance(
-            euler::RHS_trad,
-            &self.sys.0,
-            &mut self.sys.1,
-            dt,
-            &self.grid,
-            Some(&mut self.wb),
-        );
-        std::mem::swap(&mut self.sys.0, &mut self.sys.1);
-    }
-
-    pub fn vortex(&mut self, x0: f32, y0: f32) {
-        // Should parametrise such that we have radius, drop in pressure at center, etc
-        let rstar = 1.0;
-        let eps = 3.0;
-        #[allow(non_snake_case)]
-        let M = 0.5;
-
-        let p_inf = 1.0 / (euler::GAMMA * M * M);
-        let t = 0.0;
-
-        let nx = self.grid.nx();
-        let ny = self.grid.ny();
-
-        for j in 0..ny {
-            for i in 0..nx {
-                let x = self.grid.x[(j, i)];
-                let y = self.grid.y[(j, i)];
-
-                let dx = (x - x0) - t;
-                let dy = y - y0;
-                let f = (1.0 - (dx * dx + dy * dy)) / (rstar * rstar);
-
-                use euler::GAMMA;
-                use std::f32::consts::PI;
-                let u =
-                    1.0 - eps * dy / (2.0 * PI * p_inf.sqrt() * rstar * rstar) * (f / 2.0).exp();
-                let v =
-                    0.0 + eps * dx / (2.0 * PI * p_inf.sqrt() * rstar * rstar) * (f / 2.0).exp();
-                let rho = f32::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 = p_inf * rho.powf(GAMMA);
-                assert!(p > 0.0);
-                let e = p / (GAMMA - 1.0) + rho * (u * u + v * v) / 2.0;
-                assert!(e > 0.0);
-
-                self.sys.0[(0, j, i)] = rho;
-                self.sys.0[(1, j, i)] = rho * u;
-                self.sys.0[(2, j, i)] = rho * v;
-                self.sys.0[(3, j, i)] = e;
-            }
-        }
-    }
-}
-
-impl<SBP: operators::UpwindOperator> EulerSystem<SBP> {
-    pub fn advance_upwind(&mut self, dt: f32) {
-        euler::advance(
-            euler::RHS_upwind,
-            &self.sys.0,
-            &mut self.sys.1,
-            dt,
-            &self.grid,
-            Some(&mut self.wb),
-        );
-        std::mem::swap(&mut self.sys.0, &mut self.sys.1);
     }
 }