diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/composite/assemble.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/composite/assemble.rs
index 5ae8f7d..ca3d8fd 100644
--- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/composite/assemble.rs
+++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/composite/assemble.rs
@@ -57,6 +57,41 @@ pub struct Composite {
     pub rhs: Vec<f64>,
     /// The uniform coarse operator (coarse apertures under the patch too).
     pub coarse_problem: Problem,
+    /// R7-2a: the body's cut on the fine level (apertures, volumes, wall
+    /// vectors), `None` without a body.
+    pub fine_cut: Option<CutGeometry>,
+    /// R7-2a: every fine face on the patch boundary with the fine row's
+    /// interface terms (the MAC projection's gradient on that face).
+    pub iface: Vec<InterfaceFace>,
+}
+
+/// R7-2a: one fine face on the patch boundary. The fine row of `row` holds
+/// `Σ terms·p` for this face, i.e. the face's flux OUT of the fine cell is
+/// `−Σ terms·p` (and the coarse row holds its negative: conservative).
+#[derive(Debug, Clone)]
+pub struct InterfaceFace {
+    /// The fine unknown owning the face.
+    pub row: usize,
+    /// The face's axis (0 x, 1 y, 2 z) and the side of the fine cell it is
+    /// on (−1 / +1).
+    pub axis: usize,
+    pub side: isize,
+    /// The face's index in the fine grid's `axis` face array.
+    pub face: usize,
+    /// The coarse face it is part of (index in the coarse `axis` array).
+    pub coarse_face: usize,
+    pub terms: Vec<(usize, f64)>,
+}
+
+/// The index of face `(axis, c)` — the face on the minus side of cell `c`
+/// (`c[axis]` may equal the grid's extent) — in `g`'s `axis` face array.
+#[must_use]
+pub fn face_index(g: Grid, axis: usize, c: [usize; 3]) -> usize {
+    match axis {
+        0 => g.uface(c[2], c[1], c[0]),
+        1 => g.vface(c[2], c[1], c[0]),
+        _ => g.wface(c[2], c[1], c[0]),
+    }
 }
 
 const NONE: usize = usize::MAX;
@@ -184,6 +219,7 @@ impl Composite {
         let mut rows: Vec<Vec<(usize, f64)>> = vec![Vec::new(); n];
         let mut rhs = vec![0.0; n];
         let mut at_interface = vec![false; n];
+        let mut iface = Vec::new();
 
         let step = |c: [usize; 3], d: usize, s: isize| -> [isize; 3] {
             let mut o = [c[0] as isize, c[1] as isize, c[2] as isize];
@@ -371,27 +407,41 @@ impl Composite {
                             gc[2].div_euclid(2),
                         ];
                         let c_row = cid(cc);
+                        let mut terms: Vec<(usize, f64)> = Vec::with_capacity(12);
                         match spec.interface {
                             Interface::Direct => {
                                 let coef = hf / 1.5;
-                                rows[row].push((row, coef));
-                                rows[row].push((c_row, -coef));
+                                terms.push((row, coef));
+                                terms.push((c_row, -coef));
                             }
                             Interface::Octree => {
                                 let coef = hf * hf / hc;
                                 let kc = [lo[0] + i / 2, lo[1] + j / 2, lo[2] + k / 2];
                                 for sib in children(kc) {
-                                    rows[row].push((fid(sib), coef / 8.0));
+                                    terms.push((fid(sib), coef / 8.0));
                                 }
-                                rows[row].push((c_row, -coef));
+                                terms.push((c_row, -coef));
                             }
                             Interface::Quadratic => {
-                                rows[row].push((row, hf));
+                                terms.push((row, hf));
                                 for (u, w) in ghost(f, d, s) {
-                                    rows[row].push((u, -hf * w));
+                                    terms.push((u, -hf * w));
                                 }
                             }
                         }
+                        rows[row].extend_from_slice(&terms);
+                        let mut ff = f;
+                        ff[d] += usize::from(s > 0);
+                        let mut cf = [cc[0] as usize, cc[1] as usize, cc[2] as usize];
+                        cf[d] = if s < 0 { lo[d] } else { hi[d] };
+                        iface.push(InterfaceFace {
+                            row,
+                            axis: d,
+                            side: s,
+                            face: face_index(fg, d, ff),
+                            coarse_face: face_index(cg, d, cf),
+                            terms,
+                        });
                     }
                 }
             }
@@ -415,6 +465,8 @@ impl Composite {
             a,
             rhs,
             coarse_problem,
+            fine_cut,
+            iface,
         }
     }
 
diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/composite/mac.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/composite/mac.rs
new file mode 100644
index 0000000..2d22a22
--- /dev/null
+++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/composite/mac.rs
@@ -0,0 +1,472 @@
+//! R7-2a: the composite MAC projection on the coarse grid plus the nested
+//! ratio-2 patch (host prototype, nothing in the step calls it).
+//!
+//! Velocity layout: the normal velocity on every coarse face and on every
+//! fine face of the patch (the embedded3 staggered layout on each level).
+//! Ownership:
+//!
+//! - a coarse face between two uncovered coarse cells, or between one and
+//!   the domain boundary, is owned by the coarse level;
+//! - every face of the fine level is owned by the fine level, including the
+//!   fine faces on the patch boundary (the coarse–fine interface);
+//! - a coarse face on the interface, or under the patch, is SLAVED: its
+//!   flux is the sum of the fine fluxes it covers (average-down), so each
+//!   coarse cell next to the patch sees exactly the fluxes the fine cells
+//!   see (conservative across the interface).
+//!
+//! Operators, all in the integrated (flux) form of the composite pressure
+//! operator `A` of [`Composite`] (`A p` = `Σ c_f (p_i − p_nb)` per row):
+//!
+//! - `D u`: per unknown, the outward flux `Σ ± a_f |f| u_f` (fine apertures
+//!   `a_f` of the body's cut; a face with no fluid or with an inactive cell
+//!   on either side carries no flux: a Neumann wall), plus an optional wall
+//!   flux per unknown (a moving cut body);
+//! - `G p`: the face-normal gradient — `(p₊ − p₋)/h` on a coarse or fine
+//!   face, `± (p_b − p)/(h/2)` on the Dirichlet domain boundary, and on a
+//!   fine interface face the row's own interface terms (the Quadratic
+//!   ghost: `(g − p_f)/h_f`), so the coarse interface face's gradient is the
+//!   average of its four fine ones.
+//!
+//! By construction `D G p = b − A p` exactly (to round-off), with `b` the
+//! Dirichlet part of the operator's right-hand side; the projection
+//! `A p = b − D u*`, `u = u* − G p` therefore leaves `D u` equal to minus
+//! the linear solve's residual in every cell, the interface cells included.
+
+use super::assemble::face_index;
+use super::solve::{CompositeSolve, solve_bicgstab_with};
+use super::{Composite, CompositeSpec};
+use crate::solvers::incompressible::embedded3::Grid;
+use rayon::prelude::*;
+
+const NONE: usize = usize::MAX;
+
+/// What a face is on the composite grid.
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub enum FaceClass {
+    /// Coarse face between two uncovered coarse cells.
+    Coarse,
+    /// Coarse face on the (Dirichlet) domain boundary.
+    Boundary,
+    /// Coarse face on the patch boundary (slaved to its 4 fine faces).
+    CoarseInterface,
+    /// Coarse face under the patch (slaved; not part of any row).
+    Covered,
+    /// Fine face between two active fine cells (possibly cut).
+    Fine,
+    /// Fine face on the patch boundary (owns the interface flux).
+    FineInterface,
+    /// Fine face with no fluid or an inactive cell on a side (a wall).
+    Closed,
+}
+
+/// The staggered velocity of both levels (normal component per face, in
+/// each grid's u / v / w face layout).
+#[derive(Debug, Clone)]
+pub struct MacField {
+    pub coarse: [Vec<f64>; 3],
+    pub fine: [Vec<f64>; 3],
+}
+
+/// Outcome of one [`MacProjection::project`].
+#[derive(Debug, Clone, Copy)]
+pub struct ProjectionStats {
+    pub iterations: usize,
+    pub rel_residual: f64,
+    pub converged: bool,
+    /// max |D u| over the unknowns before and after.
+    pub div_before: f64,
+    pub div_after: f64,
+    pub solve_s: f64,
+    /// Divergence, gradient and average-down time.
+    pub mac_s: f64,
+}
+
+fn n_faces(g: Grid, axis: usize) -> usize {
+    (g.nx + usize::from(axis == 0))
+        * (g.ny + usize::from(axis == 1))
+        * (g.nz + usize::from(axis == 2))
+}
+
+/// `(i, j, k)` of face `idx` of `g`'s `axis` array.
+fn face_coords(g: Grid, axis: usize, idx: usize) -> [usize; 3] {
+    let nx = g.nx + usize::from(axis == 0);
+    let ny = g.ny + usize::from(axis == 1);
+    [idx % nx, (idx / nx) % ny, idx / (nx * ny)]
+}
+
+/// The composite projection: the operator, its preconditioner and the face
+/// tables.
+pub struct MacProjection {
+    pub comp: Composite,
+    pub coarse_class: [Vec<FaceClass>; 3],
+    pub fine_class: [Vec<FaceClass>; 3],
+    /// Fluid area `a_f h_f²` of every fine face (0 on a closed face).
+    pub fine_area: [Vec<f64>; 3],
+    /// The Dirichlet pressure at every boundary face (0 elsewhere).
+    pb: [Vec<f64>; 3],
+    /// The Dirichlet part of the right-hand side (the spec's source, if
+    /// any, included).
+    base_rhs: Vec<f64>,
+    /// Per unknown: its cell (coarse cells for `u < n_coarse`, else fine).
+    cell_of: Vec<usize>,
+    pre: CompositeSolve,
+}
+
+impl MacProjection {
+    /// Build the composite operator of `spec` (its `source` should be zero:
+    /// it would act as an added mass source) and the face tables.
+    #[must_use]
+    pub fn new(spec: &CompositeSpec<'_>, sweeps: usize) -> Self {
+        let comp = Composite::build(spec);
+        let (cg, fg, lo, hi) = (comp.coarse, comp.fine, comp.lo, comp.hi);
+        let hc = cg.dx;
+        let cdims = [cg.nx, cg.ny, cg.nz];
+        let fdims = [fg.nx, fg.ny, fg.nz];
+        let covered = |c: [usize; 3]| (0..3).all(|d| c[d] >= lo[d] && c[d] < hi[d]);
+        let mut coarse_class: [Vec<FaceClass>; 3] = Default::default();
+        let mut pb: [Vec<f64>; 3] = Default::default();
+        for d in 0..3 {
+            let nf = n_faces(cg, d);
+            coarse_class[d] = vec![FaceClass::Coarse; nf];
+            pb[d] = vec![0.0; nf];
+            for idx in 0..nf {
+                let c = face_coords(cg, d, idx);
+                debug_assert_eq!(face_index(cg, d, c), idx);
+                let class = if c[d] == 0 || c[d] == cdims[d] {
+                    let mut x = [
+                        (c[0] as f64 + 0.5) * hc,
+                        (c[1] as f64 + 0.5) * hc,
+                        (c[2] as f64 + 0.5) * hc,
+                    ];
+                    x[d] = c[d] as f64 * hc;
+                    pb[d][idx] = (spec.dirichlet)(x[0], x[1], x[2]);
+                    FaceClass::Boundary
+                } else {
+                    let mut m = c;
+                    m[d] -= 1;
+                    match (covered(m), covered(c)) {
+                        (true, true) => FaceClass::Covered,
+                        (false, false) => FaceClass::Coarse,
+                        _ => FaceClass::CoarseInterface,
+                    }
+                };
+                coarse_class[d][idx] = class;
+            }
+        }
+        let hf2 = fg.dx * fg.dx;
+        let fid = |c: [usize; 3]| comp.fine_id[fg.cell(c[2], c[1], c[0])];
+        let mut fine_class: [Vec<FaceClass>; 3] = Default::default();
+        let mut fine_area: [Vec<f64>; 3] = Default::default();
+        for d in 0..3 {
+            let nf = n_faces(fg, d);
+            fine_class[d] = vec![FaceClass::Closed; nf];
+            fine_area[d] = vec![0.0; nf];
+            let ap = comp.fine_cut.as_ref().map(|cut| match d {
+                0 => &cut.a_u,
+                1 => &cut.a_v,
+                _ => &cut.a_w,
+            });
+            for idx in 0..nf {
+                let c = face_coords(fg, d, idx);
+                let a = ap.map_or(1.0, |ap| ap[idx]);
+                if c[d] == 0 || c[d] == fdims[d] {
+                    fine_class[d][idx] = FaceClass::FineInterface;
+                    fine_area[d][idx] = hf2;
+                    continue;
+                }
+                let mut m = c;
+                m[d] -= 1;
+                if a > 0.0 && fid(m) != NONE && fid(c) != NONE {
+                    fine_class[d][idx] = FaceClass::Fine;
+                    fine_area[d][idx] = a * hf2;
+                }
+            }
+        }
+        let mut cell_of = vec![NONE; comp.unknowns()];
+        for (cell, &u) in comp.coarse_id.iter().enumerate() {
+            if u != NONE {
+                cell_of[u] = cell;
+            }
+        }
+        for (cell, &u) in comp.fine_id.iter().enumerate() {
+            if u != NONE {
+                cell_of[u] = cell;
+            }
+        }
+        let pre = CompositeSolve::new(&comp, sweeps);
+        let base_rhs = comp.rhs.clone();
+        Self {
+            comp,
+            coarse_class,
+            fine_class,
+            fine_area,
+            pb,
+            base_rhs,
+            cell_of,
+            pre,
+        }
+    }
+
+    /// A zero field.
+    #[must_use]
+    pub fn zeros(&self) -> MacField {
+        let (cg, fg) = (self.comp.coarse, self.comp.fine);
+        MacField {
+            coarse: std::array::from_fn(|d| vec![0.0; n_faces(cg, d)]),
+            fine: std::array::from_fn(|d| vec![0.0; n_faces(fg, d)]),
+        }
+    }
+
+    /// The centre of face `idx` of `axis` on the fine (`fine`) or coarse
+    /// level, in the coarse grid's coordinates.
+    #[must_use]
+    pub fn face_centre(&self, fine: bool, axis: usize, idx: usize) -> [f64; 3] {
+        let (g, o) = if fine {
+            let hc = self.comp.coarse.dx;
+            let lo = self.comp.lo;
+            (
+                self.comp.fine,
+                [lo[0] as f64 * hc, lo[1] as f64 * hc, lo[2] as f64 * hc],
+            )
+        } else {
+            (self.comp.coarse, [0.0; 3])
+        };
+        let c = face_coords(g, axis, idx);
+        let h = g.dx;
+        let mut x = [
+            o[0] + (c[0] as f64 + 0.5) * h,
+            o[1] + (c[1] as f64 + 0.5) * h,
+            o[2] + (c[2] as f64 + 0.5) * h,
+        ];
+        x[axis] = o[axis] + c[axis] as f64 * h;
+        x
+    }
+
+    /// The unknowns on the minus and plus side of fine face `idx` of `axis`
+    /// (`None` outside the patch or on an inactive cell).
+    #[must_use]
+    pub fn fine_face_cells(&self, axis: usize, idx: usize) -> [Option<usize>; 2] {
+        let fg = self.comp.fine;
+        let fdims = [fg.nx, fg.ny, fg.nz];
+        let c = face_coords(fg, axis, idx);
+        let id = |c: [usize; 3]| {
+            let u = self.comp.fine_id[fg.cell(c[2], c[1], c[0])];
+            (u != NONE).then_some(u)
+        };
+        let minus = (c[axis] > 0).then(|| {
+            let mut m = c;
+            m[axis] -= 1;
+            m
+        });
+        [
+            minus.and_then(id),
+            (c[axis] < fdims[axis]).then_some(c).and_then(id),
+        ]
+    }
+
+    /// Point values of `f(axis, x)` on every face (slaved coarse faces
+    /// then averaged down from the fine ones).
+    #[must_use]
+    pub fn sample(&self, f: &(dyn Fn(usize, [f64; 3]) -> f64 + Sync)) -> MacField {
+        let mut u = self.zeros();
+        for d in 0..3 {
+            u.coarse[d]
+                .par_iter_mut()
+                .enumerate()
+                .for_each(|(i, v)| *v = f(d, self.face_centre(false, d, i)));
+            u.fine[d]
+                .par_iter_mut()
+                .enumerate()
+                .for_each(|(i, v)| *v = f(d, self.face_centre(true, d, i)));
+        }
+        self.average_down(&mut u);
+        u
+    }
+
+    /// The slaved coarse faces (interface and covered): flux = the sum of
+    /// the fine fluxes under them.
+    pub fn average_down(&self, u: &mut MacField) {
+        let (cg, fg, lo) = (self.comp.coarse, self.comp.fine, self.comp.lo);
+        let inv = 1.0 / (cg.dx * cg.dx);
+        for d in 0..3 {
+            let (fine, area, class) = (&u.fine[d], &self.fine_area[d], &self.coarse_class[d]);
+            u.coarse[d].par_iter_mut().enumerate().for_each(|(idx, v)| {
+                if !matches!(class[idx], FaceClass::CoarseInterface | FaceClass::Covered) {
+                    return;
+                }
+                let c = face_coords(cg, d, idx);
+                let (t1, t2) = match d {
+                    0 => (1, 2),
+                    1 => (0, 2),
+                    _ => (0, 1),
+                };
+                let mut s = 0.0;
+                for a in 0..2 {
+                    for b in 0..2 {
+                        let mut f = [0usize; 3];
+                        f[d] = 2 * (c[d] - lo[d]);
+                        f[t1] = 2 * (c[t1] - lo[t1]) + a;
+                        f[t2] = 2 * (c[t2] - lo[t2]) + b;
+                        let fi = face_index(fg, d, f);
+                        s += area[fi] * fine[fi];
+                    }
+                }
+                *v = s * inv;
+            });
+        }
+    }
+
+    /// `out = D u (+ wall)`: the outward flux of every unknown.
+    pub fn divergence(&self, u: &MacField, wall: Option<&[f64]>, out: &mut [f64]) {
+        let (cg, fg, nc) = (self.comp.coarse, self.comp.fine, self.comp.n_coarse);
+        let ac = cg.dx * cg.dx;
+        out.par_iter_mut().enumerate().for_each(|(r, o)| {
+            let cell = self.cell_of[r];
+            let mut s = 0.0;
+            if r < nc {
+                let (k, j, i) = (cell / (cg.nx * cg.ny), (cell / cg.nx) % cg.ny, cell % cg.nx);
+                for d in 0..3 {
+                    let mut p = [i, j, k];
+                    let m = face_index(cg, d, p);
+                    p[d] += 1;
+                    let q = face_index(cg, d, p);
+                    s += ac * (u.coarse[d][q] - u.coarse[d][m]);
+                }
+            } else {
+                let (k, j, i) = (cell / (fg.nx * fg.ny), (cell / fg.nx) % fg.ny, cell % fg.nx);
+                for d in 0..3 {
+                    let mut p = [i, j, k];
+                    let m = face_index(fg, d, p);
+                    p[d] += 1;
+                    let q = face_index(fg, d, p);
+                    let a = &self.fine_area[d];
+                    s += a[q] * u.fine[d][q] - a[m] * u.fine[d][m];
+                }
+            }
+            if let Some(w) = wall {
+                s += w[r];
+            }
+            *o = s;
+        });
+    }
+
+    /// `u −= G p` on every owned face, then the slaved faces averaged down.
+    pub fn subtract_gradient(&self, p: &[f64], u: &mut MacField) {
+        let (cg, fg) = (self.comp.coarse, self.comp.fine);
+        let (hc, hf) = (cg.dx, fg.dx);
+        let cdims = [cg.nx, cg.ny, cg.nz];
+        let (cid, fid) = (&self.comp.coarse_id, &self.comp.fine_id);
+        for d in 0..3 {
+            let (class, pb) = (&self.coarse_class[d], &self.pb[d]);
+            u.coarse[d].par_iter_mut().enumerate().for_each(|(idx, v)| {
+                let c = face_coords(cg, d, idx);
+                let pc = |c: [usize; 3]| p[cid[cg.cell(c[2], c[1], c[0])]];
+                let mut m = c;
+                match class[idx] {
+                    FaceClass::Coarse => {
+                        m[d] -= 1;
+                        *v -= (pc(c) - pc(m)) / hc;
+                    }
+                    FaceClass::Boundary => {
+                        if c[d] == 0 {
+                            *v -= (pc(c) - pb[idx]) / (0.5 * hc);
+                        } else {
+                            debug_assert_eq!(c[d], cdims[d]);
+                            m[d] -= 1;
+                            *v -= (pb[idx] - pc(m)) / (0.5 * hc);
+                        }
+                    }
+                    _ => {}
+                }
+            });
+            let class = &self.fine_class[d];
+            u.fine[d].par_iter_mut().enumerate().for_each(|(idx, v)| {
+                if class[idx] != FaceClass::Fine {
+                    return;
+                }
+                let c = face_coords(fg, d, idx);
+                let mut m = c;
+                m[d] -= 1;
+                let pf = |c: [usize; 3]| p[fid[fg.cell(c[2], c[1], c[0])]];
+                *v -= (pf(c) - pf(m)) / hf;
+            });
+        }
+        let inv = 1.0 / (hf * hf);
+        for f in &self.comp.iface {
+            let t: f64 = f.terms.iter().map(|&(col, w)| w * p[col]).sum();
+            // Out of the fine cell: −t / h_f²; along +axis: × side.
+            u.fine[f.axis][f.face] -= f.side as f64 * (-t * inv);
+        }
+        self.average_down(u);
+    }
+
+    /// `G p` as a field (the Dirichlet values on the boundary faces).
+    #[must_use]
+    pub fn gradient(&self, p: &[f64]) -> MacField {
+        let mut g = self.zeros();
+        self.subtract_gradient(p, &mut g);
+        for d in 0..3 {
+            g.coarse[d].iter_mut().for_each(|v| *v = -*v);
+            g.fine[d].iter_mut().for_each(|v| *v = -*v);
+        }
+        g
+    }
+
+    /// The Dirichlet part `b` of `D G p = b − A p`.
+    #[must_use]
+    pub fn boundary_rhs(&self) -> &[f64] {
+        &self.base_rhs
+    }
+
+    /// The wall flux of a body translating at `vel` through the fine cut
+    /// cells (`vel · wall vector`), per unknown.
+    #[must_use]
+    pub fn translating_wall_flux(&self, vel: [f64; 3]) -> Vec<f64> {
+        let mut w = vec![0.0; self.comp.unknowns()];
+        if let Some(cut) = self.comp.fine_cut.as_ref() {
+            for (cell, &u) in self.comp.fine_id.iter().enumerate() {
+                if u != NONE {
+                    let a = cut.wall[cell];
+                    w[u] = vel[0] * a[0] + vel[1] * a[1] + vel[2] * a[2];
+                }
+            }
+        }
+        w
+    }
+
+    /// Project `u` onto the discretely divergence-free fields (plus the wall
+    /// flux): `A p = b − D u`, `u −= G p`. `p` is the initial guess and
+    /// returns the pressure (times dt/ρ in a step).
+    pub fn project(
+        &mut self,
+        u: &mut MacField,
+        p: &mut [f64],
+        wall: Option<&[f64]>,
+        tol: f64,
+        max_it: usize,
+    ) -> ProjectionStats {
+        let t0 = std::time::Instant::now();
+        let n = self.comp.unknowns();
+        self.average_down(u);
+        let mut div = vec![0.0; n];
+        self.divergence(u, wall, &mut div);
+        let div_before = div.iter().fold(0.0f64, |m, v| m.max(v.abs()));
+        let rhs: Vec<f64> = self.base_rhs.iter().zip(&div).map(|(b, d)| b - d).collect();
+        let mut mac_s = t0.elapsed().as_secs_f64();
+        let st = solve_bicgstab_with(&self.comp, &mut self.pre, &rhs, p, tol, max_it);
+        let t1 = std::time::Instant::now();
+        self.subtract_gradient(p, u);
+        self.divergence(u, wall, &mut div);
+        let div_after = div.iter().fold(0.0f64, |m, v| m.max(v.abs()));
+        mac_s += t1.elapsed().as_secs_f64();
+        ProjectionStats {
+            iterations: st.iterations,
+            rel_residual: st.rel_residual,
+            converged: st.converged,
+            div_before,
+            div_after,
+            solve_s: st.solve_s,
+            mac_s,
+        }
+    }
+}
diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/composite/mod.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/composite/mod.rs
index 6e49293..e3b5b35 100644
--- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/composite/mod.rs
+++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/composite/mod.rs
@@ -27,10 +27,14 @@
 //! the interface or the coarse uncovered cells.
 
 mod assemble;
+mod mac;
 mod solve;
 
-pub use assemble::{Composite, CompositeSpec, Interface, Sdf, uniform_problem};
-pub use solve::{CompositeSolve, SolveStats, solve_bicgstab};
+pub use assemble::{
+    Composite, CompositeSpec, Interface, InterfaceFace, Sdf, face_index, uniform_problem,
+};
+pub use mac::{FaceClass, MacField, MacProjection, ProjectionStats};
+pub use solve::{CompositeSolve, SolveStats, solve_bicgstab, solve_bicgstab_with};
 
 /// Compressed sparse rows.
 #[derive(Debug, Clone, Default)]
diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/composite/solve.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/composite/solve.rs
index 5bfa65b..417dc0f 100644
--- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/composite/solve.rs
+++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/composite/solve.rs
@@ -128,9 +128,25 @@ pub fn solve_bicgstab(
     let t0 = std::time::Instant::now();
     let mut pre = CompositeSolve::new(c, sweeps);
     let setup_s = t0.elapsed().as_secs_f64();
+    let mut st = solve_bicgstab_with(c, &mut pre, &c.rhs, x, tol, max_it);
+    st.setup_s = setup_s;
+    st
+}
+
+/// R7-2a: [`solve_bicgstab`] with a prepared preconditioner and an explicit
+/// right-hand side (the projection solves many right-hand sides on one
+/// operator). `setup_s` is 0.
+pub fn solve_bicgstab_with(
+    c: &Composite,
+    pre: &mut CompositeSolve,
+    b: &[f64],
+    x: &mut [f64],
+    tol: f64,
+    max_it: usize,
+) -> SolveStats {
+    let setup_s = 0.0;
     let t1 = std::time::Instant::now();
     let n = c.unknowns();
-    let b = &c.rhs;
     let bn = dot(b, b).sqrt().max(f64::MIN_POSITIVE);
     let mut r = vec![0.0; n];
     c.a.apply(x, &mut r);
diff --git a/crates/specialized/rtx-cfd/tests/embedded3_composite_projection.rs b/crates/specialized/rtx-cfd/tests/embedded3_composite_projection.rs
new file mode 100644
index 0000000..95dd9b7
--- /dev/null
+++ b/crates/specialized/rtx-cfd/tests/embedded3_composite_projection.rs
@@ -0,0 +1,1007 @@
+//! R7-2a: the composite MAC projection on a coarse grid with one nested
+//! ratio-2 patch (`embedded3::composite::MacProjection`, a host prototype
+//! nothing else calls).
+//!
+//! Gates (thresholds asserted, registered before the runs):
+//!
+//! - P1 (`projection_random`): a random field projected — `D G p = b − A p`
+//!   to ≤ 1e-12 (relative to max |A p|) for a random `p`; after the
+//!   projection max |D u| ≤ 1e-11 × max |D u*| in every class of cell (coarse
+//!   interface, fine interface, the rest); the coarse interface fluxes equal
+//!   the sums of their fine fluxes to ≤ 1e-13; the projection idempotent
+//!   (max |P P u − P u| ≤ 1e-10 × max |P u|);
+//! - P2 (`projection_mms_ladder`, ignored): `u* = u_div + ∇q` on the unit
+//!   cube, patch = the middle half, coarse n = 16/32/64: the projected
+//!   velocity against `u_div` at the face centres, split into fine interface
+//!   faces, coarse interface faces, fine interior faces, coarse faces and the
+//!   domain boundary faces; L2 and L∞ orders ≥ 1.8 on the last pair for the
+//!   four non-boundary classes (Quadratic interface; Octree / Direct
+//!   reported);
+//! - P3 (`projection_cut_sphere`, ignored): a Neumann sphere cut into the
+//!   patch (apertures on the fine faces): the div gate of P1, a translating
+//!   body (uniform `u*` + its wall flux) left unchanged, the orders of P2
+//!   away from the body (interface classes ≥ 1.8 in L2), and the composite
+//!   against the uniformly fine projection with the same cut;
+//! - P4 (`projection_cost`, ignored): wall time and unknowns against the
+//!   uniformly fine MAC projection solved by embedded3's production PCG.
+//!
+//! CSVs go to `$R7_OUT` when set.
+
+use rtx_cfd::solvers::incompressible::MultigridParameters;
+use rtx_cfd::solvers::incompressible::embedded3::Grid;
+use rtx_cfd::solvers::incompressible::embedded3::composite::{
+    CompositeSpec, FaceClass, Interface, MacField, MacProjection, Sdf, uniform_problem,
+};
+use rtx_cfd::solvers::incompressible::embedded3::poisson::{Problem, solve_pcg};
+use std::f64::consts::PI;
+use std::io::Write;
+use std::sync::Arc;
+
+type Field3 = fn(usize, [f64; 3]) -> f64;
+type Scalar = fn(f64, f64, f64) -> f64;
+
+fn zero(_x: f64, _y: f64, _z: f64) -> f64 {
+    0.0
+}
+
+fn out_file(name: &str) -> Option<std::fs::File> {
+    let dir = std::env::var("R7_OUT").ok()?;
+    std::fs::create_dir_all(&dir).ok()?;
+    std::fs::File::create(format!("{dir}/{name}")).ok()
+}
+
+fn maxabs(v: &[f64]) -> f64 {
+    v.iter().fold(0.0f64, |m, x| m.max(x.abs()))
+}
+
+struct Lcg(u64);
+
+impl Lcg {
+    fn next(&mut self) -> f64 {
+        self.0 = self
+            .0
+            .wrapping_mul(6_364_136_223_846_793_005)
+            .wrapping_add(1_442_695_040_888_963_407);
+        ((self.0 >> 11) as f64 / (1u64 << 53) as f64) * 2.0 - 1.0
+    }
+}
+
+fn projection(
+    g: Grid,
+    lo: [usize; 3],
+    hi: [usize; 3],
+    iface: Interface,
+    body: Option<Sdf>,
+    dirichlet: Scalar,
+) -> MacProjection {
+    let spec = CompositeSpec {
+        coarse: g,
+        lo,
+        hi,
+        interface: iface,
+        body,
+        source: &zero,
+        dirichlet: &dirichlet,
+    };
+    MacProjection::new(&spec, 2)
+}
+
+// ---------------------------------------------------------------- P1
+
+struct P1 {
+    identity: f64,
+    div_before: f64,
+    /// max |D u| after: coarse interface cells, fine interface cells, rest.
+    div_after: [f64; 3],
+    conservation: f64,
+    idempotence: f64,
+    iterations: usize,
+    rel: f64,
+}
+
+fn random_field(mp: &MacProjection, rng: &mut Lcg) -> MacField {
+    let mut u = mp.zeros();
+    for d in 0..3 {
+        u.coarse[d].iter_mut().for_each(|v| *v = rng.next());
+        u.fine[d].iter_mut().for_each(|v| *v = rng.next());
+    }
+    u
+}
+
+fn class_max(mp: &MacProjection, div: &[f64]) -> [f64; 3] {
+    let mut m = [0.0f64; 3];
+    for (r, v) in div.iter().enumerate() {
+        let c = if !mp.comp.at_interface[r] {
+            2
+        } else if r < mp.comp.n_coarse {
+            0
+        } else {
+            1
+        };
+        m[c] = m[c].max(v.abs());
+    }
+    m
+}
+
+/// max |coarse flux − Σ fine fluxes| over the coarse interface faces,
+/// relative to the largest such flux.
+fn conservation(mp: &MacProjection, u: &MacField) -> f64 {
+    let hc2 = mp.comp.coarse.dx * mp.comp.coarse.dx;
+    let mut sums: [std::collections::HashMap<usize, f64>; 3] = Default::default();
+    for f in &mp.comp.iface {
+        *sums[f.axis].entry(f.coarse_face).or_insert(0.0) +=
+            mp.fine_area[f.axis][f.face] * u.fine[f.axis][f.face];
+    }
+    let (mut d, mut m) = (0.0f64, 0.0f64);
+    for a in 0..3 {
+        for (&cf, &s) in &sums[a] {
+            assert_eq!(mp.coarse_class[a][cf], FaceClass::CoarseInterface);
+            d = d.max((hc2 * u.coarse[a][cf] - s).abs());
+            m = m.max(s.abs());
+        }
+    }
+    d / m.max(f64::MIN_POSITIVE)
+}
+
+fn run_p1(g: Grid, lo: [usize; 3], hi: [usize; 3], body: Option<Sdf>, seed: u64) -> P1 {
+    let mut mp = projection(g, lo, hi, Interface::Quadratic, body, zero);
+    let n = mp.comp.unknowns();
+    let mut rng = Lcg(seed);
+    // D G p = b − A p for a random p.
+    let p: Vec<f64> = (0..n).map(|_| rng.next()).collect();
+    let gp = mp.gradient(&p);
+    let mut dgp = vec![0.0; n];
+    mp.divergence(&gp, None, &mut dgp);
+    let mut ap = vec![0.0; n];
+    mp.comp.a.apply(&p, &mut ap);
+    let b = mp.boundary_rhs();
+    let mut id = 0.0f64;
+    for r in 0..n {
+        id = id.max((dgp[r] - (b[r] - ap[r])).abs());
+    }
+    let identity = id / maxabs(&ap);
+    // Project a random field, twice.
+    let mut u = random_field(&mp, &mut rng);
+    let mut p1 = vec![0.0; n];
+    let st = mp.project(&mut u, &mut p1, None, 1e-13, 400);
+    let mut div = vec![0.0; n];
+    mp.divergence(&u, None, &mut div);
+    let div_after = class_max(&mp, &div);
+    let cons = conservation(&mp, &u);
+    let mut u2 = u.clone();
+    let mut p2 = vec![0.0; n];
+    let _ = mp.project(&mut u2, &mut p2, None, 1e-13, 400);
+    let (mut dd, mut mm) = (0.0f64, 0.0f64);
+    for d in 0..3 {
+        for (a, b) in u.coarse[d].iter().zip(&u2.coarse[d]) {
+            dd = dd.max((a - b).abs());
+            mm = mm.max(a.abs());
+        }
+        for (i, (a, b)) in u.fine[d].iter().zip(&u2.fine[d]).enumerate() {
+            if mp.fine_class[d][i] != FaceClass::Closed {
+                dd = dd.max((a - b).abs());
+                mm = mm.max(a.abs());
+            }
+        }
+    }
+    P1 {
+        identity,
+        div_before: st.div_before,
+        div_after,
+        conservation: cons,
+        idempotence: dd / mm,
+        iterations: st.iterations,
+        rel: st.rel_residual,
+    }
+}
+
+fn check_p1(label: &str, r: &P1, csv: &mut Option<std::fs::File>) {
+    let rel = r.div_after.map(|v| v / r.div_before);
+    eprintln!(
+        "P1 {label}: DG=b-A {:.2e}; div before {:.3e}, after/before [c-iface {:.2e}, f-iface {:.2e}, rest {:.2e}]; conservation {:.2e}; idempotence {:.2e}; {} it rel {:.2e}",
+        r.identity,
+        r.div_before,
+        rel[0],
+        rel[1],
+        rel[2],
+        r.conservation,
+        r.idempotence,
+        r.iterations,
+        r.rel
+    );
+    if let Some(f) = csv.as_mut() {
+        writeln!(
+            f,
+            "{label},{:.3e},{:.3e},{:.3e},{:.3e},{:.3e},{:.3e},{:.3e},{},{:.3e}",
+            r.identity,
+            r.div_before,
+            rel[0],
+            rel[1],
+            rel[2],
+            r.conservation,
+            r.idempotence,
+            r.iterations,
+            r.rel
+        )
+        .ok();
+    }
+    assert!(
+        r.identity <= 1e-12,
+        "{label}: D G != b - A ({:.2e})",
+        r.identity
+    );
+    for (c, v) in rel.iter().enumerate() {
+        assert!(*v <= 1e-11, "{label}: class {c} div {v:.2e}");
+    }
+    assert!(
+        r.conservation <= 1e-13,
+        "{label}: interface not conservative"
+    );
+    assert!(
+        r.idempotence <= 1e-10,
+        "{label}: not idempotent ({:.2e})",
+        r.idempotence
+    );
+}
+
+#[test]
+fn projection_smoke() {
+    let g = Grid::cubic(12, 12, 12, 1.0 / 12.0);
+    let r = run_p1(g, [3; 3], [9; 3], None, 7);
+    check_p1("smoke n12", &r, &mut None);
+}
+
+#[test]
+fn projection_random() {
+    let mut csv = out_file("p1_random.csv");
+    if let Some(f) = csv.as_mut() {
+        writeln!(
+            f,
+            "case,dg_identity,div_before,div_c_iface_rel,div_f_iface_rel,div_rest_rel,conservation,idempotence,iters,rel_res"
+        )
+        .ok();
+    }
+    let cases: [(&str, Grid, [usize; 3], [usize; 3], Option<Sdf>); 3] = [
+        (
+            "cube24_middle",
+            Grid::cubic(24, 24, 24, 1.0 / 24.0),
+            [6; 3],
+            [18; 3],
+            None,
+        ),
+        (
+            "box20x16x12_offcentre",
+            Grid::cubic(20, 16, 12, 1.0 / 20.0),
+            [4, 3, 3],
+            [13, 11, 8],
+            None,
+        ),
+        (
+            "cube24_sphere",
+            Grid::cubic(24, 24, 24, 1.0 / 24.0),
+            [6; 3],
+            [18; 3],
+            Some(sphere()),
+        ),
+    ];
+    for (label, g, lo, hi, body) in cases {
+        let r = run_p1(g, lo, hi, body, 12345);
+        check_p1(label, &r, &mut csv);
+    }
+}
+
+// ---------------------------------------------------------------- P2 / P3
+
+/// `u_div`: each component independent of its own coordinate.
+fn udiv(d: usize, x: [f64; 3]) -> f64 {
+    match d {
+        0 => (1.3 * PI * x[1] + 0.2).cos() * (0.9 * PI * x[2] + 0.6).cos(),
+        1 => (1.1 * PI * x[2] + 0.4).cos() * (1.2 * PI * x[0] + 0.1).cos(),
+        _ => (0.8 * PI * x[0] + 0.3).cos() * (1.4 * PI * x[1] + 0.1).cos(),
+    }
+}
+
+/// `q = ½ sin(πx) sin(2πy) sin(πz)`: q = 0 and ∂²q/∂n² = 0 on the cube.
+fn q_grad(d: usize, x: [f64; 3]) -> f64 {
+    let (a, b, c) = (PI * x[0], 2.0 * PI * x[1], PI * x[2]);
+    0.5 * match d {
+        0 => PI * a.cos() * b.sin() * c.sin(),
+        1 => 2.0 * PI * a.sin() * b.cos() * c.sin(),
+        _ => PI * a.sin() * b.sin() * c.cos(),
+    }
+}
+
+fn ustar_mms(d: usize, x: [f64; 3]) -> f64 {
+    udiv(d, x) + q_grad(d, x)
+}
+
+/// A localized gradient part inside the patch (the case local refinement
+/// is for): `q + exp(−r²/s²)` about the cube's centre, s = 0.06.
+const BUMP_S: f64 = 0.06;
+
+fn q_bump(x: f64, y: f64, z: f64) -> f64 {
+    let r2 = (x - 0.5).powi(2) + (y - 0.5).powi(2) + (z - 0.5).powi(2);
+    (-r2 / (BUMP_S * BUMP_S)).exp()
+}
+
+fn ustar_bump(d: usize, x: [f64; 3]) -> f64 {
+    let g = -2.0 * (x[d] - 0.5) / (BUMP_S * BUMP_S) * q_bump(x[0], x[1], x[2]);
+    ustar_mms(d, x) + g
+}
+
+const SPH_C: [f64; 3] = [0.52, 0.49, 0.51];
+const SPH_R: f64 = 0.12;
+const SPH_A: f64 = 2.0 * PI;
+
+fn sphere() -> Sdf {
+    Arc::new(|x, y, z| {
+        ((x - SPH_C[0]).powi(2) + (y - SPH_C[1]).powi(2) + (z - SPH_C[2]).powi(2)).sqrt() - SPH_R
+    })
+}
+
+/// Potential flow past the sphere (U = 1 along x) plus a rigid swirl about
+/// the z axis through its centre: divergence-free, u·n = 0 on the sphere.
+fn udiv_sph(d: usize, x: [f64; 3]) -> f64 {
+    let r = [x[0] - SPH_C[0], x[1] - SPH_C[1], x[2] - SPH_C[2]];
+    let r2 = r[0] * r[0] + r[1] * r[1] + r[2] * r[2];
+    let rr = r2.sqrt();
+    let r5 = r2 * r2 * rr;
+    let k = 1.5 * SPH_R.powi(3);
+    let pot = match d {
+        0 => 1.0 + 0.5 * SPH_R.powi(3) / (r2 * rr) - k * r[0] * r[0] / r5,
+        _ => -k * r[0] * r[d] / r5,
+    };
+    let swirl = match d {
+        0 => -0.7 * r[1],
+        1 => 0.7 * r[0],
+        _ => 0.0,
+    };
+    pot + swirl
+}
+
+/// `q = cos(a (r − R))`: ∂q/∂n = 0 on the sphere; Dirichlet = q outside.
+fn q_sph(x: f64, y: f64, z: f64) -> f64 {
+    let r = ((x - SPH_C[0]).powi(2) + (y - SPH_C[1]).powi(2) + (z - SPH_C[2]).powi(2)).sqrt();
+    (SPH_A * (r - SPH_R)).cos()
+}
+
+fn q_sph_grad(d: usize, x: [f64; 3]) -> f64 {
+    let r = [x[0] - SPH_C[0], x[1] - SPH_C[1], x[2] - SPH_C[2]];
+    let rr = (r[0] * r[0] + r[1] * r[1] + r[2] * r[2]).sqrt();
+    -SPH_A * (SPH_A * (rr - SPH_R)).sin() * r[d] / rr
+}
+
+fn ustar_sph(d: usize, x: [f64; 3]) -> f64 {
+    udiv_sph(d, x) + q_sph_grad(d, x)
+}
+
+#[derive(Default, Clone, Copy)]
+struct Err {
+    s2: f64,
+    vol: f64,
+    max: f64,
+    count: usize,
+}
+
+impl Err {
+    fn add(&mut self, e: f64, v: f64) {
+        self.s2 += e * e * v;
+        self.vol += v;
+        self.max = self.max.max(e.abs());
+        self.count += 1;
+    }
+    fn l2(&self) -> f64 {
+        (self.s2 / self.vol.max(f64::MIN_POSITIVE)).sqrt()
+    }
+}
+
+const CLASSES: [&str; 6] = ["f_iface", "c_iface", "fine", "coarse", "boundary", "cut"];
+
+/// Face errors against `exact` by class (see `CLASSES`); a fine face is
+/// `cut` when it is partly closed or touches a cut cell.
+fn face_errors(mp: &MacProjection, u: &MacField, exact: Field3) -> [Err; 6] {
+    let mut e = [Err::default(); 6];
+    let (hc, hf) = (mp.comp.coarse.dx, mp.comp.fine.dx);
+    for d in 0..3 {
+        for (i, v) in u.coarse[d].iter().enumerate() {
+            let class = match mp.coarse_class[d][i] {
+                FaceClass::CoarseInterface => 1,
+                FaceClass::Coarse => 3,
+                FaceClass::Boundary => 4,
+                _ => continue,
+            };
+            let x = mp.face_centre(false, d, i);
+            e[class].add(v - exact(d, x), hc * hc * hc);
+        }
+        for (i, v) in u.fine[d].iter().enumerate() {
+            let class = match mp.fine_class[d][i] {
+                FaceClass::FineInterface => 0,
+                FaceClass::Fine => {
+                    let cut_cell = mp
+                        .fine_face_cells(d, i)
+                        .iter()
+                        .flatten()
+                        .any(|&c| mp.comp.cut[c]);
+                    if cut_cell || mp.fine_area[d][i] < hf * hf * (1.0 - 1e-12) {
+                        5
+                    } else {
+                        2
+                    }
+                }
+                _ => continue,
+            };
+            let x = mp.face_centre(true, d, i);
+            e[class].add(v - exact(d, x), mp.fine_area[d][i] * hf);
+        }
+    }
+    e
+}
+
+fn order(a: f64, b: f64) -> f64 {
+    (a / b).log2()
+}
+
+struct Rung {
+    n: usize,
+    err: [Err; 6],
+    div_rel: f64,
+    iterations: usize,
+    unknowns: usize,
+}
+
+fn mms_rung(
+    n: usize,
+    iface: Interface,
+    body: Option<Sdf>,
+    ustar: Field3,
+    exact: Field3,
+    dirichlet: Scalar,
+) -> (Rung, MacProjection, MacField) {
+    let g = Grid::cubic(n, n, n, 1.0 / n as f64);
+    let mut mp = projection(g, [n / 4; 3], [3 * n / 4; 3], iface, body, dirichlet);
+    let mut u = mp.sample(&ustar);
+    let mut p = vec![0.0; mp.comp.unknowns()];
+    let st = mp.project(&mut u, &mut p, None, 1e-13, 400);
+    assert!(st.converged, "n {n}: {st:?}");
+    let err = face_errors(&mp, &u, exact);
+    let r = Rung {
+        n,
+        err,
+        div_rel: st.div_after / st.div_before,
+        iterations: st.iterations,
+        unknowns: mp.comp.unknowns(),
+    };
+    (r, mp, u)
+}
+
+fn ladder_report(
+    name: &str,
+    rungs: &[Rung],
+    csv: &mut Option<std::fs::File>,
+) -> Vec<[(f64, f64); 6]> {
+    let mut orders = Vec::new();
+    for (w, r) in rungs.iter().enumerate() {
+        for (c, label) in CLASSES.iter().enumerate() {
+            let e = r.err[c];
+            if e.count == 0 {
+                continue;
+            }
+            let (o2, oi) = if w > 0 {
+                let p = rungs[w - 1].err[c];
+                (order(p.l2(), e.l2()), order(p.max, e.max))
+            } else {
+                (f64::NAN, f64::NAN)
+            };
+            eprintln!(
+                "{name} n {:3} {label:9} faces {:8} L2 {:.3e} Linf {:.3e}  orders {o2:5.2} {oi:5.2}",
+                r.n,
+                e.count,
+                e.l2(),
+                e.max
+            );
+            if let Some(f) = csv.as_mut() {
+                writeln!(
+                    f,
+                    "{name},{},{label},{},{:.6e},{:.6e},{o2:.3},{oi:.3},{:.3e},{},{}",
+                    r.n,
+                    e.count,
+                    e.l2(),
+                    e.max,
+                    r.div_rel,
+                    r.iterations,
+                    r.unknowns
+                )
+                .ok();
+            }
+        }
+        if w > 0 {
+            let mut o = [(f64::NAN, f64::NAN); 6];
+            for c in 0..6 {
+                let (p, e) = (rungs[w - 1].err[c], r.err[c]);
+                if e.count > 0 && p.count > 0 {
+                    o[c] = (order(p.l2(), e.l2()), order(p.max, e.max));
+                }
+            }
+            orders.push(o);
+        }
+    }
+    orders
+}
+
+const CSV_HEAD: &str =
+    "case,n,class,faces,l2,linf,order_l2,order_linf,div_after_rel,iters,unknowns";
+
+#[test]
+#[ignore = "P2 ladder (minutes)"]
+fn projection_mms_ladder() {
+    let mut csv = out_file("p2_mms.csv");
+    if let Some(f) = csv.as_mut() {
+        writeln!(f, "{CSV_HEAD}").ok();
+    }
+    let mut verdict = Vec::new();
+    for (name, iface) in [
+        ("quadratic", Interface::Quadratic),
+        ("octree", Interface::Octree),
+        ("direct", Interface::Direct),
+    ] {
+        let rungs: Vec<Rung> = [16, 32, 64]
+            .iter()
+            .map(|&n| mms_rung(n, iface, None, ustar_mms, udiv, zero).0)
+            .collect();
+        for r in &rungs {
+            assert!(
+                r.div_rel <= 1e-10,
+                "{name} n {}: div {:.2e}",
+                r.n,
+                r.div_rel
+            );
+        }
+        let orders = ladder_report(name, &rungs, &mut csv);
+        if iface == Interface::Quadratic {
+            let last = orders.last().expect("pair");
+            for c in 0..4 {
+                let (o2, oi) = last[c];
+                let ok = o2 >= 1.8 && oi >= 1.8;
+                eprintln!(
+                    "P2 GATE {} L2 {o2:.2} Linf {oi:.2} -> {}",
+                    CLASSES[c],
+                    if ok { "HELD" } else { "FAILED" }
+                );
+                verdict.push((CLASSES[c], ok));
+            }
+        }
+    }
+    for (c, ok) in verdict {
+        assert!(ok, "P2 order gate failed on {c}");
+    }
+}
+
+// ------------------------------------------------ uniform MAC comparator
+
+/// A uniform-grid MAC projection in the same integrated form (the
+/// comparator of P3/P4): open faces are those with a positive coefficient,
+/// the Dirichlet boundary as on the composite.
+struct UniformMac {
+    g: Grid,
+    prob: Problem,
+    /// Per axis: face area (0 closed) and the Dirichlet value on boundary
+    /// faces.
+    area: [Vec<f64>; 3],
+    pb: [Vec<f64>; 3],
+}
+
+fn ufi(g: Grid, d: usize, c: [usize; 3]) -> usize {
+    match d {
+        0 => g.uface(c[2], c[1], c[0]),
+        1 => g.vface(c[2], c[1], c[0]),
+        _ => g.wface(c[2], c[1], c[0]),
+    }
+}
+
+impl UniformMac {
+    fn new(n: usize, body: Option<&Sdf>, dirichlet: Scalar) -> Self {
+        let g = Grid::cubic(n, n, n, 1.0 / n as f64);
+        let (prob, _) = uniform_problem(g, body, &zero, &dirichlet);
+        let h = g.dx;
+        let dims = [n; 3];
+        let mut area: [Vec<f64>; 3] = Default::default();
+        let mut pb: [Vec<f64>; 3] = Default::default();
+        for d in 0..3 {
+            let nf =
+                (n + usize::from(d == 0)) * (n + usize::from(d == 1)) * (n + usize::from(d == 2));
+            area[d] = vec![0.0; nf];
+            pb[d] = vec![0.0; nf];
+        }
+        for k in 0..n {
+            for j in 0..n {
+                for i in 0..n {
+                    let idx = g.cell(k, j, i);
+                    if !prob.active[idx] {
+                        continue;
+                    }
+                    // The plus face of each axis, and the minus face on the
+                    // domain boundary.
+                    for d in 0..3 {
+                        let mut c = [i, j, k];
+                        let coef = match d {
+                            0 => prob.ae[idx],
+                            1 => prob.an[idx],
+                            _ => prob.at[idx],
+                        };
+                        let centre = |c: [usize; 3]| {
+                            let mut x = [
+                                (c[0] as f64 + 0.5) * h,
+                                (c[1] as f64 + 0.5) * h,
+                                (c[2] as f64 + 0.5) * h,
+                            ];
+                            x[d] = c[d] as f64 * h;
+                            x
+                        };
+                        if c[d] == 0 {
+                            let f = ufi(g, d, c);
+                            area[d][f] = h * h;
+                            let x = centre(c);
+                            pb[d][f] = dirichlet(x[0], x[1], x[2]);
+                        }
+                        c[d] += 1;
+                        let f = ufi(g, d, c);
+                        if c[d] == dims[d] {
+                            area[d][f] = h * h;
+                            let x = centre(c);
+                            pb[d][f] = dirichlet(x[0], x[1], x[2]);
+                        } else if coef > 0.0 {
+                            area[d][f] = coef * h;
+                        }
+                    }
+                }
+            }
+        }
+        Self { g, prob, area, pb }
+    }
+
+    fn sample(&self, f: Field3) -> [Vec<f64>; 3] {
+        let (g, h) = (self.g, self.g.dx);
+        std::array::from_fn(|d| {
+            (0..self.area[d].len())
+                .map(|idx| {
+                    let nx = g.nx + usize::from(d == 0);
+                    let ny = g.ny + usize::from(d == 1);
+                    let c = [idx % nx, (idx / nx) % ny, idx / (nx * ny)];
+                    let mut x = [
+                        (c[0] as f64 + 0.5) * h,
+                        (c[1] as f64 + 0.5) * h,
+                        (c[2] as f64 + 0.5) * h,
+                    ];
+                    x[d] = c[d] as f64 * h;
+                    f(d, x)
+                })
+                .collect()
+        })
+    }
+
+    fn divergence(&self, u: &[Vec<f64>; 3], out: &mut [f64]) {
+        let g = self.g;
+        for k in 0..g.nz {
+            for j in 0..g.ny {
+                for i in 0..g.nx {
+                    let idx = g.cell(k, j, i);
+                    if !self.prob.active[idx] {
+                        out[idx] = 0.0;
+                        continue;
+                    }
+                    let mut s = 0.0;
+                    for d in 0..3 {
+                        let mut c = [i, j, k];
+                        let m = ufi(g, d, c);
+                        c[d] += 1;
+                        let p = ufi(g, d, c);
+                        s += self.area[d][p] * u[d][p] - self.area[d][m] * u[d][m];
+                    }
+                    out[idx] = s;
+                }
+            }
+        }
+    }
+
+    fn subtract_gradient(&self, p: &[f64], u: &mut [Vec<f64>; 3]) {
+        let (g, h) = (self.g, self.g.dx);
+        let n = [g.nx, g.ny, g.nz];
+        for d in 0..3 {
+            let nx = g.nx + usize::from(d == 0);
+            let ny = g.ny + usize::from(d == 1);
+            for idx in 0..u[d].len() {
+                if self.area[d][idx] == 0.0 {
+                    continue;
+                }
+                let c = [idx % nx, (idx / nx) % ny, idx / (nx * ny)];
+                let cell = |c: [usize; 3]| g.cell(c[2], c[1], c[0]);
+                let mut m = c;
+                if c[d] == 0 {
+                    u[d][idx] -= (p[cell(c)] - self.pb[d][idx]) / (0.5 * h);
+                } else if c[d] == n[d] {
+                    m[d] -= 1;
+                    u[d][idx] -= (self.pb[d][idx] - p[cell(m)]) / (0.5 * h);
+                } else {
+                    m[d] -= 1;
+                    u[d][idx] -= (p[cell(c)] - p[cell(m)]) / h;
+                }
+            }
+        }
+    }
+
+    /// Returns (setup s, solve s, mac s, iterations, max |div| after).
+    fn project(&self, u: &mut [Vec<f64>; 3]) -> (f64, f64, f64, usize, f64) {
+        let t0 = std::time::Instant::now();
+        let nc = self.g.cells();
+        let mut div = vec![0.0; nc];
+        self.divergence(u, &mut div);
+        let mut prob = self.prob.clone();
+        for (r, d) in prob.rhs.iter_mut().zip(&div) {
+            *r -= d;
+        }
+        let mut mac_s = t0.elapsed().as_secs_f64();
+        let l1: f64 = prob.rhs.iter().map(|v| v.abs()).sum();
+        let params = MultigridParameters {
+            max_iterations: 2000,
+            ..MultigridParameters::default()
+        };
+        let mut p = vec![0.0; nc];
+        let sol = solve_pcg(&prob, &mut p, &params, 1e-12 * l1, None);
+        assert!(sol.converged, "uniform projection did not converge");
+        let t1 = std::time::Instant::now();
+        self.subtract_gradient(&p, u);
+        self.divergence(u, &mut div);
+        mac_s += t1.elapsed().as_secs_f64();
+        (
+            sol.setup_ns as f64 * 1e-9,
+            sol.iterate_ns as f64 * 1e-9,
+            mac_s,
+            sol.iterations,
+            maxabs(&div),
+        )
+    }
+}
+
+#[test]
+#[ignore = "P3 cut sphere ladder (minutes)"]
+fn projection_cut_sphere() {
+    let mut csv = out_file("p3_sphere.csv");
+    if let Some(f) = csv.as_mut() {
+        writeln!(f, "{CSV_HEAD}").ok();
+    }
+    let mut extra = out_file("p3_sphere_checks.csv");
+    if let Some(f) = extra.as_mut() {
+        writeln!(
+            f,
+            "n,div_after_rel,cut_cells,translate_div_before,translate_change,vs_fine_max_diff,vs_fine_exact_err_fine,vs_fine_exact_err_composite"
+        )
+        .ok();
+    }
+    let mut rungs = Vec::new();
+    for n in [16usize, 32, 64] {
+        let (r, mut mp, u) = mms_rung(
+            n,
+            Interface::Quadratic,
+            Some(sphere()),
+            ustar_sph,
+            udiv_sph,
+            q_sph,
+        );
+        let cut_cells = mp.comp.cut.iter().filter(|&&c| c).count();
+        // A translating body: uniform u* plus its wall flux is already
+        // divergence-free and must come out unchanged.
+        let vel = [0.3, -0.2, 0.5];
+        let wall = mp.translating_wall_flux(vel);
+        let mut ut = mp.sample(&move |d, _x| vel[d]);
+        let u0 = ut.clone();
+        let mut pt = vec![0.0; mp.comp.unknowns()];
+        let mut div0 = vec![0.0; mp.comp.unknowns()];
+        mp.divergence(&ut, Some(&wall), &mut div0);
+        let translate_div = maxabs(&div0);
+        let st = mp.project(&mut ut, &mut pt, Some(&wall), 1e-12, 400);
+        // By linearity P(u0) − u0 − P(0) = G A⁻¹ (D u0 + wall): the change
+        // beyond what the Dirichlet data (q on the boundary) drives alone.
+        let mut uz = mp.zeros();
+        let mut pz = vec![0.0; mp.comp.unknowns()];
+        let _ = mp.project(&mut uz, &mut pz, None, 1e-12, 400);
+        let mut change = 0.0f64;
+        for d in 0..3 {
+            for (i, ((a, b), z)) in ut.coarse[d]
+                .iter()
+                .zip(&u0.coarse[d])
+                .zip(&uz.coarse[d])
+                .enumerate()
+            {
+                if mp.coarse_class[d][i] != FaceClass::Covered {
+                    change = change.max((a - b - z).abs());
+                }
+            }
+            for (i, ((a, b), z)) in ut.fine[d]
+                .iter()
+                .zip(&u0.fine[d])
+                .zip(&uz.fine[d])
+                .enumerate()
+            {
+                if mp.fine_class[d][i] != FaceClass::Closed {
+                    change = change.max((a - b - z).abs());
+                }
+            }
+        }
+        assert!(st.converged);
+        // Against the uniformly fine projection with the same cut (on the
+        // fine faces of the patch, which coincide).
+        let um = UniformMac::new(2 * n, Some(&sphere()), q_sph);
+        let mut uu = um.sample(ustar_sph);
+        let _ = um.project(&mut uu);
+        let lo = mp.comp.lo;
+        let off = [2 * lo[0], 2 * lo[1], 2 * lo[2]];
+        let fg = mp.comp.fine;
+        let (mut dmax, mut efine, mut ecomp) = (0.0f64, 0.0f64, 0.0f64);
+        for d in 0..3 {
+            let nx = fg.nx + usize::from(d == 0);
+            let ny = fg.ny + usize::from(d == 1);
+            for (i, v) in u.fine[d].iter().enumerate() {
+                if mp.fine_class[d][i] != FaceClass::Fine {
+                    continue;
+                }
+                let c = [i % nx, (i / nx) % ny, i / (nx * ny)];
+                let gi = ufi(um.g, d, [c[0] + off[0], c[1] + off[1], c[2] + off[2]]);
+                let x = mp.face_centre(true, d, i);
+                let ex = udiv_sph(d, x);
+                dmax = dmax.max((v - uu[d][gi]).abs());
+                efine = efine.max((uu[d][gi] - ex).abs());
+                ecomp = ecomp.max((v - ex).abs());
+            }
+        }
+        eprintln!(
+            "P3 n {n}: div after/before {:.2e}, cut cells {cut_cells}, translating: div {translate_div:.2e} change {change:.2e}; vs uniform fine: max diff {dmax:.3e} (fine err {efine:.3e}, composite err {ecomp:.3e})",
+            r.div_rel
+        );
+        if let Some(f) = extra.as_mut() {
+            writeln!(
+                f,
+                "{n},{:.3e},{cut_cells},{translate_div:.3e},{change:.3e},{dmax:.3e},{efine:.3e},{ecomp:.3e}",
+                r.div_rel
+            )
+            .ok();
+        }
+        assert!(r.div_rel <= 1e-10, "P3 n {n}: div {:.2e}", r.div_rel);
+        assert!(
+            translate_div <= 1e-12,
+            "P3 n {n}: translating div {translate_div:.2e}"
+        );
+        assert!(change <= 1e-10, "P3 n {n}: translating change {change:.2e}");
+        rungs.push(r);
+    }
+    let orders = ladder_report("sphere", &rungs, &mut csv);
+    let last = orders.last().expect("pair");
+    let mut ok_all = true;
+    for c in [0usize, 1] {
+        let (o2, oi) = last[c];
+        let ok = o2 >= 1.8;
+        ok_all &= ok;
+        eprintln!(
+            "P3 GATE {} L2 {o2:.2} (Linf {oi:.2} reported) -> {}",
+            CLASSES[c],
+            if ok { "HELD" } else { "FAILED" }
+        );
+    }
+    for c in [2usize, 3, 4, 5] {
+        let (o2, oi) = last[c];
+        eprintln!("P3 report {} L2 {o2:.2} Linf {oi:.2}", CLASSES[c]);
+    }
+    assert!(ok_all, "P3 interface order gate failed");
+}
+
+#[test]
+#[ignore = "P4 cost (minutes)"]
+fn projection_cost() {
+    let mut csv = out_file("p4_cost.csv");
+    if let Some(f) = csv.as_mut() {
+        writeln!(
+            f,
+            "field,case,n,patch,unknowns,iters,build_s,solve_s,mac_s,total_s,l2_patch,linf_patch,div_after"
+        )
+        .ok();
+    }
+    let n = 64usize;
+    let h = 1.0 / n as f64;
+    type Case = (&'static str, Field3, Scalar);
+    let fields: [Case; 2] = [("smooth", ustar_mms, zero), ("bump", ustar_bump, q_bump)];
+    for ((field, ustar, dirichlet), (label, lo, hi)) in fields
+        .into_iter()
+        .flat_map(|f| [("half", 16usize, 48usize), ("quarter", 24, 40)].map(move |p| (f, p)))
+    {
+        let t0 = std::time::Instant::now();
+        let g = Grid::cubic(n, n, n, h);
+        let mut mp = projection(g, [lo; 3], [hi; 3], Interface::Quadratic, None, dirichlet);
+        let build_s = t0.elapsed().as_secs_f64();
+        let mut u = mp.sample(&ustar);
+        let mut p = vec![0.0; mp.comp.unknowns()];
+        let st = mp.project(&mut u, &mut p, None, 1e-10, 400);
+        assert!(st.converged);
+        let mut e = Err::default();
+        let hf = mp.comp.fine.dx;
+        for d in 0..3 {
+            for (i, v) in u.fine[d].iter().enumerate() {
+                let x = mp.face_centre(true, d, i);
+                e.add(v - udiv(d, x), hf * hf * hf);
+            }
+        }
+        let total = build_s + st.solve_s + st.mac_s;
+        eprintln!(
+            "P4 {field} composite {label}: {} unknowns, {} it, build+setup {build_s:.2} s, solve {:.2} s, mac {:.2} s, total {total:.2} s, patch L2 {:.3e} Linf {:.3e}, div {:.2e}",
+            mp.comp.unknowns(),
+            st.iterations,
+            st.solve_s,
+            st.mac_s,
+            e.l2(),
+            e.max,
+            st.div_after
+        );
+        if let Some(f) = csv.as_mut() {
+            writeln!(
+                f,
+                "{field},composite,{n},{lo}..{hi},{},{},{build_s:.3},{:.3},{:.3},{total:.3},{:.6e},{:.6e},{:.3e}",
+                mp.comp.unknowns(),
+                st.iterations,
+                st.solve_s,
+                st.mac_s,
+                e.l2(),
+                e.max,
+                st.div_after
+            )
+            .ok();
+        }
+        for m in [2 * n, n] {
+            let t0 = std::time::Instant::now();
+            let um = UniformMac::new(m, None, dirichlet);
+            let build_s = t0.elapsed().as_secs_f64();
+            let mut uu = um.sample(ustar);
+            let (setup_s, solve_s, mac_s, it, div) = um.project(&mut uu);
+            let (x0, x1) = (lo as f64 * h, hi as f64 * h);
+            let mut e = Err::default();
+            let hm = um.g.dx;
+            for d in 0..3 {
+                let nx = m + usize::from(d == 0);
+                let ny = m + usize::from(d == 1);
+                for (i, v) in uu[d].iter().enumerate() {
+                    let c = [i % nx, (i / nx) % ny, i / (nx * ny)];
+                    let mut x = [
+                        (c[0] as f64 + 0.5) * hm,
+                        (c[1] as f64 + 0.5) * hm,
+                        (c[2] as f64 + 0.5) * hm,
+                    ];
+                    x[d] = c[d] as f64 * hm;
+                    if x.iter().all(|&v| v >= x0 - 1e-12 && v <= x1 + 1e-12) {
+                        e.add(v - udiv(d, x), hm * hm * hm);
+                    }
+                }
+            }
+            let total = build_s + setup_s + solve_s + mac_s;
+            let case = if m == 2 * n {
+                "uniform-fine"
+            } else {
+                "uniform-coarse"
+            };
+            eprintln!(
+                "P4 {field} {case} n {m} (region {label}): {} cells, {it} it, build+setup {:.2} s, solve {solve_s:.2} s, mac {mac_s:.2} s, total {total:.2} s, region L2 {:.3e} Linf {:.3e}, div {div:.2e}",
+                m * m * m,
+                build_s + setup_s,
+                e.l2(),
+                e.max
+            );
+            if let Some(f) = csv.as_mut() {
+                writeln!(
+                    f,
+                    "{field},{case},{m},{lo}..{hi},{},{it},{:.3},{solve_s:.3},{mac_s:.3},{total:.3},{:.6e},{:.6e},{div:.3e}",
+                    m * m * m,
+                    build_s + setup_s,
+                    e.l2(),
+                    e.max
+                )
+                .ok();
+            }
+        }
+    }
+}