From a1acf00c5a257665b0be102a298eb2dba328dfec Mon Sep 17 00:00:00 2001
From: Magnus Ulimoen <flymagnus@gmail.com>
Date: Sat, 2 May 2020 00:22:59 +0200
Subject: [PATCH] split euler and maxwell to separate crates

---
 Cargo.toml                                    |  4 +-
 euler/Cargo.toml                              | 22 +++++
 .../euler.rs => euler/benches/bench.rs        |  0
 sbp/src/euler.rs => euler/src/lib.rs          | 81 ++++++++++---------
 {sbp => euler}/tests/convergence.rs           |  2 +-
 maxwell/Cargo.toml                            | 17 ++++
 .../maxwell.rs => maxwell/benches/bench.rs    |  2 +-
 sbp/src/maxwell.rs => maxwell/src/lib.rs      | 74 ++++++++---------
 multigrid/Cargo.toml                          |  3 +-
 multigrid/src/main.rs                         |  9 ++-
 multigrid/src/parsing.rs                      | 28 +++----
 sbp/Cargo.toml                                | 12 ---
 sbp/src/grid.rs                               | 20 ++++-
 sbp/src/lib.rs                                |  8 +-
 sbp/src/utils.rs                              | 27 +++++++
 webfront/Cargo.toml                           |  2 +
 webfront/lib.rs                               |  4 +-
 17 files changed, 202 insertions(+), 113 deletions(-)
 create mode 100644 euler/Cargo.toml
 rename sbp/benches/euler.rs => euler/benches/bench.rs (100%)
 rename sbp/src/euler.rs => euler/src/lib.rs (94%)
 rename {sbp => euler}/tests/convergence.rs (99%)
 create mode 100644 maxwell/Cargo.toml
 rename sbp/benches/maxwell.rs => maxwell/benches/bench.rs (98%)
 rename sbp/src/maxwell.rs => maxwell/src/lib.rs (88%)

diff --git a/Cargo.toml b/Cargo.toml
index d22aabb..18e79a0 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -10,6 +10,8 @@ members = [
     "sbp",
     "webfront",
     "multigrid",
+    "euler",
+    "maxwell",
 ]
 
 [profile.bench]
@@ -17,4 +19,4 @@ debug = true
 
 [patch]
 [patch.crates-io]
-hdf5 = { git = "https://github.com/mulimoen/hdf5-rust.git", branch = "feature/resizable_idx" }
+hdf5 = { git = "https://github.com/mulimoen/hdf5-rust.git", branch = "master" }
diff --git a/euler/Cargo.toml b/euler/Cargo.toml
new file mode 100644
index 0000000..c68649c
--- /dev/null
+++ b/euler/Cargo.toml
@@ -0,0 +1,22 @@
+[package]
+name = "euler"
+version = "0.1.0"
+authors = ["Magnus Ulimoen <flymagnus@gmail.com>"]
+edition = "2018"
+
+[features]
+# Internal feature flag to gate the expensive tests
+# which should be run only in release builds
+expensive_tests = []
+
+[dependencies]
+ndarray = "0.13.1"
+sbp = { path = "../sbp" }
+arrayvec = "0.5.1"
+
+[dev-dependencies]
+criterion = "0.3.1"
+
+[[bench]]
+name = "bench"
+harness = false
diff --git a/sbp/benches/euler.rs b/euler/benches/bench.rs
similarity index 100%
rename from sbp/benches/euler.rs
rename to euler/benches/bench.rs
diff --git a/sbp/src/euler.rs b/euler/src/lib.rs
similarity index 94%
rename from sbp/src/euler.rs
rename to euler/src/lib.rs
index 207968c..5f01947 100644
--- a/sbp/src/euler.rs
+++ b/euler/src/lib.rs
@@ -1,11 +1,11 @@
-use super::grid::{Grid, Metrics};
-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::*;
+use sbp::grid::{Grid, Metrics};
+use sbp::integrate;
+use sbp::operators::{InterpolationOperator, SbpOperator2d, UpwindOperator2d};
+use sbp::utils::Direction;
+use sbp::Float;
 
 pub const GAMMA: Float = 1.4;
 
@@ -101,10 +101,10 @@ impl<SBP: SbpOperator2d> System<SBP> {
     }
 
     pub fn x(&self) -> ArrayView2<Float> {
-        self.grid.0.x.view()
+        self.grid.0.x()
     }
     pub fn y(&self) -> ArrayView2<Float> {
-        self.grid.0.y.view()
+        self.grid.0.y()
     }
 
     pub fn nx(&self) -> usize {
@@ -393,7 +393,7 @@ fn h2_diff() {
     }
     let field1 = Field::new(20, 21);
 
-    use super::operators::{Upwind4, Upwind9, SBP4, SBP8};
+    use sbp::operators::{Upwind4, Upwind9, SBP4, SBP8};
 
     assert!((field0.h2_err(&field1, &Upwind4).powi(2) - 4.0).abs() < 1e-3);
     assert!((field0.h2_err(&field1, &Upwind9).powi(2) - 4.0).abs() < 1e-3);
@@ -460,7 +460,7 @@ pub fn vortex(
                 return;
             },
             Some(vortice) => {
-                use crate::consts::PI;
+                use sbp::consts::PI;
 
                 let rstar = vortice.rstar;
                 let eps = vortice.eps;
@@ -485,7 +485,7 @@ pub fn vortex(
         }
 
         for vortice in iterator {
-            use crate::consts::PI;
+            use sbp::consts::PI;
 
             let rstar = vortice.rstar;
             let eps = vortice.eps;
@@ -543,7 +543,7 @@ pub fn RHS_trad(
     azip!((out in &mut k.0,
                     eflux in &dE.0,
                     fflux in &dF.0,
-                    detj in &metrics.detj.broadcast((4, y.ny(), y.nx())).unwrap()) {
+                    detj in &metrics.detj().broadcast((4, y.ny(), y.nx())).unwrap()) {
         *out = (-eflux - fflux)/detj
     });
 
@@ -584,7 +584,7 @@ pub fn RHS_upwind(
                     fflux in &dF.0,
                     ad_xi in &ad_xi.0,
                     ad_eta in &ad_eta.0,
-                    detj in &metrics.detj.broadcast((4, y.ny(), y.nx())).unwrap()) {
+                    detj in &metrics.detj().broadcast((4, y.ny(), y.nx())).unwrap()) {
         *out = (-eflux - fflux + ad_xi + ad_eta)/detj
     });
 
@@ -611,11 +611,11 @@ fn upwind_dissipation(
         .axis_iter(ndarray::Axis(1))
         .zip(tmp0.axis_iter_mut(ndarray::Axis(1)))
         .zip(tmp1.axis_iter_mut(ndarray::Axis(1)))
-        .zip(metrics.detj.iter())
-        .zip(metrics.detj_dxi_dx.iter())
-        .zip(metrics.detj_dxi_dy.iter())
-        .zip(metrics.detj_deta_dx.iter())
-        .zip(metrics.detj_deta_dy.iter())
+        .zip(metrics.detj().iter())
+        .zip(metrics.detj_dxi_dx().iter())
+        .zip(metrics.detj_dxi_dy().iter())
+        .zip(metrics.detj_deta_dx().iter())
+        .zip(metrics.detj_deta_dy().iter())
     {
         let rho = y[0];
         assert!(rho > 0.0);
@@ -661,10 +661,10 @@ fn upwind_dissipation(
 }
 
 fn fluxes(k: (&mut Field, &mut Field), y: &Field, metrics: &Metrics) {
-    let j_dxi_dx = metrics.detj_dxi_dx.view();
-    let j_dxi_dy = metrics.detj_dxi_dy.view();
-    let j_deta_dx = metrics.detj_deta_dx.view();
-    let j_deta_dy = metrics.detj_deta_dy.view();
+    let j_dxi_dx = metrics.detj_dxi_dx();
+    let j_dxi_dy = metrics.detj_dxi_dy();
+    let j_deta_dx = metrics.detj_deta_dx();
+    let j_deta_dy = metrics.detj_deta_dy();
 
     let rho = y.rho();
     let rhou = y.rhou();
@@ -871,12 +871,17 @@ pub fn extract_boundaries<'a>(
 }
 
 /// Used for storing boundary elements
-pub type BoundaryStorage = Direction<Option<ndarray::Array2<Float>>>;
+pub struct BoundaryStorage {
+    north: Option<ndarray::Array2<Float>>,
+    south: Option<ndarray::Array2<Float>>,
+    east: Option<ndarray::Array2<Float>>,
+    west: Option<ndarray::Array2<Float>>,
+}
 
 impl BoundaryStorage {
     pub fn new(bt: &BoundaryCharacteristics, grid: &Grid) -> Self {
         Self {
-            north: match bt.north {
+            north: match bt.north() {
                 BoundaryCharacteristic::Vortex(_)
                 | BoundaryCharacteristic::Interpolate(_, _)
                 | BoundaryCharacteristic::MultiGrid(_) => {
@@ -884,7 +889,7 @@ impl BoundaryStorage {
                 }
                 _ => None,
             },
-            south: match bt.south {
+            south: match bt.south() {
                 BoundaryCharacteristic::Vortex(_)
                 | BoundaryCharacteristic::Interpolate(_, _)
                 | BoundaryCharacteristic::MultiGrid(_) => {
@@ -892,7 +897,7 @@ impl BoundaryStorage {
                 }
                 _ => None,
             },
-            east: match bt.east {
+            east: match bt.east() {
                 BoundaryCharacteristic::Vortex(_)
                 | BoundaryCharacteristic::Interpolate(_, _)
                 | BoundaryCharacteristic::MultiGrid(_) => {
@@ -900,7 +905,7 @@ impl BoundaryStorage {
                 }
                 _ => None,
             },
-            west: match bt.west {
+            west: match bt.west() {
                 BoundaryCharacteristic::Vortex(_)
                 | BoundaryCharacteristic::Interpolate(_, _)
                 | BoundaryCharacteristic::MultiGrid(_) => {
@@ -955,9 +960,9 @@ fn SAT_characteristics(
             hi,
             sign,
             tau,
-            metrics.detj.slice(slice),
-            metrics.detj_deta_dx.slice(slice),
-            metrics.detj_deta_dy.slice(slice),
+            metrics.detj().slice(slice),
+            metrics.detj_deta_dx().slice(slice),
+            metrics.detj_deta_dy().slice(slice),
         );
     }
     // South boundary
@@ -977,9 +982,9 @@ fn SAT_characteristics(
             hi,
             sign,
             tau,
-            metrics.detj.slice(slice),
-            metrics.detj_deta_dx.slice(slice),
-            metrics.detj_deta_dy.slice(slice),
+            metrics.detj().slice(slice),
+            metrics.detj_deta_dx().slice(slice),
+            metrics.detj_deta_dy().slice(slice),
         );
     }
     // West Boundary
@@ -999,9 +1004,9 @@ fn SAT_characteristics(
             hi,
             sign,
             tau,
-            metrics.detj.slice(slice),
-            metrics.detj_dxi_dx.slice(slice),
-            metrics.detj_dxi_dy.slice(slice),
+            metrics.detj().slice(slice),
+            metrics.detj_dxi_dx().slice(slice),
+            metrics.detj_dxi_dy().slice(slice),
         );
     }
     // East Boundary
@@ -1021,9 +1026,9 @@ fn SAT_characteristics(
             hi,
             sign,
             tau,
-            metrics.detj.slice(slice),
-            metrics.detj_dxi_dx.slice(slice),
-            metrics.detj_dxi_dy.slice(slice),
+            metrics.detj().slice(slice),
+            metrics.detj_dxi_dx().slice(slice),
+            metrics.detj_dxi_dy().slice(slice),
         );
     }
 }
diff --git a/sbp/tests/convergence.rs b/euler/tests/convergence.rs
similarity index 99%
rename from sbp/tests/convergence.rs
rename to euler/tests/convergence.rs
index e689669..4530553 100644
--- a/sbp/tests/convergence.rs
+++ b/euler/tests/convergence.rs
@@ -1,6 +1,6 @@
 #![cfg(feature = "expensive_tests")]
+use euler::*;
 use ndarray::prelude::*;
-use sbp::euler::*;
 use sbp::Float;
 
 fn run_with_size(size: usize, op: impl sbp::operators::UpwindOperator2d + Copy) -> Float {
diff --git a/maxwell/Cargo.toml b/maxwell/Cargo.toml
new file mode 100644
index 0000000..11f75ea
--- /dev/null
+++ b/maxwell/Cargo.toml
@@ -0,0 +1,17 @@
+[package]
+name = "maxwell"
+version = "0.1.0"
+authors = ["Magnus Ulimoen <flymagnus@gmail.com>"]
+edition = "2018"
+
+[dependencies]
+ndarray = "0.13.1"
+sbp = { path = "../sbp" }
+
+[dev-dependencies]
+criterion = "0.3.1"
+
+[[bench]]
+name = "bench"
+harness = false
+
diff --git a/sbp/benches/maxwell.rs b/maxwell/benches/bench.rs
similarity index 98%
rename from sbp/benches/maxwell.rs
rename to maxwell/benches/bench.rs
index 8dbbf99..ee48d69 100644
--- a/sbp/benches/maxwell.rs
+++ b/maxwell/benches/bench.rs
@@ -1,5 +1,5 @@
 use criterion::{black_box, criterion_group, criterion_main, Criterion};
-use sbp::maxwell::System;
+use maxwell::System;
 use sbp::operators::{SbpOperator2d, Upwind4, UpwindOperator2d, SBP4};
 use sbp::Float;
 
diff --git a/sbp/src/maxwell.rs b/maxwell/src/lib.rs
similarity index 88%
rename from sbp/src/maxwell.rs
rename to maxwell/src/lib.rs
index dcc2676..f8e0bd4 100644
--- a/sbp/src/maxwell.rs
+++ b/maxwell/src/lib.rs
@@ -1,9 +1,9 @@
-use super::grid::{Grid, Metrics};
-use super::integrate;
-use super::operators::{SbpOperator2d, UpwindOperator2d};
-use crate::Float;
 use ndarray::azip;
 use ndarray::prelude::*;
+use sbp::grid::{Grid, Metrics};
+use sbp::integrate;
+use sbp::operators::{SbpOperator2d, UpwindOperator2d};
+use sbp::Float;
 
 #[derive(Clone, Debug)]
 pub struct Field(pub(crate) Array3<Float>);
@@ -113,7 +113,7 @@ impl<SBP: SbpOperator2d> System<SBP> {
         let (ex, hz, ey) = self.sys.0.components_mut();
         ndarray::azip!(
             (ex in ex, hz in hz, ey in ey,
-             &x in &self.grid.x, &y in &self.grid.y)
+             &x in &self.grid.x(), &y in &self.grid.y())
         {
             *ex = 0.0;
             *ey = 0.0;
@@ -166,7 +166,7 @@ impl<UO: UpwindOperator2d> System<UO> {
 }
 
 fn gaussian(x: Float, x0: Float, y: Float, y0: Float) -> Float {
-    use crate::consts::PI;
+    use sbp::consts::PI;
 
     let x = x - x0;
     let y = y - y0;
@@ -211,7 +211,7 @@ fn RHS<SBP: SbpOperator2d>(
     SAT_characteristics(op, k, y, metrics, &boundaries);
 
     azip!((k in &mut k.0,
-                    &detj in &metrics.detj.broadcast((3, y.ny(), y.nx())).unwrap()) {
+                    &detj in &metrics.detj().broadcast((3, y.ny(), y.nx())).unwrap()) {
         *k /= detj;
     });
 }
@@ -237,12 +237,12 @@ fn RHS_upwind<UO: UpwindOperator2d>(
     SAT_characteristics(op, k, y, metrics, &boundaries);
 
     azip!((k in &mut k.0,
-                    &detj in &metrics.detj.broadcast((3, y.ny(), y.nx())).unwrap()) {
+                    &detj in &metrics.detj().broadcast((3, y.ny(), y.nx())).unwrap()) {
         *k /= detj;
     });
 }
 
-fn fluxes<SBP: super::operators::SbpOperator2d>(
+fn fluxes<SBP: sbp::operators::SbpOperator2d>(
     op: &SBP,
     k: &mut Field,
     y: &Field,
@@ -252,14 +252,14 @@ fn fluxes<SBP: super::operators::SbpOperator2d>(
     // ex = hz_y
     {
         ndarray::azip!((a in &mut tmp.0,
-                        &dxi_dy in &metrics.detj_dxi_dy,
+                        &dxi_dy in &metrics.detj_dxi_dy(),
                         &hz in &y.hz())
             *a = dxi_dy * hz
         );
         op.diffxi(tmp.0.view(), tmp.1.view_mut());
 
         ndarray::azip!((b in &mut tmp.2,
-                        &deta_dy in &metrics.detj_deta_dy,
+                        &deta_dy in &metrics.detj_deta_dy(),
                         &hz in &y.hz())
             *b = deta_dy * hz
         );
@@ -273,8 +273,8 @@ fn fluxes<SBP: super::operators::SbpOperator2d>(
     {
         // hz = -ey_x + ex_y
         ndarray::azip!((a in &mut tmp.0,
-                        &dxi_dx in &metrics.detj_dxi_dx,
-                        &dxi_dy in &metrics.detj_dxi_dy,
+                        &dxi_dx in &metrics.detj_dxi_dx(),
+                        &dxi_dy in &metrics.detj_dxi_dy(),
                         &ex in &y.ex(),
                         &ey in &y.ey())
             *a = dxi_dx * -ey + dxi_dy * ex
@@ -282,8 +282,8 @@ fn fluxes<SBP: super::operators::SbpOperator2d>(
         op.diffxi(tmp.0.view(), tmp.1.view_mut());
 
         ndarray::azip!((b in &mut tmp.2,
-                        &deta_dx in &metrics.detj_deta_dx,
-                        &deta_dy in &metrics.detj_deta_dy,
+                        &deta_dx in &metrics.detj_deta_dx(),
+                        &deta_dy in &metrics.detj_deta_dy(),
                         &ex in &y.ex(),
                         &ey in &y.ey())
             *b = deta_dx * -ey + deta_dy * ex
@@ -298,14 +298,14 @@ fn fluxes<SBP: super::operators::SbpOperator2d>(
     // ey = -hz_x
     {
         ndarray::azip!((a in &mut tmp.0,
-                        &dxi_dx in &metrics.detj_dxi_dx,
+                        &dxi_dx in &metrics.detj_dxi_dx(),
                         &hz in &y.hz())
             *a = dxi_dx * -hz
         );
         op.diffxi(tmp.0.view(), tmp.1.view_mut());
 
         azip!((b in &mut tmp.2,
-                        &deta_dx in &metrics.detj_deta_dx,
+                        &deta_dx in &metrics.detj_deta_dx(),
                         &hz in &y.hz())
             *b = deta_dx * -hz
         );
@@ -327,8 +327,8 @@ fn dissipation<UO: UpwindOperator2d>(
     // ex component
     {
         ndarray::azip!((a in &mut tmp.0,
-                        &kx in &metrics.detj_dxi_dx,
-                        &ky in &metrics.detj_dxi_dy,
+                        &kx in &metrics.detj_dxi_dx(),
+                        &ky in &metrics.detj_dxi_dy(),
                         &ex in &y.ex(),
                         &ey in &y.ey()) {
             let r = Float::hypot(kx, ky);
@@ -337,8 +337,8 @@ fn dissipation<UO: UpwindOperator2d>(
         op.dissxi(tmp.0.view(), tmp.1.view_mut());
 
         ndarray::azip!((b in &mut tmp.2,
-                    &kx in &metrics.detj_deta_dx,
-                    &ky in &metrics.detj_deta_dy,
+                    &kx in &metrics.detj_deta_dx(),
+                    &ky in &metrics.detj_deta_dy(),
                     &ex in &y.ex(),
                     &ey in &y.ey()) {
             let r = Float::hypot(kx, ky);
@@ -354,8 +354,8 @@ fn dissipation<UO: UpwindOperator2d>(
     // hz component
     {
         ndarray::azip!((a in &mut tmp.0,
-                        &kx in &metrics.detj_dxi_dx,
-                        &ky in &metrics.detj_dxi_dy,
+                        &kx in &metrics.detj_dxi_dx(),
+                        &ky in &metrics.detj_dxi_dy(),
                         &hz in &y.hz()) {
             let r = Float::hypot(kx, ky);
             *a = r * hz;
@@ -363,8 +363,8 @@ fn dissipation<UO: UpwindOperator2d>(
         op.dissxi(tmp.0.view(), tmp.1.view_mut());
 
         ndarray::azip!((b in &mut tmp.2,
-                        &kx in &metrics.detj_deta_dx,
-                        &ky in &metrics.detj_deta_dy,
+                        &kx in &metrics.detj_deta_dx(),
+                        &ky in &metrics.detj_deta_dy(),
                         &hz in &y.hz()) {
             let r = Float::hypot(kx, ky);
             *b = r * hz;
@@ -379,8 +379,8 @@ fn dissipation<UO: UpwindOperator2d>(
     // ey
     {
         ndarray::azip!((a in &mut tmp.0,
-                        &kx in &metrics.detj_dxi_dx,
-                        &ky in &metrics.detj_dxi_dy,
+                        &kx in &metrics.detj_dxi_dx(),
+                        &ky in &metrics.detj_dxi_dy(),
                         &ex in &y.ex(),
                         &ey in &y.ey()) {
             let r = Float::hypot(kx, ky);
@@ -389,8 +389,8 @@ fn dissipation<UO: UpwindOperator2d>(
         op.dissxi(tmp.0.view(), tmp.1.view_mut());
 
         ndarray::azip!((b in &mut tmp.2,
-                    &kx in &metrics.detj_deta_dx,
-                    &ky in &metrics.detj_deta_dy,
+                    &kx in &metrics.detj_deta_dx(),
+                    &ky in &metrics.detj_deta_dy(),
                     &ex in &y.ex(),
                     &ey in &y.ey()) {
             let r = Float::hypot(kx, ky);
@@ -462,8 +462,8 @@ fn SAT_characteristics<SBP: SbpOperator2d>(
             .into_iter()
             .zip(y.slice(s![.., .., nx - 1]).gencolumns())
             .zip(g.gencolumns())
-            .zip(metrics.detj_dxi_dx.slice(s![.., nx - 1]))
-            .zip(metrics.detj_dxi_dy.slice(s![.., nx - 1]))
+            .zip(metrics.detj_dxi_dx().slice(s![.., nx - 1]))
+            .zip(metrics.detj_dxi_dy().slice(s![.., nx - 1]))
         {
             // East boundary, positive flux
             let tau = -1.0;
@@ -500,8 +500,8 @@ fn SAT_characteristics<SBP: SbpOperator2d>(
             .into_iter()
             .zip(y.slice(s![.., .., 0]).gencolumns())
             .zip(g.gencolumns())
-            .zip(metrics.detj_dxi_dx.slice(s![.., 0]))
-            .zip(metrics.detj_dxi_dy.slice(s![.., 0]))
+            .zip(metrics.detj_dxi_dx().slice(s![.., 0]))
+            .zip(metrics.detj_dxi_dy().slice(s![.., 0]))
         {
             let tau = 1.0;
 
@@ -543,8 +543,8 @@ fn SAT_characteristics<SBP: SbpOperator2d>(
             .into_iter()
             .zip(y.slice(s![.., ny - 1, ..]).gencolumns())
             .zip(g.gencolumns())
-            .zip(metrics.detj_deta_dx.slice(s![ny - 1, ..]))
-            .zip(metrics.detj_deta_dy.slice(s![ny - 1, ..]))
+            .zip(metrics.detj_deta_dx().slice(s![ny - 1, ..]))
+            .zip(metrics.detj_deta_dy().slice(s![ny - 1, ..]))
         {
             // North boundary, positive flux
             let tau = -1.0;
@@ -580,8 +580,8 @@ fn SAT_characteristics<SBP: SbpOperator2d>(
             .into_iter()
             .zip(y.slice(s![.., 0, ..]).gencolumns())
             .zip(g.gencolumns())
-            .zip(metrics.detj_deta_dx.slice(s![0, ..]))
-            .zip(metrics.detj_deta_dy.slice(s![0, ..]))
+            .zip(metrics.detj_deta_dx().slice(s![0, ..]))
+            .zip(metrics.detj_deta_dy().slice(s![0, ..]))
         {
             // South boundary, negative flux
 
diff --git a/multigrid/Cargo.toml b/multigrid/Cargo.toml
index de81dab..e68a641 100644
--- a/multigrid/Cargo.toml
+++ b/multigrid/Cargo.toml
@@ -7,7 +7,8 @@ edition = "2018"
 
 [dependencies]
 sbp = { path = "../sbp", features = ["rayon"] }
-hdf5 = "0.6.0"
+euler = { path = "../euler" }
+hdf5 = { version = "0.6.0", features = ["static"] }
 rayon = "1.3.0"
 indicatif = "0.14.0"
 structopt = "0.3.13"
diff --git a/multigrid/src/main.rs b/multigrid/src/main.rs
index 86b325d..f00a239 100644
--- a/multigrid/src/main.rs
+++ b/multigrid/src/main.rs
@@ -173,7 +173,14 @@ fn main() {
         let max_ny = sys.grids.iter().map(|g| g.ny()).max().unwrap();
         std::cmp::max(max_nx, max_ny)
     };
-    let dt = 0.2 / (max_n as Float);
+    // Add a robust method for determining CFL, use for example the maximum speed of the initial
+    // field along with \delta x / \delta y
+    // U_max = max(rhou/u, rhov/v)
+    // This requires scaling with the determinant to obtain the "true" speed in computational
+    // space
+    // CFL = 0.2
+    // \delta t = CFL min(\delta x, \delta y) / U_max
+    let dt = 0.02 / (max_n as Float);
 
     let ntime = (integration_time / dt).round() as u64;
 
diff --git a/multigrid/src/parsing.rs b/multigrid/src/parsing.rs
index 2a398be..2eb8e2a 100644
--- a/multigrid/src/parsing.rs
+++ b/multigrid/src/parsing.rs
@@ -1,17 +1,17 @@
 use super::DiffOp;
-use crate::grid::Grid;
-use crate::Float;
 use either::*;
 use json::JsonValue;
+use sbp::grid::Grid;
 use sbp::utils::h2linspace;
+use sbp::Float;
 
 pub fn json_to_grids(
     mut jsongrids: JsonValue,
-    vortexparams: sbp::euler::VortexParameters,
+    vortexparams: euler::VortexParameters,
 ) -> (
     Vec<String>,
     Vec<sbp::grid::Grid>,
-    Vec<sbp::euler::BoundaryCharacteristics>,
+    Vec<euler::BoundaryCharacteristics>,
     Vec<DiffOp>,
 ) {
     let default = jsongrids.remove("default");
@@ -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.clone())
+                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:") {
@@ -139,13 +139,13 @@ pub fn json_to_grids(
                         Box::new(Interpolation4) as Box<dyn InterpolationOperator>,
                     )
                 };
-                sbp::euler::BoundaryCharacteristic::Interpolate(
+                euler::BoundaryCharacteristic::Interpolate(
                     names.iter().position(|other| other == grid).unwrap(),
                     int_op,
                 )
             } else if let Some(multigrid) = dir.strip_prefix("multi:") {
                 let grids = multigrid.split(':');
-                sbp::euler::BoundaryCharacteristic::MultiGrid(
+                euler::BoundaryCharacteristic::MultiGrid(
                     grids
                         .map(|g| {
                             let rparen = g.find('(').unwrap();
@@ -165,12 +165,12 @@ pub fn json_to_grids(
                         .collect::<Vec<_>>(),
                 )
             } else {
-                sbp::euler::BoundaryCharacteristic::Grid(
+                euler::BoundaryCharacteristic::Grid(
                     names.iter().position(|other| other == dir).unwrap(),
                 )
             }
         }
-        None => sbp::euler::BoundaryCharacteristic::This,
+        None => euler::BoundaryCharacteristic::This,
     };
     for name in &names {
         let bc = &jsongrids[name]["boundary_conditions"];
@@ -179,7 +179,7 @@ pub fn json_to_grids(
         let bc_e = determine_bc(bc["east"].as_str());
         let bc_w = determine_bc(bc["west"].as_str());
 
-        let bc = sbp::euler::BoundaryCharacteristics {
+        let bc = euler::BoundaryCharacteristics {
             north: bc_n,
             south: bc_s,
             east: bc_e,
@@ -349,7 +349,7 @@ fn json2grid(x: JsonValue, y: JsonValue) -> Result<Grid, String> {
     Ok(Grid::new(x, y).unwrap())
 }
 
-pub fn json_to_vortex(mut json: JsonValue) -> super::euler::VortexParameters {
+pub fn json_to_vortex(mut json: JsonValue) -> euler::VortexParameters {
     let mach = json.remove("mach").as_number().unwrap().into();
 
     // Get max length of any (potential) array
@@ -377,9 +377,9 @@ pub fn json_to_vortex(mut json: JsonValue) -> super::euler::VortexParameters {
     let rstar = into_iterator(json.remove("rstar"));
     let eps = into_iterator(json.remove("eps"));
 
-    let mut vortices = sbp::euler::ArrayVec::new();
+    let mut vortices = euler::ArrayVec::new();
     for (((x0, y0), rstar), eps) in x0.zip(y0).zip(rstar).zip(eps).take(maxlen) {
-        vortices.push(sbp::euler::Vortice { x0, y0, rstar, eps })
+        vortices.push(euler::Vortice { x0, y0, rstar, eps })
     }
 
     if !json.is_empty() {
@@ -389,5 +389,5 @@ pub fn json_to_vortex(mut json: JsonValue) -> super::euler::VortexParameters {
         }
     }
 
-    super::euler::VortexParameters { vortices, mach }
+    euler::VortexParameters { vortices, mach }
 }
diff --git a/sbp/Cargo.toml b/sbp/Cargo.toml
index 56fb68d..6b22c41 100644
--- a/sbp/Cargo.toml
+++ b/sbp/Cargo.toml
@@ -9,26 +9,14 @@ 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
-# which should be run only in release builds
-expensive_tests = []
 # Use f32 as precision, default is f64
 f32 = []
 
 [dev-dependencies]
 criterion = "0.3.1"
 
-[[bench]]
-name = "maxwell"
-harness = false
-
-[[bench]]
-name = "euler"
-harness = false
-
 [[bench]]
 name = "sbpoperators"
 harness = false
diff --git a/sbp/src/grid.rs b/sbp/src/grid.rs
index 3857486..66b3149 100644
--- a/sbp/src/grid.rs
+++ b/sbp/src/grid.rs
@@ -1,6 +1,6 @@
 use super::operators::SbpOperator2d;
 use crate::Float;
-use ndarray::Array2;
+use ndarray::{Array2, ArrayView2};
 
 #[derive(Debug, Clone)]
 pub struct Grid {
@@ -115,3 +115,21 @@ impl Metrics {
         })
     }
 }
+
+impl Metrics {
+    pub fn detj(&self) -> ArrayView2<Float> {
+        self.detj.view()
+    }
+    pub fn detj_dxi_dx(&self) -> ArrayView2<Float> {
+        self.detj_dxi_dx.view()
+    }
+    pub fn detj_dxi_dy(&self) -> ArrayView2<Float> {
+        self.detj_dxi_dy.view()
+    }
+    pub fn detj_deta_dx(&self) -> ArrayView2<Float> {
+        self.detj_deta_dx.view()
+    }
+    pub fn detj_deta_dy(&self) -> ArrayView2<Float> {
+        self.detj_deta_dy.view()
+    }
+}
diff --git a/sbp/src/lib.rs b/sbp/src/lib.rs
index 6192d39..2988f20 100644
--- a/sbp/src/lib.rs
+++ b/sbp/src/lib.rs
@@ -6,16 +6,14 @@ pub type Float = f32;
 #[cfg(not(feature = "f32"))]
 pub type Float = f64;
 
-pub(crate) mod consts {
+pub mod consts {
     #[cfg(feature = "f32")]
-    pub(crate) use std::f32::consts::*;
+    pub use std::f32::consts::*;
     #[cfg(not(feature = "f32"))]
-    pub(crate) use std::f64::consts::*;
+    pub use std::f64::consts::*;
 }
 
-pub mod euler;
 pub mod grid;
 pub mod integrate;
-pub mod maxwell;
 pub mod operators;
 pub mod utils;
diff --git a/sbp/src/utils.rs b/sbp/src/utils.rs
index ab9c227..156cb3d 100644
--- a/sbp/src/utils.rs
+++ b/sbp/src/utils.rs
@@ -7,6 +7,33 @@ pub struct Direction<T> {
     pub east: T,
 }
 
+impl<T> Direction<T> {
+    pub fn north(&self) -> &T {
+        &self.north
+    }
+    pub fn north_mut(&mut self) -> &mut T {
+        &mut self.north
+    }
+    pub fn south(&self) -> &T {
+        &self.south
+    }
+    pub fn south_mut(&mut self) -> &mut T {
+        &mut self.south
+    }
+    pub fn east(&self) -> &T {
+        &self.east
+    }
+    pub fn east_mut(&mut self) -> &mut T {
+        &mut self.east
+    }
+    pub fn west(&self) -> &T {
+        &self.west
+    }
+    pub fn west_mut(&mut self) -> &mut T {
+        &mut self.west
+    }
+}
+
 pub fn h2linspace(start: Float, end: Float, n: usize) -> ndarray::Array1<Float> {
     let h = (end - start) / (n - 2) as Float;
     ndarray::Array1::from_shape_fn(n, |i| match i {
diff --git a/webfront/Cargo.toml b/webfront/Cargo.toml
index 7532aca..ddd02d7 100644
--- a/webfront/Cargo.toml
+++ b/webfront/Cargo.toml
@@ -14,3 +14,5 @@ console_error_panic_hook = "0.1.6"
 wee_alloc = "0.4.5"
 sbp = { path = "../sbp", features = ["f32"] }
 ndarray = "0.13.0"
+euler = { path = "../euler" }
+maxwell = { path = "../maxwell" }
diff --git a/webfront/lib.rs b/webfront/lib.rs
index aa44711..71ef3d4 100644
--- a/webfront/lib.rs
+++ b/webfront/lib.rs
@@ -1,6 +1,8 @@
 use wasm_bindgen::prelude::*;
 
-use sbp::{euler, maxwell, operators};
+use euler;
+use maxwell;
+use sbp::operators;
 
 #[cfg(feature = "wee_alloc")]
 #[global_allocator]