diff --git a/sbp/benches/euler.rs b/sbp/benches/euler.rs
index 07fa2a4..2ed59b8 100644
--- a/sbp/benches/euler.rs
+++ b/sbp/benches/euler.rs
@@ -25,7 +25,7 @@ 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 = System::<Upwind4>::new(x.into_owned(), y.into_owned());
+    let mut universe = System::new(x.into_owned(), y.into_owned(), Upwind4);
     group.bench_function("advance", |b| {
         b.iter(|| {
             universe.init_with_vortex(0.0, 0.0);
@@ -33,7 +33,7 @@ fn performance_benchmark(c: &mut Criterion) {
         })
     });
 
-    let mut universe = System::<Upwind4>::new(x.into_owned(), y.into_owned());
+    let mut universe = System::new(x.into_owned(), y.into_owned(), Upwind4);
     group.bench_function("advance_upwind", |b| {
         b.iter(|| {
             universe.init_with_vortex(0.0, 0.0);
@@ -41,7 +41,7 @@ fn performance_benchmark(c: &mut Criterion) {
         })
     });
 
-    let mut universe = System::<SBP4>::new(x.into_owned(), y.into_owned());
+    let mut universe = System::new(x.into_owned(), y.into_owned(), SBP4);
     group.bench_function("advance_trad4", |b| {
         b.iter(|| {
             universe.init_with_vortex(0.0, 0.0);
diff --git a/sbp/benches/maxwell.rs b/sbp/benches/maxwell.rs
index f532e5e..70005ca 100644
--- a/sbp/benches/maxwell.rs
+++ b/sbp/benches/maxwell.rs
@@ -24,7 +24,7 @@ fn performance_benchmark(c: &mut Criterion) {
     let x = ndarray::Array2::from_shape_fn((h, w), |(_, i)| i as Float / (w - 1) as Float);
     let y = ndarray::Array2::from_shape_fn((h, w), |(j, _)| j as Float / (h - 1) as Float);
 
-    let mut universe = System::<Upwind4>::new(x.clone(), y.clone());
+    let mut universe = System::new(x.clone(), y.clone(), Upwind4);
     group.bench_function("advance", |b| {
         b.iter(|| {
             universe.set_gaussian(0.5, 0.5);
@@ -32,7 +32,7 @@ fn performance_benchmark(c: &mut Criterion) {
         })
     });
 
-    let mut universe = System::<Upwind4>::new(x.clone(), y.clone());
+    let mut universe = System::new(x.clone(), y.clone(), Upwind4);
     group.bench_function("advance_upwind", |b| {
         b.iter(|| {
             universe.set_gaussian(0.5, 0.5);
@@ -40,7 +40,7 @@ fn performance_benchmark(c: &mut Criterion) {
         })
     });
 
-    let mut universe = System::<SBP4>::new(x, y);
+    let mut universe = System::new(x, y, SBP4);
     group.bench_function("advance_trad4", |b| {
         b.iter(|| {
             universe.set_gaussian(0.5, 0.5);
diff --git a/sbp/src/euler.rs b/sbp/src/euler.rs
index 08b824a..a988159 100644
--- a/sbp/src/euler.rs
+++ b/sbp/src/euler.rs
@@ -14,15 +14,16 @@ pub struct System<SBP: SbpOperator> {
     sys: (Field, Field),
     k: [Field; 4],
     wb: WorkBuffers,
-    grid: (Grid, Metrics<SBP, SBP>),
+    grid: (Grid, Metrics),
+    op: SBP,
 }
 
 impl<SBP: SbpOperator> System<SBP> {
-    pub fn new(x: ndarray::Array2<Float>, y: ndarray::Array2<Float>) -> Self {
+    pub fn new(x: ndarray::Array2<Float>, y: ndarray::Array2<Float>, op: SBP) -> Self {
         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 metrics = grid.metrics(op, op).unwrap();
         let nx = grid.nx();
         let ny = grid.ny();
         Self {
@@ -35,6 +36,7 @@ impl<SBP: SbpOperator> System<SBP> {
                 Field::new(ny, nx),
             ],
             wb: WorkBuffers::new(ny, nx),
+            op,
         }
     }
 
@@ -45,10 +47,11 @@ impl<SBP: SbpOperator> System<SBP> {
             east: BoundaryCharacteristic::This,
             west: BoundaryCharacteristic::This,
         };
+        let op = self.op;
         let rhs_trad = |k: &mut Field, y: &Field, _time: Float, gm: &(_, _), wb: &mut _| {
             let (grid, metrics) = gm;
             let boundaries = boundary_extractor(y, grid, &bc);
-            RHS_trad(k, y, metrics, &boundaries, wb)
+            RHS_trad((op, op), k, y, metrics, &boundaries, wb)
         };
         integrate::integrate::<integrate::Rk4, _, _, _, _>(
             rhs_trad,
@@ -105,10 +108,11 @@ impl<UO: UpwindOperator> System<UO> {
             east: BoundaryCharacteristic::This,
             west: BoundaryCharacteristic::This,
         };
+        let op = self.op;
         let rhs_upwind = |k: &mut Field, y: &Field, _time: Float, gm: &(_, _), wb: &mut _| {
             let (grid, metrics) = gm;
             let boundaries = boundary_extractor(y, grid, &bc);
-            RHS_upwind(k, y, metrics, &boundaries, wb)
+            RHS_upwind((op, op), k, y, metrics, &boundaries, wb)
         };
         integrate::integrate::<integrate::Rk4, _, _, _, _>(
             rhs_upwind,
@@ -268,12 +272,14 @@ impl Field {
 
 impl Field {
     /// sqrt((self-other)^T*H*(self-other))
-    pub fn h2_err<SBP: SbpOperator>(&self, other: &Self) -> Float {
+    pub fn h2_err<SBPeta: SbpOperator, SBPxi: SbpOperator>(
+        &self,
+        (opeta, opxi): (SBPeta, SBPxi),
+        other: &Self,
+    ) -> Float {
         assert_eq!(self.nx(), other.nx());
         assert_eq!(self.ny(), other.ny());
 
-        let h = SBP::h();
-
         // Resulting structure should be
         // serialized(F0 - F1)^T (Hx kron Hy) serialized(F0 - F1)
         //
@@ -282,7 +288,7 @@ impl Field {
 
         // This chains the h block into the form [h, 1, 1, 1, rev(h)],
         // and multiplies with a factor
-        let itermaker = move |n: usize, factor: Float| {
+        let itermaker = move |h: &'static [Float], n: usize, factor: Float| {
             h.iter()
                 .copied()
                 .chain(std::iter::repeat(1.0).take(n - 2 * h.len()))
@@ -291,8 +297,9 @@ impl Field {
         };
 
         let hxiterator = itermaker(
+            opxi.h(),
             self.nx(),
-            if SBP::is_h2() {
+            if opxi.is_h2() {
                 1.0 / (self.nx() - 2) as Float
             } else {
                 1.0 / (self.nx() - 1) as Float
@@ -303,8 +310,9 @@ impl Field {
         let hxiterator = hxiterator.cycle().take(self.nx() * self.ny());
 
         let hyiterator = itermaker(
+            opeta.h(),
             self.ny(),
-            1.0 / if SBP::is_h2() {
+            1.0 / if opeta.is_h2() {
                 (self.ny() - 2) as Float
             } else {
                 (self.ny() - 1) as Float
@@ -333,10 +341,12 @@ fn h2_diff() {
     }
     let field1 = Field::new(20, 21);
 
-    assert!((field0.h2_err::<super::operators::Upwind4>(&field1).powi(2) - 4.0).abs() < 1e-3);
-    assert!((field0.h2_err::<super::operators::Upwind9>(&field1).powi(2) - 4.0).abs() < 1e-3);
-    assert!((field0.h2_err::<super::operators::SBP4>(&field1).powi(2) - 4.0).abs() < 1e-3);
-    assert!((field0.h2_err::<super::operators::SBP8>(&field1).powi(2) - 4.0).abs() < 1e-3);
+    use super::operators::{Upwind4, Upwind9, SBP4, SBP8};
+
+    assert!((field0.h2_err((Upwind4, Upwind4), &field1).powi(2) - 4.0).abs() < 1e-3);
+    assert!((field0.h2_err((Upwind9, Upwind9), &field1).powi(2) - 4.0).abs() < 1e-3);
+    assert!((field0.h2_err((SBP4, SBP4), &field1).powi(2) - 4.0).abs() < 1e-3);
+    assert!((field0.h2_err((SBP8, SBP8), &field1).powi(2) - 4.0).abs() < 1e-3);
 }
 
 #[derive(Copy, Clone, Debug)]
@@ -401,9 +411,10 @@ fn pressure(gamma: Float, rho: Float, rhou: Float, rhov: Float, e: Float) -> Flo
 
 #[allow(non_snake_case)]
 pub fn RHS_trad<SBPeta: SbpOperator, SBPxi: SbpOperator>(
+    (opeta, opxi): (SBPeta, SBPxi),
     k: &mut Field,
     y: &Field,
-    metrics: &Metrics<SBPeta, SBPxi>,
+    metrics: &Metrics,
     boundaries: &BoundaryTerms,
     tmp: &mut (Field, Field, Field, Field, Field, Field),
 ) {
@@ -413,15 +424,15 @@ pub fn RHS_trad<SBPeta: SbpOperator, SBPxi: SbpOperator>(
     let dE = &mut tmp.2;
     let dF = &mut tmp.3;
 
-    SBPxi::diffxi(ehat.rho(), dE.rho_mut());
-    SBPxi::diffxi(ehat.rhou(), dE.rhou_mut());
-    SBPxi::diffxi(ehat.rhov(), dE.rhov_mut());
-    SBPxi::diffxi(ehat.e(), dE.e_mut());
+    opxi.diffxi(ehat.rho(), dE.rho_mut());
+    opxi.diffxi(ehat.rhou(), dE.rhou_mut());
+    opxi.diffxi(ehat.rhov(), dE.rhov_mut());
+    opxi.diffxi(ehat.e(), dE.e_mut());
 
-    SBPeta::diffeta(fhat.rho(), dF.rho_mut());
-    SBPeta::diffeta(fhat.rhou(), dF.rhou_mut());
-    SBPeta::diffeta(fhat.rhov(), dF.rhov_mut());
-    SBPeta::diffeta(fhat.e(), dF.e_mut());
+    opeta.diffeta(fhat.rho(), dF.rho_mut());
+    opeta.diffeta(fhat.rhou(), dF.rhou_mut());
+    opeta.diffeta(fhat.rhov(), dF.rhov_mut());
+    opeta.diffeta(fhat.e(), dF.e_mut());
 
     azip!((out in &mut k.0,
                     eflux in &dE.0,
@@ -430,14 +441,15 @@ pub fn RHS_trad<SBPeta: SbpOperator, SBPxi: SbpOperator>(
         *out = (-eflux - fflux)/detj
     });
 
-    SAT_characteristics::<SBPeta, SBPxi>(k, y, metrics, boundaries);
+    SAT_characteristics((opeta, opxi), k, y, metrics, boundaries);
 }
 
 #[allow(non_snake_case)]
 pub fn RHS_upwind<UOeta: UpwindOperator, UOxi: UpwindOperator>(
+    (opeta, opxi): (UOeta, UOxi),
     k: &mut Field,
     y: &Field,
-    metrics: &Metrics<UOeta, UOxi>,
+    metrics: &Metrics,
     boundaries: &BoundaryTerms,
     tmp: &mut (Field, Field, Field, Field, Field, Field),
 ) {
@@ -447,19 +459,25 @@ pub fn RHS_upwind<UOeta: UpwindOperator, UOxi: UpwindOperator>(
     let dE = &mut tmp.2;
     let dF = &mut tmp.3;
 
-    UOxi::diffxi(ehat.rho(), dE.rho_mut());
-    UOxi::diffxi(ehat.rhou(), dE.rhou_mut());
-    UOxi::diffxi(ehat.rhov(), dE.rhov_mut());
-    UOxi::diffxi(ehat.e(), dE.e_mut());
+    opxi.diffxi(ehat.rho(), dE.rho_mut());
+    opxi.diffxi(ehat.rhou(), dE.rhou_mut());
+    opxi.diffxi(ehat.rhov(), dE.rhov_mut());
+    opxi.diffxi(ehat.e(), dE.e_mut());
 
-    UOeta::diffeta(fhat.rho(), dF.rho_mut());
-    UOeta::diffeta(fhat.rhou(), dF.rhou_mut());
-    UOeta::diffeta(fhat.rhov(), dF.rhov_mut());
-    UOeta::diffeta(fhat.e(), dF.e_mut());
+    opeta.diffeta(fhat.rho(), dF.rho_mut());
+    opeta.diffeta(fhat.rhou(), dF.rhou_mut());
+    opeta.diffeta(fhat.rhov(), dF.rhov_mut());
+    opeta.diffeta(fhat.e(), dF.e_mut());
 
     let ad_xi = &mut tmp.4;
     let ad_eta = &mut tmp.5;
-    upwind_dissipation::<UOeta, UOxi>((ad_xi, ad_eta), y, metrics, (&mut tmp.0, &mut tmp.1));
+    upwind_dissipation(
+        (opeta, opxi),
+        (ad_xi, ad_eta),
+        y,
+        metrics,
+        (&mut tmp.0, &mut tmp.1),
+    );
 
     azip!((out in &mut k.0,
                     eflux in &dE.0,
@@ -470,14 +488,15 @@ pub fn RHS_upwind<UOeta: UpwindOperator, UOxi: UpwindOperator>(
         *out = (-eflux - fflux + ad_xi + ad_eta)/detj
     });
 
-    SAT_characteristics::<UOeta, UOxi>(k, y, metrics, boundaries);
+    SAT_characteristics((opeta, opxi), k, y, metrics, boundaries);
 }
 
 #[allow(clippy::many_single_char_names)]
 fn upwind_dissipation<UOeta: UpwindOperator, UOxi: UpwindOperator>(
+    (opeta, opxi): (UOeta, UOxi),
     k: (&mut Field, &mut Field),
     y: &Field,
-    metrics: &Metrics<UOeta, UOxi>,
+    metrics: &Metrics,
     tmp: (&mut Field, &mut Field),
 ) {
     let n = y.nx() * y.ny();
@@ -530,22 +549,18 @@ fn upwind_dissipation<UOeta: UpwindOperator, UOxi: UpwindOperator>(
         tmp1[3] = alpha_v * e * detj;
     }
 
-    UOxi::dissxi(tmp.0.rho(), k.0.rho_mut());
-    UOxi::dissxi(tmp.0.rhou(), k.0.rhou_mut());
-    UOxi::dissxi(tmp.0.rhov(), k.0.rhov_mut());
-    UOxi::dissxi(tmp.0.e(), k.0.e_mut());
+    opxi.dissxi(tmp.0.rho(), k.0.rho_mut());
+    opxi.dissxi(tmp.0.rhou(), k.0.rhou_mut());
+    opxi.dissxi(tmp.0.rhov(), k.0.rhov_mut());
+    opxi.dissxi(tmp.0.e(), k.0.e_mut());
 
-    UOeta::disseta(tmp.1.rho(), k.1.rho_mut());
-    UOeta::disseta(tmp.1.rhou(), k.1.rhou_mut());
-    UOeta::disseta(tmp.1.rhov(), k.1.rhov_mut());
-    UOeta::disseta(tmp.1.e(), k.1.e_mut());
+    opeta.disseta(tmp.1.rho(), k.1.rho_mut());
+    opeta.disseta(tmp.1.rhou(), k.1.rhou_mut());
+    opeta.disseta(tmp.1.rhov(), k.1.rhov_mut());
+    opeta.disseta(tmp.1.e(), k.1.e_mut());
 }
 
-fn fluxes<SBP1: SbpOperator, SBP2: SbpOperator>(
-    k: (&mut Field, &mut Field),
-    y: &Field,
-    metrics: &Metrics<SBP1, SBP2>,
-) {
+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();
@@ -839,17 +854,18 @@ fn vortexify(
 #[allow(non_snake_case)]
 /// Boundary conditions (SAT)
 fn SAT_characteristics<SBPeta: SbpOperator, SBPxi: SbpOperator>(
+    (opeta, opxi): (SBPeta, SBPxi),
     k: &mut Field,
     y: &Field,
-    metrics: &Metrics<SBPeta, SBPxi>,
+    metrics: &Metrics,
     boundaries: &BoundaryTerms,
 ) {
     // North boundary
     {
-        let hi = if SBPeta::is_h2() {
-            (k.ny() - 2) as Float / SBPeta::h()[0]
+        let hi = if opeta.is_h2() {
+            (k.ny() - 2) as Float / opeta.h()[0]
         } else {
-            (k.ny() - 1) as Float / SBPeta::h()[0]
+            (k.ny() - 1) as Float / opeta.h()[0]
         };
         let sign = -1.0;
         let tau = 1.0;
@@ -868,10 +884,10 @@ fn SAT_characteristics<SBPeta: SbpOperator, SBPxi: SbpOperator>(
     }
     // South boundary
     {
-        let hi = if SBPeta::is_h2() {
-            (k.ny() - 2) as Float / SBPeta::h()[0]
+        let hi = if opeta.is_h2() {
+            (k.ny() - 2) as Float / opeta.h()[0]
         } else {
-            (k.ny() - 1) as Float / SBPeta::h()[0]
+            (k.ny() - 1) as Float / opeta.h()[0]
         };
         let sign = 1.0;
         let tau = -1.0;
@@ -890,10 +906,10 @@ fn SAT_characteristics<SBPeta: SbpOperator, SBPxi: SbpOperator>(
     }
     // West Boundary
     {
-        let hi = if SBPxi::is_h2() {
-            (k.nx() - 2) as Float / SBPxi::h()[0]
+        let hi = if opxi.is_h2() {
+            (k.nx() - 2) as Float / opxi.h()[0]
         } else {
-            (k.nx() - 1) as Float / SBPxi::h()[0]
+            (k.nx() - 1) as Float / opxi.h()[0]
         };
         let sign = 1.0;
         let tau = -1.0;
@@ -912,10 +928,10 @@ fn SAT_characteristics<SBPeta: SbpOperator, SBPxi: SbpOperator>(
     }
     // East Boundary
     {
-        let hi = if SBPxi::is_h2() {
-            (k.nx() - 2) as Float / SBPxi::h()[0]
+        let hi = if opxi.is_h2() {
+            (k.nx() - 2) as Float / opxi.h()[0]
         } else {
-            (k.nx() - 1) as Float / SBPxi::h()[0]
+            (k.nx() - 1) as Float / opxi.h()[0]
         };
         let sign = -1.0;
         let tau = 1.0;
diff --git a/sbp/src/grid.rs b/sbp/src/grid.rs
index 9db994e..c99e6a3 100644
--- a/sbp/src/grid.rs
+++ b/sbp/src/grid.rs
@@ -8,19 +8,12 @@ pub struct Grid {
 }
 
 #[derive(Debug, Clone)]
-pub struct Metrics<SBPeta, SBPxi>
-where
-    SBPeta: super::operators::SbpOperator,
-    SBPxi: super::operators::SbpOperator,
-{
+pub struct Metrics {
     pub(crate) detj: Array2<Float>,
     pub(crate) detj_dxi_dx: Array2<Float>,
     pub(crate) detj_dxi_dy: Array2<Float>,
     pub(crate) detj_deta_dx: Array2<Float>,
     pub(crate) detj_deta_dy: Array2<Float>,
-
-    operatoreta: std::marker::PhantomData<SBPeta>,
-    operatorxi: std::marker::PhantomData<SBPxi>,
 }
 
 impl Grid {
@@ -43,10 +36,12 @@ impl Grid {
         self.y.view()
     }
 
-    pub fn metrics<SBPx: super::operators::SbpOperator, SBPy: super::operators::SbpOperator>(
+    pub fn metrics<SBPeta: super::operators::SbpOperator, SBPxi: super::operators::SbpOperator>(
         &self,
-    ) -> Result<Metrics<SBPx, SBPy>, ndarray::ShapeError> {
-        Metrics::new(self)
+        opeta: SBPeta,
+        opxi: SBPxi,
+    ) -> Result<Metrics, ndarray::ShapeError> {
+        Metrics::new(self, opeta, opxi)
     }
 
     pub fn north(&self) -> (ndarray::ArrayView1<Float>, ndarray::ArrayView1<Float>) {
@@ -75,23 +70,25 @@ impl Grid {
     }
 }
 
-impl<SBPeta: super::operators::SbpOperator, SBPxi: super::operators::SbpOperator>
-    Metrics<SBPeta, SBPxi>
-{
-    fn new(grid: &Grid) -> Result<Self, ndarray::ShapeError> {
+impl Metrics {
+    fn new<SBPeta: super::operators::SbpOperator, SBPxi: super::operators::SbpOperator>(
+        grid: &Grid,
+        opeta: SBPeta,
+        opxi: SBPxi,
+    ) -> 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));
-        SBPxi::diffxi(x.view(), dx_dxi.view_mut());
+        opxi.diffxi(x.view(), dx_dxi.view_mut());
         let mut dx_deta = Array2::zeros((ny, nx));
-        SBPeta::diffeta(x.view(), dx_deta.view_mut());
+        opeta.diffeta(x.view(), dx_deta.view_mut());
         let mut dy_dxi = Array2::zeros((ny, nx));
-        SBPxi::diffxi(y.view(), dy_dxi.view_mut());
+        opxi.diffxi(y.view(), dy_dxi.view_mut());
         let mut dy_deta = Array2::zeros((ny, nx));
-        SBPeta::diffeta(y.view(), dy_deta.view_mut());
+        opeta.diffeta(y.view(), dy_deta.view_mut());
 
         let mut detj = Array2::zeros((ny, nx));
         ndarray::azip!((detj in &mut detj,
@@ -122,8 +119,6 @@ impl<SBPeta: super::operators::SbpOperator, SBPxi: super::operators::SbpOperator
             detj_dxi_dy,
             detj_deta_dx,
             detj_deta_dy,
-            operatorxi: std::marker::PhantomData,
-            operatoreta: std::marker::PhantomData,
         })
     }
 }
diff --git a/sbp/src/maxwell.rs b/sbp/src/maxwell.rs
index 3a00697..de2519c 100644
--- a/sbp/src/maxwell.rs
+++ b/sbp/src/maxwell.rs
@@ -78,19 +78,21 @@ pub struct System<SBP: SbpOperator> {
     sys: (Field, Field),
     wb: WorkBuffers,
     grid: Grid,
-    metrics: Metrics<SBP, SBP>,
+    metrics: Metrics,
+    op: SBP,
 }
 
 impl<SBP: SbpOperator> System<SBP> {
-    pub fn new(x: Array2<Float>, y: Array2<Float>) -> Self {
+    pub fn new(x: Array2<Float>, y: Array2<Float>, op: SBP) -> Self {
         assert_eq!(x.shape(), y.shape());
         let ny = x.shape()[0];
         let nx = x.shape()[1];
 
         let grid = Grid::new(x, y).unwrap();
-        let metrics = grid.metrics().unwrap();
+        let metrics = grid.metrics(op, op).unwrap();
 
         Self {
+            op,
             sys: (Field::new(ny, nx), Field::new(ny, nx)),
             grid,
             metrics,
@@ -115,16 +117,20 @@ impl<SBP: SbpOperator> System<SBP> {
     }
 
     pub fn advance(&mut self, dt: Float) {
-        fn rhs_adaptor<SBP: SbpOperator>(
-            fut: &mut Field,
-            prev: &Field,
-            _time: Float,
-            c: &(&Grid, &Metrics<SBP, SBP>),
-            m: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
-        ) {
+        let op = self.op;
+        let rhs_adaptor = move |fut: &mut Field,
+                                prev: &Field,
+                                _time: Float,
+                                c: &(&Grid, &Metrics),
+                                m: &mut (
+            Array2<Float>,
+            Array2<Float>,
+            Array2<Float>,
+            Array2<Float>,
+        )| {
             let (grid, metrics) = c;
-            RHS(fut, prev, grid, metrics, m);
-        }
+            RHS(op, fut, prev, grid, metrics, m);
+        };
         let mut _time = 0.0;
         integrate::integrate::<integrate::Rk4, _, _, _, _>(
             rhs_adaptor,
@@ -143,16 +149,20 @@ impl<SBP: SbpOperator> System<SBP> {
 impl<UO: UpwindOperator> System<UO> {
     /// Using artificial dissipation with the upwind operator
     pub fn advance_upwind(&mut self, dt: Float) {
-        fn rhs_adaptor<UO: UpwindOperator>(
-            fut: &mut Field,
-            prev: &Field,
-            _time: Float,
-            c: &(&Grid, &Metrics<UO, UO>),
-            m: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
-        ) {
+        let op = self.op;
+        let rhs_adaptor = move |fut: &mut Field,
+                                prev: &Field,
+                                _time: Float,
+                                c: &(&Grid, &Metrics),
+                                m: &mut (
+            Array2<Float>,
+            Array2<Float>,
+            Array2<Float>,
+            Array2<Float>,
+        )| {
             let (grid, metrics) = c;
-            RHS_upwind(fut, prev, grid, metrics, m);
-        }
+            RHS_upwind(op, fut, prev, grid, metrics, m);
+        };
         let mut _time = 0.0;
         integrate::integrate::<integrate::Rk4, _, _, _, _>(
             rhs_adaptor,
@@ -196,13 +206,14 @@ fn gaussian(x: Float, x0: Float, y: Float, y0: Float) -> Float {
 ///
 /// This is used both in fluxes and SAT terms
 fn RHS<SBP: SbpOperator>(
+    op: SBP,
     k: &mut Field,
     y: &Field,
     _grid: &Grid,
-    metrics: &Metrics<SBP, SBP>,
+    metrics: &Metrics,
     tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
 ) {
-    fluxes(k, y, metrics, tmp);
+    fluxes(op, k, y, metrics, tmp);
 
     let boundaries = BoundaryTerms {
         north: Boundary::This,
@@ -210,7 +221,7 @@ fn RHS<SBP: SbpOperator>(
         west: Boundary::This,
         east: Boundary::This,
     };
-    SAT_characteristics(k, y, metrics, &boundaries);
+    SAT_characteristics(op, k, y, metrics, &boundaries);
 
     azip!((k in &mut k.0,
                     &detj in &metrics.detj.broadcast((3, y.ny(), y.nx())).unwrap()) {
@@ -220,14 +231,15 @@ fn RHS<SBP: SbpOperator>(
 
 #[allow(non_snake_case)]
 fn RHS_upwind<UO: UpwindOperator>(
+    op: UO,
     k: &mut Field,
     y: &Field,
     _grid: &Grid,
-    metrics: &Metrics<UO, UO>,
+    metrics: &Metrics,
     tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
 ) {
-    fluxes(k, y, metrics, tmp);
-    dissipation(k, y, metrics, tmp);
+    fluxes(op, k, y, metrics, tmp);
+    dissipation(op, k, y, metrics, tmp);
 
     let boundaries = BoundaryTerms {
         north: Boundary::This,
@@ -235,7 +247,7 @@ fn RHS_upwind<UO: UpwindOperator>(
         west: Boundary::This,
         east: Boundary::This,
     };
-    SAT_characteristics(k, y, metrics, &boundaries);
+    SAT_characteristics(op, k, y, metrics, &boundaries);
 
     azip!((k in &mut k.0,
                     &detj in &metrics.detj.broadcast((3, y.ny(), y.nx())).unwrap()) {
@@ -244,9 +256,10 @@ fn RHS_upwind<UO: UpwindOperator>(
 }
 
 fn fluxes<SBP: SbpOperator>(
+    op: SBP,
     k: &mut Field,
     y: &Field,
-    metrics: &Metrics<SBP, SBP>,
+    metrics: &Metrics,
     tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
 ) {
     // ex = hz_y
@@ -256,14 +269,14 @@ fn fluxes<SBP: SbpOperator>(
                         &hz in &y.hz())
             *a = dxi_dy * hz
         );
-        SBP::diffxi(tmp.0.view(), tmp.1.view_mut());
+        op.diffxi(tmp.0.view(), tmp.1.view_mut());
 
         ndarray::azip!((b in &mut tmp.2,
                         &deta_dy in &metrics.detj_deta_dy,
                         &hz in &y.hz())
             *b = deta_dy * hz
         );
-        SBP::diffeta(tmp.2.view(), tmp.3.view_mut());
+        op.diffeta(tmp.2.view(), tmp.3.view_mut());
 
         ndarray::azip!((flux in &mut k.ex_mut(), &ax in &tmp.1, &by in &tmp.3)
             *flux = ax + by
@@ -279,7 +292,7 @@ fn fluxes<SBP: SbpOperator>(
                         &ey in &y.ey())
             *a = dxi_dx * -ey + dxi_dy * ex
         );
-        SBP::diffxi(tmp.0.view(), tmp.1.view_mut());
+        op.diffxi(tmp.0.view(), tmp.1.view_mut());
 
         ndarray::azip!((b in &mut tmp.2,
                         &deta_dx in &metrics.detj_deta_dx,
@@ -288,7 +301,7 @@ fn fluxes<SBP: SbpOperator>(
                         &ey in &y.ey())
             *b = deta_dx * -ey + deta_dy * ex
         );
-        SBP::diffeta(tmp.2.view(), tmp.3.view_mut());
+        op.diffeta(tmp.2.view(), tmp.3.view_mut());
 
         ndarray::azip!((flux in &mut k.hz_mut(), &ax in &tmp.1, &by in &tmp.3)
             *flux = ax + by
@@ -302,14 +315,14 @@ fn fluxes<SBP: SbpOperator>(
                         &hz in &y.hz())
             *a = dxi_dx * -hz
         );
-        SBP::diffxi(tmp.0.view(), tmp.1.view_mut());
+        op.diffxi(tmp.0.view(), tmp.1.view_mut());
 
         azip!((b in &mut tmp.2,
                         &deta_dx in &metrics.detj_deta_dx,
                         &hz in &y.hz())
             *b = deta_dx * -hz
         );
-        SBP::diffeta(tmp.2.view(), tmp.3.view_mut());
+        op.diffeta(tmp.2.view(), tmp.3.view_mut());
 
         azip!((flux in &mut k.ey_mut(), &ax in &tmp.1, &by in &tmp.3)
             *flux = ax + by
@@ -318,9 +331,10 @@ fn fluxes<SBP: SbpOperator>(
 }
 
 fn dissipation<UO: UpwindOperator>(
+    op: UO,
     k: &mut Field,
     y: &Field,
-    metrics: &Metrics<UO, UO>,
+    metrics: &Metrics,
     tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
 ) {
     // ex component
@@ -333,7 +347,7 @@ fn dissipation<UO: UpwindOperator>(
             let r = Float::hypot(kx, ky);
             *a = ky*ky/r * ex + -kx*ky/r*ey;
         });
-        UO::dissxi(tmp.0.view(), tmp.1.view_mut());
+        op.dissxi(tmp.0.view(), tmp.1.view_mut());
 
         ndarray::azip!((b in &mut tmp.2,
                     &kx in &metrics.detj_deta_dx,
@@ -343,7 +357,7 @@ fn dissipation<UO: UpwindOperator>(
             let r = Float::hypot(kx, ky);
             *b = ky*ky/r * ex + -kx*ky/r*ey;
         });
-        UO::disseta(tmp.2.view(), tmp.3.view_mut());
+        op.disseta(tmp.2.view(), tmp.3.view_mut());
 
         ndarray::azip!((flux in &mut k.ex_mut(), &ax in &tmp.1, &by in &tmp.3)
             *flux += ax + by
@@ -359,7 +373,7 @@ fn dissipation<UO: UpwindOperator>(
             let r = Float::hypot(kx, ky);
             *a = r * hz;
         });
-        UO::dissxi(tmp.0.view(), tmp.1.view_mut());
+        op.dissxi(tmp.0.view(), tmp.1.view_mut());
 
         ndarray::azip!((b in &mut tmp.2,
                         &kx in &metrics.detj_deta_dx,
@@ -368,7 +382,7 @@ fn dissipation<UO: UpwindOperator>(
             let r = Float::hypot(kx, ky);
             *b = r * hz;
         });
-        UO::disseta(tmp.2.view(), tmp.3.view_mut());
+        op.disseta(tmp.2.view(), tmp.3.view_mut());
 
         ndarray::azip!((flux in &mut k.hz_mut(), &ax in &tmp.1, &by in &tmp.3)
             *flux += ax + by
@@ -385,7 +399,7 @@ fn dissipation<UO: UpwindOperator>(
             let r = Float::hypot(kx, ky);
             *a = -kx*ky/r * ex + kx*kx/r*ey;
         });
-        UO::dissxi(tmp.0.view(), tmp.1.view_mut());
+        op.dissxi(tmp.0.view(), tmp.1.view_mut());
 
         ndarray::azip!((b in &mut tmp.2,
                     &kx in &metrics.detj_deta_dx,
@@ -395,7 +409,7 @@ fn dissipation<UO: UpwindOperator>(
             let r = Float::hypot(kx, ky);
             *b = -kx*ky/r * ex + kx*kx/r*ey;
         });
-        UO::disseta(tmp.2.view(), tmp.3.view_mut());
+        op.disseta(tmp.2.view(), tmp.3.view_mut());
 
         ndarray::azip!((flux in &mut k.hz_mut(), &ax in &tmp.1, &by in &tmp.3)
             *flux += ax + by
@@ -419,9 +433,10 @@ pub struct BoundaryTerms {
 #[allow(non_snake_case)]
 /// Boundary conditions (SAT)
 fn SAT_characteristics<SBP: SbpOperator>(
+    op: SBP,
     k: &mut Field,
     y: &Field,
-    metrics: &Metrics<SBP, SBP>,
+    metrics: &Metrics,
     boundaries: &BoundaryTerms,
 ) {
     let ny = y.ny();
@@ -449,10 +464,10 @@ fn SAT_characteristics<SBP: SbpOperator>(
             Boundary::This => y.slice(s![.., .., 0]),
         };
         // East boundary
-        let hinv = if SBP::is_h2() {
-            (nx - 2) as Float / SBP::h()[0]
+        let hinv = if op.is_h2() {
+            (nx - 2) as Float / op.h()[0]
         } else {
-            (nx - 1) as Float / SBP::h()[0]
+            (nx - 1) as Float / op.h()[0]
         };
         for ((((mut k, v), g), &kx), &ky) in k
             .slice_mut(s![.., .., nx - 1])
@@ -487,10 +502,10 @@ fn SAT_characteristics<SBP: SbpOperator>(
         let g = match boundaries.east {
             Boundary::This => y.slice(s![.., .., nx - 1]),
         };
-        let hinv = if SBP::is_h2() {
-            (nx - 2) as Float / SBP::h()[0]
+        let hinv = if op.is_h2() {
+            (nx - 2) as Float / op.h()[0]
         } else {
-            (nx - 1) as Float / SBP::h()[0]
+            (nx - 1) as Float / op.h()[0]
         };
         for ((((mut k, v), g), &kx), &ky) in k
             .slice_mut(s![.., .., 0])
@@ -530,10 +545,10 @@ fn SAT_characteristics<SBP: SbpOperator>(
         let g = match boundaries.north {
             Boundary::This => y.slice(s![.., 0, ..]),
         };
-        let hinv = if SBP::is_h2() {
-            (ny - 2) as Float / SBP::h()[0]
+        let hinv = if op.is_h2() {
+            (ny - 2) as Float / op.h()[0]
         } else {
-            (ny - 1) as Float / SBP::h()[0]
+            (ny - 1) as Float / op.h()[0]
         };
         for ((((mut k, v), g), &kx), &ky) in k
             .slice_mut(s![.., ny - 1, ..])
@@ -567,10 +582,10 @@ fn SAT_characteristics<SBP: SbpOperator>(
         let g = match boundaries.south {
             Boundary::This => y.slice(s![.., ny - 1, ..]),
         };
-        let hinv = if SBP::is_h2() {
-            (ny - 2) as Float / SBP::h()[0]
+        let hinv = if op.is_h2() {
+            (ny - 2) as Float / op.h()[0]
         } else {
-            (ny - 1) as Float / SBP::h()[0]
+            (ny - 1) as Float / op.h()[0]
         };
         for ((((mut k, v), g), &kx), &ky) in k
             .slice_mut(s![.., 0, ..])
diff --git a/sbp/src/operators.rs b/sbp/src/operators.rs
index 9f5d040..5ad2742 100644
--- a/sbp/src/operators.rs
+++ b/sbp/src/operators.rs
@@ -5,33 +5,33 @@ use crate::Float;
 
 use ndarray::{ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2};
 
-pub trait SbpOperator: Send + Sync {
-    fn diff1d(prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>);
-    fn diffxi(prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
+pub trait SbpOperator: Send + Sync + Copy {
+    fn diff1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>);
+    fn diffxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
         assert_eq!(prev.shape(), fut.shape());
         for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) {
-            Self::diff1d(r0, r1);
+            self.diff1d(r0, r1);
         }
     }
-    fn diffeta(prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
-        Self::diffxi(prev.reversed_axes(), fut.reversed_axes())
+    fn diffeta(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
+        self.diffxi(prev.reversed_axes(), fut.reversed_axes())
     }
-    fn h() -> &'static [Float];
-    fn is_h2() -> bool {
+    fn h(&self) -> &'static [Float];
+    fn is_h2(&self) -> bool {
         false
     }
 }
 
 pub trait UpwindOperator: SbpOperator {
-    fn diss1d(prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>);
-    fn dissxi(prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
+    fn diss1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>);
+    fn dissxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
         assert_eq!(prev.shape(), fut.shape());
         for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) {
-            Self::diss1d(r0, r1);
+            self.diss1d(r0, r1);
         }
     }
-    fn disseta(prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
-        Self::dissxi(prev.reversed_axes(), fut.reversed_axes())
+    fn disseta(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
+        self.dissxi(prev.reversed_axes(), fut.reversed_axes())
     }
 }
 
@@ -133,6 +133,7 @@ pub(crate) mod testing {
     }
 
     pub(crate) fn check_operator_on<SBP, F, FX, FY>(
+        op: SBP,
         n: (usize, usize),
         f: F,
         dfdx: FX,
@@ -148,11 +149,11 @@ pub(crate) mod testing {
         let x = grid_eval(n, f);
 
         y.fill(0.0);
-        SBP::diffxi(x.view(), y.view_mut());
+        op.diffxi(x.view(), y.view_mut());
         approx::assert_abs_diff_eq!(&y, &grid_eval(n, dfdx), epsilon = eps);
 
         y.fill(0.0);
-        SBP::diffeta(x.view(), y.view_mut());
+        op.diffeta(x.view(), y.view_mut());
         approx::assert_abs_diff_eq!(&y, &grid_eval(n, dfdy), epsilon = eps);
     }
 }
diff --git a/sbp/src/operators/traditional4.rs b/sbp/src/operators/traditional4.rs
index 4d4eebb..4dd9917 100644
--- a/sbp/src/operators/traditional4.rs
+++ b/sbp/src/operators/traditional4.rs
@@ -2,8 +2,8 @@ use super::SbpOperator;
 use crate::Float;
 use ndarray::{ArrayView1, ArrayViewMut1};
 
-#[derive(Debug)]
-pub struct SBP4 {}
+#[derive(Debug, Copy, Clone)]
+pub struct SBP4;
 
 impl SBP4 {
     #[rustfmt::skip]
@@ -24,7 +24,7 @@ impl SBP4 {
 }
 
 impl SbpOperator for SBP4 {
-    fn diff1d(prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
+    fn diff1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
         super::diff_op_1d(
             ndarray::arr2(Self::BLOCK).view(),
             ndarray::arr1(Self::DIAG).view(),
@@ -35,7 +35,7 @@ impl SbpOperator for SBP4 {
         )
     }
 
-    fn h() -> &'static [Float] {
+    fn h(&self) -> &'static [Float] {
         Self::HBLOCK
     }
 }
@@ -47,15 +47,24 @@ fn test_trad4() {
     let nx = 20;
     let ny = 13;
 
-    check_operator_on::<SBP4, _, _, _>((ny, nx), |x, y| x + 2.0 * y, |_, _| 1.0, |_, _| 2.0, 1e-4);
-    check_operator_on::<SBP4, _, _, _>(
+    check_operator_on(
+        SBP4,
+        (ny, nx),
+        |x, y| x + 2.0 * y,
+        |_, _| 1.0,
+        |_, _| 2.0,
+        1e-4,
+    );
+    check_operator_on(
+        SBP4,
         (ny, nx),
         |x, y| x * x + 2.0 * x * y + 3.0 * y * y,
         |x, y| 2.0 * x + 2.0 * y,
         |x, y| 2.0 * x + 6.0 * y,
         1e-3,
     );
-    check_operator_on::<SBP4, _, _, _>(
+    check_operator_on(
+        SBP4,
         (ny, nx),
         |x, y| x.powi(3) + 2.0 * x.powi(2) * y + 3.0 * x * y.powi(2) + 4.0 * y.powi(3),
         |x, y| 3.0 * x.powi(2) + 4.0 * x * y + 3.0 * y.powi(2),
diff --git a/sbp/src/operators/traditional8.rs b/sbp/src/operators/traditional8.rs
index d5f0c6f..8fe19b6 100644
--- a/sbp/src/operators/traditional8.rs
+++ b/sbp/src/operators/traditional8.rs
@@ -2,8 +2,8 @@ use super::SbpOperator;
 use crate::Float;
 use ndarray::{ArrayView1, ArrayViewMut1};
 
-#[derive(Debug)]
-pub struct SBP8 {}
+#[derive(Debug, Copy, Clone)]
+pub struct SBP8;
 
 impl SBP8 {
     #[rustfmt::skip]
@@ -28,7 +28,7 @@ impl SBP8 {
 }
 
 impl SbpOperator for SBP8 {
-    fn diff1d(prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
+    fn diff1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
         super::diff_op_1d(
             ndarray::arr2(Self::BLOCK).view(),
             ndarray::arr1(Self::DIAG).view(),
@@ -39,7 +39,7 @@ impl SbpOperator for SBP8 {
         )
     }
 
-    fn h() -> &'static [Float] {
+    fn h(&self) -> &'static [Float] {
         Self::HBLOCK
     }
 }
@@ -51,7 +51,8 @@ fn test_trad8() {
     let ny = 16;
 
     // Order one polynomial
-    check_operator_on::<SBP8, _, _, _>(
+    check_operator_on(
+        SBP8,
         (ny, nx),
         |x, y| x + 2.0 * y,
         |_x, _y| 1.0,
@@ -60,31 +61,35 @@ fn test_trad8() {
     );
 
     // Order two polynomial
-    check_operator_on::<SBP8, _, _, _>(
+    check_operator_on(
+        SBP8,
         (ny, nx),
         |x, y| x * x + 0.5 * y * y,
         |x, _y| 2.0 * x,
         |_x, y| y,
         1e-4,
     );
-    check_operator_on::<SBP8, _, _, _>((ny, nx), |x, y| x * y, |_x, y| y, |x, _y| x, 1e-4);
+    check_operator_on(SBP8, (ny, nx), |x, y| x * y, |_x, y| y, |x, _y| x, 1e-4);
 
     // Order three polynomials
-    check_operator_on::<SBP8, _, _, _>(
+    check_operator_on(
+        SBP8,
         (ny, nx),
         |x, y| x * x * x + y * y * y / 6.0,
         |x, _y| 3.0 * x * x,
         |_x, y| y * y / 2.0,
         1e-4,
     );
-    check_operator_on::<SBP8, _, _, _>(
+    check_operator_on(
+        SBP8,
         (ny, nx),
         |x, y| x * x * y + x * y * y / 2.0,
         |x, y| 2.0 * x * y + y * y / 2.0,
         |x, y| x * x + x * y,
         1e-4,
     );
-    check_operator_on::<SBP8, _, _, _>(
+    check_operator_on(
+        SBP8,
         (ny, nx),
         |x, y| x.powi(3) + 2.0 * x.powi(2) * y + 3.0 * x * y.powi(2) + 4.0 * y.powi(3),
         |x, y| 3.0 * x.powi(2) + 4.0 * x * y + 3.0 * y.powi(2),
@@ -93,7 +98,8 @@ fn test_trad8() {
     );
 
     // Order four polynomials
-    check_operator_on::<SBP8, _, _, _>(
+    check_operator_on(
+        SBP8,
         (ny, nx),
         |x, y| x.powi(4) + x.powi(3) * y + x.powi(2) * y.powi(2) + x * y.powi(3) + y.powi(4),
         |x, y| 4.0 * x.powi(3) + 3.0 * x.powi(2) * y + 2.0 * x * y.powi(2) + y.powi(3),
diff --git a/sbp/src/operators/upwind4.rs b/sbp/src/operators/upwind4.rs
index 8efa7c6..834fbeb 100644
--- a/sbp/src/operators/upwind4.rs
+++ b/sbp/src/operators/upwind4.rs
@@ -2,8 +2,8 @@ use super::{SbpOperator, UpwindOperator};
 use crate::Float;
 use ndarray::{ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2, Axis};
 
-#[derive(Debug)]
-pub struct Upwind4 {}
+#[derive(Debug, Copy, Clone)]
+pub struct Upwind4;
 
 /// Simdtype used in diff_simd_col and diff_simd_row
 #[cfg(feature = "f32")]
@@ -276,7 +276,7 @@ impl Upwind4 {
 }
 
 impl SbpOperator for Upwind4 {
-    fn diff1d(prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
+    fn diff1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
         super::diff_op_1d(
             ndarray::arr2(Self::BLOCK).view(),
             ndarray::arr1(Self::DIAG).view(),
@@ -286,7 +286,7 @@ impl SbpOperator for Upwind4 {
             fut,
         )
     }
-    fn diffxi(prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
+    fn diffxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
         assert_eq!(prev.shape(), fut.shape());
         assert!(prev.shape()[1] >= 2 * Self::BLOCK.len());
 
@@ -300,14 +300,14 @@ impl SbpOperator for Upwind4 {
             ([_, _], [_, _]) => {
                 // Fallback, work row by row
                 for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) {
-                    Self::diff1d(r0, r1);
+                    Self.diff1d(r0, r1);
                 }
             }
             _ => unreachable!("Should only be two elements in the strides vectors"),
         }
     }
 
-    fn h() -> &'static [Float] {
+    fn h(&self) -> &'static [Float] {
         Self::HBLOCK
     }
 }
@@ -326,14 +326,14 @@ fn upwind4_test() {
         target[i] = 1.0;
     }
     res.fill(0.0);
-    Upwind4::diff1d(source.view(), res.view_mut());
+    Upwind4.diff1d(source.view(), res.view_mut());
     approx::assert_abs_diff_eq!(&res, &target, epsilon = 1e-4);
     {
         let source = source.to_owned().insert_axis(ndarray::Axis(0));
         let mut res = res.to_owned().insert_axis(ndarray::Axis(0));
         let target = target.to_owned().insert_axis(ndarray::Axis(0));
         res.fill(0.0);
-        Upwind4::diffxi(source.view(), res.view_mut());
+        Upwind4.diffxi(source.view(), res.view_mut());
         approx::assert_abs_diff_eq!(&res, &target, epsilon = 1e-2);
     }
 
@@ -342,7 +342,7 @@ fn upwind4_test() {
         let target = Array2::from_shape_fn((nx, 8), |(i, _)| target[i]);
         let mut res = Array2::zeros((nx, 8));
         res.fill(0.0);
-        Upwind4::diffeta(source.view(), res.view_mut());
+        Upwind4.diffeta(source.view(), res.view_mut());
         approx::assert_abs_diff_eq!(&res.to_owned(), &target.to_owned(), epsilon = 1e-2);
     }
 
@@ -352,14 +352,14 @@ fn upwind4_test() {
         target[i] = 2.0 * x;
     }
     res.fill(0.0);
-    Upwind4::diff1d(source.view(), res.view_mut());
+    Upwind4.diff1d(source.view(), res.view_mut());
     approx::assert_abs_diff_eq!(&res, &target, epsilon = 1e-4);
     {
         let source = source.to_owned().insert_axis(ndarray::Axis(0));
         let mut res = res.to_owned().insert_axis(ndarray::Axis(0));
         let target = target.to_owned().insert_axis(ndarray::Axis(0));
         res.fill(0.0);
-        Upwind4::diffxi(source.view(), res.view_mut());
+        Upwind4.diffxi(source.view(), res.view_mut());
         approx::assert_abs_diff_eq!(&res, &target, epsilon = 1e-2);
     }
 
@@ -368,7 +368,7 @@ fn upwind4_test() {
         let target = Array2::from_shape_fn((nx, 8), |(i, _)| target[i]);
         let mut res = Array2::zeros((nx, 8));
         res.fill(0.0);
-        Upwind4::diffeta(source.view(), res.view_mut());
+        Upwind4.diffeta(source.view(), res.view_mut());
         approx::assert_abs_diff_eq!(&res.to_owned(), &target.to_owned(), epsilon = 1e-2);
     }
 
@@ -378,7 +378,7 @@ fn upwind4_test() {
         target[i] = 3.0 * x * x;
     }
     res.fill(0.0);
-    Upwind4::diff1d(source.view(), res.view_mut());
+    Upwind4.diff1d(source.view(), res.view_mut());
     approx::assert_abs_diff_eq!(&res, &target, epsilon = 1e-2);
 
     {
@@ -386,7 +386,7 @@ fn upwind4_test() {
         let mut res = res.to_owned().insert_axis(ndarray::Axis(0));
         let target = target.to_owned().insert_axis(ndarray::Axis(0));
         res.fill(0.0);
-        Upwind4::diffxi(source.view(), res.view_mut());
+        Upwind4.diffxi(source.view(), res.view_mut());
         approx::assert_abs_diff_eq!(&res, &target, epsilon = 1e-2);
     }
 
@@ -395,13 +395,13 @@ fn upwind4_test() {
         let target = Array2::from_shape_fn((nx, 8), |(i, _)| target[i]);
         let mut res = Array2::zeros((nx, 8));
         res.fill(0.0);
-        Upwind4::diffeta(source.view(), res.view_mut());
+        Upwind4.diffeta(source.view(), res.view_mut());
         approx::assert_abs_diff_eq!(&res.to_owned(), &target.to_owned(), epsilon = 1e-2);
     }
 }
 
 impl UpwindOperator for Upwind4 {
-    fn diss1d(prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
+    fn diss1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
         super::diff_op_1d(
             ndarray::arr2(Self::DISS_BLOCK).view(),
             ndarray::arr1(Self::DISS_DIAG).view(),
@@ -411,7 +411,7 @@ impl UpwindOperator for Upwind4 {
             fut,
         )
     }
-    fn dissxi(prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
+    fn dissxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
         assert_eq!(prev.shape(), fut.shape());
         assert!(prev.shape()[1] >= 2 * Self::BLOCK.len());
 
@@ -425,7 +425,7 @@ impl UpwindOperator for Upwind4 {
             ([_, _], [_, _]) => {
                 // Fallback, work row by row
                 for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) {
-                    Self::diss1d(r0, r1);
+                    Self.diss1d(r0, r1);
                 }
             }
             _ => unreachable!("Should only be two elements in the strides vectors"),
@@ -440,28 +440,32 @@ fn upwind4_test2() {
     let nx = 32;
     let ny = 16;
 
-    check_operator_on::<Upwind4, _, _, _>(
+    check_operator_on(
+        Upwind4,
         (ny, nx),
         |x, y| x + 2.0 * y,
         |_, _| 1.0,
         |_, _| 2.0,
         1e-4,
     );
-    check_operator_on::<Upwind4, _, _, _>(
+    check_operator_on(
+        Upwind4,
         (ny, nx),
         |x, y| x * x + 2.0 * x * y + 3.0 * y * y,
         |x, y| 2.0 * x + 2.0 * y,
         |x, y| 2.0 * x + 6.0 * y,
         1e-3,
     );
-    check_operator_on::<Upwind4, _, _, _>(
+    check_operator_on(
+        Upwind4,
         (ny, nx),
         |x, y| x.powi(3) + 2.0 * x.powi(2) * y + 3.0 * x * y.powi(2) + 4.0 * y.powi(3),
         |x, y| 3.0 * x.powi(2) + 4.0 * x * y + 3.0 * y.powi(2),
         |x, y| 2.0 * x.powi(2) + 6.0 * x * y + 12.0 * y.powi(2),
         1e-1,
     );
-    check_operator_on::<Upwind4, _, _, _>(
+    check_operator_on(
+        Upwind4,
         (32, 32),
         |x, y| x.powi(3) + 2.0 * x.powi(2) * y + 3.0 * x * y.powi(2) + 4.0 * y.powi(3),
         |x, y| 3.0 * x.powi(2) + 4.0 * x * y + 3.0 * y.powi(2),
diff --git a/sbp/src/operators/upwind4h2.rs b/sbp/src/operators/upwind4h2.rs
index b6857a3..d9304ba 100644
--- a/sbp/src/operators/upwind4h2.rs
+++ b/sbp/src/operators/upwind4h2.rs
@@ -2,8 +2,8 @@ use super::{SbpOperator, UpwindOperator};
 use crate::Float;
 use ndarray::{ArrayView1, ArrayViewMut1};
 
-#[derive(Debug)]
-pub struct Upwind4h2 {}
+#[derive(Debug, Copy, Clone)]
+pub struct Upwind4h2;
 
 impl Upwind4h2 {
     #[rustfmt::skip]
@@ -37,7 +37,7 @@ impl Upwind4h2 {
 }
 
 impl SbpOperator for Upwind4h2 {
-    fn diff1d(prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
+    fn diff1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
         super::diff_op_1d(
             ndarray::arr2(Self::BLOCK).view(),
             ndarray::arr1(Self::DIAG).view(),
@@ -48,10 +48,10 @@ impl SbpOperator for Upwind4h2 {
         )
     }
 
-    fn h() -> &'static [Float] {
+    fn h(&self) -> &'static [Float] {
         Self::HBLOCK
     }
-    fn is_h2() -> bool {
+    fn is_h2(&self) -> bool {
         true
     }
 }
@@ -64,25 +64,25 @@ fn upwind4h2_test() {
 
     let mut res = ndarray::Array1::zeros(nx);
 
-    Upwind4h2::diff1d(x.view(), res.view_mut());
+    Upwind4h2.diff1d(x.view(), res.view_mut());
     let ans = &x * 0.0 + 1.0;
     approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-4);
 
     res.fill(0.0);
     let y = &x * &x / 2.0;
-    Upwind4h2::diff1d(y.view(), res.view_mut());
+    Upwind4h2.diff1d(y.view(), res.view_mut());
     let ans = &x;
     approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-4);
 
     res.fill(0.0);
     let y = &x * &x * &x / 3.0;
-    Upwind4h2::diff1d(y.view(), res.view_mut());
+    Upwind4h2.diff1d(y.view(), res.view_mut());
     let ans = &x * &x;
     approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-2);
 }
 
 impl UpwindOperator for Upwind4h2 {
-    fn diss1d(prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
+    fn diss1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
         super::diff_op_1d(
             ndarray::arr2(Self::DISS_BLOCK).view(),
             ndarray::arr1(Self::DISS_DIAG).view(),
diff --git a/sbp/src/operators/upwind9.rs b/sbp/src/operators/upwind9.rs
index 18009f9..9c11143 100644
--- a/sbp/src/operators/upwind9.rs
+++ b/sbp/src/operators/upwind9.rs
@@ -2,8 +2,8 @@ use super::{SbpOperator, UpwindOperator};
 use crate::Float;
 use ndarray::{ArrayView1, ArrayViewMut1};
 
-#[derive(Debug)]
-pub struct Upwind9 {}
+#[derive(Debug, Copy, Clone)]
+pub struct Upwind9;
 
 impl Upwind9 {
     #[rustfmt::skip]
@@ -45,7 +45,7 @@ impl Upwind9 {
 }
 
 impl SbpOperator for Upwind9 {
-    fn diff1d(prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
+    fn diff1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
         super::diff_op_1d(
             ndarray::arr2(Self::BLOCK).view(),
             ndarray::arr1(Self::DIAG).view(),
@@ -56,13 +56,13 @@ impl SbpOperator for Upwind9 {
         )
     }
 
-    fn h() -> &'static [Float] {
+    fn h(&self) -> &'static [Float] {
         Self::HBLOCK
     }
 }
 
 impl UpwindOperator for Upwind9 {
-    fn diss1d(prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
+    fn diss1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
         super::diff_op_1d(
             ndarray::arr2(Self::DISS_BLOCK).view(),
             ndarray::arr1(Self::DISS_DIAG).view(),
@@ -81,7 +81,8 @@ fn test_upwind9() {
     let ny = 16;
 
     // Order one polynomial
-    check_operator_on::<Upwind9, _, _, _>(
+    check_operator_on(
+        Upwind9,
         (ny, nx),
         |x, y| x + 2.0 * y,
         |_x, _y| 1.0,
@@ -90,31 +91,35 @@ fn test_upwind9() {
     );
 
     // Order two polynomial
-    check_operator_on::<Upwind9, _, _, _>(
+    check_operator_on(
+        Upwind9,
         (ny, nx),
         |x, y| x * x + 0.5 * y * y,
         |x, _y| 2.0 * x,
         |_x, y| y,
         1e-4,
     );
-    check_operator_on::<Upwind9, _, _, _>((ny, nx), |x, y| x * y, |_x, y| y, |x, _y| x, 1e-4);
+    check_operator_on(Upwind9, (ny, nx), |x, y| x * y, |_x, y| y, |x, _y| x, 1e-4);
 
     // Order three polynomials
-    check_operator_on::<Upwind9, _, _, _>(
+    check_operator_on(
+        Upwind9,
         (ny, nx),
         |x, y| x * x * x + y * y * y / 6.0,
         |x, _y| 3.0 * x * x,
         |_x, y| y * y / 2.0,
         1e-4,
     );
-    check_operator_on::<Upwind9, _, _, _>(
+    check_operator_on(
+        Upwind9,
         (ny, nx),
         |x, y| x * x * y + x * y * y / 2.0,
         |x, y| 2.0 * x * y + y * y / 2.0,
         |x, y| x * x + x * y,
         1e-4,
     );
-    check_operator_on::<Upwind9, _, _, _>(
+    check_operator_on(
+        Upwind9,
         (ny, nx),
         |x, y| x.powi(3) + 2.0 * x.powi(2) * y + 3.0 * x * y.powi(2) + 4.0 * y.powi(3),
         |x, y| 3.0 * x.powi(2) + 4.0 * x * y + 3.0 * y.powi(2),
@@ -123,7 +128,8 @@ fn test_upwind9() {
     );
 
     // Order four polynomials
-    check_operator_on::<Upwind9, _, _, _>(
+    check_operator_on(
+        Upwind9,
         (ny, nx),
         |x, y| x.powi(4) + x.powi(3) * y + x.powi(2) * y.powi(2) + x * y.powi(3) + y.powi(4),
         |x, y| 4.0 * x.powi(3) + 3.0 * x.powi(2) * y + 2.0 * x * y.powi(2) + y.powi(3),
diff --git a/sbp/src/operators/upwind9h2.rs b/sbp/src/operators/upwind9h2.rs
index 06defed..b65296d 100644
--- a/sbp/src/operators/upwind9h2.rs
+++ b/sbp/src/operators/upwind9h2.rs
@@ -2,8 +2,8 @@ use super::{SbpOperator, UpwindOperator};
 use crate::Float;
 use ndarray::{ArrayView1, ArrayViewMut1};
 
-#[derive(Debug)]
-pub struct Upwind9h2 {}
+#[derive(Debug, Copy, Clone)]
+pub struct Upwind9h2;
 
 impl Upwind9h2 {
     #[rustfmt::skip]
@@ -45,7 +45,7 @@ impl Upwind9h2 {
 }
 
 impl SbpOperator for Upwind9h2 {
-    fn diff1d(prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
+    fn diff1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
         super::diff_op_1d(
             ndarray::arr2(Self::BLOCK).view(),
             ndarray::arr1(Self::DIAG).view(),
@@ -56,10 +56,10 @@ impl SbpOperator for Upwind9h2 {
         )
     }
 
-    fn h() -> &'static [Float] {
+    fn h(&self) -> &'static [Float] {
         Self::HBLOCK
     }
-    fn is_h2() -> bool {
+    fn is_h2(&self) -> bool {
         true
     }
 }
@@ -72,25 +72,25 @@ fn upwind9h2_test() {
 
     let mut res = ndarray::Array1::zeros(nx);
 
-    Upwind9h2::diff1d(x.view(), res.view_mut());
+    Upwind9h2.diff1d(x.view(), res.view_mut());
     let ans = &x * 0.0 + 1.0;
     approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-4);
 
     res.fill(0.0);
     let y = &x * &x / 2.0;
-    Upwind9h2::diff1d(y.view(), res.view_mut());
+    Upwind9h2.diff1d(y.view(), res.view_mut());
     let ans = &x;
     approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-4);
 
     res.fill(0.0);
     let y = &x * &x * &x / 3.0;
-    Upwind9h2::diff1d(y.view(), res.view_mut());
+    Upwind9h2.diff1d(y.view(), res.view_mut());
     let ans = &x * &x;
     approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-2);
 }
 
 impl UpwindOperator for Upwind9h2 {
-    fn diss1d(prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
+    fn diss1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
         super::diff_op_1d(
             ndarray::arr2(Self::DISS_BLOCK).view(),
             ndarray::arr1(Self::DISS_DIAG).view(),
diff --git a/webfront/lib.rs b/webfront/lib.rs
index 13f3a9a..aa44711 100644
--- a/webfront/lib.rs
+++ b/webfront/lib.rs
@@ -20,7 +20,7 @@ impl MaxwellUniverse {
     pub fn new(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(maxwell::System::new(x, y))
+        Self(maxwell::System::new(x, y, operators::Upwind4))
     }
 
     pub fn init(&mut self, x0: f32, y0: f32) {
@@ -53,7 +53,7 @@ pub struct EulerUniverse(euler::System<operators::Upwind4>);
 
 impl EulerUniverse {
     pub fn new(x: ndarray::Array2<f32>, y: ndarray::Array2<f32>) -> Self {
-        Self(euler::System::new(x, y))
+        Self(euler::System::new(x, y, operators::Upwind4))
     }
 }
 
@@ -63,7 +63,7 @@ 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(euler::System::new(x, y))
+        Self(euler::System::new(x, y, operators::Upwind4))
     }
 
     pub fn init(&mut self, x0: f32, y0: f32) {