diff --git a/multigrid/src/main.rs b/multigrid/src/main.rs
index f3b2ba2..19f9986 100644
--- a/multigrid/src/main.rs
+++ b/multigrid/src/main.rs
@@ -16,11 +16,11 @@ struct System {
 }
 
 enum Metrics {
-    Upwind4(grid::Metrics<operators::Upwind4>),
-    Upwind9(grid::Metrics<operators::Upwind9>),
-    Upwind4h2(grid::Metrics<operators::Upwind4h2>),
-    Trad4(grid::Metrics<operators::SBP4>),
-    Trad8(grid::Metrics<operators::SBP8>),
+    Upwind4(grid::Metrics<operators::Upwind4, operators::Upwind4>),
+    Upwind9(grid::Metrics<operators::Upwind9, operators::Upwind9>),
+    Upwind4h2(grid::Metrics<operators::Upwind4h2, operators::Upwind4h2>),
+    Trad4(grid::Metrics<operators::SBP4, operators::SBP4>),
+    Trad8(grid::Metrics<operators::SBP8, operators::SBP8>),
 }
 
 impl System {
@@ -42,11 +42,20 @@ impl System {
         let metrics = grids
             .iter()
             .map(|g| match operator {
-                "upwind4" => Metrics::Upwind4(g.metrics::<operators::Upwind4>().unwrap()),
-                "upwind9" => Metrics::Upwind9(g.metrics::<operators::Upwind9>().unwrap()),
-                "upwind4h2" => Metrics::Upwind4h2(g.metrics::<operators::Upwind4h2>().unwrap()),
-                "trad4" => Metrics::Trad4(g.metrics::<operators::SBP4>().unwrap()),
-                "trad8" => Metrics::Trad8(g.metrics::<operators::SBP8>().unwrap()),
+                "upwind4" => Metrics::Upwind4(
+                    g.metrics::<operators::Upwind4, operators::Upwind4>()
+                        .unwrap(),
+                ),
+                "upwind9" => Metrics::Upwind9(
+                    g.metrics::<operators::Upwind9, operators::Upwind9>()
+                        .unwrap(),
+                ),
+                "upwind4h2" => Metrics::Upwind4h2(
+                    g.metrics::<operators::Upwind4h2, operators::Upwind4h2>()
+                        .unwrap(),
+                ),
+                "trad4" => Metrics::Trad4(g.metrics::<operators::SBP4, operators::SBP4>().unwrap()),
+                "trad8" => Metrics::Trad8(g.metrics::<operators::SBP8, operators::SBP8>().unwrap()),
                 op => panic!("operator {} not known", op),
             })
             .collect::<Vec<_>>();
diff --git a/sbp/src/euler.rs b/sbp/src/euler.rs
index 9a85e08..dbf95ed 100644
--- a/sbp/src/euler.rs
+++ b/sbp/src/euler.rs
@@ -14,7 +14,7 @@ pub struct System<SBP: SbpOperator> {
     sys: (Field, Field),
     k: [Field; 4],
     wb: WorkBuffers,
-    grid: (Grid, Metrics<SBP>),
+    grid: (Grid, Metrics<SBP, SBP>),
 }
 
 impl<SBP: SbpOperator> System<SBP> {
@@ -403,7 +403,7 @@ fn pressure(gamma: Float, rho: Float, rhou: Float, rhov: Float, e: Float) -> Flo
 pub fn RHS_trad<SBP: SbpOperator>(
     k: &mut Field,
     y: &Field,
-    metrics: &Metrics<SBP>,
+    metrics: &Metrics<SBP, SBP>,
     boundaries: &BoundaryTerms,
     tmp: &mut (Field, Field, Field, Field, Field, Field),
 ) {
@@ -437,7 +437,7 @@ pub fn RHS_trad<SBP: SbpOperator>(
 pub fn RHS_upwind<UO: UpwindOperator>(
     k: &mut Field,
     y: &Field,
-    metrics: &Metrics<UO>,
+    metrics: &Metrics<UO, UO>,
     boundaries: &BoundaryTerms,
     tmp: &mut (Field, Field, Field, Field, Field, Field),
 ) {
@@ -477,7 +477,7 @@ pub fn RHS_upwind<UO: UpwindOperator>(
 fn upwind_dissipation<UO: UpwindOperator>(
     k: (&mut Field, &mut Field),
     y: &Field,
-    metrics: &Metrics<UO>,
+    metrics: &Metrics<UO, UO>,
     tmp: (&mut Field, &mut Field),
 ) {
     let n = y.nx() * y.ny();
@@ -541,7 +541,7 @@ 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, metrics: &Metrics<SBP>) {
+fn fluxes<SBP: SbpOperator>(k: (&mut Field, &mut Field), y: &Field, metrics: &Metrics<SBP, 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();
@@ -815,7 +815,7 @@ fn vortexify(
 fn SAT_characteristics<SBP: SbpOperator>(
     k: &mut Field,
     y: &Field,
-    metrics: &Metrics<SBP>,
+    metrics: &Metrics<SBP, SBP>,
     boundaries: &BoundaryTerms,
 ) {
     // North boundary
diff --git a/sbp/src/grid.rs b/sbp/src/grid.rs
index cda2470..9db994e 100644
--- a/sbp/src/grid.rs
+++ b/sbp/src/grid.rs
@@ -8,9 +8,10 @@ pub struct Grid {
 }
 
 #[derive(Debug, Clone)]
-pub struct Metrics<SBP>
+pub struct Metrics<SBPeta, SBPxi>
 where
-    SBP: super::operators::SbpOperator,
+    SBPeta: super::operators::SbpOperator,
+    SBPxi: super::operators::SbpOperator,
 {
     pub(crate) detj: Array2<Float>,
     pub(crate) detj_dxi_dx: Array2<Float>,
@@ -18,7 +19,8 @@ where
     pub(crate) detj_deta_dx: Array2<Float>,
     pub(crate) detj_deta_dy: Array2<Float>,
 
-    operator: std::marker::PhantomData<SBP>,
+    operatoreta: std::marker::PhantomData<SBPeta>,
+    operatorxi: std::marker::PhantomData<SBPxi>,
 }
 
 impl Grid {
@@ -41,9 +43,9 @@ impl Grid {
         self.y.view()
     }
 
-    pub fn metrics<SBP: super::operators::SbpOperator>(
+    pub fn metrics<SBPx: super::operators::SbpOperator, SBPy: super::operators::SbpOperator>(
         &self,
-    ) -> Result<Metrics<SBP>, ndarray::ShapeError> {
+    ) -> Result<Metrics<SBPx, SBPy>, ndarray::ShapeError> {
         Metrics::new(self)
     }
 
@@ -73,7 +75,9 @@ impl Grid {
     }
 }
 
-impl<SBP: super::operators::SbpOperator> Metrics<SBP> {
+impl<SBPeta: super::operators::SbpOperator, SBPxi: super::operators::SbpOperator>
+    Metrics<SBPeta, SBPxi>
+{
     fn new(grid: &Grid) -> Result<Self, ndarray::ShapeError> {
         let ny = grid.ny();
         let nx = grid.nx();
@@ -81,13 +85,13 @@ impl<SBP: super::operators::SbpOperator> Metrics<SBP> {
         let y = &grid.y;
 
         let mut dx_dxi = Array2::zeros((ny, nx));
-        SBP::diffxi(x.view(), dx_dxi.view_mut());
+        SBPxi::diffxi(x.view(), dx_dxi.view_mut());
         let mut dx_deta = Array2::zeros((ny, nx));
-        SBP::diffeta(x.view(), dx_deta.view_mut());
+        SBPeta::diffeta(x.view(), dx_deta.view_mut());
         let mut dy_dxi = Array2::zeros((ny, nx));
-        SBP::diffxi(y.view(), dy_dxi.view_mut());
+        SBPxi::diffxi(y.view(), dy_dxi.view_mut());
         let mut dy_deta = Array2::zeros((ny, nx));
-        SBP::diffeta(y.view(), dy_deta.view_mut());
+        SBPeta::diffeta(y.view(), dy_deta.view_mut());
 
         let mut detj = Array2::zeros((ny, nx));
         ndarray::azip!((detj in &mut detj,
@@ -118,7 +122,8 @@ impl<SBP: super::operators::SbpOperator> Metrics<SBP> {
             detj_dxi_dy,
             detj_deta_dx,
             detj_deta_dy,
-            operator: std::marker::PhantomData,
+            operatorxi: std::marker::PhantomData,
+            operatoreta: std::marker::PhantomData,
         })
     }
 }
diff --git a/sbp/src/maxwell.rs b/sbp/src/maxwell.rs
index 0dacf35..3a00697 100644
--- a/sbp/src/maxwell.rs
+++ b/sbp/src/maxwell.rs
@@ -78,7 +78,7 @@ pub struct System<SBP: SbpOperator> {
     sys: (Field, Field),
     wb: WorkBuffers,
     grid: Grid,
-    metrics: Metrics<SBP>,
+    metrics: Metrics<SBP, SBP>,
 }
 
 impl<SBP: SbpOperator> System<SBP> {
@@ -119,7 +119,7 @@ impl<SBP: SbpOperator> System<SBP> {
             fut: &mut Field,
             prev: &Field,
             _time: Float,
-            c: &(&Grid, &Metrics<SBP>),
+            c: &(&Grid, &Metrics<SBP, SBP>),
             m: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
         ) {
             let (grid, metrics) = c;
@@ -147,7 +147,7 @@ impl<UO: UpwindOperator> System<UO> {
             fut: &mut Field,
             prev: &Field,
             _time: Float,
-            c: &(&Grid, &Metrics<UO>),
+            c: &(&Grid, &Metrics<UO, UO>),
             m: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
         ) {
             let (grid, metrics) = c;
@@ -199,7 +199,7 @@ fn RHS<SBP: SbpOperator>(
     k: &mut Field,
     y: &Field,
     _grid: &Grid,
-    metrics: &Metrics<SBP>,
+    metrics: &Metrics<SBP, SBP>,
     tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
 ) {
     fluxes(k, y, metrics, tmp);
@@ -223,7 +223,7 @@ fn RHS_upwind<UO: UpwindOperator>(
     k: &mut Field,
     y: &Field,
     _grid: &Grid,
-    metrics: &Metrics<UO>,
+    metrics: &Metrics<UO, UO>,
     tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
 ) {
     fluxes(k, y, metrics, tmp);
@@ -246,7 +246,7 @@ fn RHS_upwind<UO: UpwindOperator>(
 fn fluxes<SBP: SbpOperator>(
     k: &mut Field,
     y: &Field,
-    metrics: &Metrics<SBP>,
+    metrics: &Metrics<SBP, SBP>,
     tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
 ) {
     // ex = hz_y
@@ -320,7 +320,7 @@ fn fluxes<SBP: SbpOperator>(
 fn dissipation<UO: UpwindOperator>(
     k: &mut Field,
     y: &Field,
-    metrics: &Metrics<UO>,
+    metrics: &Metrics<UO, UO>,
     tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
 ) {
     // ex component
@@ -421,7 +421,7 @@ pub struct BoundaryTerms {
 fn SAT_characteristics<SBP: SbpOperator>(
     k: &mut Field,
     y: &Field,
-    metrics: &Metrics<SBP>,
+    metrics: &Metrics<SBP, SBP>,
     boundaries: &BoundaryTerms,
 ) {
     let ny = y.ny();