make System generic over SbpOperator

src/benches/bench.rs

@@ -1,30 +1,48 @@
 use criterion::{black_box, criterion_group, criterion_main, Criterion};
-use maxwell::Universe;
+use maxwell::operators::{SbpOperator, Upwind4, UpwindOperator, SBP4};
+use maxwell::System;
 
-fn advance_system(universe: &mut Universe, n: usize) {
+fn advance_system<SBP: SbpOperator>(universe: &mut System<SBP>, n: usize) {
     for _ in 0..n {
         universe.advance(0.01);
     }
 }
 
+fn advance_system_upwind<UO: UpwindOperator>(universe: &mut System<UO>, n: usize) {
+    for _ in 0..n {
+        universe.advance_upwind(0.01);
+    }
+}
+
 fn performance_benchmark(c: &mut Criterion) {
     let mut group = c.benchmark_group("System");
-    group.sample_size(15);
+    group.sample_size(25);
 
     let w = 40;
     let h = 26;
     let x = ndarray::Array2::from_shape_fn((h, w), |(_, i)| i as f32 / (w - 1) as f32);
     let y = ndarray::Array2::from_shape_fn((h, w), |(j, _)| j as f32 / (h - 1) as f32);
-    let mut universe = Universe::new(
-        w as u32,
-        h as u32,
-        x.as_slice().unwrap(),
-        y.as_slice().unwrap(),
-    );
 
+    let mut universe = System::<Upwind4>::new(w, h, x.as_slice().unwrap(), y.as_slice().unwrap());
     group.bench_function("advance", |b| {
         b.iter(|| {
-            universe.init(0.0, 0.0);
+            universe.set_gaussian(0.5, 0.5);
+            advance_system(&mut universe, black_box(20))
+        })
+    });
+
+    let mut universe = System::<Upwind4>::new(w, h, x.as_slice().unwrap(), y.as_slice().unwrap());
+    group.bench_function("advance_upwind", |b| {
+        b.iter(|| {
+            universe.set_gaussian(0.5, 0.5);
+            advance_system_upwind(&mut universe, black_box(20))
+        })
+    });
+
+    let mut universe = System::<SBP4>::new(w, h, x.as_slice().unwrap(), y.as_slice().unwrap());
+    group.bench_function("advance_trad4", |b| {
+        b.iter(|| {
+            universe.set_gaussian(0.5, 0.5);
             advance_system(&mut universe, black_box(20))
         })
     });

src/grid.rs

@@ -16,9 +16,7 @@ where
 }
 
 impl<SBP: super::operators::SbpOperator> Grid<SBP> {
-    pub fn new(width: u32, height: u32, x: &[f32], y: &[f32]) -> Result<Self, ndarray::ShapeError> {
-        let nx = width as usize;
-        let ny = height as usize;
+    pub fn new(nx: usize, ny: usize, x: &[f32], y: &[f32]) -> Result<Self, ndarray::ShapeError> {
         let x = Array2::from_shape_vec((ny, nx), x.to_vec())?;
         let y = Array2::from_shape_vec((ny, nx), y.to_vec())?;
 

src/lib.rs

@@ -2,7 +2,7 @@ use wasm_bindgen::prelude::*;
 
 mod grid;
 mod maxwell;
-mod operators;
+pub mod operators;
 pub use crate::maxwell::{Field, WorkBuffers};
 pub(crate) use grid::Grid;
 
@@ -17,33 +17,62 @@ pub fn set_panic_hook() {
 }
 
 #[wasm_bindgen]
-pub struct Universe {
-    sys: (Field, Field),
-    wb: WorkBuffers,
-    grid: Grid<operators::Upwind4>,
-}
+pub struct Universe(System<operators::Upwind4>);
 
 #[wasm_bindgen]
 impl Universe {
     #[wasm_bindgen(constructor)]
-    pub fn new(width: u32, height: u32, x: &[f32], y: &[f32]) -> Self {
-        assert_eq!((width * height) as usize, x.len());
-        assert_eq!((width * height) as usize, y.len());
+    pub fn new(width: usize, height: usize, x: &[f32], y: &[f32]) -> Self {
+        Self(System::new(width as usize, height as usize, x, y))
+    }
+
+    pub fn init(&mut self, x0: f32, y0: f32) {
+        self.0.set_gaussian(x0, y0);
+    }
+
+    pub fn advance(&mut self, dt: f32) {
+        self.0.advance(dt)
+    }
+
+    pub fn advance_upwind(&mut self, dt: f32) {
+        self.0.advance_upwind(dt)
+    }
+
+    pub fn get_ex_ptr(&self) -> *const u8 {
+        self.0.sys.0.ex().as_ptr() as *const u8
+    }
+
+    pub fn get_ey_ptr(&self) -> *const u8 {
+        self.0.sys.0.ey().as_ptr() as *const u8
+    }
+
+    pub fn get_hz_ptr(&self) -> *const u8 {
+        self.0.sys.0.hz().as_ptr() as *const u8
+    }
+}
+
+pub struct System<SBP: operators::SbpOperator> {
+    sys: (Field, Field),
+    wb: WorkBuffers,
+    grid: Grid<SBP>,
+}
+
+impl<SBP: operators::SbpOperator> System<SBP> {
+    pub fn new(width: usize, height: usize, x: &[f32], y: &[f32]) -> Self {
+        assert_eq!((width * height), x.len());
+        assert_eq!((width * height), y.len());
 
         let grid = Grid::new(width, height, x, y).expect(
             "Could not create grid. Different number of elements compared to width*height?",
         );
         Self {
-            sys: (
-                Field::new(width as usize, height as usize),
-                Field::new(width as usize, height as usize),
-            ),
+            sys: (Field::new(width, height), Field::new(width, height)),
             grid,
-            wb: WorkBuffers::new(width as usize, height as usize),
+            wb: WorkBuffers::new(width, height),
         }
     }
 
-    fn set_gaussian(&mut self, x0: f32, y0: f32) {
+    pub fn set_gaussian(&mut self, x0: f32, y0: f32) {
         let (ex, hz, ey) = self.sys.0.components_mut();
         ndarray::azip!(
             (ex in ex, hz in hz, ey in ey,
@@ -55,22 +84,6 @@ impl Universe {
         });
     }
 
-    pub fn init(&mut self, x0: f32, y0: f32) {
-        self.set_gaussian(x0, y0);
-    }
-
-    /// Using artificial dissipation with the upwind operator
-    pub fn advance_upwind(&mut self, dt: f32) {
-        maxwell::advance_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);
-    }
-
     pub fn advance(&mut self, dt: f32) {
         maxwell::advance(
             &self.sys.0,
@@ -81,17 +94,19 @@ impl Universe {
         );
         std::mem::swap(&mut self.sys.0, &mut self.sys.1);
     }
-
-    pub fn get_ex_ptr(&self) -> *const u8 {
-        self.sys.0.ex().as_ptr() as *const u8
 }
 
-    pub fn get_ey_ptr(&self) -> *const u8 {
-        self.sys.0.ey().as_ptr() as *const u8
-    }
-
-    pub fn get_hz_ptr(&self) -> *const u8 {
-        self.sys.0.hz().as_ptr() as *const u8
+impl<UO: operators::UpwindOperator> System<UO> {
+    /// Using artificial dissipation with the upwind operator
+    pub fn advance_upwind(&mut self, dt: f32) {
+        maxwell::advance_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);
     }
 }