diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/loads.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/loads.rs
new file mode 100644
index 0000000..4f4a473
--- /dev/null
+++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/loads.rs
@@ -0,0 +1,406 @@
+//! The two load routes of the 2D solver in 3D. Route A, surface stress
+//! reconstruction: at each surface sample the pressure is linearly
+//! extrapolated to the wall from probes at `h` and `2h` along the normal;
+//! the wall-normal derivatives of the normal and two tangential velocity
+//! components from the same probes with the surface velocity at the wall
+//! (quadratic fit); the tangential derivatives of the normal velocity from
+//! the surface-velocity function; `σ·n = (−p + 2μ ∂ₙuₙ) n + μ Σₜ (∂ₙuₜ +
+//! ∂ₜuₙ) t`. Route B, a momentum balance over a box of whole cells
+//! enclosing the body — it reads no near-wall value; the two are unrelated
+//! readings of one solution.
+
+use super::body::Body;
+use super::field::Field;
+use super::wall::{FaceKind, Mask, linear_fit, stencil_nodes, z_planes};
+
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct SurfaceForce {
+    pub f: [f64; 3],
+    pub samples: usize,
+    pub skipped: usize,
+}
+
+fn cross(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
+    [
+        a[1] * b[2] - a[2] * b[1],
+        a[2] * b[0] - a[0] * b[2],
+        a[0] * b[1] - a[1] * b[0],
+    ]
+}
+
+fn dot(a: [f64; 3], b: [f64; 3]) -> f64 {
+    a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
+}
+
+/// Two unit tangents orthonormal to `n`.
+fn tangents(n: [f64; 3]) -> ([f64; 3], [f64; 3]) {
+    let a = if n[0].abs() < 0.9 {
+        [1.0, 0.0, 0.0]
+    } else {
+        [0.0, 1.0, 0.0]
+    };
+    let t1 = cross(n, a);
+    let l = dot(t1, t1).sqrt();
+    let t1 = [t1[0] / l, t1[1] / l, t1[2] / l];
+    (t1, cross(n, t1))
+}
+
+impl Mask {
+    /// Pressure at a point from the cell centres: trilinear when the
+    /// surrounding cells are fluid, else the least-squares fit through the
+    /// fluid ones; `None` if degenerate.
+    pub fn pressure_at(&self, p: &[f64], x: f64, y: f64, z: f64) -> Option<f64> {
+        let g = self.grid();
+        let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz);
+        let (gx, gy, gz) = (x / dx - 0.5, y / dy - 0.5, z / dz - 0.5);
+        let i0 = gx.floor().clamp(0.0, (nx - 2) as f64) as usize;
+        let j0 = gy.floor().clamp(0.0, (ny - 2) as f64) as usize;
+        let fx = (gx - i0 as f64).clamp(0.0, 1.0);
+        let fy = (gy - j0 as f64).clamp(0.0, 1.0);
+        let (k0, k1, fz) = z_planes(gz, nz, self.periodic_z());
+        // A query on a single weighted plane lives on it (the fit has no
+        // z variation to determine).
+        let z = if nz == 1 || fz == 0.0 {
+            (k0 as f64 + 0.5) * dz
+        } else if fz == 1.0 {
+            (k1 as f64 + 0.5) * dz
+        } else {
+            z
+        };
+        let dirs = |m: usize, f: f64| {
+            if m <= 1 {
+                vec![(0usize, 1.0)]
+            } else {
+                vec![(0, 1.0 - f), (1, f)]
+            }
+        };
+        let kz = |dk: usize| if dk == 0 { k0 } else { k1 };
+        let mut nodes = Vec::with_capacity(8);
+        for (dk, wk) in dirs(nz, fz) {
+            for (dj, wj) in dirs(ny, fy) {
+                for (di, wi) in dirs(nx, fx) {
+                    nodes.push((g.cell(kz(dk), j0 + dj, i0 + di), wi * wj * wk));
+                }
+            }
+        }
+        if nodes.iter().all(|&(idx, _)| self.is_fluid_cell(idx)) {
+            return Some(nodes.iter().map(|&(idx, w)| w * p[idx]).sum());
+        }
+        // Weighted by the z-direction weight (a z-invariant field then fits
+        // as the 2D four-cell fit); the z weight of node `(dk, dj, di)` is
+        // `wk`, recovered from the trilinear product.
+        let zw = |wk: f64| wk;
+        let mut pts: Vec<(f64, f64, f64, f64, f64)> = Vec::new();
+        let mut it = nodes.iter();
+        for (dk, wk) in dirs(nz, fz) {
+            for _ in 0..dirs(ny, fy).len() * dirs(nx, fx).len() {
+                let &(idx, _) = it.next().expect("node");
+                let _ = dk;
+                if self.is_fluid_cell(idx) {
+                    let (k, j, i) = g.kji(idx);
+                    pts.push((
+                        (i as f64 + 0.5) * dx,
+                        (j as f64 + 0.5) * dy,
+                        (k as f64 + 0.5) * dz,
+                        p[idx],
+                        zw(wk),
+                    ));
+                }
+            }
+        }
+        linear_fit(&pts, (x, y, z))
+    }
+
+    /// Velocity at a point: trilinear over a component's nodes when all are
+    /// fluid faces, else the fit through the fluid ones and the point's own
+    /// boundary intercept with its surface velocity.
+    pub fn velocity_at(
+        &self,
+        body: &Body,
+        f: &Field,
+        x: f64,
+        y: f64,
+        z: f64,
+        t: f64,
+    ) -> Option<[f64; 3]> {
+        let g = self.grid();
+        let eps = 1e-6 * g.dx.min(g.dy).min(g.dz);
+        let s = body.phi(x, y, z, t);
+        let (n1, n2, n3) = body.normal(x, y, z, t, eps);
+        let foot = (x - s * n1, y - s * n2, z - s * n3);
+        let vf = body.surface_velocity(foot.0, foot.1, foot.2, t);
+        let vel_foot = [vf.0, vf.1, vf.2];
+        let mut out = [0.0; 3];
+        for c in 0..3 {
+            // A query on a single weighted plane of this component's nodes
+            // lives on it (the fit has no z variation to determine).
+            let (planes, offset) = if c == 2 {
+                (if self.periodic_z() { g.nz } else { g.nz + 1 }, 0.0)
+            } else {
+                (g.nz, 0.5)
+            };
+            let gz = z / g.dz - offset;
+            let (k0, k1, fz) = z_planes(gz, planes, self.periodic_z());
+            let zq = if planes <= 1 || fz == 0.0 {
+                (k0 as f64 + offset) * g.dz
+            } else if fz == 1.0 {
+                (k1 as f64 + offset) * g.dz
+            } else {
+                z
+            };
+            let foot_c = (foot.0, foot.1, if zq == z { foot.2 } else { zq });
+            let nodes = stencil_nodes((x, y, zq), c, g, self.periodic_z(), |_| None);
+            let values: &[f64] = [&f.u, &f.v, &f.w][c];
+            let fluid = |idx: usize| match c {
+                0 => self.u_kind(idx) == FaceKind::Fluid,
+                1 => self.v_kind(idx) == FaceKind::Fluid,
+                _ => self.w_kind(idx) == FaceKind::Fluid,
+            };
+            if nodes.iter().all(|n| fluid(n.idx)) {
+                out[c] = nodes.iter().map(|n| n.weight * values[n.idx]).sum();
+            } else {
+                let mut pts: Vec<(f64, f64, f64, f64, f64)> = nodes
+                    .iter()
+                    .filter(|n| fluid(n.idx))
+                    .map(|n| (n.x, n.y, n.z, values[n.idx], n.zw))
+                    .collect();
+                pts.push((foot_c.0, foot_c.1, foot_c.2, vel_foot[c], 1.0));
+                out[c] = linear_fit(&pts, (x, y, zq))?;
+            }
+        }
+        Some(out)
+    }
+
+    /// The reconstructed traction at one surface point with outward normal `n`.
+    pub fn traction_at(
+        &self,
+        body: &Body,
+        f: &Field,
+        mu: f64,
+        t: f64,
+        x: [f64; 3],
+        n: [f64; 3],
+    ) -> Option<[f64; 3]> {
+        let g = self.grid();
+        let h = g.dx.min(g.dy).min(g.dz);
+        let (d1, d2) = (h, 2.0 * h);
+        let at = |d: f64| [x[0] + d * n[0], x[1] + d * n[1], x[2] + d * n[2]];
+        let (x1, x2) = (at(d1), at(d2));
+        let p1 = self.pressure_at(&f.p, x1[0], x1[1], x1[2])?;
+        let p2 = self.pressure_at(&f.p, x2[0], x2[1], x2[2])?;
+        let u1 = self.velocity_at(body, f, x1[0], x1[1], x1[2], t)?;
+        let u2 = self.velocity_at(body, f, x2[0], x2[1], x2[2], t)?;
+        let p_wall = p1 + (p1 - p2) * d1 / (d2 - d1);
+        let (t1, t2) = tangents(n);
+        let us = body.surface_velocity(x[0], x[1], x[2], t);
+        let us = [us.0, us.1, us.2];
+        let wall_gradient =
+            |f1: f64, f2: f64| (f1 * d2 * d2 - f2 * d1 * d1) / (d1 * d2 * (d2 - d1));
+        let dn =
+            |dir: [f64; 3]| wall_gradient(dot(u1, dir) - dot(us, dir), dot(u2, dir) - dot(us, dir));
+        let eps = 1e-6 * h;
+        let dt_un = |dir: [f64; 3]| {
+            let p = body.surface_velocity(
+                x[0] + eps * dir[0],
+                x[1] + eps * dir[1],
+                x[2] + eps * dir[2],
+                t,
+            );
+            let m = body.surface_velocity(
+                x[0] - eps * dir[0],
+                x[1] - eps * dir[1],
+                x[2] - eps * dir[2],
+                t,
+            );
+            ((p.0 - m.0) * n[0] + (p.1 - m.1) * n[1] + (p.2 - m.2) * n[2]) / (2.0 * eps)
+        };
+        let traction_n = -p_wall + 2.0 * mu * dn(n);
+        let traction_t1 = mu * (dn(t1) + dt_un(t1));
+        let traction_t2 = mu * (dn(t2) + dt_un(t2));
+        Some([
+            traction_n * n[0] + traction_t1 * t1[0] + traction_t2 * t2[0],
+            traction_n * n[1] + traction_t1 * t1[1] + traction_t2 * t2[1],
+            traction_n * n[2] + traction_t1 * t1[2] + traction_t2 * t2[2],
+        ])
+    }
+
+    /// Route A: the surface integral of the reconstructed traction over the
+    /// body's samples at spacing `ds`.
+    pub fn surface_force(&self, body: &Body, f: &Field, mu: f64, t: f64, ds: f64) -> SurfaceForce {
+        let samples = body.surface_samples(ds);
+        let mut force = [0.0; 3];
+        let mut skipped = 0;
+        for s in &samples {
+            match self.traction_at(body, f, mu, t, [s.x, s.y, s.z], [s.nx, s.ny, s.nz]) {
+                Some(tr) => {
+                    for c in 0..3 {
+                        force[c] += tr[c] * s.area;
+                    }
+                }
+                None => skipped += 1,
+            }
+        }
+        SurfaceForce {
+            f: force,
+            samples: samples.len(),
+            skipped,
+        }
+    }
+
+    /// Route B: the momentum balance over the box of whole cells
+    /// `[i0, i1) × [j0, j1) × [k0, k1)` (in the fluid on its boundary):
+    /// `F = Σ_outer (σ·n − ρ u (u·n)) A − d/dt ∫ ρ u dV + ∫ f dV`.
+    #[allow(clippy::too_many_arguments)]
+    pub fn control_volume_force(
+        &self,
+        f: &Field,
+        dt: f64,
+        rho: f64,
+        mu: f64,
+        source: Option<&dyn Fn(f64, f64, f64) -> (f64, f64, f64)>,
+        (i0, i1, j0, j1, k0, k1): (usize, usize, usize, usize, usize, usize),
+    ) -> [f64; 3] {
+        let g = self.grid();
+        let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz);
+        let (u, v, w, p) = (&f.u, &f.v, &f.w, &f.p);
+        let mut force = [0.0; 3];
+        let uc =
+            |k: usize, j: usize, i: usize| 0.5 * (u[g.uface(k, j, i)] + u[g.uface(k, j, i + 1)]);
+        let vc =
+            |k: usize, j: usize, i: usize| 0.5 * (v[g.vface(k, j, i)] + v[g.vface(k, j + 1, i)]);
+        let wc =
+            |k: usize, j: usize, i: usize| 0.5 * (w[g.wface(k, j, i)] + w[g.wface(k + 1, j, i)]);
+        // Central inside, one-sided at the domain edge.
+        let dd = |prev: Option<f64>, here: f64, next: Option<f64>, h: f64| match (prev, next) {
+            (Some(a), Some(b)) => (b - a) / (2.0 * h),
+            (None, Some(b)) => (b - here) / h,
+            (Some(a), None) => (here - a) / h,
+            (None, None) => 0.0,
+        };
+        // x faces: u lives there.
+        for k in k0..k1 {
+            for j in j0..j1 {
+                for (i, sign) in [(i1, 1.0), (i0, -1.0)] {
+                    let un = u[g.uface(k, j, i)];
+                    let p_f = 0.5 * (p[g.cell(k, j, i - 1)] + p[g.cell(k, j, i)]);
+                    let dudx = (u[g.uface(k, j, i + 1)] - u[g.uface(k, j, i - 1)]) / (2.0 * dx);
+                    let dvdx = (vc(k, j, i) - vc(k, j, i - 1)) / dx;
+                    let dwdx = (wc(k, j, i) - wc(k, j, i - 1)) / dx;
+                    let dudy = dd(
+                        (j > 0).then(|| u[g.uface(k, j - 1, i)]),
+                        un,
+                        (j + 1 < ny).then(|| u[g.uface(k, j + 1, i)]),
+                        dy,
+                    );
+                    let dudz = dd(
+                        (k > 0).then(|| u[g.uface(k - 1, j, i)]),
+                        un,
+                        (k + 1 < nz).then(|| u[g.uface(k + 1, j, i)]),
+                        dz,
+                    );
+                    let v_f = 0.5 * (vc(k, j, i - 1) + vc(k, j, i));
+                    let w_f = 0.5 * (wc(k, j, i - 1) + wc(k, j, i));
+                    let a = dy * dz;
+                    force[0] += sign * ((-p_f + 2.0 * mu * dudx) - rho * un * un) * a;
+                    force[1] += sign * (mu * (dudy + dvdx) - rho * v_f * un) * a;
+                    force[2] += sign * (mu * (dudz + dwdx) - rho * w_f * un) * a;
+                }
+            }
+        }
+        // y faces: v lives there.
+        for k in k0..k1 {
+            for i in i0..i1 {
+                for (j, sign) in [(j1, 1.0), (j0, -1.0)] {
+                    let vn = v[g.vface(k, j, i)];
+                    let p_f = 0.5 * (p[g.cell(k, j - 1, i)] + p[g.cell(k, j, i)]);
+                    let dvdy = (v[g.vface(k, j + 1, i)] - v[g.vface(k, j - 1, i)]) / (2.0 * dy);
+                    let dudy = (uc(k, j, i) - uc(k, j - 1, i)) / dy;
+                    let dwdy = (wc(k, j, i) - wc(k, j - 1, i)) / dy;
+                    let dvdx = dd(
+                        (i > 0).then(|| v[g.vface(k, j, i - 1)]),
+                        vn,
+                        (i + 1 < nx).then(|| v[g.vface(k, j, i + 1)]),
+                        dx,
+                    );
+                    let dvdz = dd(
+                        (k > 0).then(|| v[g.vface(k - 1, j, i)]),
+                        vn,
+                        (k + 1 < nz).then(|| v[g.vface(k + 1, j, i)]),
+                        dz,
+                    );
+                    let u_f = 0.5 * (uc(k, j - 1, i) + uc(k, j, i));
+                    let w_f = 0.5 * (wc(k, j - 1, i) + wc(k, j, i));
+                    let a = dx * dz;
+                    force[0] += sign * (mu * (dudy + dvdx) - rho * u_f * vn) * a;
+                    force[1] += sign * ((-p_f + 2.0 * mu * dvdy) - rho * vn * vn) * a;
+                    force[2] += sign * (mu * (dvdz + dwdy) - rho * w_f * vn) * a;
+                }
+            }
+        }
+        // z faces: w lives there. A box spanning the whole z range has its
+        // z faces on the domain's z sides: no momentum flux and no shear on
+        // a slip wall, and the two faces cancel on a periodic pair — skipped.
+        let full_span = k0 == 0 && k1 == nz;
+        for j in (j0..j1).filter(|_| !full_span) {
+            for i in i0..i1 {
+                for (k, sign) in [(k1, 1.0), (k0, -1.0)] {
+                    let wn = w[g.wface(k, j, i)];
+                    let p_f = 0.5 * (p[g.cell(k - 1, j, i)] + p[g.cell(k, j, i)]);
+                    let dwdz = (w[g.wface(k + 1, j, i)] - w[g.wface(k - 1, j, i)]) / (2.0 * dz);
+                    let dudz = (uc(k, j, i) - uc(k - 1, j, i)) / dz;
+                    let dvdz = (vc(k, j, i) - vc(k - 1, j, i)) / dz;
+                    let dwdx = dd(
+                        (i > 0).then(|| w[g.wface(k, j, i - 1)]),
+                        wn,
+                        (i + 1 < nx).then(|| w[g.wface(k, j, i + 1)]),
+                        dx,
+                    );
+                    let dwdy = dd(
+                        (j > 0).then(|| w[g.wface(k, j - 1, i)]),
+                        wn,
+                        (j + 1 < ny).then(|| w[g.wface(k, j + 1, i)]),
+                        dy,
+                    );
+                    let u_f = 0.5 * (uc(k - 1, j, i) + uc(k, j, i));
+                    let v_f = 0.5 * (vc(k - 1, j, i) + vc(k, j, i));
+                    let a = dx * dy;
+                    force[0] += sign * (mu * (dudz + dwdx) - rho * u_f * wn) * a;
+                    force[1] += sign * (mu * (dvdz + dwdy) - rho * v_f * wn) * a;
+                    force[2] += sign * ((-p_f + 2.0 * mu * dwdz) - rho * wn * wn) * a;
+                }
+            }
+        }
+        // Unsteady term and source over the fluid cells of the box.
+        let dv = dx * dy * dz;
+        for k in k0..k1 {
+            for j in j0..j1 {
+                for i in i0..i1 {
+                    let idx = g.cell(k, j, i);
+                    if !self.is_fluid_cell(idx) {
+                        continue;
+                    }
+                    let (fu0, fu1) = (g.uface(k, j, i), g.uface(k, j, i + 1));
+                    let (fv0, fv1) = (g.vface(k, j, i), g.vface(k, j + 1, i));
+                    let (fw0, fw1) = (g.wface(k, j, i), g.wface(k + 1, j, i));
+                    let du = 0.5 * ((u[fu0] - f.u_old[fu0]) + (u[fu1] - f.u_old[fu1]));
+                    let dvv = 0.5 * ((v[fv0] - f.v_old[fv0]) + (v[fv1] - f.v_old[fv1]));
+                    let dw = 0.5 * ((w[fw0] - f.w_old[fw0]) + (w[fw1] - f.w_old[fw1]));
+                    force[0] -= rho * du / dt * dv;
+                    force[1] -= rho * dvv / dt * dv;
+                    force[2] -= rho * dw / dt * dv;
+                    if let Some(s) = source {
+                        let (sx, sy, sz) = s(
+                            (i as f64 + 0.5) * dx,
+                            (j as f64 + 0.5) * dy,
+                            (k as f64 + 0.5) * dz,
+                        );
+                        force[0] += sx * dv;
+                        force[1] += sy * dv;
+                        force[2] += sz * dv;
+                    }
+                }
+            }
+        }
+        force
+    }
+}
diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs
index a0f9f97..0b1227f 100644
--- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs
+++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs
@@ -10,11 +10,15 @@ pub mod body;
 pub mod cut;
 pub mod field;
 pub mod grid;
+pub mod loads;
 pub mod poisson;
 pub mod step;
+pub mod wall;
 
 pub use body::{Body, SurfaceSample};
 pub use cut::CutGeometry;
 pub use field::Field;
 pub use grid::Grid;
+pub use loads::SurfaceForce;
 pub use step::{Boundaries, Fluid, Parameters, Side, Solver, StepResult};
+pub use wall::{FaceKind, Mask, WallScheme};
diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs
index 9a176c6..1c1e5ca 100644
--- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs
+++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs
@@ -10,8 +10,10 @@ mod predictor;
 mod projection;
 
 use super::Grid;
+use super::body::Body;
 use super::field::Field;
 use super::poisson::{PcgCache, Problem, solve_pcg_cached};
+use super::wall::{FaceKind, Mask, WallScheme};
 use crate::solvers::incompressible::poisson::{
     MgPrecision, MgSmoother, MultigridParameters, PoissonSolution,
 };
@@ -40,7 +42,7 @@ pub struct Boundaries {
 }
 
 impl Boundaries {
-    fn any_outlet(self) -> bool {
+    pub(crate) fn any_outlet(self) -> bool {
         [self.x0, self.x1, self.y0, self.y1, self.z0, self.z1].contains(&Side::PressureOutlet)
     }
 
@@ -69,6 +71,8 @@ pub struct Parameters {
     pub convection_scheme: ConvectionScheme,
     /// Relative part of the pressure solve's inner stop (the 2D 1e-2).
     pub inner_stop_factor: f64,
+    /// The wall treatment of an embedded body.
+    pub wall_scheme: WallScheme,
 }
 
 impl Default for Parameters {
@@ -81,6 +85,7 @@ impl Default for Parameters {
             poisson_precision: MgPrecision::F64,
             convection_scheme: ConvectionScheme::Upwind,
             inner_stop_factor: 1e-2,
+            wall_scheme: WallScheme::GhostBinary,
         }
     }
 }
@@ -102,6 +107,9 @@ pub struct Solver {
     pub params: Parameters,
     pub(super) momentum_source: Option<Vec3Fn>,
     boundary_velocity: Option<Vec3Fn>,
+    body: Option<Body>,
+    mask: Option<Mask>,
+    last_ghost_correction: f64,
     pcg_cache: PcgCache,
     time: f64,
     initialized: bool,
@@ -126,6 +134,9 @@ impl Solver {
             params,
             momentum_source: None,
             boundary_velocity: None,
+            body: None,
+            mask: None,
+            last_ghost_correction: 0.0,
             pcg_cache: PcgCache::default(),
             time: 0.0,
             initialized: false,
@@ -147,6 +158,28 @@ impl Solver {
         self.boundary_velocity = Some(Box::new(f));
     }
 
+    /// A static embedded body (the mask is built at initialisation).
+    pub fn set_body(&mut self, body: Body) {
+        self.body = Some(body);
+        self.mask = None;
+    }
+
+    #[must_use]
+    pub fn body(&self) -> Option<&Body> {
+        self.body.as_ref()
+    }
+
+    #[must_use]
+    pub fn mask(&self) -> Option<&Mask> {
+        self.mask.as_ref()
+    }
+
+    /// The last step's ghost compatibility correction.
+    #[must_use]
+    pub fn ghost_correction(&self) -> f64 {
+        self.last_ghost_correction
+    }
+
     #[must_use]
     pub fn time(&self) -> f64 {
         self.time
@@ -167,22 +200,30 @@ impl Solver {
             .map_or((0.0, 0.0, 0.0), |f| f(x, y, z, t))
     }
 
-    // The fluid predicates: everything is fluid until the wall arrives.
+    // The fluid predicates (everything is fluid without a body).
     #[inline]
-    fn u_is_fluid(&self, _k: usize, _j: usize, _i: usize) -> bool {
-        true
+    pub(super) fn u_is_fluid(&self, k: usize, j: usize, i: usize) -> bool {
+        self.mask
+            .as_ref()
+            .is_none_or(|m| m.u_kind(m.grid().uface(k, j, i)) == FaceKind::Fluid)
     }
     #[inline]
-    fn v_is_fluid(&self, _k: usize, _j: usize, _i: usize) -> bool {
-        true
+    pub(super) fn v_is_fluid(&self, k: usize, j: usize, i: usize) -> bool {
+        self.mask
+            .as_ref()
+            .is_none_or(|m| m.v_kind(m.grid().vface(k, j, i)) == FaceKind::Fluid)
     }
     #[inline]
-    fn w_is_fluid(&self, _k: usize, _j: usize, _i: usize) -> bool {
-        true
+    pub(super) fn w_is_fluid(&self, k: usize, j: usize, i: usize) -> bool {
+        self.mask
+            .as_ref()
+            .is_none_or(|m| m.w_kind(m.grid().wface(k, j, i)) == FaceKind::Fluid)
     }
     #[inline]
-    fn cell_is_fluid(&self, _k: usize, _j: usize, _i: usize) -> bool {
-        true
+    pub(super) fn cell_is_fluid(&self, k: usize, j: usize, i: usize) -> bool {
+        self.mask
+            .as_ref()
+            .is_none_or(|m| m.is_fluid_cell(m.grid().cell(k, j, i)))
     }
 
     pub(super) fn upwind(face_velocity: f64, upstream: f64, downstream: f64) -> f64 {
@@ -318,10 +359,27 @@ impl Solver {
         field.copy_to_starred();
     }
 
-    /// Stamp the `t = time` boundary data (lazily called by the first step).
+    /// Build the mask (with a body) and stamp the `t = time` boundary data
+    /// and ghost values (lazily called by the first step).
     pub fn initialize(&mut self, field: &mut Field) {
         let t = self.time;
+        if let Some(body) = &self.body {
+            if self.mask.is_none() {
+                assert_eq!(
+                    self.params.wall_scheme,
+                    WallScheme::GhostBinary,
+                    "item 10 brings CutCell"
+                );
+                self.mask = Some(
+                    Mask::build(body, field.grid, t, self.params.boundaries)
+                        .expect("embedded mask"),
+                );
+            }
+        }
         self.apply_boundary_normals(field, t);
+        if let (Some(body), Some(mask)) = (&self.body, &self.mask) {
+            mask.impose(body, &mut field.u, &mut field.v, &mut field.w, t);
+        }
         self.initialized = true;
     }
 
@@ -354,6 +412,11 @@ impl Solver {
             }
             field.copy_to_starred();
         }
+        // Ghost faces follow the corrected field (the next step's stencil data).
+        if let (Some(body), Some(mask)) = (&self.body, &self.mask) {
+            self.last_ghost_correction =
+                mask.impose(body, &mut field.u, &mut field.v, &mut field.w, t_new);
+        }
         self.time = t_new;
         StepResult {
             converged: final_residual < self.params.tolerance,
diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/projection.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/projection.rs
index f9bf472..860c8c2 100644
--- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/projection.rs
+++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/projection.rs
@@ -93,7 +93,11 @@ impl Solver {
     /// 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))
+        (!self.params.boundaries.any_outlet()).then(|| {
+            self.mask
+                .as_ref()
+                .map_or(g.cell(0, 1, 1), super::super::wall::Mask::anchor)
+        })
     }
 
     /// The inner stop of a projection from the source scale (the 2D rule).
diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs
new file mode 100644
index 0000000..fbeda60
--- /dev/null
+++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs
@@ -0,0 +1,655 @@
+//! The binary ghost wall — the 2D `EmbeddedMask` on the 3D grid. Cells are
+//! fluid where φ(centre) > 0; an interior face between two cells that are
+//! not both fluid is a ghost (within 1.5 h of the surface) or solid; a
+//! ghost face is prescribed the least-squares linear reconstruction
+//! through the fluid nodes of the trilinear stencil around the probe (one
+//! cell beyond the face's mirror image along the normal) plus the foot of
+//! the normal with its surface velocity — exact for linear fields — with
+//! the profile along the normal as the fallback; the ghost faces bounding
+//! fluid cells share a flux compatibility correction.
+
+use super::Grid;
+use super::body::Body;
+use super::step::{Boundaries, Side};
+
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub enum FaceKind {
+    Fluid,
+    Ghost,
+    Solid,
+}
+
+/// The wall treatment of an embedded body.
+#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
+pub enum WallScheme {
+    /// The 2D binary mask ported (this module).
+    #[default]
+    GhostBinary,
+    /// The apertured cut-cell wall (item 10).
+    CutCell,
+}
+
+#[derive(Debug, Clone, Copy)]
+pub(crate) struct StencilNode {
+    pub(crate) idx: usize,
+    pub(crate) x: f64,
+    pub(crate) y: f64,
+    pub(crate) z: f64,
+    /// The trilinear weight (the fallback profile's interpolation).
+    pub(crate) weight: f64,
+    /// The z-direction weight alone: the least-squares weight of the node,
+    /// so that a z-invariant field fits exactly as the 2D four-node fit
+    /// (the two planes share the in-plane nodes' unit weight).
+    pub(crate) zw: f64,
+    pub(crate) fallback: Option<f64>,
+}
+
+#[derive(Debug, Clone)]
+struct Ghost {
+    idx: usize,
+    x: f64,
+    y: f64,
+    z: f64,
+    foot: (f64, f64, f64),
+    u_surface: f64,
+    s_face: f64,
+    s_probe: f64,
+    nodes: Vec<StencilNode>,
+    /// Outward-from-fluid sign for the compatibility correction (0 when no
+    /// fluid cell is adjacent).
+    flux_sign: f64,
+}
+
+#[derive(Clone)]
+pub struct Mask {
+    grid: Grid,
+    periodic_z: bool,
+    cell_fluid: Vec<bool>,
+    u_kind: Vec<FaceKind>,
+    v_kind: Vec<FaceKind>,
+    w_kind: Vec<FaceKind>,
+    u_ghosts: Vec<Ghost>,
+    v_ghosts: Vec<Ghost>,
+    w_ghosts: Vec<Ghost>,
+    anchor: usize,
+    fluid_cells: usize,
+}
+
+/// The z lattice position of a query: the lower plane index, the upper
+/// plane index and the weight of the upper plane. Periodic z wraps; a wall
+/// clamps (and a query outside the lattice takes the nearest plane).
+pub(crate) fn z_planes(gz: f64, planes: usize, periodic: bool) -> (usize, usize, f64) {
+    if planes <= 1 {
+        return (0, 0, 0.0);
+    }
+    if periodic {
+        let n = planes as f64;
+        let w = gz.rem_euclid(n);
+        let k0 = w.floor() as usize % planes;
+        ((k0) % planes, (k0 + 1) % planes, w - w.floor())
+    } else {
+        let k0 = gz.floor().clamp(0.0, (planes - 2) as f64) as usize;
+        (k0, k0 + 1, (gz - k0 as f64).clamp(0.0, 1.0))
+    }
+}
+
+/// The nodes of velocity component `c` (0 u, 1 v, 2 w) around `point`
+/// with trilinear weights (a single plane of nodes in a direction
+/// collapses to one node), clamped in x and y, wrapped or clamped in z;
+/// `fallback(idx)` supplies the value of a node that is not a fluid face.
+pub(crate) fn stencil_nodes(
+    point: (f64, f64, f64),
+    c: usize,
+    g: Grid,
+    periodic_z: bool,
+    fallback: impl Fn(usize) -> Option<f64>,
+) -> Vec<StencilNode> {
+    let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz);
+    // w faces on a periodic grid: the plane k = nz is the plane 0.
+    let (gx, gy, gz, max_i, max_j, planes_k) = match c {
+        0 => (
+            point.0 / dx,
+            point.1 / dy - 0.5,
+            point.2 / dz - 0.5,
+            nx,
+            ny - 1,
+            nz,
+        ),
+        1 => (
+            point.0 / dx - 0.5,
+            point.1 / dy,
+            point.2 / dz - 0.5,
+            nx - 1,
+            ny,
+            nz,
+        ),
+        _ => (
+            point.0 / dx - 0.5,
+            point.1 / dy - 0.5,
+            point.2 / dz,
+            nx - 1,
+            ny - 1,
+            if periodic_z { nz } else { nz + 1 },
+        ),
+    };
+    let i0 = gx.floor().clamp(0.0, (max_i.max(1) - 1) as f64) as usize;
+    let j0 = gy.floor().clamp(0.0, (max_j.max(1) - 1) as f64) as usize;
+    let fx = (gx - i0 as f64).clamp(0.0, 1.0);
+    let fy = (gy - j0 as f64).clamp(0.0, 1.0);
+    let (k0, k1, fz) = z_planes(gz, planes_k, periodic_z);
+    let pos = |k: usize, j: usize, i: usize| match c {
+        0 => (i as f64 * dx, (j as f64 + 0.5) * dy, (k as f64 + 0.5) * dz),
+        1 => ((i as f64 + 0.5) * dx, j as f64 * dy, (k as f64 + 0.5) * dz),
+        _ => ((i as f64 + 0.5) * dx, (j as f64 + 0.5) * dy, k as f64 * dz),
+    };
+    let index = |k: usize, j: usize, i: usize| match c {
+        0 => g.uface(k, j, i),
+        1 => g.vface(k, j, i),
+        _ => g.wface(k, j, i),
+    };
+    let dirs = |m: usize, f: f64| {
+        if m <= 1 {
+            vec![(0usize, 1.0)]
+        } else {
+            vec![(0, 1.0 - f), (1, f)]
+        }
+    };
+    let mut out = Vec::with_capacity(8);
+    let kz = |dk: usize| if dk == 0 { k0 } else { k1 };
+    for (dk, wk) in dirs(planes_k, fz) {
+        for (dj, wj) in dirs(max_j, fy) {
+            for (di, wi) in dirs(max_i, fx) {
+                let (k, j, i) = (kz(dk), j0 + dj, i0 + di);
+                let (x, y, z) = pos(k, j, i);
+                let idx = index(k, j, i);
+                out.push(StencilNode {
+                    idx,
+                    x,
+                    y,
+                    z,
+                    weight: wi * wj * wk,
+                    zw: wk,
+                    fallback: fallback(idx),
+                });
+            }
+        }
+    }
+    out
+}
+
+/// Weighted least-squares fit of `a + b (x−x0) + c (y−y0) + d (z−z0)`
+/// through `pts` (each `(x, y, z, value, weight)`; zero-weight points drop
+/// out) evaluated at `at`; directions in which every weighted point has
+/// the same coordinate are dropped. `None` if the normal matrix is
+/// singular relative to its scale.
+pub(crate) fn linear_fit(pts_w: &[(f64, f64, f64, f64, f64)], at: (f64, f64, f64)) -> Option<f64> {
+    let pts: Vec<(f64, f64, f64, f64, f64)> = pts_w.iter().copied().filter(|p| p.4 > 0.0).collect();
+    let spread = |f: &dyn Fn(&(f64, f64, f64, f64, f64)) -> f64| {
+        let (lo, hi) = pts
+            .iter()
+            .fold((f64::INFINITY, f64::NEG_INFINITY), |(lo, hi), p| {
+                (lo.min(f(p)), hi.max(f(p)))
+            });
+        hi - lo > 1e-12 * (hi.abs() + lo.abs() + 1e-300)
+    };
+    let on = [true, spread(&|p| p.0), spread(&|p| p.1), spread(&|p| p.2)];
+    let cols: Vec<usize> = (0..4).filter(|&c| on[c]).collect();
+    let m = cols.len();
+    if pts.len() < m {
+        return None;
+    }
+    let mut a = vec![vec![0.0; m + 1]; m];
+    for p in &pts {
+        let full = [1.0, p.0 - at.0, p.1 - at.1, p.2 - at.2];
+        let r: Vec<f64> = cols.iter().map(|&c| full[c]).collect();
+        for i in 0..m {
+            for j in 0..m {
+                a[i][j] += p.4 * r[i] * r[j];
+            }
+            a[i][m] += p.4 * r[i] * p.3;
+        }
+    }
+    let scale: f64 = (0..m).map(|i| a[i][i]).product();
+    if scale <= 0.0 {
+        return None;
+    }
+    let mut det = 1.0;
+    for col in 0..m {
+        let piv = (col..m)
+            .max_by(|&p, &q| a[p][col].abs().partial_cmp(&a[q][col].abs()).unwrap())
+            .unwrap();
+        a.swap(col, piv);
+        let d = a[col][col];
+        det *= d;
+        if d.abs() <= 1e-14 * scale.powf(1.0 / m as f64) {
+            return None;
+        }
+        for r in col + 1..m {
+            let f = a[r][col] / d;
+            for c in col..=m {
+                a[r][c] -= f * a[col][c];
+            }
+        }
+    }
+    if det.abs() <= 1e-10 * scale {
+        return None;
+    }
+    let mut x = vec![0.0; m];
+    for i in (0..m).rev() {
+        let mut s = a[i][m];
+        for j in i + 1..m {
+            s -= a[i][j] * x[j];
+        }
+        x[i] = s / a[i][i];
+    }
+    Some(x[0])
+}
+
+impl Ghost {
+    fn reconstruct(&self, values: &[f64]) -> f64 {
+        let mut pts: Vec<(f64, f64, f64, f64, f64)> = self
+            .nodes
+            .iter()
+            .filter(|n| n.fallback.is_none())
+            .map(|n| (n.x, n.y, n.z, values[n.idx], n.zw))
+            .collect();
+        pts.push((self.foot.0, self.foot.1, self.foot.2, self.u_surface, 1.0));
+        if let Some(val) = linear_fit(&pts, (self.x, self.y, self.z)) {
+            return val;
+        }
+        let mut probe = 0.0;
+        for n in &self.nodes {
+            probe += n.weight * n.fallback.unwrap_or_else(|| values[n.idx]);
+        }
+        self.u_surface + (probe - self.u_surface) * (self.s_face / self.s_probe)
+    }
+}
+
+impl Mask {
+    /// Classify the grid against `body` at `t`. A solid cell on a domain
+    /// side is refused unless that side is a Velocity, SlipWall or Periodic
+    /// side (an outlet may not be blocked).
+    pub fn build(body: &Body, g: Grid, t: f64, b: Boundaries) -> Result<Self, String> {
+        let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz);
+        let periodic = b.z0 == Side::Periodic;
+        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 cell_fluid = vec![true; g.cells()];
+        let mut fluid_cells = 0;
+        let mut anchor = None;
+        let allowed = |side: Side| matches!(side, Side::Velocity | Side::Periodic | Side::SlipWall);
+        for k in 0..nz {
+            for j in 0..ny {
+                for i in 0..nx {
+                    let fluid = body.phi(xc(i), yc(j), zc(k), t) > 0.0;
+                    let idx = g.cell(k, j, i);
+                    cell_fluid[idx] = fluid;
+                    if fluid {
+                        fluid_cells += 1;
+                        if anchor.is_none() {
+                            anchor = Some(idx);
+                        }
+                    } else {
+                        let touches = (i == 0 && !allowed(b.x0))
+                            || (i + 1 == nx && !allowed(b.x1))
+                            || (j == 0 && !allowed(b.y0))
+                            || (j + 1 == ny && !allowed(b.y1))
+                            || (k == 0 && !allowed(b.z0))
+                            || (k + 1 == nz && !allowed(b.z1));
+                        if touches {
+                            return Err(format!(
+                                "embedded body reaches a domain side that is not a Velocity/Periodic side at cell ({k}, {j}, {i})"
+                            ));
+                        }
+                    }
+                }
+            }
+        }
+        let Some(anchor) = anchor else {
+            return Err("embedded body covers the whole domain".into());
+        };
+        let h_min = dx.min(dy).min(dz);
+        let reach = 1.5 * h_min;
+        let eps = 1e-6 * h_min;
+        let is_fluid = |k: usize, j: usize, i: usize| cell_fluid[g.cell(k, j, i)];
+        let mut u_kind = vec![FaceKind::Fluid; g.n_ufaces()];
+        let mut v_kind = vec![FaceKind::Fluid; g.n_vfaces()];
+        let mut w_kind = vec![FaceKind::Fluid; g.n_wfaces()];
+        let kind_of = |phi: f64| {
+            if phi > -reach {
+                FaceKind::Ghost
+            } else {
+                FaceKind::Solid
+            }
+        };
+        for k in 0..nz {
+            for j in 0..ny {
+                for i in 1..nx {
+                    if !(is_fluid(k, j, i - 1) && is_fluid(k, j, i)) {
+                        u_kind[g.uface(k, j, i)] =
+                            kind_of(body.phi(i as f64 * dx, yc(j), zc(k), t));
+                    }
+                }
+            }
+            for j in 1..ny {
+                for i in 0..nx {
+                    if !(is_fluid(k, j - 1, i) && is_fluid(k, j, i)) {
+                        v_kind[g.vface(k, j, i)] =
+                            kind_of(body.phi(xc(i), j as f64 * dy, zc(k), t));
+                    }
+                }
+            }
+        }
+        let w_range = if periodic { 0..nz } else { 1..nz };
+        for k in w_range.clone() {
+            let below = if k > 0 { k - 1 } else { nz - 1 };
+            for j in 0..ny {
+                for i in 0..nx {
+                    if !(is_fluid(below, j, i) && is_fluid(k, j, i)) {
+                        w_kind[g.wface(k, j, i)] =
+                            kind_of(body.phi(xc(i), yc(j), k as f64 * dz, t));
+                    }
+                }
+            }
+        }
+        if periodic {
+            for j in 0..ny {
+                for i in 0..nx {
+                    w_kind[g.wface(nz, j, i)] = w_kind[g.wface(0, j, i)];
+                }
+            }
+        }
+        let face_pos = |c: usize, k: usize, j: usize, i: usize| match c {
+            0 => (i as f64 * dx, yc(j), zc(k)),
+            1 => (xc(i), j as f64 * dy, zc(k)),
+            _ => (xc(i), yc(j), k as f64 * dz),
+        };
+        let kji_of = |c: usize, idx: usize| -> (usize, usize, usize) {
+            match c {
+                0 => (idx / ((nx + 1) * ny), (idx / (nx + 1)) % ny, idx % (nx + 1)),
+                1 => (idx / (nx * (ny + 1)), (idx / nx) % (ny + 1), idx % nx),
+                _ => (idx / (nx * ny), (idx / nx) % ny, idx % nx),
+            }
+        };
+        let build_ghost = |c: usize,
+                           kinds: &[FaceKind],
+                           idx: usize,
+                           (x, y, z): (f64, f64, f64),
+                           flux_sign: f64|
+         -> Ghost {
+            let s_face = body.phi(x, y, z, t);
+            let (n1, n2, n3) = body.normal(x, y, z, t, eps);
+            let foot = (x - s_face * n1, y - s_face * n2, z - s_face * n3);
+            let s_probe = s_face.abs() + h_min;
+            let probe = (
+                foot.0 + s_probe * n1,
+                foot.1 + s_probe * n2,
+                foot.2 + s_probe * n3,
+            );
+            let vel = body.surface_velocity(foot.0, foot.1, foot.2, t);
+            let u_surface = [vel.0, vel.1, vel.2][c];
+            let nodes = stencil_nodes(probe, c, g, periodic, |nidx| {
+                if kinds[nidx] == FaceKind::Fluid {
+                    None
+                } else {
+                    let (kk, jj, ii) = kji_of(c, nidx);
+                    let (px, py, pz) = face_pos(c, kk, jj, ii);
+                    let s = body.phi(px, py, pz, t);
+                    let (m1, m2, m3) = body.normal(px, py, pz, t, eps);
+                    let f = body.surface_velocity(px - s * m1, py - s * m2, pz - s * m3, t);
+                    Some([f.0, f.1, f.2][c])
+                }
+            });
+            Ghost {
+                idx,
+                x,
+                y,
+                z,
+                foot,
+                u_surface,
+                s_face,
+                s_probe,
+                nodes,
+                flux_sign,
+            }
+        };
+        let mut u_ghosts = Vec::new();
+        let mut v_ghosts = Vec::new();
+        let mut w_ghosts = Vec::new();
+        for k in 0..nz {
+            for j in 0..ny {
+                for i in 1..nx {
+                    let idx = g.uface(k, j, i);
+                    if u_kind[idx] == FaceKind::Ghost {
+                        let sign = if is_fluid(k, j, i - 1) {
+                            1.0
+                        } else if is_fluid(k, j, i) {
+                            -1.0
+                        } else {
+                            0.0
+                        };
+                        u_ghosts.push(build_ghost(0, &u_kind, idx, face_pos(0, k, j, i), sign));
+                    }
+                }
+            }
+            for j in 1..ny {
+                for i in 0..nx {
+                    let idx = g.vface(k, j, i);
+                    if v_kind[idx] == FaceKind::Ghost {
+                        let sign = if is_fluid(k, j - 1, i) {
+                            1.0
+                        } else if is_fluid(k, j, i) {
+                            -1.0
+                        } else {
+                            0.0
+                        };
+                        v_ghosts.push(build_ghost(1, &v_kind, idx, face_pos(1, k, j, i), sign));
+                    }
+                }
+            }
+        }
+        for k in w_range {
+            let below = if k > 0 { k - 1 } else { nz - 1 };
+            for j in 0..ny {
+                for i in 0..nx {
+                    let idx = g.wface(k, j, i);
+                    if w_kind[idx] == FaceKind::Ghost {
+                        let sign = if is_fluid(below, j, i) {
+                            1.0
+                        } else if is_fluid(k, j, i) {
+                            -1.0
+                        } else {
+                            0.0
+                        };
+                        w_ghosts.push(build_ghost(2, &w_kind, idx, face_pos(2, k, j, i), sign));
+                    }
+                }
+            }
+        }
+        Ok(Self {
+            grid: g,
+            periodic_z: periodic,
+            cell_fluid,
+            u_kind,
+            v_kind,
+            w_kind,
+            u_ghosts,
+            v_ghosts,
+            w_ghosts,
+            anchor,
+            fluid_cells,
+        })
+    }
+
+    #[inline]
+    #[must_use]
+    pub fn is_fluid_cell(&self, idx: usize) -> bool {
+        self.cell_fluid[idx]
+    }
+    #[inline]
+    #[must_use]
+    pub fn u_kind(&self, idx: usize) -> FaceKind {
+        self.u_kind[idx]
+    }
+    #[inline]
+    #[must_use]
+    pub fn v_kind(&self, idx: usize) -> FaceKind {
+        self.v_kind[idx]
+    }
+    #[inline]
+    #[must_use]
+    pub fn w_kind(&self, idx: usize) -> FaceKind {
+        self.w_kind[idx]
+    }
+    #[must_use]
+    pub fn anchor(&self) -> usize {
+        self.anchor
+    }
+    #[must_use]
+    pub fn fluid_cells(&self) -> usize {
+        self.fluid_cells
+    }
+    #[must_use]
+    pub fn ghost_faces(&self) -> usize {
+        self.u_ghosts.len() + self.v_ghosts.len() + self.w_ghosts.len()
+    }
+    #[must_use]
+    pub fn grid(&self) -> Grid {
+        self.grid
+    }
+    #[must_use]
+    pub fn periodic_z(&self) -> bool {
+        self.periodic_z
+    }
+
+    /// Impose the wall on `(u, v, w)` from the same field.
+    pub fn impose(&self, body: &Body, u: &mut [f64], v: &mut [f64], w: &mut [f64], t: f64) -> f64 {
+        let (us, vs, ws) = (u.to_vec(), v.to_vec(), w.to_vec());
+        self.impose_from(body, &us, &vs, &ws, u, v, w, t)
+    }
+
+    /// Solid faces: the surface velocity; ghost faces: the reconstruction
+    /// from the SOURCE field, minus the shared flux compatibility
+    /// correction over the flux-carrying ghosts. Returns the correction.
+    #[allow(clippy::too_many_arguments)]
+    pub fn impose_from(
+        &self,
+        body: &Body,
+        u_src: &[f64],
+        v_src: &[f64],
+        w_src: &[f64],
+        u: &mut [f64],
+        v: &mut [f64],
+        w: &mut [f64],
+        t: f64,
+    ) -> f64 {
+        let g = self.grid;
+        let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz);
+        for k in 0..nz {
+            for j in 0..ny {
+                for i in 1..nx {
+                    let idx = g.uface(k, j, i);
+                    if self.u_kind[idx] == FaceKind::Solid {
+                        u[idx] = body
+                            .surface_velocity(
+                                i as f64 * dx,
+                                (j as f64 + 0.5) * dy,
+                                (k as f64 + 0.5) * dz,
+                                t,
+                            )
+                            .0;
+                    }
+                }
+            }
+            for j in 1..ny {
+                for i in 0..nx {
+                    let idx = g.vface(k, j, i);
+                    if self.v_kind[idx] == FaceKind::Solid {
+                        v[idx] = body
+                            .surface_velocity(
+                                (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 {
+                    let idx = g.wface(k, j, i);
+                    if self.w_kind[idx] == FaceKind::Solid {
+                        w[idx] = body
+                            .surface_velocity(
+                                (i as f64 + 0.5) * dx,
+                                (j as f64 + 0.5) * dy,
+                                k as f64 * dz,
+                                t,
+                            )
+                            .2;
+                    }
+                }
+            }
+        }
+        let u_vals: Vec<f64> = self
+            .u_ghosts
+            .iter()
+            .map(|gh| gh.reconstruct(u_src))
+            .collect();
+        let v_vals: Vec<f64> = self
+            .v_ghosts
+            .iter()
+            .map(|gh| gh.reconstruct(v_src))
+            .collect();
+        let w_vals: Vec<f64> = self
+            .w_ghosts
+            .iter()
+            .map(|gh| gh.reconstruct(w_src))
+            .collect();
+        let (au, av, aw) = (dy * dz, dx * dz, dx * dy);
+        let mut net = 0.0;
+        let mut area = 0.0;
+        for (gh, &val) in self.u_ghosts.iter().zip(&u_vals) {
+            if gh.flux_sign != 0.0 {
+                net += gh.flux_sign * val * au;
+                area += au;
+            }
+        }
+        for (gh, &val) in self.v_ghosts.iter().zip(&v_vals) {
+            if gh.flux_sign != 0.0 {
+                net += gh.flux_sign * val * av;
+                area += av;
+            }
+        }
+        for (gh, &val) in self.w_ghosts.iter().zip(&w_vals) {
+            if gh.flux_sign != 0.0 {
+                net += gh.flux_sign * val * aw;
+                area += aw;
+            }
+        }
+        let correction = if area > 0.0 { net / area } else { 0.0 };
+        for (gh, &val) in self.u_ghosts.iter().zip(&u_vals) {
+            u[gh.idx] = val - gh.flux_sign * correction;
+        }
+        for (gh, &val) in self.v_ghosts.iter().zip(&v_vals) {
+            v[gh.idx] = val - gh.flux_sign * correction;
+        }
+        for (gh, &val) in self.w_ghosts.iter().zip(&w_vals) {
+            w[gh.idx] = val - gh.flux_sign * correction;
+        }
+        // The periodic seam: the w face at k = nz is the face at k = 0.
+        for j in 0..ny {
+            for i in 0..nx {
+                let (f0, fn_) = (g.wface(0, j, i), g.wface(nz, j, i));
+                if self.w_kind[f0] != FaceKind::Fluid && self.w_kind[fn_] == self.w_kind[f0] {
+                    w[fn_] = w[f0];
+                }
+            }
+        }
+        correction
+    }
+}
diff --git a/crates/specialized/rtx-cfd/tests/embedded3_cylinder_identity.rs b/crates/specialized/rtx-cfd/tests/embedded3_cylinder_identity.rs
new file mode 100644
index 0000000..a2fc196
--- /dev/null
+++ b/crates/specialized/rtx-cfd/tests/embedded3_cylinder_identity.rs
@@ -0,0 +1,351 @@
+//! embedded3 gate 9b: Turek–Hron CFD1 (the cylinder with the rigid flag,
+//! Re 20) on the 3D solver at ny = 41 — the 2D geometry extruded, at nz = 1
+//! (dz = 1, z slip) and nz = 4 periodic: the settled control-volume drag
+//! 15.6156 and surface drag 15.7126 of the 2D embedded record to
+//! `rel < 5e-4` (printed-digit identity across the regimes).
+
+use rtx_cfd::solvers::incompressible::embedded3::{
+    Body, Boundaries, Field, Fluid, Grid, Parameters, Side, Solver,
+};
+use rtx_cfd::solvers::incompressible::{EmbeddedBody, MgSmoother};
+
+const L: f64 = 2.5;
+const H: f64 = 0.41;
+const RHO: f64 = 1000.0;
+const NU: f64 = 1e-3;
+const U_MEAN: f64 = 0.2;
+const SOR_DRAG_CV: f64 = 15.6156;
+const SOR_DRAG_SURFACE: f64 = 15.7126;
+
+fn inflow(y: f64) -> f64 {
+    1.5 * U_MEAN * y * (H - y) / (0.5 * H).powi(2)
+}
+
+fn body2() -> EmbeddedBody {
+    EmbeddedBody::union(
+        EmbeddedBody::circle(0.2, 0.2, 0.05),
+        EmbeddedBody::rectangle(0.20, 0.19, 0.6, 0.21),
+    )
+}
+
+fn run(ny: usize, nz: usize, dz: f64, periodic: bool) -> (f64, f64, usize, usize) {
+    let h = H / ny as f64;
+    let nx = (L / h).round() as usize;
+    let mu = RHO * NU;
+    let u_peak = 1.5 * 1.5 * U_MEAN;
+    let dt = 0.25 / (2.0 * u_peak / h + 4.0 * NU / (h * h));
+    let z = if periodic {
+        Side::Periodic
+    } else {
+        Side::SlipWall
+    };
+    let mut solver = Solver::new(
+        Fluid {
+            density: RHO,
+            viscosity: mu,
+            reference_velocity: U_MEAN,
+            reference_length: 0.1,
+        },
+        Parameters {
+            corrector_steps: 2,
+            tolerance: 1e-7,
+            boundaries: Boundaries {
+                x1: Side::PressureOutlet,
+                z0: z,
+                z1: z,
+                ..Boundaries::default()
+            },
+            poisson_smoother: MgSmoother::Lexicographic,
+            ..Parameters::default()
+        },
+    );
+    solver.set_boundary_velocity(|x, y, _z, _t| {
+        if x <= 0.0 {
+            (inflow(y), 0.0, 0.0)
+        } else {
+            (0.0, 0.0, 0.0)
+        }
+    });
+    let lz = nz as f64 * dz;
+    solver.set_body(Body::extruded(body2(), lz));
+    let g = Grid {
+        nx,
+        ny,
+        nz,
+        dx: h,
+        dy: h,
+        dz,
+    };
+    let mut f = Field::new(g);
+    for k in 0..nz {
+        for j in 0..ny {
+            let u0 = inflow((j as f64 + 0.5) * h);
+            for i in 0..=nx {
+                f.u[g.uface(k, j, i)] = u0;
+            }
+        }
+    }
+    solver.initialize(&mut f);
+    let cv = (
+        (0.10 / h).round() as usize,
+        (0.75 / h).round() as usize,
+        (0.05 / h).round() as usize,
+        (0.36 / h).round() as usize,
+        0,
+        nz,
+    );
+    let flow_through = L / U_MEAN;
+    let min_steps = (flow_through / dt).ceil() as usize;
+    let mut history: Vec<f64> = Vec::new();
+    let mut steps = 0;
+    loop {
+        solver.advance(&mut f, dt);
+        steps += 1;
+        if steps % 50 == 0 {
+            let fx = solver
+                .mask()
+                .unwrap()
+                .control_volume_force(&f, dt, RHO, mu, None, cv)[0]
+                / lz;
+            history.push(fx);
+            let umax = f.u.iter().fold(0.0_f64, |m, v| m.max(v.abs()));
+            assert!(umax.is_finite(), "non-finite at step {steps}");
+            if steps >= min_steps && history.len() > 4 {
+                let now = history[history.len() - 1];
+                let then = history[history.len() - 5];
+                if ((now - then) / now).abs() < 1e-4 {
+                    break;
+                }
+            }
+        }
+        assert!(steps < 400_000, "did not settle");
+    }
+    let mask = solver.mask().unwrap();
+    let surface = mask.surface_force(solver.body().unwrap(), &f, mu, solver.time(), 0.5 * h);
+    let drag_cv = mask.control_volume_force(&f, dt, RHO, mu, None, cv)[0] / lz;
+    (drag_cv, surface.f[0] / lz, surface.skipped, steps)
+}
+
+#[test]
+fn cfd1_at_ny_41_reproduces_the_two_d_record() {
+    let ny = 41;
+    let h = H / ny as f64;
+    for (nz, dz, periodic) in [(1usize, 1.0, false), (4, h, true)] {
+        let (cv, surface, skipped, steps) = run(ny, nz, dz, periodic);
+        let rel_cv = ((cv - SOR_DRAG_CV) / SOR_DRAG_CV).abs();
+        let rel_s = ((surface - SOR_DRAG_SURFACE) / SOR_DRAG_SURFACE).abs();
+        println!(
+            "  ny 41 nz {nz} periodic {periodic}: {steps} steps; CV drag {cv:.4} (record 15.6156, rel {rel_cv:.2e}); surface drag {surface:.4} (record 15.7126, rel {rel_s:.2e}, skipped {skipped})"
+        );
+        assert!(rel_cv < 5e-4, "CV drag {cv:.4} vs the record 15.6156");
+        assert!(
+            rel_s < 5e-4,
+            "surface drag {surface:.4} vs the record 15.7126"
+        );
+    }
+}
+
+/// Diagnostic: which probes fail on the skipped surface samples, and the
+/// surface force per z level, at nz 4 periodic after 200 steps.
+#[test]
+#[ignore = "diagnostic: skipped surface samples and per-level force on CFD1 at nz 4"]
+fn skipped_samples_diagnostic() {
+    let ny = 41;
+    let h = H / ny as f64;
+    let nx = (L / h).round() as usize;
+    let mu = RHO * NU;
+    let u_peak = 1.5 * 1.5 * U_MEAN;
+    let dt = 0.25 / (2.0 * u_peak / h + 4.0 * NU / (h * h));
+    let (nz, dz) = (4usize, h);
+    let lz = nz as f64 * dz;
+    let mut solver = Solver::new(
+        Fluid {
+            density: RHO,
+            viscosity: mu,
+            reference_velocity: U_MEAN,
+            reference_length: 0.1,
+        },
+        Parameters {
+            corrector_steps: 2,
+            tolerance: 1e-7,
+            boundaries: Boundaries {
+                x1: Side::PressureOutlet,
+                z0: Side::Periodic,
+                z1: Side::Periodic,
+                ..Boundaries::default()
+            },
+            poisson_smoother: MgSmoother::Lexicographic,
+            ..Parameters::default()
+        },
+    );
+    solver.set_boundary_velocity(|x, y, _z, _t| {
+        if x <= 0.0 {
+            (inflow(y), 0.0, 0.0)
+        } else {
+            (0.0, 0.0, 0.0)
+        }
+    });
+    solver.set_body(Body::extruded(body2(), lz));
+    let g = Grid {
+        nx,
+        ny,
+        nz,
+        dx: h,
+        dy: h,
+        dz,
+    };
+    let mut f = Field::new(g);
+    for k in 0..nz {
+        for j in 0..ny {
+            let u0 = inflow((j as f64 + 0.5) * h);
+            for i in 0..=nx {
+                f.u[g.uface(k, j, i)] = u0;
+            }
+        }
+    }
+    solver.initialize(&mut f);
+    for _ in 0..200 {
+        solver.advance(&mut f, dt);
+    }
+    let mask = solver.mask().unwrap();
+    let body = solver.body().unwrap();
+    let samples = body.surface_samples(0.5 * h);
+    let mut by_z: std::collections::BTreeMap<i64, (usize, usize, f64)> =
+        std::collections::BTreeMap::new();
+    let mut shown = 0;
+    for s in &samples {
+        let n = [s.nx, s.ny, s.nz];
+        let key = (s.z * 1e4).round() as i64;
+        let e = by_z.entry(key).or_insert((0, 0, 0.0));
+        e.0 += 1;
+        match mask.traction_at(body, &f, mu, solver.time(), [s.x, s.y, s.z], n) {
+            Some(tr) => e.2 += tr[0] * s.area,
+            None => {
+                e.1 += 1;
+                if shown < 6 {
+                    shown += 1;
+                    let at = |d: f64| [s.x + d * n[0], s.y + d * n[1], s.z + d * n[2]];
+                    let (x1, x2) = (at(h), at(2.0 * h));
+                    println!(
+                        "  skipped ({:.4}, {:.4}, {:.4}) n ({:.2}, {:.2}): p1 {} p2 {} u1 {} u2 {}",
+                        s.x,
+                        s.y,
+                        s.z,
+                        s.nx,
+                        s.ny,
+                        mask.pressure_at(&f.p, x1[0], x1[1], x1[2]).is_some(),
+                        mask.pressure_at(&f.p, x2[0], x2[1], x2[2]).is_some(),
+                        mask.velocity_at(body, &f, x1[0], x1[1], x1[2], 0.0)
+                            .is_some(),
+                        mask.velocity_at(body, &f, x2[0], x2[1], x2[2], 0.0)
+                            .is_some()
+                    );
+                }
+            }
+        }
+    }
+    for (z, (n, sk, fx)) in &by_z {
+        println!(
+            "  z {:.4}: {n} samples, {sk} skipped, drag contribution per unit depth {:.4}",
+            *z as f64 / 1e4,
+            fx / (lz / by_z.len() as f64)
+        );
+    }
+}
+
+/// Diagnostic: is the periodic nz 4 solution z-invariant, and does its
+/// plane 0 equal the nz 1 solution, after 200 steps from the same start?
+#[test]
+#[ignore = "diagnostic: plane symmetry of CFD1 at nz 4 periodic"]
+fn plane_symmetry_diagnostic() {
+    let ny = 41;
+    let h = H / ny as f64;
+    let nx = (L / h).round() as usize;
+    let mu = RHO * NU;
+    let u_peak = 1.5 * 1.5 * U_MEAN;
+    let dt = 0.25 / (2.0 * u_peak / h + 4.0 * NU / (h * h));
+    let mk = |nz: usize, dz: f64, z: Side| {
+        let mut s = Solver::new(
+            Fluid {
+                density: RHO,
+                viscosity: mu,
+                reference_velocity: U_MEAN,
+                reference_length: 0.1,
+            },
+            Parameters {
+                corrector_steps: 2,
+                tolerance: 1e-7,
+                boundaries: Boundaries {
+                    x1: Side::PressureOutlet,
+                    z0: z,
+                    z1: z,
+                    ..Boundaries::default()
+                },
+                poisson_smoother: MgSmoother::Lexicographic,
+                ..Parameters::default()
+            },
+        );
+        s.set_boundary_velocity(|x, y, _z, _t| {
+            if x <= 0.0 {
+                (inflow(y), 0.0, 0.0)
+            } else {
+                (0.0, 0.0, 0.0)
+            }
+        });
+        s.set_body(Body::extruded(body2(), nz as f64 * dz));
+        let g = Grid {
+            nx,
+            ny,
+            nz,
+            dx: h,
+            dy: h,
+            dz,
+        };
+        let mut f = Field::new(g);
+        for k in 0..nz {
+            for j in 0..ny {
+                let u0 = inflow((j as f64 + 0.5) * h);
+                for i in 0..=nx {
+                    f.u[g.uface(k, j, i)] = u0;
+                }
+            }
+        }
+        s.initialize(&mut f);
+        (s, f, g)
+    };
+    let (mut s1, mut f1, g1) = mk(1, 1.0, Side::SlipWall);
+    let (mut s4, mut f4, g4) = mk(4, h, Side::Periodic);
+    println!(
+        "  ghost faces: nz 1 {} / nz 4 {} (per plane {})",
+        s1.mask().unwrap().ghost_faces(),
+        s4.mask().unwrap().ghost_faces(),
+        s4.mask().unwrap().ghost_faces() / 4
+    );
+    for step in 1..=200 {
+        s1.advance(&mut f1, dt);
+        s4.advance(&mut f4, dt);
+        if [1, 2, 10, 50, 200].contains(&step) {
+            let plane = |f: &Field, g: &Grid, k: usize| {
+                f.u[k * g.ny * (g.nx + 1)..(k + 1) * g.ny * (g.nx + 1)].to_vec()
+            };
+            let p0 = plane(&f4, &g4, 0);
+            let mut zinv = 0.0_f64;
+            for k in 1..4 {
+                for (a, b) in plane(&f4, &g4, k).iter().zip(&p0) {
+                    zinv = zinv.max((a - b).abs());
+                }
+            }
+            let p1 = plane(&f1, &g1, 0);
+            let vs1 = p0
+                .iter()
+                .zip(&p1)
+                .fold(0.0_f64, |m, (a, b)| m.max((a - b).abs()));
+            let wmax = f4.w.iter().fold(0.0_f64, |m, v| m.max(v.abs()));
+            println!(
+                "  step {step}: nz 4 planes within {zinv:.3e}; plane 0 vs nz 1 {vs1:.3e}; max |w| {wmax:.3e}; ghost corr nz1 {:.3e} / nz4 {:.3e}",
+                s1.ghost_correction(),
+                s4.ghost_correction()
+            );
+        }
+    }
+}
diff --git a/crates/specialized/rtx-cfd/tests/embedded3_embedded_mms.rs b/crates/specialized/rtx-cfd/tests/embedded3_embedded_mms.rs
new file mode 100644
index 0000000..6582a9e
--- /dev/null
+++ b/crates/specialized/rtx-cfd/tests/embedded3_embedded_mms.rs
@@ -0,0 +1,301 @@
+//! embedded3 gate 9a: the manufactured solution with an embedded sphere
+//! (centre (0.6, 0.45, 0.5), r 0.2, off-centre so the exact force is not
+//! zero by symmetry) carrying the exact field as its surface velocity, on
+//! the binary ghost wall. The velocity error falls at the scheme's order,
+//! every fluid cell is divergence-free, the compatibility correction
+//! shrinks, and both load routes converge to the exact surface integral of
+//! the manufactured stress (the control-volume route measures F − M with M
+//! the momentum flux through the porous manufactured surface).
+
+use rtx_cfd::solvers::incompressible::ConvectionScheme;
+use rtx_cfd::solvers::incompressible::embedded3::{Body, Field, Fluid, Grid, Parameters, Solver};
+use std::f64::consts::PI;
+
+const RHO: f64 = 1.0;
+const MU: f64 = 0.05;
+const C: (f64, f64, f64) = (0.6, 0.45, 0.5);
+const R: f64 = 0.2;
+
+fn u3(x: f64, y: f64, z: f64) -> f64 {
+    (PI * x).sin() * (PI * y).cos() * (PI * z).cos()
+}
+fn v3(x: f64, y: f64, z: f64) -> f64 {
+    (PI * x).cos() * (PI * y).sin() * (PI * z).cos()
+}
+fn w3(x: f64, y: f64, z: f64) -> f64 {
+    -2.0 * (PI * x).cos() * (PI * y).cos() * (PI * z).sin()
+}
+fn p3(x: f64, y: f64, z: f64) -> f64 {
+    (PI * x).sin() * (PI * y).sin() * (PI * z).sin()
+}
+/// The velocity gradient ∂u_i/∂x_j and the pressure gradient.
+fn grads(x: f64, y: f64, z: f64) -> ([[f64; 3]; 3], [f64; 3]) {
+    let (sx, cx) = (PI * x).sin_cos();
+    let (sy, cy) = (PI * y).sin_cos();
+    let (sz, cz) = (PI * z).sin_cos();
+    (
+        [
+            [PI * cx * cy * cz, -PI * sx * sy * cz, -PI * sx * cy * sz],
+            [-PI * sx * sy * cz, PI * cx * cy * cz, -PI * cx * sy * sz],
+            [
+                2.0 * PI * sx * cy * sz,
+                2.0 * PI * cx * sy * sz,
+                -2.0 * PI * cx * cy * cz,
+            ],
+        ],
+        [PI * cx * sy * sz, PI * sx * cy * sz, PI * sx * sy * cz],
+    )
+}
+fn source3(x: f64, y: f64, z: f64) -> (f64, f64, f64) {
+    let (g, gp) = grads(x, y, z);
+    let u = [u3(x, y, z), v3(x, y, z), w3(x, y, z)];
+    let lap = -3.0 * PI * PI;
+    let conv = |i: usize| u[0] * g[i][0] + u[1] * g[i][1] + u[2] * g[i][2];
+    (
+        RHO * conv(0) + gp[0] - MU * lap * u[0],
+        RHO * conv(1) + gp[1] - MU * lap * u[1],
+        RHO * conv(2) + gp[2] - MU * lap * u[2],
+    )
+}
+fn boundary3(x: f64, y: f64, z: f64) -> (f64, f64, f64) {
+    let u = if x <= 0.0 || x >= 1.0 {
+        0.0
+    } else {
+        u3(x, y, z)
+    };
+    let v = if y <= 0.0 || y >= 1.0 {
+        0.0
+    } else {
+        v3(x, y, z)
+    };
+    let w = if z <= 0.0 || z >= 1.0 {
+        0.0
+    } else {
+        w3(x, y, z)
+    };
+    (u, v, w)
+}
+
+/// Exact force `∮ (−p I + μ(∇u + ∇uᵀ)) n dA` and momentum flux `∮ ρ u (u·n) dA`
+/// over the sphere by a fine Fibonacci quadrature.
+fn exact_force_and_flux() -> ([f64; 3], [f64; 3]) {
+    let n = 200_000;
+    let golden = PI * (3.0 - 5.0_f64.sqrt());
+    let (mut f, mut m) = ([0.0; 3], [0.0; 3]);
+    let da = 4.0 * PI * R * R / n as f64;
+    for k in 0..n {
+        let zz = 1.0 - 2.0 * (k as f64 + 0.5) / n as f64;
+        let rr = (1.0 - zz * zz).sqrt();
+        let th = golden * k as f64;
+        let nrm = [rr * th.cos(), rr * th.sin(), zz];
+        let (x, y, z) = (C.0 + R * nrm[0], C.1 + R * nrm[1], C.2 + R * nrm[2]);
+        let (g, _) = grads(x, y, z);
+        let p = p3(x, y, z);
+        let u = [u3(x, y, z), v3(x, y, z), w3(x, y, z)];
+        let un = u[0] * nrm[0] + u[1] * nrm[1] + u[2] * nrm[2];
+        for i in 0..3 {
+            let mut t = -p * nrm[i];
+            for j in 0..3 {
+                t += MU * (g[i][j] + g[j][i]) * nrm[j];
+            }
+            f[i] += t * da;
+            m[i] += RHO * u[i] * un * da;
+        }
+    }
+    (f, m)
+}
+
+struct Measurement {
+    l2_velocity: f64,
+    max_div: f64,
+    ghost_correction: f64,
+    force_surface: [f64; 3],
+    skipped: usize,
+    force_cv: [f64; 3],
+}
+
+fn measure(n: usize) -> Measurement {
+    let h = 1.0 / n as f64;
+    let dt = 0.4 * (h * h / (4.0 * MU / RHO)).min(h);
+    let mut solver = Solver::new(
+        Fluid {
+            density: RHO,
+            viscosity: MU,
+            reference_velocity: 1.0,
+            reference_length: 1.0,
+        },
+        Parameters {
+            corrector_steps: 2,
+            tolerance: 1e-8,
+            convection_scheme: ConvectionScheme::Upwind,
+            ..Parameters::default()
+        },
+    );
+    solver.set_momentum_source(|x, y, z, _t| source3(x, y, z));
+    solver.set_boundary_velocity(|x, y, z, _t| boundary3(x, y, z));
+    solver.set_body(
+        Body::sphere(|_t| C, R)
+            .with_surface_velocity(|x, y, z, _t| (u3(x, y, z), v3(x, y, z), w3(x, y, z))),
+    );
+    let g = Grid::cubic(n, n, n, h);
+    let mut f = Field::new(g);
+    solver.initialize(&mut f);
+    for _ in 0..200_000 {
+        let (bu, bv, bw) = (f.u.clone(), f.v.clone(), f.w.clone());
+        solver.advance(&mut f, dt);
+        let mut change = 0.0_f64;
+        for (a, b) in
+            f.u.iter()
+                .zip(&bu)
+                .chain(f.v.iter().zip(&bv))
+                .chain(f.w.iter().zip(&bw))
+        {
+            change = change.max((a - b).abs());
+        }
+        if change / dt < 1e-6 {
+            break;
+        }
+    }
+    let mask = solver.mask().expect("mask");
+    use rtx_cfd::solvers::incompressible::embedded3::FaceKind;
+    let (mut sq, mut vol) = (0.0, 0.0);
+    let dv = h * h * h;
+    for k in 0..n {
+        for j in 0..n {
+            for i in 1..n {
+                if mask.u_kind(g.uface(k, j, i)) == FaceKind::Fluid {
+                    let e = f.u[g.uface(k, j, i)]
+                        - u3(i as f64 * h, (j as f64 + 0.5) * h, (k as f64 + 0.5) * h);
+                    sq += e * e * dv;
+                    vol += dv;
+                }
+            }
+        }
+        for j in 1..n {
+            for i in 0..n {
+                if mask.v_kind(g.vface(k, j, i)) == FaceKind::Fluid {
+                    let e = f.v[g.vface(k, j, i)]
+                        - v3((i as f64 + 0.5) * h, j as f64 * h, (k as f64 + 0.5) * h);
+                    sq += e * e * dv;
+                    vol += dv;
+                }
+            }
+        }
+    }
+    for k in 1..n {
+        for j in 0..n {
+            for i in 0..n {
+                if mask.w_kind(g.wface(k, j, i)) == FaceKind::Fluid {
+                    let e = f.w[g.wface(k, j, i)]
+                        - w3((i as f64 + 0.5) * h, (j as f64 + 0.5) * h, k as f64 * h);
+                    sq += e * e * dv;
+                    vol += dv;
+                }
+            }
+        }
+    }
+    let mut max_div = 0.0_f64;
+    for k in 0..n {
+        for j in 0..n {
+            for i in 0..n {
+                if mask.is_fluid_cell(g.cell(k, j, i)) {
+                    let div = (f.u[g.uface(k, j, i + 1)] - f.u[g.uface(k, j, i)]) / h
+                        + (f.v[g.vface(k, j + 1, i)] - f.v[g.vface(k, j, i)]) / h
+                        + (f.w[g.wface(k + 1, j, i)] - f.w[g.wface(k, j, i)]) / h;
+                    max_div = max_div.max(div.abs());
+                }
+            }
+        }
+    }
+    let body = solver.body().expect("body");
+    let surface = mask.surface_force(body, &f, MU, solver.time(), 0.5 * h);
+    let (i0, i1) = (n / 8, n - n / 8);
+    let src = |x: f64, y: f64, z: f64| source3(x, y, z);
+    let force_cv = mask.control_volume_force(&f, dt, RHO, MU, Some(&src), (i0, i1, i0, i1, i0, i1));
+    Measurement {
+        l2_velocity: (sq / vol).sqrt(),
+        max_div,
+        ghost_correction: solver.ghost_correction().abs(),
+        force_surface: surface.f,
+        skipped: surface.skipped,
+        force_cv,
+    }
+}
+
+fn norm(a: [f64; 3]) -> f64 {
+    (a[0] * a[0] + a[1] * a[1] + a[2] * a[2]).sqrt()
+}
+
+fn ladder(resolutions: &[usize]) {
+    let (fe, m) = exact_force_and_flux();
+    let f_scale = norm(fe);
+    let fcv = [fe[0] - m[0], fe[1] - m[1], fe[2] - m[2]];
+    println!(
+        "  exact force {fe:.5?}; momentum flux {m:.5?}; the control-volume route measures {fcv:.5?}"
+    );
+    let ms: Vec<Measurement> = resolutions.iter().map(|&n| measure(n)).collect();
+    let errors: Vec<f64> = ms.iter().map(|x| x.l2_velocity).collect();
+    let mut se = Vec::new();
+    let mut ce = Vec::new();
+    for (k, (mm, &n)) in ms.iter().zip(resolutions).enumerate() {
+        let rate = if k == 0 {
+            "    -".to_string()
+        } else {
+            format!("{:5.2}", (errors[k - 1] / errors[k]).log2())
+        };
+        let s = norm([
+            mm.force_surface[0] - fe[0],
+            mm.force_surface[1] - fe[1],
+            mm.force_surface[2] - fe[2],
+        ]) / f_scale;
+        let c = norm([
+            mm.force_cv[0] - fcv[0],
+            mm.force_cv[1] - fcv[1],
+            mm.force_cv[2] - fcv[2],
+        ]) / f_scale;
+        println!(
+            "  n = {n:3}  L2 u {:.4e} (order {rate})  max div {:.2e}  ghost corr {:.2e}  F_surface {:.4?} rel {s:.3e} (skipped {})  F_cv {:.4?} rel {c:.3e}",
+            mm.l2_velocity,
+            mm.max_div,
+            mm.ghost_correction,
+            mm.force_surface,
+            mm.skipped,
+            mm.force_cv
+        );
+        se.push(s);
+        ce.push(c);
+    }
+    assert!(
+        errors.windows(2).all(|w| w[1] < w[0]),
+        "errors not monotone {errors:?}"
+    );
+    for w in errors.windows(2) {
+        let rate = (w[0] / w[1]).log2();
+        assert!(
+            rate > 0.75 && rate < 2.3,
+            "order {rate:.3} outside [0.75, 2.3]"
+        );
+    }
+    for mm in &ms {
+        assert!(mm.max_div < 1e-5, "max div {:.3e}", mm.max_div);
+    }
+    assert!(
+        se.windows(2).all(|w| w[1] < w[0]),
+        "surface-route error not falling {se:?}"
+    );
+    assert!(
+        ce.windows(2).all(|w| w[1] < w[0]),
+        "control-volume-route error not falling {ce:?}"
+    );
+}
+
+#[test]
+fn embedded_sphere_recovers_the_manufactured_solution() {
+    ladder(&[12, 24]);
+}
+
+#[test]
+#[ignore = "the three-rung ladder to n = 48 (minutes on the host)"]
+fn embedded_sphere_three_rungs() {
+    ladder(&[12, 24, 48]);
+}