rtx-cfd embedded3 items 4–6: field.rs + step/{mod, predictor, projection} (364/700/419 lines) — the host PISO step re-laid from the verified three_d code; gates HELD: MMS + Poiseuille marches value-identical to the 2D embedded solver at nz=1 over 200 steps; 3D MMS orders 0.88 upwind / 1.61 TVD, div ≤ 5e-9; Beltrami 1.08 / 1.25 with face-averaged data, div−mean ≤ 7e-8; Poiseuille |u−û| ≤ 8e-10 at nz 1 and periodic nz 4

Co-Authored-By: Claude Fable 5.1 <[email protected]>

crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/field.rs

@@ -0,0 +1,98 @@
+//! The 3D staggered field: u on `(nx + 1) × ny × nz` faces, v on
+//! `nx × (ny + 1) × nz`, w on `nx × ny × (nz + 1)`, p on the cells; flat
+//! storage (see [`super::Grid`] for the index conventions).
+
+use super::Grid;
+
+#[derive(Debug, Clone)]
+pub struct Field {
+    pub grid: Grid,
+    pub u: Vec<f64>,
+    pub v: Vec<f64>,
+    pub w: Vec<f64>,
+    pub p: Vec<f64>,
+    pub u_old: Vec<f64>,
+    pub v_old: Vec<f64>,
+    pub w_old: Vec<f64>,
+    pub u_star: Vec<f64>,
+    pub v_star: Vec<f64>,
+    pub w_star: Vec<f64>,
+    pub p_prime: Vec<f64>,
+    /// The continuity source of the projection.
+    pub sp: Vec<f64>,
+}
+
+impl Field {
+    #[must_use]
+    pub fn new(grid: Grid) -> Self {
+        let nu = (grid.nx + 1) * grid.ny * grid.nz;
+        let nv = grid.nx * (grid.ny + 1) * grid.nz;
+        let nw = grid.nx * grid.ny * (grid.nz + 1);
+        let nc = grid.cells();
+        Self {
+            grid,
+            u: vec![0.0; nu],
+            v: vec![0.0; nv],
+            w: vec![0.0; nw],
+            p: vec![0.0; nc],
+            u_old: vec![0.0; nu],
+            v_old: vec![0.0; nv],
+            w_old: vec![0.0; nw],
+            u_star: vec![0.0; nu],
+            v_star: vec![0.0; nv],
+            w_star: vec![0.0; nw],
+            p_prime: vec![0.0; nc],
+            sp: vec![0.0; nc],
+        }
+    }
+
+    pub fn update_old_values(&mut self) {
+        self.u_old.copy_from_slice(&self.u);
+        self.v_old.copy_from_slice(&self.v);
+        self.w_old.copy_from_slice(&self.w);
+    }
+
+    pub fn copy_to_starred(&mut self) {
+        self.u_star.copy_from_slice(&self.u);
+        self.v_star.copy_from_slice(&self.v);
+        self.w_star.copy_from_slice(&self.w);
+    }
+
+    /// `max |∇·u|` over the cells (per unit volume).
+    #[must_use]
+    pub fn max_divergence(&self) -> f64 {
+        let g = self.grid;
+        let mut worst = 0.0_f64;
+        for k in 0..g.nz {
+            for j in 0..g.ny {
+                for i in 0..g.nx {
+                    let div = (self.u[g.uface(k, j, i + 1)] - self.u[g.uface(k, j, i)]) / g.dx
+                        + (self.v[g.vface(k, j + 1, i)] - self.v[g.vface(k, j, i)]) / g.dy
+                        + (self.w[g.wface(k + 1, j, i)] - self.w[g.wface(k, j, i)]) / g.dz;
+                    worst = worst.max(div.abs());
+                }
+            }
+        }
+        worst
+    }
+
+    /// Kinetic energy `½ Σ (u² + v² + w²) dV` over the cells (face values
+    /// averaged to the cell).
+    #[must_use]
+    pub fn kinetic_energy(&self, rho: f64) -> f64 {
+        let g = self.grid;
+        let dv = g.dx * g.dy * g.dz;
+        let mut e = 0.0;
+        for k in 0..g.nz {
+            for j in 0..g.ny {
+                for i in 0..g.nx {
+                    let uc = 0.5 * (self.u[g.uface(k, j, i)] + self.u[g.uface(k, j, i + 1)]);
+                    let vc = 0.5 * (self.v[g.vface(k, j, i)] + self.v[g.vface(k, j + 1, i)]);
+                    let wc = 0.5 * (self.w[g.wface(k, j, i)] + self.w[g.wface(k + 1, j, i)]);
+                    e += 0.5 * rho * (uc * uc + vc * vc + wc * wc) * dv;
+                }
+            }
+        }
+        e
+    }
+}

crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs

@@ -6,7 +6,11 @@
 //! Layout: cells `(k, j, i)` row-major, `cell = (k·ny + j)·nx + i`;
 //! u faces on `(nx + 1)·ny·nz`, v on `nx·(ny + 1)·nz`, w on `nx·ny·(nz + 1)`.
 
+pub mod field;
 pub mod grid;
 pub mod poisson;
+pub mod step;
 
+pub use field::Field;
 pub use grid::Grid;
+pub use step::{Boundaries, Fluid, Parameters, Side, Solver, StepResult};

crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs

@@ -0,0 +1,363 @@
+//! The host reference of the PISO step: the 2D embedded solver's predictor
+//! and projection transcribed expression for expression, the z terms
+//! APPENDED after the 2D expression so that at `nz = 1` (z sides slip,
+//! `dz = 1`) every number is the 2D solver's. The fluid predicates are the
+//! wall's hooks (item 9).
+
+mod predictor;
+mod projection;
+
+use super::Grid;
+use super::field::Field;
+use super::poisson::{PcgCache, Problem, solve_pcg_cached};
+use crate::solvers::incompressible::poisson::{
+    MgPrecision, MgSmoother, MultigridParameters, PoissonSolution,
+};
+use crate::solvers::incompressible::simple::ConvectionScheme;
+
+/// Boundary type of one domain side (the 2D `SideBoundary` semantics, plus
+/// `Periodic`, allowed on the z pair only in Stage 1).
+#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
+pub enum Side {
+    #[default]
+    Velocity,
+    SlipWall,
+    PressureOutlet,
+    Periodic,
+}
+
+/// The six sides: `x0` at `x = 0`, `x1` at `x = nx dx`, and so on.
+#[derive(Debug, Clone, Copy, Default)]
+pub struct Boundaries {
+    pub x0: Side,
+    pub x1: Side,
+    pub y0: Side,
+    pub y1: Side,
+    pub z0: Side,
+    pub z1: Side,
+}
+
+impl Boundaries {
+    fn any_outlet(self) -> bool {
+        [self.x0, self.x1, self.y0, self.y1, self.z0, self.z1].contains(&Side::PressureOutlet)
+    }
+
+    pub(super) fn periodic_z(self) -> bool {
+        self.z0 == Side::Periodic
+    }
+}
+
+/// The fluid and the reference scales (the 2D `CfdConfig` fields used).
+#[derive(Debug, Clone, Copy)]
+pub struct Fluid {
+    pub density: f64,
+    pub viscosity: f64,
+    pub reference_velocity: f64,
+    pub reference_length: f64,
+}
+
+#[derive(Debug, Clone)]
+pub struct Parameters {
+    pub corrector_steps: usize,
+    /// Tolerance on the normalised mass imbalance after correction.
+    pub tolerance: f64,
+    pub boundaries: Boundaries,
+    pub poisson_smoother: MgSmoother,
+    pub poisson_precision: MgPrecision,
+    pub convection_scheme: ConvectionScheme,
+    /// Relative part of the pressure solve's inner stop (the 2D 1e-2).
+    pub inner_stop_factor: f64,
+}
+
+impl Default for Parameters {
+    fn default() -> Self {
+        Self {
+            corrector_steps: 2,
+            tolerance: 1e-6,
+            boundaries: Boundaries::default(),
+            poisson_smoother: MgSmoother::Lexicographic,
+            poisson_precision: MgPrecision::F64,
+            convection_scheme: ConvectionScheme::Upwind,
+            inner_stop_factor: 1e-2,
+        }
+    }
+}
+
+/// One step's outcome.
+#[derive(Debug, Clone, Copy)]
+pub struct StepResult {
+    pub converged: bool,
+    pub corrector_steps_performed: usize,
+    pub final_residual: f64,
+    /// CG iterations summed over the step's projections.
+    pub poisson_iterations: usize,
+}
+
+type Vec3Fn = Box<dyn Fn(f64, f64, f64, f64) -> (f64, f64, f64) + Send + Sync>;
+
+pub struct Solver {
+    pub fluid: Fluid,
+    pub params: Parameters,
+    pub(super) momentum_source: Option<Vec3Fn>,
+    boundary_velocity: Option<Vec3Fn>,
+    pcg_cache: PcgCache,
+    time: f64,
+    initialized: bool,
+    /// `(setup ns, iterate ns, solves, CG iterations)` summed.
+    poisson_profile: (u64, u64, u64, u64),
+}
+
+impl Solver {
+    #[must_use]
+    pub fn new(fluid: Fluid, params: Parameters) -> Self {
+        let b = params.boundaries;
+        assert_eq!(
+            b.z0 == Side::Periodic,
+            b.z1 == Side::Periodic,
+            "periodic z needs both z sides periodic"
+        );
+        for side in [b.x0, b.x1, b.y0, b.y1] {
+            assert_ne!(side, Side::Periodic, "Stage 1: periodic only in z");
+        }
+        Self {
+            fluid,
+            params,
+            momentum_source: None,
+            boundary_velocity: None,
+            pcg_cache: PcgCache::default(),
+            time: 0.0,
+            initialized: false,
+            poisson_profile: (0, 0, 0, 0),
+        }
+    }
+
+    pub fn set_momentum_source<F>(&mut self, f: F)
+    where
+        F: Fn(f64, f64, f64, f64) -> (f64, f64, f64) + Send + Sync + 'static,
+    {
+        self.momentum_source = Some(Box::new(f));
+    }
+
+    pub fn set_boundary_velocity<F>(&mut self, f: F)
+    where
+        F: Fn(f64, f64, f64, f64) -> (f64, f64, f64) + Send + Sync + 'static,
+    {
+        self.boundary_velocity = Some(Box::new(f));
+    }
+
+    #[must_use]
+    pub fn time(&self) -> f64 {
+        self.time
+    }
+
+    pub fn set_time(&mut self, t: f64) {
+        self.time = t;
+    }
+
+    #[must_use]
+    pub fn poisson_profile(&self) -> (u64, u64, u64, u64) {
+        self.poisson_profile
+    }
+
+    pub(super) fn boundary(&self, x: f64, y: f64, z: f64, t: f64) -> (f64, f64, f64) {
+        self.boundary_velocity
+            .as_ref()
+            .map_or((0.0, 0.0, 0.0), |f| f(x, y, z, t))
+    }
+
+    // The fluid predicates: everything is fluid until the wall arrives.
+    #[inline]
+    fn u_is_fluid(&self, _k: usize, _j: usize, _i: usize) -> bool {
+        true
+    }
+    #[inline]
+    fn v_is_fluid(&self, _k: usize, _j: usize, _i: usize) -> bool {
+        true
+    }
+    #[inline]
+    fn w_is_fluid(&self, _k: usize, _j: usize, _i: usize) -> bool {
+        true
+    }
+    #[inline]
+    fn cell_is_fluid(&self, _k: usize, _j: usize, _i: usize) -> bool {
+        true
+    }
+
+    pub(super) fn upwind(face_velocity: f64, upstream: f64, downstream: f64) -> f64 {
+        if face_velocity >= 0.0 {
+            upstream
+        } else {
+            downstream
+        }
+    }
+
+    /// Stamp the prescribed normal velocities on the six sides at `t`.
+    pub fn apply_boundary_normals(&self, field: &mut Field, t: f64) {
+        let g = field.grid;
+        let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz);
+        let b = self.params.boundaries;
+        let outlet = Side::PressureOutlet;
+        for k in 0..nz {
+            let z = (k as f64 + 0.5) * dz;
+            for j in 0..ny {
+                let y = (j as f64 + 0.5) * dy;
+                if b.x0 != outlet {
+                    field.u[g.uface(k, j, 0)] = self.boundary(0.0, y, z, t).0;
+                }
+                if b.x1 != outlet {
+                    field.u[g.uface(k, j, nx)] = self.boundary(nx as f64 * dx, y, z, t).0;
+                }
+            }
+            for i in 0..nx {
+                let x = (i as f64 + 0.5) * dx;
+                if b.y0 != outlet {
+                    field.v[g.vface(k, 0, i)] = self.boundary(x, 0.0, z, t).1;
+                }
+                if b.y1 != outlet {
+                    field.v[g.vface(k, ny, i)] = self.boundary(x, ny as f64 * dy, z, t).1;
+                }
+            }
+        }
+        if !b.periodic_z() {
+            for j in 0..ny {
+                let y = (j as f64 + 0.5) * dy;
+                for i in 0..nx {
+                    let x = (i as f64 + 0.5) * dx;
+                    if b.z0 != outlet {
+                        field.w[g.wface(0, j, i)] = self.boundary(x, y, 0.0, t).2;
+                    }
+                    if b.z1 != outlet {
+                        field.w[g.wface(nz, j, i)] = self.boundary(x, y, nz as f64 * dz, t).2;
+                    }
+                }
+            }
+        }
+    }
+
+    /// The explicit predictor on the fluid faces; outlet faces zero-gradient.
+    fn momentum_predictor(&self, field: &mut Field, dt: f64, t_old: f64) {
+        let g = field.grid;
+        let (nx, ny, nz) = (g.nx, g.ny, g.nz);
+        let b = self.params.boundaries;
+        let periodic = b.periodic_z();
+        for k in 0..nz {
+            for j in 0..ny {
+                for i in 1..nx {
+                    if !self.u_is_fluid(k, j, i) {
+                        continue;
+                    }
+                    let rhs = self.u_rhs(field, k, j, i, t_old);
+                    let f = g.uface(k, j, i);
+                    field.u[f] = field.u_old[f] + dt * rhs;
+                }
+            }
+        }
+        for k in 0..nz {
+            for j in 1..ny {
+                for i in 0..nx {
+                    if !self.v_is_fluid(k, j, i) {
+                        continue;
+                    }
+                    let rhs = self.v_rhs(field, k, j, i, t_old);
+                    let f = g.vface(k, j, i);
+                    field.v[f] = field.v_old[f] + dt * rhs;
+                }
+            }
+        }
+        let k_range = if periodic { 0..nz } else { 1..nz };
+        for k in k_range {
+            for j in 0..ny {
+                for i in 0..nx {
+                    if !self.w_is_fluid(k, j, i) {
+                        continue;
+                    }
+                    let rhs = self.w_rhs(field, k, j, i, t_old);
+                    let f = g.wface(k, j, i);
+                    field.w[f] = field.w_old[f] + dt * rhs;
+                }
+            }
+        }
+        if periodic {
+            for j in 0..ny {
+                for i in 0..nx {
+                    field.w[g.wface(nz, j, i)] = field.w[g.wface(0, j, i)];
+                }
+            }
+        }
+        let outlet = Side::PressureOutlet;
+        for k in 0..nz {
+            for j in 0..ny {
+                if b.x0 == outlet {
+                    field.u[g.uface(k, j, 0)] = field.u[g.uface(k, j, 1)];
+                }
+                if b.x1 == outlet {
+                    field.u[g.uface(k, j, nx)] = field.u[g.uface(k, j, nx - 1)];
+                }
+            }
+            for i in 0..nx {
+                if b.y0 == outlet {
+                    field.v[g.vface(k, 0, i)] = field.v[g.vface(k, 1, i)];
+                }
+                if b.y1 == outlet {
+                    field.v[g.vface(k, ny, i)] = field.v[g.vface(k, ny - 1, i)];
+                }
+            }
+        }
+        for j in 0..ny {
+            for i in 0..nx {
+                if b.z0 == outlet {
+                    field.w[g.wface(0, j, i)] = field.w[g.wface(1, j, i)];
+                }
+                if b.z1 == outlet {
+                    field.w[g.wface(nz, j, i)] = field.w[g.wface(nz - 1, j, i)];
+                }
+            }
+        }
+        field.copy_to_starred();
+    }
+
+    /// Stamp the `t = time` boundary data (lazily called by the first step).
+    pub fn initialize(&mut self, field: &mut Field) {
+        let t = self.time;
+        self.apply_boundary_normals(field, t);
+        self.initialized = true;
+    }
+
+    /// One step of `dt`: predictor, correctors, clock.
+    pub fn advance(&mut self, field: &mut Field, dt: f64) -> StepResult {
+        assert!(
+            dt > 0.0 && dt.is_finite(),
+            "time step must be positive and finite, got {dt}"
+        );
+        if !self.initialized {
+            self.initialize(field);
+        }
+        let t_old = self.time;
+        let t_new = t_old + dt;
+        field.update_old_values();
+        self.momentum_predictor(field, dt, t_old);
+        self.apply_boundary_normals(field, t_new);
+        field.copy_to_starred();
+        let mut total = 0;
+        let mut final_residual = f64::INFINITY;
+        let mut poisson_iterations = 0;
+        for corrector in 0..self.params.corrector_steps.max(1) {
+            let sol = self.solve_correction(field, dt, corrector == 0);
+            poisson_iterations += sol.iterations;
+            let mass_residual = self.apply_correction(field, dt);
+            final_residual = mass_residual;
+            total += 1;
+            if mass_residual < self.params.tolerance {
+                break;
+            }
+            field.copy_to_starred();
+        }
+        self.time = t_new;
+        StepResult {
+            converged: final_residual < self.params.tolerance,
+            corrector_steps_performed: total,
+            final_residual,
+            poisson_iterations,
+        }
+    }
+}

crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/predictor.rs

@@ -0,0 +1,679 @@
+//! The three explicit predictors of the 3D PISO step (`Solver`): the
+//! 2D `u_rhs`/`v_rhs` expression for expression with the z terms appended,
+//! and `w_rhs` as the v pattern turned along z. Split from `piso_host.rs`
+//! for the file-size rule; `impl Solver` continues here.
+
+use super::{Side, Solver};
+use crate::solvers::incompressible::embedded3::field::Field;
+use crate::solvers::incompressible::simple::ConvectionScheme;
+
+impl Solver {
+    /// Plane above `k` (wrapping when periodic).
+    #[inline]
+    fn k_up(&self, k: usize, nz: usize) -> Option<usize> {
+        if k + 1 < nz {
+            Some(k + 1)
+        } else if self.params.boundaries.periodic_z() {
+            Some(0)
+        } else {
+            None
+        }
+    }
+
+    #[inline]
+    fn k_down(&self, k: usize, nz: usize) -> Option<usize> {
+        if k > 0 {
+            Some(k - 1)
+        } else if self.params.boundaries.periodic_z() {
+            Some(nz - 1)
+        } else {
+            None
+        }
+    }
+
+    /// The predictor's right-hand side on the u face `(k, j, i)`, `i = 1..nx`:
+    /// the 2D `u_rhs` expression for expression, then `− conv_z + diff_z`.
+    #[allow(clippy::too_many_lines)]
+    pub(crate) fn u_rhs(&self, field: &Field, k: usize, j: usize, i: usize, t_old: f64) -> f64 {
+        let g = field.grid;
+        let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz);
+        let rho = self.fluid.density;
+        let nu = self.fluid.viscosity / rho;
+        let b = self.params.boundaries;
+        let velocity = Side::Velocity;
+        let uo = &field.u_old;
+        let vo = &field.v_old;
+        let wo = &field.w_old;
+        let uf = |kk: usize, jj: usize, ii: usize| g.uface(kk, jj, ii);
+        let vf = |kk: usize, jj: usize, ii: usize| g.vface(kk, jj, ii);
+        let wf = |kk: usize, jj: usize, ii: usize| g.wface(kk, jj, ii);
+        let zc = (k as f64 + 0.5) * dz;
+        let u_p = uo[uf(k, j, i)];
+
+        let ue_face = 0.5 * (uo[uf(k, j, i)] + uo[uf(k, j, i + 1)]);
+        let uw_face = 0.5 * (uo[uf(k, j, i - 1)] + uo[uf(k, j, i)]);
+
+        let south_is_wall = j == 0;
+        let north_is_wall = j + 1 == ny;
+
+        let vn_face = 0.5 * (vo[vf(k, j + 1, i - 1)] + vo[vf(k, j + 1, i)]);
+        let vs_face = 0.5 * (vo[vf(k, j, i - 1)] + vo[vf(k, j, i)]);
+
+        let beyond_north = if b.y1 == velocity {
+            self.boundary(i as f64 * dx, ny as f64 * dy, zc, t_old).0
+        } else {
+            u_p
+        };
+        let beyond_south = if b.y0 == velocity {
+            self.boundary(i as f64 * dx, 0.0, zc, t_old).0
+        } else {
+            u_p
+        };
+
+        let conv_x = (ue_face * Self::upwind(ue_face, uo[uf(k, j, i)], uo[uf(k, j, i + 1)])
+            - uw_face * Self::upwind(uw_face, uo[uf(k, j, i - 1)], uo[uf(k, j, i)]))
+            / dx;
+        let conv_y = (vn_face
+            * if north_is_wall {
+                Self::upwind(vn_face, u_p, beyond_north)
+            } else {
+                Self::upwind(vn_face, uo[uf(k, j, i)], uo[uf(k, j + 1, i)])
+            }
+            - vs_face
+                * if south_is_wall {
+                    Self::upwind(vs_face, beyond_south, u_p)
+                } else {
+                    Self::upwind(vs_face, uo[uf(k, j - 1, i)], uo[uf(k, j, i)])
+                })
+            / dy;
+
+        let scheme = self.params.convection_scheme;
+        let mut conv_x = conv_x;
+        let mut conv_y = conv_y;
+        if scheme != ConvectionScheme::Upwind {
+            let delta_e = if ue_face >= 0.0 {
+                scheme.face_correction(
+                    Some(uo[uf(k, j, i - 1)]),
+                    uo[uf(k, j, i)],
+                    uo[uf(k, j, i + 1)],
+                )
+            } else {
+                let far = (i + 2 <= nx).then(|| uo[uf(k, j, i + 2)]);
+                scheme.face_correction(far, uo[uf(k, j, i + 1)], uo[uf(k, j, i)])
+            };
+            let delta_w = if uw_face >= 0.0 {
+                let far = (i >= 2).then(|| uo[uf(k, j, i - 2)]);
+                scheme.face_correction(far, uo[uf(k, j, i - 1)], uo[uf(k, j, i)])
+            } else {
+                scheme.face_correction(
+                    Some(uo[uf(k, j, i + 1)]),
+                    uo[uf(k, j, i)],
+                    uo[uf(k, j, i - 1)],
+                )
+            };
+            let delta_n = if north_is_wall {
+                0.0
+            } else if vn_face >= 0.0 {
+                let far = (j >= 1).then(|| uo[uf(k, j - 1, i)]);
+                scheme.face_correction(far, uo[uf(k, j, i)], uo[uf(k, j + 1, i)])
+            } else {
+                let far = (j + 2 < ny).then(|| uo[uf(k, j + 2, i)]);
+                scheme.face_correction(far, uo[uf(k, j + 1, i)], uo[uf(k, j, i)])
+            };
+            let delta_s = if south_is_wall {
+                0.0
+            } else if vs_face >= 0.0 {
+                let far = (j >= 2).then(|| uo[uf(k, j - 2, i)]);
+                scheme.face_correction(far, uo[uf(k, j - 1, i)], uo[uf(k, j, i)])
+            } else {
+                let far = (j + 1 < ny).then(|| uo[uf(k, j + 1, i)]);
+                scheme.face_correction(far, uo[uf(k, j, i)], uo[uf(k, j - 1, i)])
+            };
+            conv_x += (ue_face * delta_e - uw_face * delta_w) / dx;
+            conv_y += (vn_face * delta_n - vs_face * delta_s) / dy;
+        }
+
+        let diff_x = nu * (uo[uf(k, j, i + 1)] - 2.0 * u_p + uo[uf(k, j, i - 1)]) / (dx * dx);
+
+        let flux_north = if north_is_wall {
+            if b.y1 == velocity {
+                let u_wall = self.boundary(i as f64 * dx, ny as f64 * dy, zc, t_old).0;
+                nu * (u_wall - u_p) / (0.5 * dy)
+            } else {
+                0.0
+            }
+        } else {
+            nu * (uo[uf(k, j + 1, i)] - u_p) / dy
+        };
+        let flux_south = if south_is_wall {
+            if b.y0 == velocity {
+                let u_wall = self.boundary(i as f64 * dx, 0.0, zc, t_old).0;
+                nu * (u_p - u_wall) / (0.5 * dy)
+            } else {
+                0.0
+            }
+        } else {
+            nu * (u_p - uo[uf(k, j - 1, i)]) / dy
+        };
+        let diff_y = (flux_north - flux_south) / dy;
+
+        let pressure_gradient =
+            -(field.p[g.cell(k, j, i)] - field.p[g.cell(k, j, i - 1)]) / (rho * dx);
+
+        let body_force = self.momentum_source.as_ref().map_or(0.0, |f| {
+            f(i as f64 * dx, (j as f64 + 0.5) * dy, zc, t_old).0 / rho
+        });
+
+        let rhs_2d = -conv_x - conv_y + diff_x + diff_y + pressure_gradient + body_force;
+
+        // --- the z terms, the y pattern turned along k ---
+        let ku = self.k_up(k, nz);
+        let kd = self.k_down(k, nz);
+        let top_is_wall = ku.is_none();
+        let bottom_is_wall = kd.is_none();
+        // The w faces above/below the u face: on top of the cells west and
+        // east of it (face k + 1 of cell k is face index k + 1; periodic:
+        // the face at k = nz equals the face at 0).
+        let wt_face = 0.5 * (wo[wf(k + 1, j, i - 1)] + wo[wf(k + 1, j, i)]);
+        let wb_face = 0.5 * (wo[wf(k, j, i - 1)] + wo[wf(k, j, i)]);
+        let beyond_top = if b.z1 == velocity {
+            self.boundary(i as f64 * dx, (j as f64 + 0.5) * dy, nz as f64 * dz, t_old)
+                .0
+        } else {
+            u_p
+        };
+        let beyond_bottom = if b.z0 == velocity {
+            self.boundary(i as f64 * dx, (j as f64 + 0.5) * dy, 0.0, t_old)
+                .0
+        } else {
+            u_p
+        };
+        let u_up = ku.map(|kk| uo[uf(kk, j, i)]);
+        let u_dn = kd.map(|kk| uo[uf(kk, j, i)]);
+        let mut conv_z = (wt_face
+            * match u_up {
+                Some(un) => Self::upwind(wt_face, u_p, un),
+                None => Self::upwind(wt_face, u_p, beyond_top),
+            }
+            - wb_face
+                * match u_dn {
+                    Some(ud) => Self::upwind(wb_face, ud, u_p),
+                    None => Self::upwind(wb_face, beyond_bottom, u_p),
+                })
+            / dz;
+        if scheme != ConvectionScheme::Upwind {
+            let far_up2 = ku
+                .and_then(|kk| self.k_up(kk, nz))
+                .map(|kk| uo[uf(kk, j, i)]);
+            let far_dn2 = kd
+                .and_then(|kk| self.k_down(kk, nz))
+                .map(|kk| uo[uf(kk, j, i)]);
+            let delta_t = if top_is_wall {
+                0.0
+            } else if wt_face >= 0.0 {
+                scheme.face_correction(u_dn, u_p, u_up.unwrap_or(u_p))
+            } else {
+                scheme.face_correction(far_up2, u_up.unwrap_or(u_p), u_p)
+            };
+            let delta_b = if bottom_is_wall {
+                0.0
+            } else if wb_face >= 0.0 {
+                scheme.face_correction(far_dn2, u_dn.unwrap_or(u_p), u_p)
+            } else {
+                scheme.face_correction(u_up, u_p, u_dn.unwrap_or(u_p))
+            };
+            conv_z += (wt_face * delta_t - wb_face * delta_b) / dz;
+        }
+        let flux_top = match u_up {
+            Some(un) => nu * (un - u_p) / dz,
+            None => {
+                if b.z1 == velocity {
+                    nu * (beyond_top - u_p) / (0.5 * dz)
+                } else {
+                    0.0
+                }
+            }
+        };
+        let flux_bottom = match u_dn {
+            Some(ud) => nu * (u_p - ud) / dz,
+            None => {
+                if b.z0 == velocity {
+                    nu * (u_p - beyond_bottom) / (0.5 * dz)
+                } else {
+                    0.0
+                }
+            }
+        };
+        let diff_z = (flux_top - flux_bottom) / dz;
+
+        rhs_2d - conv_z + diff_z
+    }
+
+    /// The v face `(k, j, i)`, `j = 1..ny`: the 2D `v_rhs` then the z terms.
+    #[allow(clippy::too_many_lines)]
+    pub(crate) fn v_rhs(&self, field: &Field, k: usize, j: usize, i: usize, t_old: f64) -> f64 {
+        let g = field.grid;
+        let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz);
+        let rho = self.fluid.density;
+        let nu = self.fluid.viscosity / rho;
+        let b = self.params.boundaries;
+        let velocity = Side::Velocity;
+        let uo = &field.u_old;
+        let vo = &field.v_old;
+        let wo = &field.w_old;
+        let uf = |kk: usize, jj: usize, ii: usize| g.uface(kk, jj, ii);
+        let vf = |kk: usize, jj: usize, ii: usize| g.vface(kk, jj, ii);
+        let wf = |kk: usize, jj: usize, ii: usize| g.wface(kk, jj, ii);
+        let zc = (k as f64 + 0.5) * dz;
+        let v_p = vo[vf(k, j, i)];
+
+        let vn_face = 0.5 * (vo[vf(k, j, i)] + vo[vf(k, j + 1, i)]);
+        let vs_face = 0.5 * (vo[vf(k, j - 1, i)] + vo[vf(k, j, i)]);
+
+        let west_is_wall = i == 0;
+        let east_is_wall = i + 1 == nx;
+
+        let ue_face = 0.5 * (uo[uf(k, j - 1, i + 1)] + uo[uf(k, j, i + 1)]);
+        let uw_face = 0.5 * (uo[uf(k, j - 1, i)] + uo[uf(k, j, i)]);
+        let beyond_east = if b.x1 == velocity {
+            self.boundary(nx as f64 * dx, j as f64 * dy, zc, t_old).1
+        } else {
+            v_p
+        };
+        let beyond_west = if b.x0 == velocity {
+            self.boundary(0.0, j as f64 * dy, zc, t_old).1
+        } else {
+            v_p
+        };
+
+        let conv_y = (vn_face * Self::upwind(vn_face, vo[vf(k, j, i)], vo[vf(k, j + 1, i)])
+            - vs_face * Self::upwind(vs_face, vo[vf(k, j - 1, i)], vo[vf(k, j, i)]))
+            / dy;
+        let conv_x = (ue_face
+            * if east_is_wall {
+                Self::upwind(ue_face, v_p, beyond_east)
+            } else {
+                Self::upwind(ue_face, vo[vf(k, j, i)], vo[vf(k, j, i + 1)])
+            }
+            - uw_face
+                * if west_is_wall {
+                    Self::upwind(uw_face, beyond_west, v_p)
+                } else {
+                    Self::upwind(uw_face, vo[vf(k, j, i - 1)], vo[vf(k, j, i)])
+                })
+            / dx;
+
+        let scheme = self.params.convection_scheme;
+        let mut conv_x = conv_x;
+        let mut conv_y = conv_y;
+        if scheme != ConvectionScheme::Upwind {
+            let delta_n = if vn_face >= 0.0 {
+                scheme.face_correction(
+                    Some(vo[vf(k, j - 1, i)]),
+                    vo[vf(k, j, i)],
+                    vo[vf(k, j + 1, i)],
+                )
+            } else {
+                let far = (j + 2 <= ny).then(|| vo[vf(k, j + 2, i)]);
+                scheme.face_correction(far, vo[vf(k, j + 1, i)], vo[vf(k, j, i)])
+            };
+            let delta_s = if vs_face >= 0.0 {
+                let far = (j >= 2).then(|| vo[vf(k, j - 2, i)]);
+                scheme.face_correction(far, vo[vf(k, j - 1, i)], vo[vf(k, j, i)])
+            } else {
+                scheme.face_correction(
+                    Some(vo[vf(k, j + 1, i)]),
+                    vo[vf(k, j, i)],
+                    vo[vf(k, j - 1, i)],
+                )
+            };
+            let delta_e = if east_is_wall {
+                0.0
+            } else if ue_face >= 0.0 {
+                let far = (i >= 1).then(|| vo[vf(k, j, i - 1)]);
+                scheme.face_correction(far, vo[vf(k, j, i)], vo[vf(k, j, i + 1)])
+            } else {
+                let far = (i + 2 < nx).then(|| vo[vf(k, j, i + 2)]);
+                scheme.face_correction(far, vo[vf(k, j, i + 1)], vo[vf(k, j, i)])
+            };
+            let delta_w = if west_is_wall {
+                0.0
+            } else if uw_face >= 0.0 {
+                let far = (i >= 2).then(|| vo[vf(k, j, i - 2)]);
+                scheme.face_correction(far, vo[vf(k, j, i - 1)], vo[vf(k, j, i)])
+            } else {
+                let far = (i + 1 < nx).then(|| vo[vf(k, j, i + 1)]);
+                scheme.face_correction(far, vo[vf(k, j, i)], vo[vf(k, j, i - 1)])
+            };
+            conv_y += (vn_face * delta_n - vs_face * delta_s) / dy;
+            conv_x += (ue_face * delta_e - uw_face * delta_w) / dx;
+        }
+
+        let diff_y = nu * (vo[vf(k, j + 1, i)] - 2.0 * v_p + vo[vf(k, j - 1, i)]) / (dy * dy);
+
+        let flux_east = if east_is_wall {
+            if b.x1 == velocity {
+                let v_wall = self.boundary(nx as f64 * dx, j as f64 * dy, zc, t_old).1;
+                nu * (v_wall - v_p) / (0.5 * dx)
+            } else {
+                0.0
+            }
+        } else {
+            nu * (vo[vf(k, j, i + 1)] - v_p) / dx
+        };
+        let flux_west = if west_is_wall {
+            if b.x0 == velocity {
+                let v_wall = self.boundary(0.0, j as f64 * dy, zc, t_old).1;
+                nu * (v_p - v_wall) / (0.5 * dx)
+            } else {
+                0.0
+            }
+        } else {
+            nu * (v_p - vo[vf(k, j, i - 1)]) / dx
+        };
+        let diff_x = (flux_east - flux_west) / dx;
+
+        let pressure_gradient =
+            -(field.p[g.cell(k, j, i)] - field.p[g.cell(k, j - 1, i)]) / (rho * dy);
+
+        let body_force = self.momentum_source.as_ref().map_or(0.0, |f| {
+            f((i as f64 + 0.5) * dx, j as f64 * dy, zc, t_old).1 / rho
+        });
+
+        let rhs_2d = -conv_x - conv_y + diff_x + diff_y + pressure_gradient + body_force;
+
+        // --- z terms ---
+        let ku = self.k_up(k, nz);
+        let kd = self.k_down(k, nz);
+        let top_is_wall = ku.is_none();
+        let bottom_is_wall = kd.is_none();
+        let wt_face = 0.5 * (wo[wf(k + 1, j - 1, i)] + wo[wf(k + 1, j, i)]);
+        let wb_face = 0.5 * (wo[wf(k, j - 1, i)] + wo[wf(k, j, i)]);
+        let beyond_top = if b.z1 == velocity {
+            self.boundary((i as f64 + 0.5) * dx, j as f64 * dy, nz as f64 * dz, t_old)
+                .1
+        } else {
+            v_p
+        };
+        let beyond_bottom = if b.z0 == velocity {
+            self.boundary((i as f64 + 0.5) * dx, j as f64 * dy, 0.0, t_old)
+                .1
+        } else {
+            v_p
+        };
+        let v_up = ku.map(|kk| vo[vf(kk, j, i)]);
+        let v_dn = kd.map(|kk| vo[vf(kk, j, i)]);
+        let mut conv_z = (wt_face
+            * match v_up {
+                Some(vn) => Self::upwind(wt_face, v_p, vn),
+                None => Self::upwind(wt_face, v_p, beyond_top),
+            }
+            - wb_face
+                * match v_dn {
+                    Some(vd) => Self::upwind(wb_face, vd, v_p),
+                    None => Self::upwind(wb_face, beyond_bottom, v_p),
+                })
+            / dz;
+        if scheme != ConvectionScheme::Upwind {
+            let far_up2 = ku
+                .and_then(|kk| self.k_up(kk, nz))
+                .map(|kk| vo[vf(kk, j, i)]);
+            let far_dn2 = kd
+                .and_then(|kk| self.k_down(kk, nz))
+                .map(|kk| vo[vf(kk, j, i)]);
+            let delta_t = if top_is_wall {
+                0.0
+            } else if wt_face >= 0.0 {
+                scheme.face_correction(v_dn, v_p, v_up.unwrap_or(v_p))
+            } else {
+                scheme.face_correction(far_up2, v_up.unwrap_or(v_p), v_p)
+            };
+            let delta_b = if bottom_is_wall {
+                0.0
+            } else if wb_face >= 0.0 {
+                scheme.face_correction(far_dn2, v_dn.unwrap_or(v_p), v_p)
+            } else {
+                scheme.face_correction(v_up, v_p, v_dn.unwrap_or(v_p))
+            };
+            conv_z += (wt_face * delta_t - wb_face * delta_b) / dz;
+        }
+        let flux_top = match v_up {
+            Some(vn) => nu * (vn - v_p) / dz,
+            None => {
+                if b.z1 == velocity {
+                    nu * (beyond_top - v_p) / (0.5 * dz)
+                } else {
+                    0.0
+                }
+            }
+        };
+        let flux_bottom = match v_dn {
+            Some(vd) => nu * (v_p - vd) / dz,
+            None => {
+                if b.z0 == velocity {
+                    nu * (v_p - beyond_bottom) / (0.5 * dz)
+                } else {
+                    0.0
+                }
+            }
+        };
+        let diff_z = (flux_top - flux_bottom) / dz;
+
+        rhs_2d - conv_z + diff_z
+    }
+
+    /// The w face `(k, j, i)` between cells `k − 1` (wrapping when periodic)
+    /// and `k`: the v pattern with z as its own direction and x, y transverse.
+    #[allow(clippy::too_many_lines)]
+    pub(crate) fn w_rhs(&self, field: &Field, k: usize, j: usize, i: usize, t_old: f64) -> f64 {
+        let g = field.grid;
+        let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz);
+        let rho = self.fluid.density;
+        let nu = self.fluid.viscosity / rho;
+        let b = self.params.boundaries;
+        let velocity = Side::Velocity;
+        let periodic = b.periodic_z();
+        let uo = &field.u_old;
+        let vo = &field.v_old;
+        let wo = &field.w_old;
+        let uf = |kk: usize, jj: usize, ii: usize| g.uface(kk, jj, ii);
+        let vf = |kk: usize, jj: usize, ii: usize| g.vface(kk, jj, ii);
+        let wf = |kk: usize, jj: usize, ii: usize| g.wface(kk, jj, ii);
+        // The cells below and above this face (k = 0 only when periodic).
+        let k_below = if k > 0 { k - 1 } else { nz - 1 };
+        let k_above = k % nz;
+        // Own-direction neighbours: the faces k − 1 and k + 1 (wrapping).
+        let w_dn_idx = if k > 0 {
+            wf(k - 1, j, i)
+        } else {
+            wf(nz - 1, j, i)
+        };
+        let w_up_idx = if k + 1 == nz && periodic {
+            wf(0, j, i)
+        } else {
+            wf(k + 1, j, i)
+        };
+        let zf = k as f64 * dz;
+        let w_p = wo[wf(k, j, i)];
+
+        let wt_face = 0.5 * (wo[wf(k, j, i)] + wo[w_up_idx]);
+        let wb_face = 0.5 * (wo[w_dn_idx] + wo[wf(k, j, i)]);
+
+        let west_is_wall = i == 0;
+        let east_is_wall = i + 1 == nx;
+        let south_is_wall = j == 0;
+        let north_is_wall = j + 1 == ny;
+
+        let ue_face = 0.5 * (uo[uf(k_below, j, i + 1)] + uo[uf(k_above, j, i + 1)]);
+        let uw_face = 0.5 * (uo[uf(k_below, j, i)] + uo[uf(k_above, j, i)]);
+        let vn_face = 0.5 * (vo[vf(k_below, j + 1, i)] + vo[vf(k_above, j + 1, i)]);
+        let vs_face = 0.5 * (vo[vf(k_below, j, i)] + vo[vf(k_above, j, i)]);
+        let beyond_east = if b.x1 == velocity {
+            self.boundary(nx as f64 * dx, (j as f64 + 0.5) * dy, zf, t_old)
+                .2
+        } else {
+            w_p
+        };
+        let beyond_west = if b.x0 == velocity {
+            self.boundary(0.0, (j as f64 + 0.5) * dy, zf, t_old).2
+        } else {
+            w_p
+        };
+        let beyond_north = if b.y1 == velocity {
+            self.boundary((i as f64 + 0.5) * dx, ny as f64 * dy, zf, t_old)
+                .2
+        } else {
+            w_p
+        };
+        let beyond_south = if b.y0 == velocity {
+            self.boundary((i as f64 + 0.5) * dx, 0.0, zf, t_old).2
+        } else {
+            w_p
+        };
+
+        let conv_z = (wt_face * Self::upwind(wt_face, w_p, wo[w_up_idx])
+            - wb_face * Self::upwind(wb_face, wo[w_dn_idx], w_p))
+            / dz;
+        let conv_x = (ue_face
+            * if east_is_wall {
+                Self::upwind(ue_face, w_p, beyond_east)
+            } else {
+                Self::upwind(ue_face, w_p, wo[wf(k, j, i + 1)])
+            }
+            - uw_face
+                * if west_is_wall {
+                    Self::upwind(uw_face, beyond_west, w_p)
+                } else {
+                    Self::upwind(uw_face, wo[wf(k, j, i - 1)], w_p)
+                })
+            / dx;
+        let conv_y = (vn_face
+            * if north_is_wall {
+                Self::upwind(vn_face, w_p, beyond_north)
+            } else {
+                Self::upwind(vn_face, w_p, wo[wf(k, j + 1, i)])
+            }
+            - vs_face
+                * if south_is_wall {
+                    Self::upwind(vs_face, beyond_south, w_p)
+                } else {
+                    Self::upwind(vs_face, wo[wf(k, j - 1, i)], w_p)
+                })
+            / dy;
+
+        let scheme = self.params.convection_scheme;
+        let mut conv_x = conv_x;
+        let mut conv_y = conv_y;
+        let mut conv_z = conv_z;
+        if scheme != ConvectionScheme::Upwind {
+            // Own direction: far nodes two faces away (wrapping when periodic).
+            let far_up2 = if periodic {
+                Some(wo[wf((k + 2) % nz, j, i)])
+            } else {
+                (k + 2 <= nz).then(|| wo[wf(k + 2, j, i)])
+            };
+            let far_dn2 = if periodic {
+                Some(wo[wf((k + nz - 2) % nz, j, i)])
+            } else {
+                (k >= 2).then(|| wo[wf(k - 2, j, i)])
+            };
+            let delta_t = if wt_face >= 0.0 {
+                scheme.face_correction(Some(wo[w_dn_idx]), w_p, wo[w_up_idx])
+            } else {
+                scheme.face_correction(far_up2, wo[w_up_idx], w_p)
+            };
+            let delta_b = if wb_face >= 0.0 {
+                scheme.face_correction(far_dn2, wo[w_dn_idx], w_p)
+            } else {
+                scheme.face_correction(Some(wo[w_up_idx]), w_p, wo[w_dn_idx])
+            };
+            let delta_e = if east_is_wall {
+                0.0
+            } else if ue_face >= 0.0 {
+                let far = (i >= 1).then(|| wo[wf(k, j, i - 1)]);
+                scheme.face_correction(far, w_p, wo[wf(k, j, i + 1)])
+            } else {
+                let far = (i + 2 < nx).then(|| wo[wf(k, j, i + 2)]);
+                scheme.face_correction(far, wo[wf(k, j, i + 1)], w_p)
+            };
+            let delta_w = if west_is_wall {
+                0.0
+            } else if uw_face >= 0.0 {
+                let far = (i >= 2).then(|| wo[wf(k, j, i - 2)]);
+                scheme.face_correction(far, wo[wf(k, j, i - 1)], w_p)
+            } else {
+                let far = (i + 1 < nx).then(|| wo[wf(k, j, i + 1)]);
+                scheme.face_correction(far, w_p, wo[wf(k, j, i - 1)])
+            };
+            let delta_n = if north_is_wall {
+                0.0
+            } else if vn_face >= 0.0 {
+                let far = (j >= 1).then(|| wo[wf(k, j - 1, i)]);
+                scheme.face_correction(far, w_p, wo[wf(k, j + 1, i)])
+            } else {
+                let far = (j + 2 < ny).then(|| wo[wf(k, j + 2, i)]);
+                scheme.face_correction(far, wo[wf(k, j + 1, i)], w_p)
+            };
+            let delta_s = if south_is_wall {
+                0.0
+            } else if vs_face >= 0.0 {
+                let far = (j >= 2).then(|| wo[wf(k, j - 2, i)]);
+                scheme.face_correction(far, wo[wf(k, j - 1, i)], w_p)
+            } else {
+                let far = (j + 1 < ny).then(|| wo[wf(k, j + 1, i)]);
+                scheme.face_correction(far, w_p, wo[wf(k, j - 1, i)])
+            };
+            conv_z += (wt_face * delta_t - wb_face * delta_b) / dz;
+            conv_x += (ue_face * delta_e - uw_face * delta_w) / dx;
+            conv_y += (vn_face * delta_n - vs_face * delta_s) / dy;
+        }
+
+        let diff_z = nu * (wo[w_up_idx] - 2.0 * w_p + wo[w_dn_idx]) / (dz * dz);
+        let flux_east = if east_is_wall {
+            if b.x1 == velocity {
+                nu * (beyond_east - w_p) / (0.5 * dx)
+            } else {
+                0.0
+            }
+        } else {
+            nu * (wo[wf(k, j, i + 1)] - w_p) / dx
+        };
+        let flux_west = if west_is_wall {
+            if b.x0 == velocity {
+                nu * (w_p - beyond_west) / (0.5 * dx)
+            } else {
+                0.0
+            }
+        } else {
+            nu * (w_p - wo[wf(k, j, i - 1)]) / dx
+        };
+        let diff_x = (flux_east - flux_west) / dx;
+        let flux_north = if north_is_wall {
+            if b.y1 == velocity {
+                nu * (beyond_north - w_p) / (0.5 * dy)
+            } else {
+                0.0
+            }
+        } else {
+            nu * (wo[wf(k, j + 1, i)] - w_p) / dy
+        };
+        let flux_south = if south_is_wall {
+            if b.y0 == velocity {
+                nu * (w_p - beyond_south) / (0.5 * dy)
+            } else {
+                0.0
+            }
+        } else {
+            nu * (w_p - wo[wf(k, j - 1, i)]) / dy
+        };
+        let diff_y = (flux_north - flux_south) / dy;
+
+        let pressure_gradient =
+            -(field.p[g.cell(k_above, j, i)] - field.p[g.cell(k_below, j, i)]) / (rho * dz);
+        let body_force = self.momentum_source.as_ref().map_or(0.0, |f| {
+            f((i as f64 + 0.5) * dx, (j as f64 + 0.5) * dy, zf, t_old).2 / rho
+        });
+
+        -conv_x - conv_y - conv_z + diff_x + diff_y + diff_z + pressure_gradient + body_force
+    }
+}

crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/projection.rs

@@ -0,0 +1,418 @@
+//! The projection: the pressure-correction operator, the boundary and
+//! source tables, the solve and apply halves of one corrector.
+
+use super::{Boundaries, Side, Solver};
+use crate::solvers::incompressible::embedded3::Grid;
+use crate::solvers::incompressible::embedded3::field::Field;
+use crate::solvers::incompressible::embedded3::poisson::{Problem, solve_pcg_cached};
+use crate::solvers::incompressible::poisson::{MultigridParameters, PoissonSolution};
+
+impl Solver {
+    /// The pressure-correction OPERATOR (the 2D `poisson_problem` with the z
+    /// faces): coefficient `dt · A / δ` across every fluid interior face,
+    /// zero across prescribed ones, the outlet's Dirichlet half a cell out
+    /// as `extra_diag`; the right-hand side left zero.
+    pub(crate) fn poisson_operator(&self, g: Grid, dt: f64) -> Problem {
+        let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz);
+        let b = self.params.boundaries;
+        let outlet = Side::PressureOutlet;
+        let periodic = b.periodic_z();
+        let ae_interior = dt * (dy * dz) / dx;
+        let an_interior = dt * (dx * dz) / dy;
+        let at_interior = dt * (dx * dy) / dz;
+        let ae_outlet = dt * (dy * dz) / (0.5 * dx);
+        let an_outlet = dt * (dx * dz) / (0.5 * dy);
+        let at_outlet = dt * (dx * dy) / (0.5 * dz);
+        let mut problem = Problem::new(nx, ny, nz);
+        problem.periodic_z = periodic;
+        for k in 0..nz {
+            for j in 0..ny {
+                for i in 0..nx {
+                    let idx = g.cell(k, j, i);
+                    if !self.cell_is_fluid(k, j, i) {
+                        problem.active[idx] = false;
+                        continue;
+                    }
+                    let mut extra = 0.0;
+                    if i + 1 == nx {
+                        if b.x1 == outlet {
+                            extra += ae_outlet;
+                        }
+                    } else if self.u_is_fluid(k, j, i + 1) {
+                        problem.ae[idx] = ae_interior;
+                    }
+                    if i == 0 {
+                        if b.x0 == outlet {
+                            extra += ae_outlet;
+                        }
+                    } else if self.u_is_fluid(k, j, i) {
+                        problem.aw[idx] = ae_interior;
+                    }
+                    if j + 1 == ny {
+                        if b.y1 == outlet {
+                            extra += an_outlet;
+                        }
+                    } else if self.v_is_fluid(k, j + 1, i) {
+                        problem.an[idx] = an_interior;
+                    }
+                    if j == 0 {
+                        if b.y0 == outlet {
+                            extra += an_outlet;
+                        }
+                    } else if self.v_is_fluid(k, j, i) {
+                        problem.as_[idx] = an_interior;
+                    }
+                    if k + 1 == nz && !periodic {
+                        if b.z1 == outlet {
+                            extra += at_outlet;
+                        }
+                    } else if self.w_is_fluid((k + 1) % nz, j, i) {
+                        problem.at[idx] = at_interior;
+                    }
+                    if k == 0 && !periodic {
+                        if b.z0 == outlet {
+                            extra += at_outlet;
+                        }
+                    } else if self.w_is_fluid(k, j, i) {
+                        problem.ab[idx] = at_interior;
+                    }
+                    problem.extra_diag[idx] = extra;
+                }
+            }
+        }
+        problem
+    }
+
+    /// The operator with `field.sp` as the right-hand side.
+    pub(crate) fn poisson_problem(&self, field: &Field, dt: f64) -> Problem {
+        let mut problem = self.poisson_operator(field.grid, dt);
+        problem.rhs.copy_from_slice(&field.sp);
+        problem
+    }
+
+    /// The anchor cell of a pure-Neumann projection (the first fluid cell,
+    /// the 2D `(1, 1)` at `k = 0`), or `None` with an outlet.
+    pub(crate) fn anchor_cell(&self, g: Grid) -> Option<usize> {
+        (!self.params.boundaries.any_outlet()).then_some(g.cell(0, 1, 1))
+    }
+
+    /// The inner stop of a projection from the source scale (the 2D rule).
+    pub(crate) fn inner_stop(&self, g: Grid, source_scale: f64) -> f64 {
+        (self.params.inner_stop_factor * source_scale)
+            .max(0.1 * self.params.tolerance * self.reference_flux(g))
+            + 1e-14
+    }
+
+    /// The boundary function on the six sides at every point a predictor or
+    /// the stamping reads (`[side][component]`, see the device driver):
+    /// x sides: u at `(k, j)`, v at `(k, j = 0..=ny)`, w at `(k = 0..=nz, j)`;
+    /// y sides: u at `(k, i = 0..=nx)`, v at `(k, i)`, w at `(k = 0..=nz, i)`;
+    /// z sides: u at `(j, i = 0..=nx)`, v at `(j = 0..=ny, i)`, w at `(j, i)`.
+    #[allow(clippy::type_complexity)]
+    pub fn boundary_tables(&self, g: Grid, t: f64) -> [[Vec<f64>; 3]; 6] {
+        let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz);
+        let xc = |i: usize| (i as f64 + 0.5) * dx;
+        let yc = |j: usize| (j as f64 + 0.5) * dy;
+        let zc = |k: usize| (k as f64 + 0.5) * dz;
+        let mut out: [[Vec<f64>; 3]; 6] = Default::default();
+        for (s, x) in [(0usize, 0.0), (1, nx as f64 * dx)] {
+            let mut u = vec![0.0; ny * nz];
+            let mut v = vec![0.0; (ny + 1) * nz];
+            let mut w = vec![0.0; ny * (nz + 1)];
+            for k in 0..nz {
+                for j in 0..ny {
+                    u[k * ny + j] = self.boundary(x, yc(j), zc(k), t).0;
+                }
+                for j in 0..=ny {
+                    v[k * (ny + 1) + j] = self.boundary(x, j as f64 * dy, zc(k), t).1;
+                }
+            }
+            for k in 0..=nz {
+                for j in 0..ny {
+                    w[k * ny + j] = self.boundary(x, yc(j), k as f64 * dz, t).2;
+                }
+            }
+            out[s] = [u, v, w];
+        }
+        for (s, y) in [(2usize, 0.0), (3, ny as f64 * dy)] {
+            let mut u = vec![0.0; (nx + 1) * nz];
+            let mut v = vec![0.0; nx * nz];
+            let mut w = vec![0.0; nx * (nz + 1)];
+            for k in 0..nz {
+                for i in 0..=nx {
+                    u[k * (nx + 1) + i] = self.boundary(i as f64 * dx, y, zc(k), t).0;
+                }
+                for i in 0..nx {
+                    v[k * nx + i] = self.boundary(xc(i), y, zc(k), t).1;
+                }
+            }
+            for k in 0..=nz {
+                for i in 0..nx {
+                    w[k * nx + i] = self.boundary(xc(i), y, k as f64 * dz, t).2;
+                }
+            }
+            out[s] = [u, v, w];
+        }
+        for (s, z) in [(4usize, 0.0), (5, nz as f64 * dz)] {
+            let mut u = vec![0.0; (nx + 1) * ny];
+            let mut v = vec![0.0; nx * (ny + 1)];
+            let mut w = vec![0.0; nx * ny];
+            for j in 0..ny {
+                for i in 0..=nx {
+                    u[j * (nx + 1) + i] = self.boundary(i as f64 * dx, yc(j), z, t).0;
+                }
+                for i in 0..nx {
+                    w[j * nx + i] = self.boundary(xc(i), yc(j), z, t).2;
+                }
+            }
+            for j in 0..=ny {
+                for i in 0..nx {
+                    v[j * nx + i] = self.boundary(xc(i), j as f64 * dy, z, t).1;
+                }
+            }
+            out[s] = [u, v, w];
+        }
+        out
+    }
+
+    /// The momentum source at every u, v, w face at time `t`.
+    pub fn source_tables(&self, g: Grid, t: f64) -> (Vec<f64>, Vec<f64>, Vec<f64>) {
+        let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz);
+        let mut su = vec![0.0; (nx + 1) * ny * nz];
+        let mut sv = vec![0.0; nx * (ny + 1) * nz];
+        let mut sw = vec![0.0; nx * ny * (nz + 1)];
+        if let Some(f) = self.momentum_source.as_ref() {
+            for k in 0..nz {
+                for j in 0..ny {
+                    for i in 0..=nx {
+                        su[g.uface(k, j, i)] = f(
+                            i as f64 * dx,
+                            (j as f64 + 0.5) * dy,
+                            (k as f64 + 0.5) * dz,
+                            t,
+                        )
+                        .0;
+                    }
+                }
+            }
+            for k in 0..nz {
+                for j in 0..=ny {
+                    for i in 0..nx {
+                        sv[g.vface(k, j, i)] = f(
+                            (i as f64 + 0.5) * dx,
+                            j as f64 * dy,
+                            (k as f64 + 0.5) * dz,
+                            t,
+                        )
+                        .1;
+                    }
+                }
+            }
+            for k in 0..=nz {
+                for j in 0..ny {
+                    for i in 0..nx {
+                        sw[g.wface(k, j, i)] = f(
+                            (i as f64 + 0.5) * dx,
+                            (j as f64 + 0.5) * dy,
+                            k as f64 * dz,
+                            t,
+                        )
+                        .2;
+                    }
+                }
+            }
+        }
+        (su, sv, sw)
+    }
+
+    pub(crate) fn reference_flux(&self, g: Grid) -> f64 {
+        self.fluid.density
+            * self.fluid.reference_velocity
+            * self.fluid.reference_length
+            * (g.nz as f64 * g.dz)
+    }
+
+    /// The solve half of a projection: `sp` from `u*`, then `p'`.
+    pub(crate) fn solve_correction(
+        &mut self,
+        field: &mut Field,
+        dt: f64,
+        warm_start: bool,
+    ) -> PoissonSolution {
+        let g = field.grid;
+        let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz);
+        let rho = self.fluid.density;
+        let mut source_scale = 0.0;
+        for k in 0..nz {
+            for j in 0..ny {
+                for i in 0..nx {
+                    let idx = g.cell(k, j, i);
+                    if !self.cell_is_fluid(k, j, i) {
+                        field.sp[idx] = 0.0;
+                        continue;
+                    }
+                    let divergence_flux = rho
+                        * ((field.u_star[g.uface(k, j, i + 1)] - field.u_star[g.uface(k, j, i)])
+                            * (dy * dz)
+                            + (field.v_star[g.vface(k, j + 1, i)]
+                                - field.v_star[g.vface(k, j, i)])
+                                * (dx * dz)
+                            + (field.w_star[g.wface(k + 1, j, i)]
+                                - field.w_star[g.wface(k, j, i)])
+                                * (dx * dy));
+                    field.sp[idx] = -divergence_flux;
+                    source_scale += divergence_flux.abs();
+                }
+            }
+        }
+        let inner_stop = self.inner_stop(g, source_scale);
+        let problem = self.poisson_problem(field, dt);
+        let mut p_prime = vec![0.0; g.cells()];
+        if warm_start {
+            for k in 0..nz {
+                for j in 0..ny {
+                    for i in 0..nx {
+                        if self.cell_is_fluid(k, j, i) {
+                            let idx = g.cell(k, j, i);
+                            p_prime[idx] = field.p_prime[idx];
+                        }
+                    }
+                }
+            }
+        }
+        let anchor_cell = self.anchor_cell(g);
+        let params = MultigridParameters {
+            precision: self.params.poisson_precision,
+            smoother: self.params.poisson_smoother,
+            ..MultigridParameters::default()
+        };
+        let solution = solve_pcg_cached(
+            &problem,
+            &mut p_prime,
+            &params,
+            inner_stop,
+            anchor_cell,
+            &mut self.pcg_cache,
+        );
+        let (s0, i0, c0, k0) = self.poisson_profile;
+        self.poisson_profile = (
+            s0 + solution.setup_ns,
+            i0 + solution.iterate_ns,
+            c0 + 1,
+            k0 + solution.iterations as u64,
+        );
+        field.p_prime.copy_from_slice(&p_prime);
+        solution
+    }
+
+    /// The apply half: correct the fluid faces from `u*`, add `p'` to `p`,
+    /// return the normalised mass imbalance.
+    pub(crate) fn apply_correction(&self, field: &mut Field, dt: f64) -> f64 {
+        let g = field.grid;
+        let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz);
+        let rho = self.fluid.density;
+        let b = self.params.boundaries;
+        let outlet = Side::PressureOutlet;
+        let periodic = b.periodic_z();
+        let pp = &field.p_prime;
+        for k in 0..nz {
+            for j in 0..ny {
+                for i in 1..nx {
+                    if self.u_is_fluid(k, j, i) {
+                        let dp_dx = (pp[g.cell(k, j, i)] - pp[g.cell(k, j, i - 1)]) / dx;
+                        let f = g.uface(k, j, i);
+                        field.u[f] = field.u_star[f] - (dt / rho) * dp_dx;
+                    }
+                }
+                if b.x0 == outlet {
+                    let dp_dx = (pp[g.cell(k, j, 0)] - 0.0) / (0.5 * dx);
+                    let f = g.uface(k, j, 0);
+                    field.u[f] = field.u_star[f] - (dt / rho) * dp_dx;
+                }
+                if b.x1 == outlet {
+                    let dp_dx = (0.0 - pp[g.cell(k, j, nx - 1)]) / (0.5 * dx);
+                    let f = g.uface(k, j, nx);
+                    field.u[f] = field.u_star[f] - (dt / rho) * dp_dx;
+                }
+            }
+            for i in 0..nx {
+                for j in 1..ny {
+                    if self.v_is_fluid(k, j, i) {
+                        let dp_dy = (pp[g.cell(k, j, i)] - pp[g.cell(k, j - 1, i)]) / dy;
+                        let f = g.vface(k, j, i);
+                        field.v[f] = field.v_star[f] - (dt / rho) * dp_dy;
+                    }
+                }
+                if b.y0 == outlet {
+                    let dp_dy = (pp[g.cell(k, 0, i)] - 0.0) / (0.5 * dy);
+                    let f = g.vface(k, 0, i);
+                    field.v[f] = field.v_star[f] - (dt / rho) * dp_dy;
+                }
+                if b.y1 == outlet {
+                    let dp_dy = (0.0 - pp[g.cell(k, ny - 1, i)]) / (0.5 * dy);
+                    let f = g.vface(k, ny, i);
+                    field.v[f] = field.v_star[f] - (dt / rho) * dp_dy;
+                }
+            }
+        }
+        for j in 0..ny {
+            for i in 0..nx {
+                let k_range = if periodic { 0..nz } else { 1..nz };
+                for k in k_range {
+                    if self.w_is_fluid(k, j, i) {
+                        let below = if k > 0 { k - 1 } else { nz - 1 };
+                        let dp_dz = (pp[g.cell(k, j, i)] - pp[g.cell(below, j, i)]) / dz;
+                        let f = g.wface(k, j, i);
+                        field.w[f] = field.w_star[f] - (dt / rho) * dp_dz;
+                    }
+                }
+                if periodic {
+                    field.w[g.wface(nz, j, i)] = field.w[g.wface(0, j, i)];
+                }
+                if b.z0 == outlet {
+                    let dp_dz = (pp[g.cell(0, j, i)] - 0.0) / (0.5 * dz);
+                    let f = g.wface(0, j, i);
+                    field.w[f] = field.w_star[f] - (dt / rho) * dp_dz;
+                }
+                if b.z1 == outlet {
+                    let dp_dz = (0.0 - pp[g.cell(nz - 1, j, i)]) / (0.5 * dz);
+                    let f = g.wface(nz, j, i);
+                    field.w[f] = field.w_star[f] - (dt / rho) * dp_dz;
+                }
+            }
+        }
+        for k in 0..nz {
+            for j in 0..ny {
+                for i in 0..nx {
+                    if self.cell_is_fluid(k, j, i) {
+                        let idx = g.cell(k, j, i);
+                        field.p[idx] += pp[idx];
+                    }
+                }
+            }
+        }
+        let mut mass_imbalance = 0.0;
+        for k in 0..nz {
+            for j in 0..ny {
+                for i in 0..nx {
+                    if !self.cell_is_fluid(k, j, i) {
+                        continue;
+                    }
+                    let divergence_flux = rho
+                        * ((field.u[g.uface(k, j, i + 1)] - field.u[g.uface(k, j, i)]) * (dy * dz)
+                            + (field.v[g.vface(k, j + 1, i)] - field.v[g.vface(k, j, i)])
+                                * (dx * dz)
+                            + (field.w[g.wface(k + 1, j, i)] - field.w[g.wface(k, j, i)])
+                                * (dx * dy));
+                    mass_imbalance += divergence_flux.abs();
+                }
+            }
+        }
+        let reference_flux = self.reference_flux(g);
+        if reference_flux > 0.0 {
+            mass_imbalance / reference_flux
+        } else {
+            mass_imbalance
+        }
+    }
+}