split euler and maxwell to separate crates

maxwell/Cargo.toml

@@ -0,0 +1,17 @@
+[package]
+name = "maxwell"
+version = "0.1.0"
+authors = ["Magnus Ulimoen <flymagnus@gmail.com>"]
+edition = "2018"
+
+[dependencies]
+ndarray = "0.13.1"
+sbp = { path = "../sbp" }
+
+[dev-dependencies]
+criterion = "0.3.1"
+
+[[bench]]
+name = "bench"
+harness = false
+

maxwell/benches/bench.rs

@@ -0,0 +1,55 @@
+use criterion::{black_box, criterion_group, criterion_main, Criterion};
+use maxwell::System;
+use sbp::operators::{SbpOperator2d, Upwind4, UpwindOperator2d, SBP4};
+use sbp::Float;
+
+fn advance_system<SBP: SbpOperator2d>(universe: &mut System<SBP>, n: usize) {
+    for _ in 0..n {
+        universe.advance(0.01);
+    }
+}
+
+fn advance_system_upwind<UO: UpwindOperator2d>(universe: &mut System<UO>, n: usize) {
+    for _ in 0..n {
+        universe.advance_upwind(0.01);
+    }
+}
+
+fn performance_benchmark(c: &mut Criterion) {
+    let mut group = c.benchmark_group("MaxwellSystem");
+    group.sample_size(25);
+
+    let w = 40;
+    let h = 26;
+    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::new(x.clone(), y.clone(), Upwind4);
+    group.bench_function("advance", |b| {
+        b.iter(|| {
+            universe.set_gaussian(0.5, 0.5);
+            advance_system(&mut universe, black_box(20))
+        })
+    });
+
+    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);
+            advance_system_upwind(&mut universe, black_box(20))
+        })
+    });
+
+    let mut universe = System::new(x, y, SBP4);
+    group.bench_function("advance_trad4", |b| {
+        b.iter(|| {
+            universe.set_gaussian(0.5, 0.5);
+            advance_system(&mut universe, black_box(20))
+        })
+    });
+
+    group.finish();
+}
+
+criterion_group!(benches, performance_benchmark);
+criterion_main!(benches);

maxwell/src/lib.rs

@@ -0,0 +1,628 @@
+use ndarray::azip;
+use ndarray::prelude::*;
+use sbp::grid::{Grid, Metrics};
+use sbp::integrate;
+use sbp::operators::{SbpOperator2d, UpwindOperator2d};
+use sbp::Float;
+
+#[derive(Clone, Debug)]
+pub struct Field(pub(crate) Array3<Float>);
+
+impl std::ops::Deref for Field {
+    type Target = Array3<Float>;
+    fn deref(&self) -> &Self::Target {
+        &self.0
+    }
+}
+
+impl std::ops::DerefMut for Field {
+    fn deref_mut(&mut self) -> &mut Self::Target {
+        &mut self.0
+    }
+}
+
+impl<'a> std::convert::From<&'a Field> for ArrayView3<'a, Float> {
+    fn from(f: &'a Field) -> Self {
+        f.0.view()
+    }
+}
+impl<'a> std::convert::From<&'a mut Field> for ArrayViewMut3<'a, Float> {
+    fn from(f: &'a mut Field) -> Self {
+        f.0.view_mut()
+    }
+}
+
+impl Field {
+    pub fn new(height: usize, width: usize) -> Self {
+        let field = Array3::zeros((3, height, width));
+
+        Self(field)
+    }
+
+    pub fn nx(&self) -> usize {
+        self.0.shape()[2]
+    }
+    pub fn ny(&self) -> usize {
+        self.0.shape()[1]
+    }
+
+    pub fn ex(&self) -> ArrayView2<Float> {
+        self.slice(s![0, .., ..])
+    }
+    pub fn hz(&self) -> ArrayView2<Float> {
+        self.slice(s![1, .., ..])
+    }
+    pub fn ey(&self) -> ArrayView2<Float> {
+        self.slice(s![2, .., ..])
+    }
+
+    pub fn ex_mut(&mut self) -> ArrayViewMut2<Float> {
+        self.slice_mut(s![0, .., ..])
+    }
+    pub fn hz_mut(&mut self) -> ArrayViewMut2<Float> {
+        self.slice_mut(s![1, .., ..])
+    }
+    pub fn ey_mut(&mut self) -> ArrayViewMut2<Float> {
+        self.slice_mut(s![2, .., ..])
+    }
+
+    pub fn components_mut(
+        &mut self,
+    ) -> (
+        ArrayViewMut2<Float>,
+        ArrayViewMut2<Float>,
+        ArrayViewMut2<Float>,
+    ) {
+        self.0
+            .multi_slice_mut((s![0, .., ..], s![1, .., ..], s![2, .., ..]))
+    }
+}
+
+#[derive(Debug, Clone)]
+pub struct System<SBP: SbpOperator2d> {
+    sys: (Field, Field),
+    wb: WorkBuffers,
+    grid: Grid,
+    metrics: Metrics,
+    op: SBP,
+}
+
+impl<SBP: SbpOperator2d> System<SBP> {
+    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(&op).unwrap();
+
+        Self {
+            op,
+            sys: (Field::new(ny, nx), Field::new(ny, nx)),
+            grid,
+            metrics,
+            wb: WorkBuffers::new(ny, nx),
+        }
+    }
+
+    pub fn field(&self) -> &Field {
+        &self.sys.0
+    }
+
+    pub fn set_gaussian(&mut self, x0: Float, y0: Float) {
+        let (ex, hz, ey) = self.sys.0.components_mut();
+        ndarray::azip!(
+            (ex in ex, hz in hz, ey in ey,
+             &x in &self.grid.x(), &y in &self.grid.y())
+        {
+            *ex = 0.0;
+            *ey = 0.0;
+            *hz = gaussian(x, x0, y, y0)/32.0;
+        });
+    }
+
+    pub fn advance(&mut self, dt: Float) {
+        let op = &self.op;
+        let grid = &self.grid;
+        let metrics = &self.metrics;
+        let wb = &mut self.wb.tmp;
+        let rhs_adaptor = move |fut: &mut Field, prev: &Field, _time: Float| {
+            RHS(op, fut, prev, grid, metrics, wb);
+        };
+        let mut _time = 0.0;
+        integrate::integrate::<integrate::Rk4, _, _>(
+            rhs_adaptor,
+            &self.sys.0,
+            &mut self.sys.1,
+            &mut _time,
+            dt,
+            &mut self.wb.k,
+        );
+        std::mem::swap(&mut self.sys.0, &mut self.sys.1);
+    }
+}
+
+impl<UO: UpwindOperator2d> System<UO> {
+    /// Using artificial dissipation with the upwind operator
+    pub fn advance_upwind(&mut self, dt: Float) {
+        let op = &self.op;
+        let grid = &self.grid;
+        let metrics = &self.metrics;
+        let wb = &mut self.wb.tmp;
+        let rhs_adaptor = move |fut: &mut Field, prev: &Field, _time: Float| {
+            RHS_upwind(op, fut, prev, grid, metrics, wb);
+        };
+        let mut _time = 0.0;
+        integrate::integrate::<integrate::Rk4, _, _>(
+            rhs_adaptor,
+            &self.sys.0,
+            &mut self.sys.1,
+            &mut _time,
+            dt,
+            &mut self.wb.k,
+        );
+        std::mem::swap(&mut self.sys.0, &mut self.sys.1);
+    }
+}
+
+fn gaussian(x: Float, x0: Float, y: Float, y0: Float) -> Float {
+    use sbp::consts::PI;
+
+    let x = x - x0;
+    let y = y - y0;
+
+    let sigma = 0.05;
+
+    1.0 / (2.0 * PI * sigma * sigma) * (-(x * x + y * y) / (2.0 * sigma * sigma)).exp()
+}
+
+#[allow(non_snake_case)]
+/// Solving (Au)_x + (Bu)_y
+/// with:
+///        A               B
+///  [ 0,  0,  0]    [ 0,  1,  0]
+///  [ 0,  0, -1]    [ 1,  0,  0]
+///  [ 0, -1,  0]    [ 0,  0,  0]
+///
+/// This flux is rotated by the grid metrics
+/// (Au)_x + (Bu)_y = 1/J [
+///          (J xi_x Au)_xi + (J eta_x Au)_eta
+///          (J xi_y Bu)_xi + (J eta_y Bu)_eta
+///      ]
+/// where J is the grid determinant
+///
+/// This is used both in fluxes and SAT terms
+fn RHS<SBP: SbpOperator2d>(
+    op: &SBP,
+    k: &mut Field,
+    y: &Field,
+    _grid: &Grid,
+    metrics: &Metrics,
+    tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
+) {
+    fluxes(op, k, y, metrics, tmp);
+
+    let boundaries = BoundaryTerms {
+        north: Boundary::This,
+        south: Boundary::This,
+        west: Boundary::This,
+        east: Boundary::This,
+    };
+    SAT_characteristics(op, k, y, metrics, &boundaries);
+
+    azip!((k in &mut k.0,
+                    &detj in &metrics.detj().broadcast((3, y.ny(), y.nx())).unwrap()) {
+        *k /= detj;
+    });
+}
+
+#[allow(non_snake_case)]
+fn RHS_upwind<UO: UpwindOperator2d>(
+    op: &UO,
+    k: &mut Field,
+    y: &Field,
+    _grid: &Grid,
+    metrics: &Metrics,
+    tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
+) {
+    fluxes(op, k, y, metrics, tmp);
+    dissipation(op, k, y, metrics, tmp);
+
+    let boundaries = BoundaryTerms {
+        north: Boundary::This,
+        south: Boundary::This,
+        west: Boundary::This,
+        east: Boundary::This,
+    };
+    SAT_characteristics(op, k, y, metrics, &boundaries);
+
+    azip!((k in &mut k.0,
+                    &detj in &metrics.detj().broadcast((3, y.ny(), y.nx())).unwrap()) {
+        *k /= detj;
+    });
+}
+
+fn fluxes<SBP: sbp::operators::SbpOperator2d>(
+    op: &SBP,
+    k: &mut Field,
+    y: &Field,
+    metrics: &Metrics,
+    tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
+) {
+    // ex = hz_y
+    {
+        ndarray::azip!((a in &mut tmp.0,
+                        &dxi_dy in &metrics.detj_dxi_dy(),
+                        &hz in &y.hz())
+            *a = dxi_dy * hz
+        );
+        op.diffxi(tmp.0.view(), tmp.1.view_mut());
+
+        ndarray::azip!((b in &mut tmp.2,
+                        &deta_dy in &metrics.detj_deta_dy(),
+                        &hz in &y.hz())
+            *b = deta_dy * hz
+        );
+        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
+        );
+    }
+
+    {
+        // hz = -ey_x + ex_y
+        ndarray::azip!((a in &mut tmp.0,
+                        &dxi_dx in &metrics.detj_dxi_dx(),
+                        &dxi_dy in &metrics.detj_dxi_dy(),
+                        &ex in &y.ex(),
+                        &ey in &y.ey())
+            *a = dxi_dx * -ey + dxi_dy * ex
+        );
+        op.diffxi(tmp.0.view(), tmp.1.view_mut());
+
+        ndarray::azip!((b in &mut tmp.2,
+                        &deta_dx in &metrics.detj_deta_dx(),
+                        &deta_dy in &metrics.detj_deta_dy(),
+                        &ex in &y.ex(),
+                        &ey in &y.ey())
+            *b = deta_dx * -ey + deta_dy * ex
+        );
+        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
+        );
+    }
+
+    // ey = -hz_x
+    {
+        ndarray::azip!((a in &mut tmp.0,
+                        &dxi_dx in &metrics.detj_dxi_dx(),
+                        &hz in &y.hz())
+            *a = dxi_dx * -hz
+        );
+        op.diffxi(tmp.0.view(), tmp.1.view_mut());
+
+        azip!((b in &mut tmp.2,
+                        &deta_dx in &metrics.detj_deta_dx(),
+                        &hz in &y.hz())
+            *b = deta_dx * -hz
+        );
+        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
+        );
+    }
+}
+
+fn dissipation<UO: UpwindOperator2d>(
+    op: &UO,
+    k: &mut Field,
+    y: &Field,
+    metrics: &Metrics,
+    tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
+) {
+    // ex component
+    {
+        ndarray::azip!((a in &mut tmp.0,
+                        &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);
+            *a = ky*ky/r * ex + -kx*ky/r*ey;
+        });
+        op.dissxi(tmp.0.view(), tmp.1.view_mut());
+
+        ndarray::azip!((b in &mut tmp.2,
+                    &kx in &metrics.detj_deta_dx(),
+                    &ky in &metrics.detj_deta_dy(),
+                    &ex in &y.ex(),
+                    &ey in &y.ey()) {
+            let r = Float::hypot(kx, ky);
+            *b = ky*ky/r * ex + -kx*ky/r*ey;
+        });
+        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
+        );
+    }
+
+    // hz component
+    {
+        ndarray::azip!((a in &mut tmp.0,
+                        &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;
+        });
+        op.dissxi(tmp.0.view(), tmp.1.view_mut());
+
+        ndarray::azip!((b in &mut tmp.2,
+                        &kx in &metrics.detj_deta_dx(),
+                        &ky in &metrics.detj_deta_dy(),
+                        &hz in &y.hz()) {
+            let r = Float::hypot(kx, ky);
+            *b = r * hz;
+        });
+        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
+        );
+    }
+
+    // ey
+    {
+        ndarray::azip!((a in &mut tmp.0,
+                        &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);
+            *a = -kx*ky/r * ex + kx*kx/r*ey;
+        });
+        op.dissxi(tmp.0.view(), tmp.1.view_mut());
+
+        ndarray::azip!((b in &mut tmp.2,
+                    &kx in &metrics.detj_deta_dx(),
+                    &ky in &metrics.detj_deta_dy(),
+                    &ex in &y.ex(),
+                    &ey in &y.ey()) {
+            let r = Float::hypot(kx, ky);
+            *b = -kx*ky/r * ex + kx*kx/r*ey;
+        });
+        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
+        );
+    }
+}
+
+#[derive(Clone, Debug)]
+pub enum Boundary {
+    This,
+}
+
+#[derive(Clone, Debug)]
+pub struct BoundaryTerms {
+    pub north: Boundary,
+    pub south: Boundary,
+    pub east: Boundary,
+    pub west: Boundary,
+}
+
+#[allow(non_snake_case)]
+/// Boundary conditions (SAT)
+fn SAT_characteristics<SBP: SbpOperator2d>(
+    op: &SBP,
+    k: &mut Field,
+    y: &Field,
+    metrics: &Metrics,
+    boundaries: &BoundaryTerms,
+) {
+    let ny = y.ny();
+    let nx = y.nx();
+
+    fn positive_flux(kx: Float, ky: Float) -> [[Float; 3]; 3] {
+        let r = (kx * kx + ky * ky).sqrt();
+        [
+            [ky * ky / r / 2.0, ky / 2.0, -kx * ky / r / 2.0],
+            [ky / 2.0, r / 2.0, -kx / 2.0],
+            [-kx * ky / r / 2.0, -kx / 2.0, kx * kx / r / 2.0],
+        ]
+    }
+    fn negative_flux(kx: Float, ky: Float) -> [[Float; 3]; 3] {
+        let r = (kx * kx + ky * ky).sqrt();
+        [
+            [-ky * ky / r / 2.0, ky / 2.0, kx * ky / r / 2.0],
+            [ky / 2.0, -r / 2.0, -kx / 2.0],
+            [kx * ky / r / 2.0, -kx / 2.0, -kx * kx / r / 2.0],
+        ]
+    }
+
+    {
+        let g = match boundaries.east {
+            Boundary::This => y.slice(s![.., .., 0]),
+        };
+        // East boundary
+        let hinv = if op.is_h2xi() {
+            (nx - 2) as Float / op.hxi()[0]
+        } else {
+            (nx - 1) as Float / op.hxi()[0]
+        };
+        for ((((mut k, v), g), &kx), &ky) in k
+            .slice_mut(s![.., .., nx - 1])
+            .gencolumns_mut()
+            .into_iter()
+            .zip(y.slice(s![.., .., nx - 1]).gencolumns())
+            .zip(g.gencolumns())
+            .zip(metrics.detj_dxi_dx().slice(s![.., nx - 1]))
+            .zip(metrics.detj_dxi_dy().slice(s![.., nx - 1]))
+        {
+            // East boundary, positive flux
+            let tau = -1.0;
+
+            let v = (v[0], v[1], v[2]);
+            let g = (g[0], g[1], g[2]);
+
+            let plus = positive_flux(kx, ky);
+
+            k[0] += tau
+                * hinv
+                * (plus[0][0] * (v.0 - g.0) + plus[0][1] * (v.1 - g.1) + plus[0][2] * (v.2 - g.2));
+            k[1] += tau
+                * hinv
+                * (plus[1][0] * (v.0 - g.0) + plus[1][1] * (v.1 - g.1) + plus[1][2] * (v.2 - g.2));
+            k[2] += tau
+                * hinv
+                * (plus[2][0] * (v.0 - g.0) + plus[2][1] * (v.1 - g.1) + plus[2][2] * (v.2 - g.2));
+        }
+    }
+    {
+        // West boundary, negative flux
+        let g = match boundaries.east {
+            Boundary::This => y.slice(s![.., .., nx - 1]),
+        };
+        let hinv = if op.is_h2xi() {
+            (nx - 2) as Float / op.hxi()[0]
+        } else {
+            (nx - 1) as Float / op.hxi()[0]
+        };
+        for ((((mut k, v), g), &kx), &ky) in k
+            .slice_mut(s![.., .., 0])
+            .gencolumns_mut()
+            .into_iter()
+            .zip(y.slice(s![.., .., 0]).gencolumns())
+            .zip(g.gencolumns())
+            .zip(metrics.detj_dxi_dx().slice(s![.., 0]))
+            .zip(metrics.detj_dxi_dy().slice(s![.., 0]))
+        {
+            let tau = 1.0;
+
+            let v = (v[0], v[1], v[2]);
+            let g = (g[0], g[1], g[2]);
+
+            let minus = negative_flux(kx, ky);
+
+            k[0] += tau
+                * hinv
+                * (minus[0][0] * (v.0 - g.0)
+                    + minus[0][1] * (v.1 - g.1)
+                    + minus[0][2] * (v.2 - g.2));
+            k[1] += tau
+                * hinv
+                * (minus[1][0] * (v.0 - g.0)
+                    + minus[1][1] * (v.1 - g.1)
+                    + minus[1][2] * (v.2 - g.2));
+            k[2] += tau
+                * hinv
+                * (minus[2][0] * (v.0 - g.0)
+                    + minus[2][1] * (v.1 - g.1)
+                    + minus[2][2] * (v.2 - g.2));
+        }
+    }
+
+    {
+        let g = match boundaries.north {
+            Boundary::This => y.slice(s![.., 0, ..]),
+        };
+        let hinv = if op.is_h2eta() {
+            (ny - 2) as Float / op.heta()[0]
+        } else {
+            (ny - 1) as Float / op.heta()[0]
+        };
+        for ((((mut k, v), g), &kx), &ky) in k
+            .slice_mut(s![.., ny - 1, ..])
+            .gencolumns_mut()
+            .into_iter()
+            .zip(y.slice(s![.., ny - 1, ..]).gencolumns())
+            .zip(g.gencolumns())
+            .zip(metrics.detj_deta_dx().slice(s![ny - 1, ..]))
+            .zip(metrics.detj_deta_dy().slice(s![ny - 1, ..]))
+        {
+            // North boundary, positive flux
+            let tau = -1.0;
+            let v = (v[0], v[1], v[2]);
+            let g = (g[0], g[1], g[2]);
+
+            let plus = positive_flux(kx, ky);
+
+            k[0] += tau
+                * hinv
+                * (plus[0][0] * (v.0 - g.0) + plus[0][1] * (v.1 - g.1) + plus[0][2] * (v.2 - g.2));
+            k[1] += tau
+                * hinv
+                * (plus[1][0] * (v.0 - g.0) + plus[1][1] * (v.1 - g.1) + plus[1][2] * (v.2 - g.2));
+            k[2] += tau
+                * hinv
+                * (plus[2][0] * (v.0 - g.0) + plus[2][1] * (v.1 - g.1) + plus[2][2] * (v.2 - g.2));
+        }
+    }
+
+    {
+        let g = match boundaries.south {
+            Boundary::This => y.slice(s![.., ny - 1, ..]),
+        };
+        let hinv = if op.is_h2eta() {
+            (ny - 2) as Float / op.heta()[0]
+        } else {
+            (ny - 1) as Float / op.heta()[0]
+        };
+        for ((((mut k, v), g), &kx), &ky) in k
+            .slice_mut(s![.., 0, ..])
+            .gencolumns_mut()
+            .into_iter()
+            .zip(y.slice(s![.., 0, ..]).gencolumns())
+            .zip(g.gencolumns())
+            .zip(metrics.detj_deta_dx().slice(s![0, ..]))
+            .zip(metrics.detj_deta_dy().slice(s![0, ..]))
+        {
+            // South boundary, negative flux
+
+            let tau = 1.0;
+            let v = (v[0], v[1], v[2]);
+            let g = (g[0], g[1], g[2]);
+
+            let minus = negative_flux(kx, ky);
+
+            k[0] += tau
+                * hinv
+                * (minus[0][0] * (v.0 - g.0)
+                    + minus[0][1] * (v.1 - g.1)
+                    + minus[0][2] * (v.2 - g.2));
+            k[1] += tau
+                * hinv
+                * (minus[1][0] * (v.0 - g.0)
+                    + minus[1][1] * (v.1 - g.1)
+                    + minus[1][2] * (v.2 - g.2));
+            k[2] += tau
+                * hinv
+                * (minus[2][0] * (v.0 - g.0)
+                    + minus[2][1] * (v.1 - g.1)
+                    + minus[2][2] * (v.2 - g.2));
+        }
+    }
+}
+
+#[derive(Clone, Debug)]
+pub struct WorkBuffers {
+    k: [Field; 4],
+    tmp: (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
+}
+
+impl WorkBuffers {
+    pub fn new(ny: usize, nx: usize) -> Self {
+        let arr2 = Array2::zeros((ny, nx));
+        let arr3 = Field::new(ny, nx);
+        Self {
+            k: [arr3.clone(), arr3.clone(), arr3.clone(), arr3],
+            tmp: (arr2.clone(), arr2.clone(), arr2.clone(), arr2),
+        }
+    }
+}