diff --git a/sbp/examples/multigrid.rs b/sbp/examples/multigrid.rs
index 42f508d..71c54ce 100644
--- a/sbp/examples/multigrid.rs
+++ b/sbp/examples/multigrid.rs
@@ -14,12 +14,13 @@ struct System<T: operators::UpwindOperator> {
         euler::Field,
     )>,
     k: [Vec<euler::Field>; 4],
-    grids: Vec<grid::Grid<T>>,
+    grids: Vec<grid::Grid>,
+    metrics: Vec<grid::Metrics<T>>,
     bt: Vec<euler::BoundaryCharacteristics>,
 }
 
 impl<T: operators::UpwindOperator> System<T> {
-    fn new(grids: Vec<grid::Grid<T>>, bt: Vec<euler::BoundaryCharacteristics>) -> Self {
+    fn new(grids: Vec<grid::Grid>, bt: Vec<euler::BoundaryCharacteristics>) -> Self {
         let fnow = grids
             .iter()
             .map(|g| euler::Field::new(g.ny(), g.nx()))
@@ -33,6 +34,7 @@ impl<T: operators::UpwindOperator> System<T> {
             })
             .collect();
         let k = [fnow.clone(), fnow.clone(), fnow.clone(), fnow.clone()];
+        let metrics = grids.iter().map(|g| g.metrics().unwrap()).collect();
 
         Self {
             fnow,
@@ -40,6 +42,7 @@ impl<T: operators::UpwindOperator> System<T> {
             k,
             wb,
             grids,
+            metrics,
             bt,
         }
     }
@@ -142,14 +145,14 @@ impl<T: operators::UpwindOperator> System<T> {
                 })
                 .collect::<Vec<_>>();
 
-            for ((((prev, fut), grid), wb), bt) in fields
+            for ((((prev, fut), metrics), wb), bt) in fields
                 .iter()
                 .zip(fnext)
-                .zip(&self.grids)
+                .zip(&self.metrics)
                 .zip(&mut self.wb)
                 .zip(bt)
             {
-                euler::RHS_upwind(fut, prev, grid, &bt, wb)
+                euler::RHS_upwind(fut, prev, metrics, &bt, wb)
             }
         }
     }
@@ -257,14 +260,14 @@ impl<T: operators::UpwindOperator> System<T> {
                         },
                     })
                     .collect::<Vec<_>>();
-                for ((((prev, fut), grid), wb), bt) in fields
+                for ((((prev, fut), metrics), wb), bt) in fields
                     .iter()
                     .zip(&mut self.k[i])
-                    .zip(&self.grids)
+                    .zip(&self.metrics)
                     .zip(&mut self.wb)
                     .zip(bt)
                 {
-                    s.spawn(move |_| euler::RHS_upwind(fut, prev, grid, &bt, wb));
+                    s.spawn(move |_| euler::RHS_upwind(fut, prev, metrics, &bt, wb));
                 }
             });
         }
@@ -307,15 +310,13 @@ fn main() {
             west: determine_bc(grid.dirw.as_ref()),
         });
     }
-    let mut grids: Vec<grid::Grid<operators::Upwind4>> = Vec::with_capacity(jgrids.len());
-    for grid in jgrids {
-        grids.push(grid::Grid::new(grid.x, grid.y).unwrap());
-    }
+    let grids = jgrids.into_iter().map(|egrid| egrid.grid).collect();
+
     let integration_time: f64 = json["integration_time"].as_number().unwrap().into();
 
     let vortexparams = utils::json_to_vortex(json["vortex"].clone());
 
-    let mut sys = System::new(grids, bt);
+    let mut sys = System::<sbp::operators::Upwind4>::new(grids, bt);
     sys.vortex(0.0, vortexparams);
 
     let max_n = {
diff --git a/sbp/src/euler.rs b/sbp/src/euler.rs
index 6c03857..3132f71 100644
--- a/sbp/src/euler.rs
+++ b/sbp/src/euler.rs
@@ -1,4 +1,4 @@
-use super::grid::Grid;
+use super::grid::{Grid, Metrics};
 use super::integrate;
 use super::operators::{SbpOperator, UpwindOperator};
 use super::Float;
@@ -13,7 +13,7 @@ pub const GAMMA: Float = 1.4;
 pub struct System<SBP: SbpOperator> {
     sys: (Field, Field),
     wb: WorkBuffers,
-    grid: Grid<SBP>,
+    grid: (Grid, Metrics<SBP>),
 }
 
 impl<SBP: SbpOperator> System<SBP> {
@@ -21,11 +21,12 @@ impl<SBP: SbpOperator> System<SBP> {
         let grid = Grid::new(x, y).expect(
             "Could not create grid. Different number of elements compared to width*height?",
         );
+        let metrics = grid.metrics().unwrap();
         let nx = grid.nx();
         let ny = grid.ny();
         Self {
             sys: (Field::new(ny, nx), Field::new(ny, nx)),
-            grid,
+            grid: (grid, metrics),
             wb: WorkBuffers::new(ny, nx),
         }
     }
@@ -37,16 +38,17 @@ impl<SBP: SbpOperator> System<SBP> {
             east: BoundaryCharacteristic::This,
             west: BoundaryCharacteristic::This,
         };
-        let rhs_trad = |k: &mut Field, y: &Field, grid: &Grid<_>, wb: &mut _| {
+        let rhs_trad = |k: &mut Field, y: &Field, grid: &Grid, metrics: &Metrics<_>, wb: &mut _| {
             let boundaries = boundary_extractor(y, grid, &bc);
-            RHS_trad(k, y, grid, &boundaries, wb)
+            RHS_trad(k, y, metrics, &boundaries, wb)
         };
         integrate::rk4(
             rhs_trad,
             &self.sys.0,
             &mut self.sys.1,
             dt,
-            &self.grid,
+            &self.grid.0,
+            &self.grid.1,
             &mut self.wb.k,
             &mut self.wb.tmp,
         );
@@ -56,7 +58,7 @@ impl<SBP: SbpOperator> System<SBP> {
     pub fn vortex(&mut self, t: Float, vortex_parameters: VortexParameters) {
         self.sys
             .0
-            .vortex(self.grid.x.view(), self.grid.y.view(), t, vortex_parameters);
+            .vortex(self.grid.0.x(), self.grid.0.y(), t, vortex_parameters);
     }
 
     #[allow(clippy::many_single_char_names)]
@@ -70,12 +72,9 @@ impl<SBP: SbpOperator> System<SBP> {
             mach: 0.5,
         };
 
-        self.sys.0.vortex(
-            self.grid.x.view(),
-            self.grid.y.view(),
-            0.0,
-            vortex_parameters,
-        )
+        self.sys
+            .0
+            .vortex(self.grid.0.x(), self.grid.0.y(), 0.0, vortex_parameters)
     }
 
     pub fn field(&self) -> &Field {
@@ -83,10 +82,10 @@ impl<SBP: SbpOperator> System<SBP> {
     }
 
     pub fn x(&self) -> ArrayView2<Float> {
-        self.grid.x.view()
+        self.grid.0.x.view()
     }
     pub fn y(&self) -> ArrayView2<Float> {
-        self.grid.y.view()
+        self.grid.0.y.view()
     }
 }
 
@@ -98,16 +97,18 @@ impl<UO: UpwindOperator> System<UO> {
             east: BoundaryCharacteristic::This,
             west: BoundaryCharacteristic::This,
         };
-        let rhs_upwind = |k: &mut Field, y: &Field, grid: &Grid<_>, wb: &mut _| {
-            let boundaries = boundary_extractor(y, grid, &bc);
-            RHS_upwind(k, y, grid, &boundaries, wb)
-        };
+        let rhs_upwind =
+            |k: &mut Field, y: &Field, grid: &Grid, metrics: &Metrics<_>, wb: &mut _| {
+                let boundaries = boundary_extractor(y, grid, &bc);
+                RHS_upwind(k, y, metrics, &boundaries, wb)
+            };
         integrate::rk4(
             rhs_upwind,
             &self.sys.0,
             &mut self.sys.1,
             dt,
-            &self.grid,
+            &self.grid.0,
+            &self.grid.1,
             &mut self.wb.k,
             &mut self.wb.tmp,
         );
@@ -382,13 +383,13 @@ fn pressure(gamma: Float, rho: Float, rhou: Float, rhov: Float, e: Float) -> Flo
 pub(crate) fn RHS_trad<SBP: SbpOperator>(
     k: &mut Field,
     y: &Field,
-    grid: &Grid<SBP>,
+    metrics: &Metrics<SBP>,
     boundaries: &BoundaryTerms,
     tmp: &mut (Field, Field, Field, Field, Field, Field),
 ) {
     let ehat = &mut tmp.0;
     let fhat = &mut tmp.1;
-    fluxes((ehat, fhat), y, grid);
+    fluxes((ehat, fhat), y, metrics);
     let dE = &mut tmp.2;
     let dF = &mut tmp.3;
 
@@ -405,24 +406,24 @@ pub(crate) fn RHS_trad<SBP: SbpOperator>(
     azip!((out in &mut k.0,
                     eflux in &dE.0,
                     fflux in &dF.0,
-                    detj in &grid.detj.broadcast((4, y.ny(), y.nx())).unwrap()) {
+                    detj in &metrics.detj.broadcast((4, y.ny(), y.nx())).unwrap()) {
         *out = (-eflux - fflux)/detj
     });
 
-    SAT_characteristics(k, y, grid, boundaries);
+    SAT_characteristics(k, y, metrics, boundaries);
 }
 
 #[allow(non_snake_case)]
 pub fn RHS_upwind<UO: UpwindOperator>(
     k: &mut Field,
     y: &Field,
-    grid: &Grid<UO>,
+    metrics: &Metrics<UO>,
     boundaries: &BoundaryTerms,
     tmp: &mut (Field, Field, Field, Field, Field, Field),
 ) {
     let ehat = &mut tmp.0;
     let fhat = &mut tmp.1;
-    fluxes((ehat, fhat), y, grid);
+    fluxes((ehat, fhat), y, metrics);
     let dE = &mut tmp.2;
     let dF = &mut tmp.3;
 
@@ -438,25 +439,25 @@ pub fn RHS_upwind<UO: UpwindOperator>(
 
     let ad_xi = &mut tmp.4;
     let ad_eta = &mut tmp.5;
-    upwind_dissipation((ad_xi, ad_eta), y, grid, (&mut tmp.0, &mut tmp.1));
+    upwind_dissipation((ad_xi, ad_eta), y, metrics, (&mut tmp.0, &mut tmp.1));
 
     azip!((out in &mut k.0,
                     eflux in &dE.0,
                     fflux in &dF.0,
                     ad_xi in &ad_xi.0,
                     ad_eta in &ad_eta.0,
-                    detj in &grid.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
     });
 
-    SAT_characteristics(k, y, grid, boundaries);
+    SAT_characteristics(k, y, metrics, boundaries);
 }
 
 #[allow(clippy::many_single_char_names)]
 fn upwind_dissipation<UO: UpwindOperator>(
     k: (&mut Field, &mut Field),
     y: &Field,
-    grid: &Grid<UO>,
+    metrics: &Metrics<UO>,
     tmp: (&mut Field, &mut Field),
 ) {
     let n = y.nx() * y.ny();
@@ -471,11 +472,11 @@ fn upwind_dissipation<UO: UpwindOperator>(
         .axis_iter(ndarray::Axis(1))
         .zip(tmp0.axis_iter_mut(ndarray::Axis(1)))
         .zip(tmp1.axis_iter_mut(ndarray::Axis(1)))
-        .zip(grid.detj.iter())
-        .zip(grid.detj_dxi_dx.iter())
-        .zip(grid.detj_dxi_dy.iter())
-        .zip(grid.detj_deta_dx.iter())
-        .zip(grid.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);
@@ -520,11 +521,11 @@ fn upwind_dissipation<UO: UpwindOperator>(
     UO::disseta(tmp.1.e(), k.1.e_mut());
 }
 
-fn fluxes<SBP: SbpOperator>(k: (&mut Field, &mut Field), y: &Field, grid: &Grid<SBP>) {
-    let j_dxi_dx = grid.detj_dxi_dx.view();
-    let j_dxi_dy = grid.detj_dxi_dy.view();
-    let j_deta_dx = grid.detj_deta_dx.view();
-    let j_deta_dy = grid.detj_deta_dy.view();
+fn fluxes<SBP: SbpOperator>(k: (&mut Field, &mut Field), y: &Field, metrics: &Metrics<SBP>) {
+    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 rho = y.rho();
     let rhou = y.rhou();
@@ -590,9 +591,9 @@ pub struct BoundaryCharacteristics {
     pub west: BoundaryCharacteristic,
 }
 
-fn boundary_extractor<'a, SBP: SbpOperator>(
+fn boundary_extractor<'a>(
     field: &'a Field,
-    _grid: &Grid<SBP>,
+    _grid: &Grid,
     bc: &BoundaryCharacteristics,
 ) -> BoundaryTerms<'a> {
     BoundaryTerms {
@@ -624,7 +625,7 @@ fn boundary_extractor<'a, SBP: SbpOperator>(
 fn SAT_characteristics<SBP: SbpOperator>(
     k: &mut Field,
     y: &Field,
-    grid: &Grid<SBP>,
+    metrics: &Metrics<SBP>,
     boundaries: &BoundaryTerms,
 ) {
     // North boundary
@@ -640,9 +641,9 @@ fn SAT_characteristics<SBP: SbpOperator>(
             hi,
             sign,
             tau,
-            grid.detj.slice(slice),
-            grid.detj_deta_dx.slice(slice),
-            grid.detj_deta_dy.slice(slice),
+            metrics.detj.slice(slice),
+            metrics.detj_deta_dx.slice(slice),
+            metrics.detj_deta_dy.slice(slice),
         );
     }
     // South boundary
@@ -658,9 +659,9 @@ fn SAT_characteristics<SBP: SbpOperator>(
             hi,
             sign,
             tau,
-            grid.detj.slice(slice),
-            grid.detj_deta_dx.slice(slice),
-            grid.detj_deta_dy.slice(slice),
+            metrics.detj.slice(slice),
+            metrics.detj_deta_dx.slice(slice),
+            metrics.detj_deta_dy.slice(slice),
         );
     }
     // West Boundary
@@ -676,9 +677,9 @@ fn SAT_characteristics<SBP: SbpOperator>(
             hi,
             sign,
             tau,
-            grid.detj.slice(slice),
-            grid.detj_dxi_dx.slice(slice),
-            grid.detj_dxi_dy.slice(slice),
+            metrics.detj.slice(slice),
+            metrics.detj_dxi_dx.slice(slice),
+            metrics.detj_dxi_dy.slice(slice),
         );
     }
     // East Boundary
@@ -694,9 +695,9 @@ fn SAT_characteristics<SBP: SbpOperator>(
             hi,
             sign,
             tau,
-            grid.detj.slice(slice),
-            grid.detj_dxi_dx.slice(slice),
-            grid.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/src/grid.rs b/sbp/src/grid.rs
index ed50d43..0cbce71 100644
--- a/sbp/src/grid.rs
+++ b/sbp/src/grid.rs
@@ -2,13 +2,16 @@ use crate::Float;
 use ndarray::Array2;
 
 #[derive(Debug, Clone)]
-pub struct Grid<SBP>
+pub struct Grid {
+    pub(crate) x: Array2<Float>,
+    pub(crate) y: Array2<Float>,
+}
+
+#[derive(Debug, Clone)]
+pub struct Metrics<SBP>
 where
     SBP: super::operators::SbpOperator,
 {
-    pub(crate) x: Array2<Float>,
-    pub(crate) y: Array2<Float>,
-
     pub(crate) detj: Array2<Float>,
     pub(crate) detj_dxi_dx: Array2<Float>,
     pub(crate) detj_dxi_dy: Array2<Float>,
@@ -18,11 +21,39 @@ where
     operator: std::marker::PhantomData<SBP>,
 }
 
-impl<SBP: super::operators::SbpOperator> Grid<SBP> {
+impl Grid {
     pub fn new(x: Array2<Float>, y: Array2<Float>) -> Result<Self, ndarray::ShapeError> {
         assert_eq!(x.shape(), y.shape());
-        let ny = x.shape()[0];
-        let nx = x.shape()[1];
+
+        Ok(Self { x, y })
+    }
+    pub fn nx(&self) -> usize {
+        self.x.shape()[1]
+    }
+    pub fn ny(&self) -> usize {
+        self.x.shape()[0]
+    }
+
+    pub fn x(&self) -> ndarray::ArrayView2<Float> {
+        self.x.view()
+    }
+    pub fn y(&self) -> ndarray::ArrayView2<Float> {
+        self.y.view()
+    }
+
+    pub fn metrics<SBP: super::operators::SbpOperator>(
+        &self,
+    ) -> Result<Metrics<SBP>, ndarray::ShapeError> {
+        Metrics::new(self)
+    }
+}
+
+impl<SBP: super::operators::SbpOperator> Metrics<SBP> {
+    fn new(grid: &Grid) -> Result<Self, ndarray::ShapeError> {
+        let ny = grid.ny();
+        let nx = grid.nx();
+        let x = &grid.x;
+        let y = &grid.y;
 
         let mut dx_dxi = Array2::zeros((ny, nx));
         SBP::diffxi(x.view(), dx_dxi.view_mut());
@@ -57,8 +88,6 @@ impl<SBP: super::operators::SbpOperator> Grid<SBP> {
         let detj_deta_dy = dx_dxi;
 
         Ok(Self {
-            x,
-            y,
             detj,
             detj_dxi_dx,
             detj_dxi_dy,
@@ -67,17 +96,4 @@ impl<SBP: super::operators::SbpOperator> Grid<SBP> {
             operator: std::marker::PhantomData,
         })
     }
-    pub fn nx(&self) -> usize {
-        self.x.shape()[1]
-    }
-    pub fn ny(&self) -> usize {
-        self.x.shape()[0]
-    }
-
-    pub fn x(&self) -> ndarray::ArrayView2<Float> {
-        self.x.view()
-    }
-    pub fn y(&self) -> ndarray::ArrayView2<Float> {
-        self.y.view()
-    }
 }
diff --git a/sbp/src/integrate.rs b/sbp/src/integrate.rs
index 7969561..7cce642 100644
--- a/sbp/src/integrate.rs
+++ b/sbp/src/integrate.rs
@@ -1,4 +1,4 @@
-use super::grid::Grid;
+use super::grid::{Grid, Metrics};
 use super::operators::SbpOperator;
 use super::Float;
 use ndarray::{Array3, Zip};
@@ -8,13 +8,14 @@ pub(crate) fn rk4<'a, F: 'a, SBP, RHS, WB>(
     prev: &F,
     fut: &mut F,
     dt: Float,
-    grid: &Grid<SBP>,
+    grid: &Grid,
+    metrics: &Metrics<SBP>,
     k: &mut [F; 4],
     mut wb: &mut WB,
 ) where
     F: std::ops::Deref<Target = Array3<Float>> + std::ops::DerefMut<Target = Array3<Float>>,
     SBP: SbpOperator,
-    RHS: Fn(&mut F, &F, &Grid<SBP>, &mut WB),
+    RHS: Fn(&mut F, &F, &Grid, &Metrics<SBP>, &mut WB),
 {
     assert_eq!(prev.shape(), fut.shape());
 
@@ -48,6 +49,6 @@ pub(crate) fn rk4<'a, F: 'a, SBP, RHS, WB>(
             }
         };
 
-        rhs(&mut k[i], &fut, grid, &mut wb);
+        rhs(&mut k[i], &fut, grid, metrics, &mut wb);
     }
 }
diff --git a/sbp/src/maxwell.rs b/sbp/src/maxwell.rs
index 21e3a06..1e7c287 100644
--- a/sbp/src/maxwell.rs
+++ b/sbp/src/maxwell.rs
@@ -1,4 +1,4 @@
-use super::grid::Grid;
+use super::grid::{Grid, Metrics};
 use super::integrate;
 use super::operators::{SbpOperator, UpwindOperator};
 use crate::Float;
@@ -77,7 +77,8 @@ impl Field {
 pub struct System<SBP: SbpOperator> {
     sys: (Field, Field),
     wb: WorkBuffers,
-    grid: Grid<SBP>,
+    grid: Grid,
+    metrics: Metrics<SBP>,
 }
 
 impl<SBP: SbpOperator> System<SBP> {
@@ -87,10 +88,12 @@ impl<SBP: SbpOperator> System<SBP> {
         let nx = x.shape()[1];
 
         let grid = Grid::new(x, y).unwrap();
+        let metrics = grid.metrics().unwrap();
 
         Self {
             sys: (Field::new(ny, nx), Field::new(ny, nx)),
             grid,
+            metrics,
             wb: WorkBuffers::new(ny, nx),
         }
     }
@@ -118,6 +121,7 @@ impl<SBP: SbpOperator> System<SBP> {
             &mut self.sys.1,
             dt,
             &self.grid,
+            &self.metrics,
             &mut self.wb.k,
             &mut self.wb.tmp,
         );
@@ -134,6 +138,7 @@ impl<UO: UpwindOperator> System<UO> {
             &mut self.sys.1,
             dt,
             &self.grid,
+            &self.metrics,
             &mut self.wb.k,
             &mut self.wb.tmp,
         );
@@ -171,10 +176,11 @@ fn gaussian(x: Float, x0: Float, y: Float, y0: Float) -> Float {
 fn RHS<SBP: SbpOperator>(
     k: &mut Field,
     y: &Field,
-    grid: &Grid<SBP>,
+    _grid: &Grid,
+    metrics: &Metrics<SBP>,
     tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
 ) {
-    fluxes(k, y, grid, tmp);
+    fluxes(k, y, metrics, tmp);
 
     let boundaries = BoundaryTerms {
         north: Boundary::This,
@@ -182,10 +188,10 @@ fn RHS<SBP: SbpOperator>(
         west: Boundary::This,
         east: Boundary::This,
     };
-    SAT_characteristics(k, y, grid, &boundaries);
+    SAT_characteristics(k, y, metrics, &boundaries);
 
     azip!((k in &mut k.0,
-                    &detj in &grid.detj.broadcast((3, y.ny(), y.nx())).unwrap()) {
+                    &detj in &metrics.detj.broadcast((3, y.ny(), y.nx())).unwrap()) {
         *k /= detj;
     });
 }
@@ -194,11 +200,12 @@ fn RHS<SBP: SbpOperator>(
 fn RHS_upwind<UO: UpwindOperator>(
     k: &mut Field,
     y: &Field,
-    grid: &Grid<UO>,
+    _grid: &Grid,
+    metrics: &Metrics<UO>,
     tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
 ) {
-    fluxes(k, y, grid, tmp);
-    dissipation(k, y, grid, tmp);
+    fluxes(k, y, metrics, tmp);
+    dissipation(k, y, metrics, tmp);
 
     let boundaries = BoundaryTerms {
         north: Boundary::This,
@@ -206,10 +213,10 @@ fn RHS_upwind<UO: UpwindOperator>(
         west: Boundary::This,
         east: Boundary::This,
     };
-    SAT_characteristics(k, y, grid, &boundaries);
+    SAT_characteristics(k, y, metrics, &boundaries);
 
     azip!((k in &mut k.0,
-                    &detj in &grid.detj.broadcast((3, y.ny(), y.nx())).unwrap()) {
+                    &detj in &metrics.detj.broadcast((3, y.ny(), y.nx())).unwrap()) {
         *k /= detj;
     });
 }
@@ -217,20 +224,20 @@ fn RHS_upwind<UO: UpwindOperator>(
 fn fluxes<SBP: SbpOperator>(
     k: &mut Field,
     y: &Field,
-    grid: &Grid<SBP>,
+    metrics: &Metrics<SBP>,
     tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
 ) {
     // ex = hz_y
     {
         ndarray::azip!((a in &mut tmp.0,
-                        &dxi_dy in &grid.detj_dxi_dy,
+                        &dxi_dy in &metrics.detj_dxi_dy,
                         &hz in &y.hz())
             *a = dxi_dy * hz
         );
         SBP::diffxi(tmp.0.view(), tmp.1.view_mut());
 
         ndarray::azip!((b in &mut tmp.2,
-                        &deta_dy in &grid.detj_deta_dy,
+                        &deta_dy in &metrics.detj_deta_dy,
                         &hz in &y.hz())
             *b = deta_dy * hz
         );
@@ -244,8 +251,8 @@ fn fluxes<SBP: SbpOperator>(
     {
         // hz = -ey_x + ex_y
         ndarray::azip!((a in &mut tmp.0,
-                        &dxi_dx in &grid.detj_dxi_dx,
-                        &dxi_dy in &grid.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
@@ -253,8 +260,8 @@ fn fluxes<SBP: SbpOperator>(
         SBP::diffxi(tmp.0.view(), tmp.1.view_mut());
 
         ndarray::azip!((b in &mut tmp.2,
-                        &deta_dx in &grid.detj_deta_dx,
-                        &deta_dy in &grid.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
@@ -269,14 +276,14 @@ fn fluxes<SBP: SbpOperator>(
     // ey = -hz_x
     {
         ndarray::azip!((a in &mut tmp.0,
-                        &dxi_dx in &grid.detj_dxi_dx,
+                        &dxi_dx in &metrics.detj_dxi_dx,
                         &hz in &y.hz())
             *a = dxi_dx * -hz
         );
         SBP::diffxi(tmp.0.view(), tmp.1.view_mut());
 
         azip!((b in &mut tmp.2,
-                        &deta_dx in &grid.detj_deta_dx,
+                        &deta_dx in &metrics.detj_deta_dx,
                         &hz in &y.hz())
             *b = deta_dx * -hz
         );
@@ -291,14 +298,14 @@ fn fluxes<SBP: SbpOperator>(
 fn dissipation<UO: UpwindOperator>(
     k: &mut Field,
     y: &Field,
-    grid: &Grid<UO>,
+    metrics: &Metrics<UO>,
     tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
 ) {
     // ex component
     {
         ndarray::azip!((a in &mut tmp.0,
-                        &kx in &grid.detj_dxi_dx,
-                        &ky in &grid.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);
@@ -307,8 +314,8 @@ fn dissipation<UO: UpwindOperator>(
         UO::dissxi(tmp.0.view(), tmp.1.view_mut());
 
         ndarray::azip!((b in &mut tmp.2,
-                    &kx in &grid.detj_deta_dx,
-                    &ky in &grid.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);
@@ -324,8 +331,8 @@ fn dissipation<UO: UpwindOperator>(
     // hz component
     {
         ndarray::azip!((a in &mut tmp.0,
-                        &kx in &grid.detj_dxi_dx,
-                        &ky in &grid.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;
@@ -333,8 +340,8 @@ fn dissipation<UO: UpwindOperator>(
         UO::dissxi(tmp.0.view(), tmp.1.view_mut());
 
         ndarray::azip!((b in &mut tmp.2,
-                        &kx in &grid.detj_deta_dx,
-                        &ky in &grid.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;
@@ -349,8 +356,8 @@ fn dissipation<UO: UpwindOperator>(
     // ey
     {
         ndarray::azip!((a in &mut tmp.0,
-                        &kx in &grid.detj_dxi_dx,
-                        &ky in &grid.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);
@@ -359,8 +366,8 @@ fn dissipation<UO: UpwindOperator>(
         UO::dissxi(tmp.0.view(), tmp.1.view_mut());
 
         ndarray::azip!((b in &mut tmp.2,
-                    &kx in &grid.detj_deta_dx,
-                    &ky in &grid.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);
@@ -392,7 +399,7 @@ pub struct BoundaryTerms {
 fn SAT_characteristics<SBP: SbpOperator>(
     k: &mut Field,
     y: &Field,
-    grid: &Grid<SBP>,
+    metrics: &Metrics<SBP>,
     boundaries: &BoundaryTerms,
 ) {
     let ny = y.ny();
@@ -427,8 +434,8 @@ fn SAT_characteristics<SBP: SbpOperator>(
             .into_iter()
             .zip(y.slice(s![.., .., nx - 1]).gencolumns())
             .zip(g.gencolumns())
-            .zip(grid.detj_dxi_dx.slice(s![.., nx - 1]))
-            .zip(grid.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;
@@ -461,8 +468,8 @@ fn SAT_characteristics<SBP: SbpOperator>(
             .into_iter()
             .zip(y.slice(s![.., .., 0]).gencolumns())
             .zip(g.gencolumns())
-            .zip(grid.detj_dxi_dx.slice(s![.., 0]))
-            .zip(grid.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;
 
@@ -500,8 +507,8 @@ fn SAT_characteristics<SBP: SbpOperator>(
             .into_iter()
             .zip(y.slice(s![.., ny - 1, ..]).gencolumns())
             .zip(g.gencolumns())
-            .zip(grid.detj_deta_dx.slice(s![ny - 1, ..]))
-            .zip(grid.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;
@@ -533,8 +540,8 @@ fn SAT_characteristics<SBP: SbpOperator>(
             .into_iter()
             .zip(y.slice(s![.., 0, ..]).gencolumns())
             .zip(g.gencolumns())
-            .zip(grid.detj_deta_dx.slice(s![0, ..]))
-            .zip(grid.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/sbp/src/utils.rs b/sbp/src/utils.rs
index aa2d6c7..dd84fa4 100644
--- a/sbp/src/utils.rs
+++ b/sbp/src/utils.rs
@@ -1,10 +1,10 @@
+use crate::grid::Grid;
 use crate::Float;
 use json::JsonValue;
 
 #[derive(Debug, Clone)]
-pub struct SimpleGrid {
-    pub x: ndarray::Array2<Float>,
-    pub y: ndarray::Array2<Float>,
+pub struct ExtendedGrid {
+    pub grid: Grid,
     pub name: Option<String>,
     pub dire: Option<String>,
     pub dirw: Option<String>,
@@ -28,8 +28,8 @@ pub struct SimpleGrid {
 /// Optional parameters:
 /// * name (for relating boundaries)
 /// * dir{e,w,n,s} (for boundary terms)
-pub fn json_to_grids(json: JsonValue) -> Result<Vec<SimpleGrid>, String> {
-    fn json_to_grid(mut grid: JsonValue) -> Result<SimpleGrid, String> {
+pub fn json_to_grids(json: JsonValue) -> Result<Vec<ExtendedGrid>, String> {
+    fn json_to_grid(mut grid: JsonValue) -> Result<ExtendedGrid, String> {
         #[derive(Debug)]
         enum ArrayForm {
             /// Only know the one dimension, will broadcast to
@@ -173,9 +173,8 @@ pub fn json_to_grids(json: JsonValue) -> Result<Vec<SimpleGrid>, String> {
             }
         }
 
-        Ok(SimpleGrid {
-            x,
-            y,
+        Ok(ExtendedGrid {
+            grid: Grid::new(x, y).unwrap(),
             name,
             dire,
             dirw,
@@ -195,56 +194,73 @@ pub fn json_to_grids(json: JsonValue) -> Result<Vec<SimpleGrid>, String> {
 
 #[test]
 fn parse_linspace() {
-    let grids =
-        json_to_grids(r#"[{"name": "main", "x": "linspace:0:10:20", "y": "linspace:0:10:21"}]"#)
-            .unwrap();
+    let grids = json_to_grids(
+        json::parse(r#"[{"name": "main", "x": "linspace:0:10:20", "y": "linspace:0:10:21"}]"#)
+            .unwrap(),
+    )
+    .unwrap();
     assert_eq!(grids.len(), 1);
-    assert_eq!(grids[0].x.shape(), [21, 20]);
-    assert_eq!(grids[0].y.shape(), [21, 20]);
+    assert_eq!(grids[0].grid.x.shape(), [21, 20]);
+    assert_eq!(grids[0].grid.y.shape(), [21, 20]);
     assert_eq!(grids[0].name.as_ref().unwrap(), "main");
-    let grids =
-        json_to_grids(r#"{"name": "main", "x": "linspace:0:10:20", "y": "linspace:0:10:21"}"#)
-            .unwrap();
+    let grids = json_to_grids(
+        json::parse(r#"{"name": "main", "x": "linspace:0:10:20", "y": "linspace:0:10:21"}"#)
+            .unwrap(),
+    )
+    .unwrap();
     assert_eq!(grids.len(), 1);
-    assert_eq!(grids[0].x.shape(), [21, 20]);
-    assert_eq!(grids[0].y.shape(), [21, 20]);
+    assert_eq!(grids[0].grid.x.shape(), [21, 20]);
+    assert_eq!(grids[0].grid.y.shape(), [21, 20]);
     assert_eq!(grids[0].name.as_ref().unwrap(), "main");
 }
 
 #[test]
 fn parse_1d() {
-    let grids = json_to_grids(r#"{"x": [1, 2, 3, 4, 5.1, 3], "y": [1, 2]}"#).unwrap();
+    let grids =
+        json_to_grids(json::parse(r#"{"x": [1, 2, 3, 4, 5.1, 3], "y": [1, 2]}"#).unwrap()).unwrap();
     assert_eq!(grids.len(), 1);
     let grid = &grids[0];
-    assert_eq!(grid.x.shape(), &[2, 6]);
-    assert_eq!(grid.x.shape(), grid.y.shape());
+    assert_eq!(grid.grid.x.shape(), &[2, 6]);
+    assert_eq!(grid.grid.x.shape(), grid.grid.y.shape());
 }
 
 #[test]
 fn parse_2d() {
-    let grids = json_to_grids(r#"{"x": [[1, 2], [3, 4], [5.1, 3]], "y": [1, 2, 3]}"#).unwrap();
-    assert_eq!(grids.len(), 1);
-    let grid = &grids[0];
-    assert_eq!(grid.x.shape(), &[3, 2]);
-    assert_eq!(grid.x.shape(), grid.y.shape());
-    json_to_grids(r#"{"x": [[1, 2], [3, 4], [5.1, 3], [1]], "y": [1, 2, 3]}"#).unwrap_err();
-    json_to_grids(r#"{"y": [[1, 2], [3, 4], [5.1, 3], [1]], "x": [1, 2, 3]}"#).unwrap_err();
     let grids =
-        json_to_grids(r#"{"x": [[1, 2], [3, 4], [5.1, 3]], "y": [[1, 2], [3, 4], [5, 6]]}"#)
+        json_to_grids(json::parse(r#"{"x": [[1, 2], [3, 4], [5.1, 3]], "y": [1, 2, 3]}"#).unwrap())
             .unwrap();
     assert_eq!(grids.len(), 1);
-    json_to_grids(r#"{"x": [[1, 2], [3, 4], [5.1, 3]], "y": [[1, 2], [3, 4], [5]]}"#).unwrap_err();
+    let grid = &grids[0];
+    assert_eq!(grid.grid.x.shape(), &[3, 2]);
+    assert_eq!(grid.grid.x.shape(), grid.grid.y.shape());
+    json_to_grids(
+        json::parse(r#"{"x": [[1, 2], [3, 4], [5.1, 3], [1]], "y": [1, 2, 3]}"#).unwrap(),
+    )
+    .unwrap_err();
+    json_to_grids(
+        json::parse(r#"{"y": [[1, 2], [3, 4], [5.1, 3], [1]], "x": [1, 2, 3]}"#).unwrap(),
+    )
+    .unwrap_err();
+    let grids = json_to_grids(
+        json::parse(r#"{"x": [[1, 2], [3, 4], [5.1, 3]], "y": [[1, 2], [3, 4], [5, 6]]}"#).unwrap(),
+    )
+    .unwrap();
+    assert_eq!(grids.len(), 1);
+    json_to_grids(
+        json::parse(r#"{"x": [[1, 2], [3, 4], [5.1, 3]], "y": [[1, 2], [3, 4], [5]]}"#).unwrap(),
+    )
+    .unwrap_err();
 }
 
 #[test]
 fn parse_err() {
-    json_to_grids(r#"[{"#).unwrap_err();
-    json_to_grids(r#"{}"#).unwrap_err();
-    json_to_grids(r#"0.45"#).unwrap_err();
-    json_to_grids(r#"{"x": "linspace", "y": [0.1, 0.2]}"#).unwrap_err();
-    json_to_grids(r#"{"x": "linspace:::", "y": [0.1, 0.2]}"#).unwrap_err();
-    json_to_grids(r#"{"x": "linspace:1.2:3.1:412.2", "y": [0.1, 0.2]}"#).unwrap_err();
-    json_to_grids(r#"{"x": [-2, -3, "dfd"], "y": [0.1, 0.2]}"#).unwrap_err();
+    json_to_grids(json::parse(r#"{}"#).unwrap()).unwrap_err();
+    json_to_grids(json::parse(r#"0.45"#).unwrap()).unwrap_err();
+    json_to_grids(json::parse(r#"{"x": "linspace", "y": [0.1, 0.2]}"#).unwrap()).unwrap_err();
+    json_to_grids(json::parse(r#"{"x": "linspace:::", "y": [0.1, 0.2]}"#).unwrap()).unwrap_err();
+    json_to_grids(json::parse(r#"{"x": "linspace:1.2:3.1:412.2", "y": [0.1, 0.2]}"#).unwrap())
+        .unwrap_err();
+    json_to_grids(json::parse(r#"{"x": [-2, -3, "dfd"], "y": [0.1, 0.2]}"#).unwrap()).unwrap_err();
 }
 
 pub fn json_to_vortex(mut json: JsonValue) -> super::euler::VortexParameters {