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
new file mode 100644
index 0000000..5ae8f7d
--- /dev/null
+++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/composite/assemble.rs
@@ -0,0 +1,519 @@
+//! The composite operator: coarse–coarse, fine–fine and coarse–fine faces,
+//! the Dirichlet outer boundary, the cut apertures of an optional body on
+//! the fine level, and the uniform coarse [`Problem`] the preconditioner's
+//! coarse correction solves.
+
+use super::Csr;
+use crate::solvers::incompressible::embedded3::poisson::Problem;
+use crate::solvers::incompressible::embedded3::{Body, CutGeometry, Grid};
+
+/// The coarse–fine face flux (see the module docs).
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub enum Interface {
+    Direct,
+    Octree,
+    Quadratic,
+}
+
+/// What to build. The coarse grid must be cubic; the patch `lo ≤ (i, j, k)
+/// < hi` (coarse indices) keeps two coarse cells from the domain boundary.
+pub struct CompositeSpec<'a> {
+    pub coarse: Grid,
+    pub lo: [usize; 3],
+    pub hi: [usize; 3],
+    pub interface: Interface,
+    /// A body cutting the fine level (φ > 0 fluid), Neumann on its wall.
+    pub body: Option<Sdf>,
+    /// `f = −Δp` (integrated by the midpoint rule over the fluid volume).
+    pub source: &'a (dyn Fn(f64, f64, f64) -> f64 + Sync),
+    /// The Dirichlet value on the outer boundary.
+    pub dirichlet: &'a (dyn Fn(f64, f64, f64) -> f64 + Sync),
+}
+
+/// The assembled composite problem.
+pub struct Composite {
+    pub coarse: Grid,
+    pub fine: Grid,
+    pub lo: [usize; 3],
+    pub hi: [usize; 3],
+    pub interface: Interface,
+    /// Per coarse cell: its unknown, or `usize::MAX` (covered / solid).
+    pub coarse_id: Vec<usize>,
+    /// Per fine cell: its unknown, or `usize::MAX` (solid).
+    pub fine_id: Vec<usize>,
+    /// Unknowns `0..n_coarse` are coarse, the rest fine.
+    pub n_coarse: usize,
+    /// Per unknown: the coarse cell holding it (itself, or the fine
+    /// cell's parent) — the restriction / prolongation map.
+    pub parent: Vec<usize>,
+    /// Per unknown: the cell centre and the fluid volume.
+    pub centre: Vec<[f64; 3]>,
+    pub volume: Vec<f64>,
+    /// Per unknown: the cell owns a coarse–fine face.
+    pub at_interface: Vec<bool>,
+    /// Per unknown: the cell is cut by the body (0 < fluid fraction < 1).
+    pub cut: Vec<bool>,
+    pub a: Csr,
+    pub rhs: Vec<f64>,
+    /// The uniform coarse operator (coarse apertures under the patch too).
+    pub coarse_problem: Problem,
+}
+
+const NONE: usize = usize::MAX;
+/// (axis, sign) of the six faces: −x, +x, −y, +y, −z, +z.
+const DIRS: [(usize, isize); 6] = [(0, -1), (0, 1), (1, -1), (1, 1), (2, -1), (2, 1)];
+
+/// A signed distance (φ > 0 fluid) shared by the coarse and fine builds.
+pub type Sdf = std::sync::Arc<dyn Fn(f64, f64, f64) -> f64 + Send + Sync>;
+
+/// The body of `phi` on a grid starting at `origin`.
+fn body_at(phi: &Sdf, origin: [f64; 3]) -> Body {
+    let phi = phi.clone();
+    Body::from_sdf(move |x, y, z, _t| phi(x + origin[0], y + origin[1], z + origin[2]))
+}
+
+/// The fluid aperture of face `dir` of cell `c = (i, j, k)` of `cut`.
+fn aperture(cut: Option<&CutGeometry>, g: Grid, c: [usize; 3], dir: usize) -> f64 {
+    let Some(cut) = cut else { return 1.0 };
+    let [i, j, k] = c;
+    match dir {
+        0 => cut.a_u[g.uface(k, j, i)],
+        1 => cut.a_u[g.uface(k, j, i + 1)],
+        2 => cut.a_v[g.vface(k, j, i)],
+        3 => cut.a_v[g.vface(k, j + 1, i)],
+        4 => cut.a_w[g.wface(k, j, i)],
+        _ => cut.a_w[g.wface(k + 1, j, i)],
+    }
+}
+
+fn fraction(cut: Option<&CutGeometry>, g: Grid, c: [usize; 3]) -> f64 {
+    cut.map_or(1.0, |cut| cut.vol[g.cell(c[2], c[1], c[0])])
+}
+
+fn active(cut: Option<&CutGeometry>, g: Grid, c: [usize; 3]) -> bool {
+    fraction(cut, g, c) > 0.0 && (0..6).any(|d| aperture(cut, g, c, d) > 0.0)
+}
+
+impl Composite {
+    #[must_use]
+    pub fn build(spec: &CompositeSpec<'_>) -> Self {
+        let cg = spec.coarse;
+        let (lo, hi) = (spec.lo, spec.hi);
+        let hc = cg.dx;
+        assert!(
+            (cg.dy - hc).abs() < 1e-14 * hc && (cg.dz - hc).abs() < 1e-14 * hc,
+            "cubic cells"
+        );
+        let dims = [cg.nx, cg.ny, cg.nz];
+        for d in 0..3 {
+            assert!(
+                lo[d] >= 2 && hi[d] + 2 <= dims[d] && hi[d] > lo[d],
+                "patch margin"
+            );
+        }
+        let hf = 0.5 * hc;
+        let fg = Grid::cubic(
+            2 * (hi[0] - lo[0]),
+            2 * (hi[1] - lo[1]),
+            2 * (hi[2] - lo[2]),
+            hf,
+        );
+        let fdims = [fg.nx, fg.ny, fg.nz];
+        let origin = [lo[0] as f64 * hc, lo[1] as f64 * hc, lo[2] as f64 * hc];
+        let fine_cut = spec
+            .body
+            .as_ref()
+            .map(|b| CutGeometry::build(&body_at(b, origin), fg, 0.0));
+        let coarse_cut = spec
+            .body
+            .as_ref()
+            .map(|b| CutGeometry::build(&body_at(b, [0.0; 3]), cg, 0.0));
+        let (fcut, ccut) = (fine_cut.as_ref(), coarse_cut.as_ref());
+        let covered = |c: [usize; 3]| (0..3).all(|d| c[d] >= lo[d] && c[d] < hi[d]);
+
+        // Unknowns.
+        let mut coarse_id = vec![NONE; cg.cells()];
+        let mut fine_id = vec![NONE; fg.cells()];
+        let (mut parent, mut centre, mut volume, mut cut_flag) =
+            (Vec::new(), Vec::new(), Vec::new(), Vec::new());
+        for k in 0..cg.nz {
+            for j in 0..cg.ny {
+                for i in 0..cg.nx {
+                    let c = [i, j, k];
+                    if covered(c) {
+                        continue;
+                    }
+                    assert!(
+                        fraction(ccut, cg, c) > 1.0 - 1e-12,
+                        "the body must stay inside the patch"
+                    );
+                    coarse_id[cg.cell(k, j, i)] = centre.len();
+                    parent.push(cg.cell(k, j, i));
+                    centre.push([
+                        (i as f64 + 0.5) * hc,
+                        (j as f64 + 0.5) * hc,
+                        (k as f64 + 0.5) * hc,
+                    ]);
+                    volume.push(hc * hc * hc);
+                    cut_flag.push(false);
+                }
+            }
+        }
+        let n_coarse = centre.len();
+        for k in 0..fg.nz {
+            for j in 0..fg.ny {
+                for i in 0..fg.nx {
+                    let c = [i, j, k];
+                    if !active(fcut, fg, c) {
+                        continue;
+                    }
+                    let frac = fraction(fcut, fg, c);
+                    fine_id[fg.cell(k, j, i)] = centre.len();
+                    parent.push(cg.cell(lo[2] + k / 2, lo[1] + j / 2, lo[0] + i / 2));
+                    centre.push([
+                        origin[0] + (i as f64 + 0.5) * hf,
+                        origin[1] + (j as f64 + 0.5) * hf,
+                        origin[2] + (k as f64 + 0.5) * hf,
+                    ]);
+                    volume.push(frac * hf * hf * hf);
+                    cut_flag.push(frac < 1.0 - 1e-9);
+                }
+            }
+        }
+        let n = centre.len();
+        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 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];
+            o[d] += s;
+            o
+        };
+        let cid = |c: [isize; 3]| coarse_id[cg.cell(c[2] as usize, c[1] as usize, c[0] as usize)];
+        let fid = |c: [usize; 3]| fine_id[fg.cell(c[2], c[1], c[0])];
+        // The fine children of covered coarse cell `k` (coarse coords).
+        let children = |kc: [usize; 3]| -> Vec<[usize; 3]> {
+            let b = [
+                2 * (kc[0] - lo[0]),
+                2 * (kc[1] - lo[1]),
+                2 * (kc[2] - lo[2]),
+            ];
+            let mut v = Vec::with_capacity(8);
+            for dk in 0..2 {
+                for dj in 0..2 {
+                    for di in 0..2 {
+                        v.push([b[0] + di, b[1] + dj, b[2] + dk]);
+                    }
+                }
+            }
+            v
+        };
+        // The quadratic ghost of fine cell `f` across its face (d, s): the
+        // weights on unknowns.
+        let ghost = |f: [usize; 3], d: usize, s: isize| -> Vec<(usize, f64)> {
+            let mut w = Vec::with_capacity(12);
+            let mut fin = f;
+            fin[d] = (f[d] as isize - s) as usize;
+            w.push((fid(fin), -0.2));
+            w.push((fid(f), 2.0 / 3.0));
+            // The coarse cell holding the ghost point.
+            let gf = [
+                (2 * lo[0] + f[0]) as isize,
+                (2 * lo[1] + f[1]) as isize,
+                (2 * lo[2] + f[2]) as isize,
+            ];
+            let mut gc = gf;
+            gc[d] += s;
+            let cc = [
+                gc[0].div_euclid(2),
+                gc[1].div_euclid(2),
+                gc[2].div_euclid(2),
+            ];
+            let (t1, t2) = match d {
+                0 => (1, 2),
+                1 => (0, 2),
+                _ => (0, 1),
+            };
+            let off = |t: usize| if gf[t] % 2 == 0 { -0.25 } else { 0.25 };
+            let (a, b) = (off(t1), off(t2));
+            let wn = 8.0 / 15.0;
+            let at = |dt1: isize, dt2: isize| {
+                let mut c = cc;
+                c[t1] += dt1;
+                c[t2] += dt2;
+                let id = cid(c);
+                assert!(
+                    id != NONE,
+                    "interpolation stencil on a covered/solid coarse cell"
+                );
+                id
+            };
+            w.push((at(0, 0), wn * (1.0 - a * a - b * b)));
+            w.push((at(-1, 0), wn * (0.5 * a * a - 0.5 * a)));
+            w.push((at(1, 0), wn * (0.5 * a * a + 0.5 * a)));
+            w.push((at(0, -1), wn * (0.5 * b * b - 0.5 * b)));
+            w.push((at(0, 1), wn * (0.5 * b * b + 0.5 * b)));
+            let x = wn * 0.25 * a * b;
+            w.push((at(1, 1), x));
+            w.push((at(-1, -1), x));
+            w.push((at(1, -1), -x));
+            w.push((at(-1, 1), -x));
+            w
+        };
+
+        // Coarse rows.
+        for k in 0..cg.nz {
+            for j in 0..cg.ny {
+                for i in 0..cg.nx {
+                    let row = coarse_id[cg.cell(k, j, i)];
+                    if row == NONE {
+                        continue;
+                    }
+                    let c = [i, j, k];
+                    let x = centre[row];
+                    rhs[row] += volume[row] * (spec.source)(x[0], x[1], x[2]);
+                    for &(d, s) in &DIRS {
+                        let nb = step(c, d, s);
+                        if nb[d] < 0 || nb[d] >= dims[d] as isize {
+                            let coef = 2.0 * hc;
+                            let mut fx = x;
+                            fx[d] += s as f64 * 0.5 * hc;
+                            rows[row].push((row, coef));
+                            rhs[row] += coef * (spec.dirichlet)(fx[0], fx[1], fx[2]);
+                            continue;
+                        }
+                        let nbu = [nb[0] as usize, nb[1] as usize, nb[2] as usize];
+                        if !covered(nbu) {
+                            let other = cid(nb);
+                            if other != NONE {
+                                rows[row].push((row, hc));
+                                rows[row].push((other, -hc));
+                            }
+                            continue;
+                        }
+                        at_interface[row] = true;
+                        let kids = children(nbu);
+                        // The four children on the face (normal index nearest C).
+                        let near: Vec<[usize; 3]> = kids
+                            .iter()
+                            .copied()
+                            .filter(|f| {
+                                let want = 2 * (nbu[d] - lo[d]) + usize::from(s < 0);
+                                f[d] == want
+                            })
+                            .collect();
+                        match spec.interface {
+                            Interface::Direct => {
+                                let coef = hf / 1.5;
+                                for f in near {
+                                    rows[row].push((row, coef));
+                                    rows[row].push((fid(f), -coef));
+                                }
+                            }
+                            Interface::Octree => {
+                                rows[row].push((row, hc));
+                                for f in kids {
+                                    rows[row].push((fid(f), -hc / 8.0));
+                                }
+                            }
+                            Interface::Quadratic => {
+                                for f in near {
+                                    // Out of the fine cell through its face −s.
+                                    rows[row].push((fid(f), -hf));
+                                    for (u, w) in ghost(f, d, -s) {
+                                        rows[row].push((u, hf * w));
+                                    }
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+        // Fine rows.
+        for k in 0..fg.nz {
+            for j in 0..fg.ny {
+                for i in 0..fg.nx {
+                    let row = fine_id[fg.cell(k, j, i)];
+                    if row == NONE {
+                        continue;
+                    }
+                    let f = [i, j, k];
+                    let x = centre[row];
+                    rhs[row] += volume[row] * (spec.source)(x[0], x[1], x[2]);
+                    for (dir, &(d, s)) in DIRS.iter().enumerate() {
+                        let nb = step(f, d, s);
+                        if nb[d] >= 0 && nb[d] < fdims[d] as isize {
+                            let nbu = [nb[0] as usize, nb[1] as usize, nb[2] as usize];
+                            let a = aperture(fcut, fg, f, dir);
+                            let other = fid(nbu);
+                            if a > 0.0 && other != NONE {
+                                rows[row].push((row, a * hf));
+                                rows[row].push((other, -a * hf));
+                            }
+                            continue;
+                        }
+                        at_interface[row] = true;
+                        assert!(
+                            aperture(fcut, fg, f, dir) > 1.0 - 1e-12,
+                            "the body must not reach the interface"
+                        );
+                        let mut gc = [
+                            (2 * lo[0] + f[0]) as isize,
+                            (2 * lo[1] + f[1]) as isize,
+                            (2 * lo[2] + f[2]) as isize,
+                        ];
+                        gc[d] += s;
+                        let cc = [
+                            gc[0].div_euclid(2),
+                            gc[1].div_euclid(2),
+                            gc[2].div_euclid(2),
+                        ];
+                        let c_row = cid(cc);
+                        match spec.interface {
+                            Interface::Direct => {
+                                let coef = hf / 1.5;
+                                rows[row].push((row, coef));
+                                rows[row].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));
+                                }
+                                rows[row].push((c_row, -coef));
+                            }
+                            Interface::Quadratic => {
+                                rows[row].push((row, hf));
+                                for (u, w) in ghost(f, d, s) {
+                                    rows[row].push((u, -hf * w));
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+        let a = Csr::from_rows(rows);
+        let coarse_problem = Self::coarse_operator(cg, ccut);
+        Self {
+            coarse: cg,
+            fine: fg,
+            lo,
+            hi,
+            interface: spec.interface,
+            coarse_id,
+            fine_id,
+            n_coarse,
+            parent,
+            centre,
+            volume,
+            at_interface,
+            cut: cut_flag,
+            a,
+            rhs,
+            coarse_problem,
+        }
+    }
+
+    /// The uniform coarse operator of the same problem (the coarse cut
+    /// under the patch), Dirichlet outside: the FAC coarse correction.
+    fn coarse_operator(cg: Grid, ccut: Option<&CutGeometry>) -> Problem {
+        let h = cg.dx;
+        let mut p = Problem::new(cg.nx, cg.ny, cg.nz);
+        let dims = [cg.nx, cg.ny, cg.nz];
+        for k in 0..cg.nz {
+            for j in 0..cg.ny {
+                for i in 0..cg.nx {
+                    let idx = cg.cell(k, j, i);
+                    p.active[idx] = active(ccut, cg, [i, j, k]);
+                }
+            }
+        }
+        for k in 0..cg.nz {
+            for j in 0..cg.ny {
+                for i in 0..cg.nx {
+                    let idx = cg.cell(k, j, i);
+                    if !p.active[idx] {
+                        continue;
+                    }
+                    let c = [i, j, k];
+                    for (dir, &(d, s)) in DIRS.iter().enumerate() {
+                        let mut nb = [i as isize, j as isize, k as isize];
+                        nb[d] += s;
+                        if nb[d] < 0 || nb[d] >= dims[d] as isize {
+                            p.extra_diag[idx] += 2.0 * h;
+                            continue;
+                        }
+                        let other = cg.cell(nb[2] as usize, nb[1] as usize, nb[0] as usize);
+                        let coef = if p.active[other] {
+                            aperture(ccut, cg, c, dir) * h
+                        } else {
+                            0.0
+                        };
+                        match dir {
+                            0 => p.aw[idx] = coef,
+                            1 => p.ae[idx] = coef,
+                            2 => p.as_[idx] = coef,
+                            3 => p.an[idx] = coef,
+                            4 => p.ab[idx] = coef,
+                            _ => p.at[idx] = coef,
+                        }
+                    }
+                }
+            }
+        }
+        p
+    }
+
+    #[must_use]
+    pub fn unknowns(&self) -> usize {
+        self.centre.len()
+    }
+}
+
+/// A uniform grid's operator and right-hand side in the same form (the
+/// reference / cost comparator), optionally cut by `body`; Dirichlet
+/// outside. Returns the problem and the per-cell fluid fraction.
+#[must_use]
+pub fn uniform_problem(
+    g: Grid,
+    body: Option<&Sdf>,
+    source: &(dyn Fn(f64, f64, f64) -> f64 + Sync),
+    dirichlet: &(dyn Fn(f64, f64, f64) -> f64 + Sync),
+) -> (Problem, Vec<f64>) {
+    let cut = body.map(|b| CutGeometry::build(&body_at(b, [0.0; 3]), g, 0.0));
+    let mut p = Composite::coarse_operator(g, cut.as_ref());
+    let h = g.dx;
+    let dims = [g.nx, g.ny, g.nz];
+    let mut frac = vec![0.0; g.cells()];
+    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 !p.active[idx] {
+                    continue;
+                }
+                let fr = fraction(cut.as_ref(), g, [i, j, k]);
+                frac[idx] = fr;
+                let x = [
+                    (i as f64 + 0.5) * h,
+                    (j as f64 + 0.5) * h,
+                    (k as f64 + 0.5) * h,
+                ];
+                p.rhs[idx] = fr * h * h * h * source(x[0], x[1], x[2]);
+                for &(d, s) in &DIRS {
+                    let c = [i, j, k][d] as isize + s;
+                    if c < 0 || c >= dims[d] as isize {
+                        let mut fx = x;
+                        fx[d] += s as f64 * 0.5 * h;
+                        p.rhs[idx] += 2.0 * h * dirichlet(fx[0], fx[1], fx[2]);
+                    }
+                }
+            }
+        }
+    }
+    (p, frac)
+}
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
new file mode 100644
index 0000000..6e49293
--- /dev/null
+++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/composite/mod.rs
@@ -0,0 +1,103 @@
+//! R7 phase 1 (omni-cortex roadmap R7, A3-ii): the composite pressure
+//! Poisson problem on a uniform coarse grid with ONE nested ratio-2 box
+//! patch, cell-centred, in the integrated (finite-volume) form of
+//! [`super::poisson::Problem`]: `Σ_faces c_f (p_i − p_nb) = rhs_i`.
+//!
+//! A host prototype, not wired into the solver (nothing in the step calls
+//! it): it exists to measure what the coarse–fine interface costs in order
+//! and in work before the patch is carried into the predictor and onto the
+//! device. Three interface fluxes are built so the choice is a measurement:
+//!
+//! - [`Interface::Direct`]: each fine sub-face couples its fine cell to the
+//!   coarse cell by a two-point flux over the centre distance `1.5 h_f`
+//!   (symmetric, conservative, not consistent tangentially);
+//! - [`Interface::Octree`]: 2604.18886 eq. 13 — the coarse face's flux
+//!   `(p_C − mean of the covered coarse cell's 8 children) / h_c`, shared
+//!   equally by the 4 fine sub-faces (conservative, non-symmetric);
+//! - [`Interface::Quadratic`]: the fine ghost by quadratic interpolation
+//!   (tangential quadratic with the cross term on the coarse layer, then
+//!   quadratic in the normal direction through two fine cells), the coarse
+//!   face's flux = the sum of its four fine fluxes (refluxing; conservative,
+//!   non-symmetric) — the AMReX / Martin–Cartwright construction.
+//!
+//! Layout: coarse cells `(k, j, i)` of [`Grid`]; the patch covers the
+//! coarse box `lo ≤ (i, j, k) < hi` and is refined by 2 in every direction.
+//! Unknowns: the uncovered active coarse cells first, then the active fine
+//! cells. A body (φ > 0 fluid) may cut the fine level; it must not reach
+//! the interface or the coarse uncovered cells.
+
+mod assemble;
+mod solve;
+
+pub use assemble::{Composite, CompositeSpec, Interface, Sdf, uniform_problem};
+pub use solve::{CompositeSolve, SolveStats, solve_bicgstab};
+
+/// Compressed sparse rows.
+#[derive(Debug, Clone, Default)]
+pub struct Csr {
+    pub start: Vec<usize>,
+    pub col: Vec<usize>,
+    pub val: Vec<f64>,
+}
+
+impl Csr {
+    /// Rows from per-row `(column, value)` lists; duplicates are summed in
+    /// the order they appear, zero entries dropped.
+    #[must_use]
+    pub fn from_rows(rows: Vec<Vec<(usize, f64)>>) -> Self {
+        let mut start = Vec::with_capacity(rows.len() + 1);
+        let mut col = Vec::new();
+        let mut val = Vec::new();
+        start.push(0);
+        for mut row in rows {
+            row.sort_by_key(|e| e.0);
+            let mut last: Option<usize> = None;
+            for (c, v) in row {
+                if last == Some(c) {
+                    *val.last_mut().expect("entry") += v;
+                } else {
+                    col.push(c);
+                    val.push(v);
+                    last = Some(c);
+                }
+            }
+            start.push(col.len());
+        }
+        Self { start, col, val }
+    }
+
+    #[must_use]
+    pub fn rows(&self) -> usize {
+        self.start.len() - 1
+    }
+
+    pub fn apply(&self, x: &[f64], y: &mut [f64]) {
+        use rayon::prelude::*;
+        y.par_iter_mut().enumerate().for_each(|(r, yr)| {
+            let mut s = 0.0;
+            for e in self.start[r]..self.start[r + 1] {
+                s += self.val[e] * x[self.col[e]];
+            }
+            *yr = s;
+        });
+    }
+
+    /// The largest `|a_ij − a_ji|` relative to the largest |a_ij|.
+    #[must_use]
+    pub fn asymmetry(&self) -> f64 {
+        let mut map = std::collections::HashMap::with_capacity(self.val.len());
+        let mut amax: f64 = 0.0;
+        for r in 0..self.rows() {
+            for e in self.start[r]..self.start[r + 1] {
+                map.insert((r, self.col[e]), self.val[e]);
+                amax = amax.max(self.val[e].abs());
+            }
+        }
+        let mut d: f64 = 0.0;
+        for (&(r, c), &v) in &map {
+            let t = map.get(&(c, r)).copied().unwrap_or(0.0);
+            d = d.max((v - t).abs());
+        }
+        d / amax.max(f64::MIN_POSITIVE)
+    }
+}
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
new file mode 100644
index 0000000..5bfa65b
--- /dev/null
+++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/composite/solve.rs
@@ -0,0 +1,195 @@
+//! The composite solve: BiCGStab (the Quadratic and Octree operators are
+//! not symmetric) right-preconditioned by one FAC cycle (McCormick's fast
+//! adaptive composite grid method in correction form):
+//!
+//! 1. `ν` forward Gauss–Seidel sweeps of the composite operator from zero;
+//! 2. the composite residual, restricted to the whole coarse grid — an
+//!    uncovered coarse cell keeps its own, a covered coarse cell receives
+//!    the SUM of its children's (integrated form);
+//! 3. one V-cycle of the uniform coarse operator (embedded3's aggregation
+//!    [`Hierarchy`], unchanged) on that residual;
+//! 4. the correction added to the coarse unknowns and injected (piecewise
+//!    constant) into the fine ones;
+//! 5. `ν` backward Gauss–Seidel sweeps.
+//!
+//! Correction form over the COMPOSITE residual is the linear counterpart of
+//! 2604.18886's FAS right-hand side `b = β R r + A_c u*` (Algorithm 4): the
+//! interface flux enters the coarse problem exactly once, through the
+//! composite residual, so nothing is double-counted at the T-junctions.
+
+use super::Composite;
+use crate::solvers::incompressible::MultigridParameters;
+use crate::solvers::incompressible::embedded3::poisson::Hierarchy;
+
+/// The prepared preconditioner.
+pub struct CompositeSolve {
+    hier: Hierarchy<f64>,
+    diag: Vec<f64>,
+    sweeps: usize,
+    rc: Vec<f64>,
+    ec: Vec<f64>,
+    res: Vec<f64>,
+}
+
+/// Outcome of [`solve_bicgstab`].
+#[derive(Debug, Clone, Copy)]
+pub struct SolveStats {
+    pub iterations: usize,
+    /// `‖b − A x‖₂ / ‖b‖₂` at exit (recomputed).
+    pub rel_residual: f64,
+    pub converged: bool,
+    pub setup_s: f64,
+    pub solve_s: f64,
+}
+
+impl CompositeSolve {
+    #[must_use]
+    pub fn new(c: &Composite, sweeps: usize) -> Self {
+        let params = MultigridParameters::default();
+        let hier = Hierarchy::<f64>::build(&c.coarse_problem, &params);
+        let n = c.unknowns();
+        let mut diag = vec![0.0; n];
+        for (r, d) in diag.iter_mut().enumerate() {
+            for e in c.a.start[r]..c.a.start[r + 1] {
+                if c.a.col[e] == r {
+                    *d = c.a.val[e];
+                }
+            }
+            assert!(*d > 0.0, "row {r} has no positive diagonal");
+        }
+        let nc = c.coarse.cells();
+        Self {
+            hier,
+            diag,
+            sweeps: sweeps.max(1),
+            rc: vec![0.0; nc],
+            ec: vec![0.0; nc],
+            res: vec![0.0; n],
+        }
+    }
+
+    fn gs_row(&self, c: &Composite, b: &[f64], x: &mut [f64], r: usize) {
+        let mut s = b[r];
+        for e in c.a.start[r]..c.a.start[r + 1] {
+            let col = c.a.col[e];
+            if col != r {
+                s -= c.a.val[e] * x[col];
+            }
+        }
+        x[r] = s / self.diag[r];
+    }
+
+    /// `z = M⁻¹ r` (one FAC cycle from zero).
+    pub fn precondition(&mut self, c: &Composite, r: &[f64], z: &mut [f64]) {
+        let n = c.unknowns();
+        z.iter_mut().for_each(|v| *v = 0.0);
+        for _ in 0..self.sweeps {
+            for row in 0..n {
+                self.gs_row(c, r, z, row);
+            }
+        }
+        c.a.apply(z, &mut self.res);
+        self.rc.iter_mut().for_each(|v| *v = 0.0);
+        for u in 0..n {
+            self.rc[c.parent[u]] += r[u] - self.res[u];
+        }
+        for (idx, v) in self.rc.iter_mut().enumerate() {
+            if !c.coarse_problem.active[idx] {
+                *v = 0.0;
+            }
+        }
+        self.hier.apply_preconditioner(&self.rc, &mut self.ec);
+        for u in 0..n {
+            let p = c.parent[u];
+            if c.coarse_problem.active[p] {
+                z[u] += self.ec[p];
+            }
+        }
+        for _ in 0..self.sweeps {
+            for row in (0..n).rev() {
+                self.gs_row(c, r, z, row);
+            }
+        }
+    }
+}
+
+fn dot(a: &[f64], b: &[f64]) -> f64 {
+    a.iter().zip(b).map(|(x, y)| x * y).sum()
+}
+
+/// Solve `A x = rhs` to `‖r‖₂ ≤ tol ‖rhs‖₂` from the given `x`.
+pub fn solve_bicgstab(
+    c: &Composite,
+    x: &mut [f64],
+    tol: f64,
+    max_it: usize,
+    sweeps: usize,
+) -> SolveStats {
+    let t0 = std::time::Instant::now();
+    let mut pre = CompositeSolve::new(c, sweeps);
+    let setup_s = t0.elapsed().as_secs_f64();
+    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);
+    for i in 0..n {
+        r[i] = b[i] - r[i];
+    }
+    let rhat = r.clone();
+    let (mut rho, mut alpha, mut omega) = (1.0, 1.0, 1.0);
+    let mut v = vec![0.0; n];
+    let mut p = vec![0.0; n];
+    let mut ph = vec![0.0; n];
+    let mut sh = vec![0.0; n];
+    let mut s = vec![0.0; n];
+    let mut t = vec![0.0; n];
+    let mut it = 0;
+    let mut rel = dot(&r, &r).sqrt() / bn;
+    while rel > tol && it < max_it {
+        it += 1;
+        let rho_n = dot(&rhat, &r);
+        let beta = (rho_n / rho) * (alpha / omega);
+        for i in 0..n {
+            p[i] = r[i] + beta * (p[i] - omega * v[i]);
+        }
+        pre.precondition(c, &p, &mut ph);
+        c.a.apply(&ph, &mut v);
+        alpha = rho_n / dot(&rhat, &v);
+        for i in 0..n {
+            s[i] = r[i] - alpha * v[i];
+        }
+        if dot(&s, &s).sqrt() / bn <= tol {
+            for i in 0..n {
+                x[i] += alpha * ph[i];
+            }
+            r.copy_from_slice(&s);
+            break;
+        }
+        pre.precondition(c, &s, &mut sh);
+        c.a.apply(&sh, &mut t);
+        omega = dot(&t, &s) / dot(&t, &t);
+        for i in 0..n {
+            x[i] += alpha * ph[i] + omega * sh[i];
+            r[i] = s[i] - omega * t[i];
+        }
+        rho = rho_n;
+        rel = dot(&r, &r).sqrt() / bn;
+    }
+    // The true residual.
+    c.a.apply(x, &mut r);
+    let mut rr = 0.0;
+    for i in 0..n {
+        let d = b[i] - r[i];
+        rr += d * d;
+    }
+    let rel = rr.sqrt() / bn;
+    SolveStats {
+        iterations: it,
+        rel_residual: rel,
+        converged: rel <= 10.0 * tol,
+        setup_s,
+        solve_s: t1.elapsed().as_secs_f64(),
+    }
+}
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 5f55ba2..0e539b1 100644
--- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs
+++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs
@@ -8,6 +8,7 @@
 
 pub mod body;
 pub mod closure;
+pub mod composite;
 pub mod cut;
 pub mod cutwall;
 pub mod exchange;
diff --git a/crates/specialized/rtx-cfd/tests/embedded3_composite_poisson.rs b/crates/specialized/rtx-cfd/tests/embedded3_composite_poisson.rs
new file mode 100644
index 0000000..677873d
--- /dev/null
+++ b/crates/specialized/rtx-cfd/tests/embedded3_composite_poisson.rs
@@ -0,0 +1,733 @@
+//! R7 phase 1: the composite Poisson problem with one nested ratio-2 patch
+//! (`embedded3::composite`, a host prototype nothing else calls).
+//!
+//! - P1 (`composite_mms_ladder`, ignored): manufactured solution on the unit
+//!   cube, patch = the middle half, n = 16/32/64 coarse; errors split into
+//!   the interface cells (owners of a coarse–fine face), the rest of the
+//!   patch, and the rest of the coarse grid, for the three interface fluxes
+//!   and the uniform coarse / fine grids;
+//! - P2 (`composite_cut_sphere`, ignored): a Neumann sphere inside the patch
+//!   cut by embedded3's `CutGeometry` on the fine level, against the
+//!   uniformly fine grid with the same cut;
+//! - P3 (`composite_cost`, ignored): unknowns and wall time against the
+//!   uniformly fine grid solved by embedded3's production PCG.
+//!
+//! Output CSVs go to `$R7_OUT` when set. The default (non-ignored) test is a
+//! seconds-scale smoke of the assembly and the solve.
+
+use rtx_cfd::solvers::incompressible::MultigridParameters;
+use rtx_cfd::solvers::incompressible::embedded3::Grid;
+use rtx_cfd::solvers::incompressible::embedded3::composite::{
+    Composite, CompositeSpec, Interface, Sdf, solve_bicgstab, uniform_problem,
+};
+use rtx_cfd::solvers::incompressible::embedded3::poisson::solve_pcg;
+use std::f64::consts::PI;
+use std::io::Write;
+use std::sync::Arc;
+
+fn mms_u(x: f64, y: f64, z: f64) -> f64 {
+    (1.3 * PI * x + 0.2).cos() * (1.1 * PI * y + 0.4).cos() * (0.9 * PI * z + 0.6).cos()
+}
+
+fn mms_f(x: f64, y: f64, z: f64) -> f64 {
+    PI * PI * (1.69 + 1.21 + 0.81) * mms_u(x, y, z)
+}
+
+/// A localized feature inside the patch on top of the smooth field: the
+/// case local refinement is for (a Gaussian of width 0.06 at the centre).
+const BUMP_S: f64 = 0.06;
+
+fn bump_u(x: f64, y: f64, z: f64) -> f64 {
+    let r2 = (x - 0.5).powi(2) + (y - 0.5).powi(2) + (z - 0.5).powi(2);
+    mms_u(x, y, z) + (-r2 / (BUMP_S * BUMP_S)).exp()
+}
+
+fn bump_f(x: f64, y: f64, z: f64) -> f64 {
+    let r2 = (x - 0.5).powi(2) + (y - 0.5).powi(2) + (z - 0.5).powi(2);
+    let s2 = BUMP_S * BUMP_S;
+    mms_f(x, y, z) - (4.0 * r2 / (s2 * s2) - 6.0 / s2) * (-r2 / s2).exp()
+}
+
+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 sph_r(x: f64, y: f64, z: f64) -> f64 {
+    ((x - SPH_C[0]).powi(2) + (y - SPH_C[1]).powi(2) + (z - SPH_C[2]).powi(2)).sqrt()
+}
+
+/// `cos(a (r − R))`: zero normal derivative on the sphere.
+fn sph_u(x: f64, y: f64, z: f64) -> f64 {
+    (SPH_A * (sph_r(x, y, z) - SPH_R)).cos()
+}
+
+fn sph_f(x: f64, y: f64, z: f64) -> f64 {
+    let r = sph_r(x, y, z);
+    let q = SPH_A * (r - SPH_R);
+    SPH_A * SPH_A * q.cos() + 2.0 * SPH_A * q.sin() / r
+}
+
+fn sphere() -> Sdf {
+    Arc::new(|x, y, z| sph_r(x, y, z) - SPH_R)
+}
+
+/// Volume-weighted RMS and max of the error over a subset.
+#[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 {
+        if self.vol > 0.0 {
+            (self.s2 / self.vol).sqrt()
+        } else {
+            f64::NAN
+        }
+    }
+}
+
+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()
+}
+
+struct CompositeRun {
+    unknowns: usize,
+    fine_unknowns: usize,
+    iterations: usize,
+    rel: f64,
+    setup_s: f64,
+    solve_s: f64,
+    asym: f64,
+    /// interface, patch interior (fine, not interface), coarse (not
+    /// interface), cut cells, all.
+    err: [Err; 5],
+}
+
+fn run_composite(
+    n: usize,
+    lo: usize,
+    hi: usize,
+    iface: Interface,
+    body: Option<Sdf>,
+    exact: fn(f64, f64, f64) -> f64,
+    source: fn(f64, f64, f64) -> f64,
+    measure_asym: bool,
+) -> CompositeRun {
+    let g = Grid::cubic(n, n, n, 1.0 / n as f64);
+    let spec = CompositeSpec {
+        coarse: g,
+        lo: [lo; 3],
+        hi: [hi; 3],
+        interface: iface,
+        body,
+        source: &source,
+        dirichlet: &exact,
+    };
+    let t = std::time::Instant::now();
+    let c = Composite::build(&spec);
+    let build_s = t.elapsed().as_secs_f64();
+    let mut x = vec![0.0; c.unknowns()];
+    let st = solve_bicgstab(&c, &mut x, 1e-11, 400, 2);
+    assert!(st.converged, "composite solve did not converge: {st:?}");
+    let mut err = [Err::default(); 5];
+    for u in 0..c.unknowns() {
+        let p = c.centre[u];
+        let e = x[u] - exact(p[0], p[1], p[2]);
+        let v = c.volume[u];
+        let fine = u >= c.n_coarse;
+        let class = if c.at_interface[u] {
+            0
+        } else if fine {
+            1
+        } else {
+            2
+        };
+        err[class].add(e, v);
+        if c.cut[u] {
+            err[3].add(e, v);
+        }
+        err[4].add(e, v);
+    }
+    CompositeRun {
+        unknowns: c.unknowns(),
+        fine_unknowns: c.unknowns() - c.n_coarse,
+        iterations: st.iterations,
+        rel: st.rel_residual,
+        setup_s: build_s + st.setup_s,
+        solve_s: st.solve_s,
+        asym: if measure_asym {
+            c.a.asymmetry()
+        } else {
+            f64::NAN
+        },
+        err,
+    }
+}
+
+struct UniformRun {
+    cells: usize,
+    iterations: usize,
+    setup_s: f64,
+    solve_s: f64,
+    /// inside the patch region [lo, hi)·h_c, outside it, cut cells, all.
+    err: [Err; 4],
+}
+
+fn run_uniform(
+    n: usize,
+    region: (f64, f64),
+    body: Option<&Sdf>,
+    exact: fn(f64, f64, f64) -> f64,
+    source: fn(f64, f64, f64) -> f64,
+) -> UniformRun {
+    let g = Grid::cubic(n, n, n, 1.0 / n as f64);
+    let t = std::time::Instant::now();
+    let (prob, frac) = uniform_problem(g, body, &source, &exact);
+    let build_s = t.elapsed().as_secs_f64();
+    let mut p = vec![0.0; g.cells()];
+    let l1: f64 = prob.rhs.iter().map(|v| v.abs()).sum();
+    let params = MultigridParameters {
+        max_iterations: 2000,
+        ..MultigridParameters::default()
+    };
+    let sol = solve_pcg(&prob, &mut p, &params, 1e-11 * l1, None);
+    assert!(sol.converged, "uniform solve did not converge");
+    let h = g.dx;
+    let mut err = [Err::default(); 4];
+    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;
+                }
+                let x = [
+                    (i as f64 + 0.5) * h,
+                    (j as f64 + 0.5) * h,
+                    (k as f64 + 0.5) * h,
+                ];
+                let e = p[idx] - exact(x[0], x[1], x[2]);
+                let v = frac[idx] * h * h * h;
+                let inside = x.iter().all(|&c| c > region.0 && c < region.1);
+                err[usize::from(!inside)].add(e, v);
+                if frac[idx] < 1.0 - 1e-9 {
+                    err[2].add(e, v);
+                }
+                err[3].add(e, v);
+            }
+        }
+    }
+    UniformRun {
+        cells: prob.active.iter().filter(|&&a| a).count(),
+        iterations: sol.iterations,
+        setup_s: build_s + sol.setup_ns as f64 * 1e-9,
+        solve_s: sol.iterate_ns as f64 * 1e-9,
+        err,
+    }
+}
+
+fn order(a: f64, b: f64) -> f64 {
+    (a / b).log2()
+}
+
+#[test]
+fn composite_smoke() {
+    for iface in [Interface::Direct, Interface::Octree, Interface::Quadratic] {
+        let r = run_composite(12, 3, 9, iface, None, mms_u, mms_f, true);
+        eprintln!(
+            "{iface:?}: {} unknowns, {} it, rel {:.2e}, asym {:.2e}, L2 iface {:.3e} all {:.3e}",
+            r.unknowns,
+            r.iterations,
+            r.rel,
+            r.asym,
+            r.err[0].l2(),
+            r.err[4].l2()
+        );
+        assert!(r.err[4].l2() < 5e-3);
+        if iface == Interface::Direct {
+            assert!(r.asym < 1e-14);
+        }
+    }
+}
+
+#[test]
+#[ignore = "P1 ladder (minutes)"]
+fn composite_mms_ladder() {
+    ladder("p1_mms.csv", mms_u, mms_f);
+}
+
+#[test]
+#[ignore = "P1b ladder with a localized bump (minutes)"]
+fn composite_bump_ladder() {
+    ladder("p1b_bump.csv", bump_u, bump_f);
+}
+
+fn ladder(name: &str, exact: fn(f64, f64, f64) -> f64, source: fn(f64, f64, f64) -> f64) {
+    let mut csv = out_file(name);
+    if let Some(f) = csv.as_mut() {
+        writeln!(
+            f,
+            "scheme,n,unknowns,iters,rel,l2_iface,linf_iface,l2_patch,linf_patch,l2_coarse,linf_coarse,l2_all,linf_all,setup_s,solve_s"
+        )
+        .ok();
+    }
+    let ns = [16usize, 32, 64];
+    let mut table: Vec<(String, Vec<[f64; 8]>)> = Vec::new();
+    for iface in [Interface::Direct, Interface::Octree, Interface::Quadratic] {
+        let mut rows = Vec::new();
+        for &n in &ns {
+            let r = run_composite(n, n / 4, 3 * n / 4, iface, None, exact, source, false);
+            let e = r.err;
+            eprintln!(
+                "{iface:?} n {n}: {} unk, {} it, rel {:.1e}; iface L2 {:.3e} Linf {:.3e} | patch L2 {:.3e} Linf {:.3e} | coarse L2 {:.3e} Linf {:.3e} | all L2 {:.3e} | {:.2}+{:.2} s",
+                r.unknowns,
+                r.iterations,
+                r.rel,
+                e[0].l2(),
+                e[0].max,
+                e[1].l2(),
+                e[1].max,
+                e[2].l2(),
+                e[2].max,
+                e[4].l2(),
+                r.setup_s,
+                r.solve_s
+            );
+            if let Some(f) = csv.as_mut() {
+                writeln!(
+                    f,
+                    "{iface:?},{n},{},{},{:.3e},{:.6e},{:.6e},{:.6e},{:.6e},{:.6e},{:.6e},{:.6e},{:.6e},{:.3},{:.3}",
+                    r.unknowns,
+                    r.iterations,
+                    r.rel,
+                    e[0].l2(),
+                    e[0].max,
+                    e[1].l2(),
+                    e[1].max,
+                    e[2].l2(),
+                    e[2].max,
+                    e[4].l2(),
+                    e[4].max,
+                    r.setup_s,
+                    r.solve_s
+                )
+                .ok();
+            }
+            rows.push([
+                e[0].l2(),
+                e[0].max,
+                e[1].l2(),
+                e[1].max,
+                e[2].l2(),
+                e[2].max,
+                e[4].l2(),
+                e[4].max,
+            ]);
+        }
+        table.push((format!("{iface:?}"), rows));
+    }
+    for (label, n_of) in [("uniform-coarse", 1usize), ("uniform-fine", 2)] {
+        for &n in &ns {
+            let r = run_uniform(n_of * n, (0.25, 0.75), None, exact, source);
+            let e = r.err;
+            eprintln!(
+                "{label} n {}: {} cells, {} it; patch-region L2 {:.3e} Linf {:.3e} | outside L2 {:.3e} Linf {:.3e} | all L2 {:.3e} | {:.2}+{:.2} s",
+                n_of * n,
+                r.cells,
+                r.iterations,
+                e[0].l2(),
+                e[0].max,
+                e[1].l2(),
+                e[1].max,
+                e[3].l2(),
+                r.setup_s,
+                r.solve_s
+            );
+            if let Some(f) = csv.as_mut() {
+                writeln!(
+                    f,
+                    "{label},{},{},{},0,nan,nan,{:.6e},{:.6e},{:.6e},{:.6e},{:.6e},{:.6e},{:.3},{:.3}",
+                    n_of * n,
+                    r.cells,
+                    r.iterations,
+                    e[0].l2(),
+                    e[0].max,
+                    e[1].l2(),
+                    e[1].max,
+                    e[3].l2(),
+                    e[3].max,
+                    r.setup_s,
+                    r.solve_s
+                )
+                .ok();
+            }
+        }
+    }
+    eprintln!("orders (16→32, 32→64): iface L2 / Linf | patch L2 | coarse L2 | all L2 / Linf");
+    for (label, rows) in &table {
+        let o = |c: usize| (order(rows[0][c], rows[1][c]), order(rows[1][c], rows[2][c]));
+        eprintln!(
+            "{label}: iface {:.2},{:.2} / {:.2},{:.2} | patch {:.2},{:.2} | coarse {:.2},{:.2} | all {:.2},{:.2} / {:.2},{:.2}",
+            o(0).0,
+            o(0).1,
+            o(1).0,
+            o(1).1,
+            o(2).0,
+            o(2).1,
+            o(4).0,
+            o(4).1,
+            o(6).0,
+            o(6).1,
+            o(7).0,
+            o(7).1
+        );
+    }
+}
+
+#[test]
+#[ignore = "P2 cut sphere in the patch (minutes)"]
+fn composite_cut_sphere() {
+    let mut csv = out_file("p2_sphere.csv");
+    if let Some(f) = csv.as_mut() {
+        writeln!(
+            f,
+            "scheme,n,unknowns,iters,l2_iface,linf_iface,l2_patch,linf_patch,l2_coarse,l2_cut,linf_cut,l2_all,linf_all"
+        )
+        .ok();
+    }
+    let body = sphere();
+    for &n in &[16usize, 32, 64] {
+        for iface in [Interface::Direct, Interface::Octree, Interface::Quadratic] {
+            let r = run_composite(
+                n,
+                n / 4,
+                3 * n / 4,
+                iface,
+                Some(body.clone()),
+                sph_u,
+                sph_f,
+                false,
+            );
+            let e = r.err;
+            eprintln!(
+                "{iface:?} n {n}: {} unk ({} fine), {} it; iface L2 {:.3e} | patch L2 {:.3e} Linf {:.3e} | coarse L2 {:.3e} | cut L2 {:.3e} Linf {:.3e} | all L2 {:.3e} Linf {:.3e}",
+                r.unknowns,
+                r.fine_unknowns,
+                r.iterations,
+                e[0].l2(),
+                e[1].l2(),
+                e[1].max,
+                e[2].l2(),
+                e[3].l2(),
+                e[3].max,
+                e[4].l2(),
+                e[4].max
+            );
+            if let Some(f) = csv.as_mut() {
+                writeln!(
+                    f,
+                    "{iface:?},{n},{},{},{:.6e},{:.6e},{:.6e},{:.6e},{:.6e},{:.6e},{:.6e},{:.6e},{:.6e}",
+                    r.unknowns,
+                    r.iterations,
+                    e[0].l2(),
+                    e[0].max,
+                    e[1].l2(),
+                    e[1].max,
+                    e[2].l2(),
+                    e[3].l2(),
+                    e[3].max,
+                    e[4].l2(),
+                    e[4].max
+                )
+                .ok();
+            }
+        }
+        for (label, m) in [("uniform-coarse", n), ("uniform-fine", 2 * n)] {
+            let r = run_uniform(m, (0.25, 0.75), Some(&body), sph_u, sph_f);
+            let e = r.err;
+            eprintln!(
+                "{label} n {m}: {} cells, {} it; patch-region L2 {:.3e} Linf {:.3e} | outside L2 {:.3e} | cut L2 {:.3e} Linf {:.3e} | all L2 {:.3e} Linf {:.3e}",
+                r.cells,
+                r.iterations,
+                e[0].l2(),
+                e[0].max,
+                e[1].l2(),
+                e[2].l2(),
+                e[2].max,
+                e[3].l2(),
+                e[3].max
+            );
+            if let Some(f) = csv.as_mut() {
+                writeln!(
+                    f,
+                    "{label},{m},{},{},nan,nan,{:.6e},{:.6e},{:.6e},{:.6e},{:.6e},{:.6e},{:.6e}",
+                    r.cells,
+                    r.iterations,
+                    e[0].l2(),
+                    e[0].max,
+                    e[1].l2(),
+                    e[2].l2(),
+                    e[2].max,
+                    e[3].l2(),
+                    e[3].max
+                )
+                .ok();
+            }
+        }
+    }
+}
+
+#[test]
+#[ignore = "P3 cost against the uniformly fine grid (minutes)"]
+fn composite_cost() {
+    let mut csv = out_file("p3_cost.csv");
+    if let Some(f) = csv.as_mut() {
+        writeln!(
+            f,
+            "field,case,n,patch,unknowns,iters,setup_s,solve_s,l2_patch_region,linf_patch_region,l2_all"
+        )
+        .ok();
+    }
+    let n = 64usize;
+    let h = 1.0 / n as f64;
+    type Pair = (
+        &'static str,
+        fn(f64, f64, f64) -> f64,
+        fn(f64, f64, f64) -> f64,
+    );
+    let fields: [Pair; 2] = [("smooth", mms_u, mms_f), ("bump", bump_u, bump_f)];
+    for (field, exact, source) in fields {
+        for (label, lo, hi) in [("half", 16usize, 48usize), ("quarter", 24, 40)] {
+            let r = run_composite(n, lo, hi, Interface::Quadratic, None, exact, source, false);
+            // The patch region's error: interface fine cells + patch interior.
+            let mut e = Err::default();
+            let c = {
+                let g = Grid::cubic(n, n, n, h);
+                let spec = CompositeSpec {
+                    coarse: g,
+                    lo: [lo; 3],
+                    hi: [hi; 3],
+                    interface: Interface::Quadratic,
+                    body: None,
+                    source: &source,
+                    dirichlet: &exact,
+                };
+                Composite::build(&spec)
+            };
+            let mut x = vec![0.0; c.unknowns()];
+            let _ = solve_bicgstab(&c, &mut x, 1e-11, 400, 2);
+            for u in c.n_coarse..c.unknowns() {
+                let p = c.centre[u];
+                e.add(x[u] - exact(p[0], p[1], p[2]), c.volume[u]);
+            }
+            eprintln!(
+                "{field} composite {label} (n {n}, patch {lo}..{hi}): {} unknowns, {} it, setup {:.2} s, solve {:.2} s, patch L2 {:.3e} Linf {:.3e}, all L2 {:.3e}",
+                r.unknowns,
+                r.iterations,
+                r.setup_s,
+                r.solve_s,
+                e.l2(),
+                e.max,
+                r.err[4].l2()
+            );
+            if let Some(f) = csv.as_mut() {
+                writeln!(
+                    f,
+                    "{field},composite-{label},{n},{lo}..{hi},{},{},{:.3},{:.3},{:.6e},{:.6e},{:.6e}",
+                    r.unknowns,
+                    r.iterations,
+                    r.setup_s,
+                    r.solve_s,
+                    e.l2(),
+                    e.max,
+                    r.err[4].l2()
+                )
+                .ok();
+            }
+            for (ulabel, m) in [("uniform-fine", 2 * n), ("uniform-coarse", n)] {
+                let r = run_uniform(m, (lo as f64 * h, hi as f64 * h), None, exact, source);
+                eprintln!(
+                    "{field} {ulabel} n {m} (region {lo}..{hi}): {} cells, {} it, setup {:.2} s, solve {:.2} s, region L2 {:.3e} Linf {:.3e}, all L2 {:.3e}",
+                    r.cells,
+                    r.iterations,
+                    r.setup_s,
+                    r.solve_s,
+                    r.err[0].l2(),
+                    r.err[0].max,
+                    r.err[3].l2()
+                );
+                if let Some(f) = csv.as_mut() {
+                    writeln!(
+                        f,
+                        "{field},{ulabel},{m},{lo}..{hi},{},{},{:.3},{:.3},{:.6e},{:.6e},{:.6e}",
+                        r.cells,
+                        r.iterations,
+                        r.setup_s,
+                        r.solve_s,
+                        r.err[0].l2(),
+                        r.err[0].max,
+                        r.err[3].l2()
+                    )
+                    .ok();
+                }
+            }
+        }
+    }
+}
+
+fn quad_u(x: f64, y: f64, z: f64) -> f64 {
+    x * x + 2.0 * y * y - 1.5 * z * z + 0.7 * x * y - 0.4 * y * z + 0.9 * x * z + 0.3 * x
+}
+
+fn quad_f(_x: f64, _y: f64, _z: f64) -> f64 {
+    -(2.0 + 4.0 - 3.0)
+}
+
+/// Local consistency: the exact quadratic's residual on the rows away from
+/// the outer boundary, per interface flux (the Quadratic ghost is exact for
+/// quadratics, so its interface rows must be at round-off).
+#[test]
+fn composite_quadratic_consistency() {
+    let n = 12;
+    let g = Grid::cubic(n, n, n, 1.0 / n as f64);
+    for iface in [Interface::Direct, Interface::Octree, Interface::Quadratic] {
+        let spec = CompositeSpec {
+            coarse: g,
+            lo: [3; 3],
+            hi: [9; 3],
+            interface: iface,
+            body: None,
+            source: &quad_f,
+            dirichlet: &quad_u,
+        };
+        let c = Composite::build(&spec);
+        let u: Vec<f64> = c.centre.iter().map(|p| quad_u(p[0], p[1], p[2])).collect();
+        let mut au = vec![0.0; u.len()];
+        c.a.apply(&u, &mut au);
+        let h = g.dx;
+        let (mut m_if, mut m_in) = (0.0f64, 0.0f64);
+        for r in 0..u.len() {
+            let p = c.centre[r];
+            if p.iter().any(|&x| x < h || x > 1.0 - h) {
+                continue;
+            }
+            // Relative to the row's volume source |f| V.
+            let res = (c.rhs[r] - au[r]).abs() / (3.0 * c.volume[r]);
+            if c.at_interface[r] {
+                m_if = m_if.max(res);
+            } else {
+                m_in = m_in.max(res);
+            }
+        }
+        eprintln!("{iface:?}: max |truncation| / |f V|: interface {m_if:.3e}, interior {m_in:.3e}");
+        assert!(m_in < 1e-10);
+        if iface == Interface::Quadratic {
+            assert!(m_if < 1e-10, "quadratic ghost not exact on a quadratic");
+        }
+    }
+}
+
+/// P2, direct: the composite solution on the patch's fine cells against
+/// the uniformly fine solution on the SAME cells with the same cut (the
+/// fine cut geometry of the patch vs the uniform grid's, compared too).
+#[test]
+#[ignore = "P2 composite vs uniform fine on the patch cells (minutes)"]
+fn composite_sphere_vs_fine() {
+    let mut csv = out_file("p2_vs_fine.csv");
+    if let Some(f) = csv.as_mut() {
+        writeln!(
+            f,
+            "scheme,n,max_dfrac,l2_diff_patch,linf_diff_patch,l2_diff_cut,linf_diff_cut,l2_err_fine_patch,linf_err_fine_patch"
+        )
+        .ok();
+    }
+    let body = sphere();
+    for &n in &[16usize, 32, 64] {
+        let m = 2 * n;
+        let gf = Grid::cubic(m, m, m, 1.0 / m as f64);
+        let (prob, frac) = uniform_problem(gf, Some(&body), &sph_f, &sph_u);
+        let mut pf = vec![0.0; gf.cells()];
+        let l1: f64 = prob.rhs.iter().map(|v| v.abs()).sum();
+        let params = MultigridParameters {
+            max_iterations: 2000,
+            ..MultigridParameters::default()
+        };
+        assert!(solve_pcg(&prob, &mut pf, &params, 1e-12 * l1, None).converged);
+        for iface in [Interface::Octree, Interface::Quadratic] {
+            let g = Grid::cubic(n, n, n, 1.0 / n as f64);
+            let spec = CompositeSpec {
+                coarse: g,
+                lo: [n / 4; 3],
+                hi: [3 * n / 4; 3],
+                interface: iface,
+                body: Some(body.clone()),
+                source: &sph_f,
+                dirichlet: &sph_u,
+            };
+            let c = Composite::build(&spec);
+            let mut x = vec![0.0; c.unknowns()];
+            assert!(solve_bicgstab(&c, &mut x, 1e-12, 400, 2).converged);
+            let (mut d_all, mut d_cut, mut e_fine) =
+                (Err::default(), Err::default(), Err::default());
+            let mut max_dfrac: f64 = 0.0;
+            let off = n / 2; // the patch's first fine index on the uniform fine grid
+            for k in 0..c.fine.nz {
+                for j in 0..c.fine.ny {
+                    for i in 0..c.fine.nx {
+                        let u = c.fine_id[c.fine.cell(k, j, i)];
+                        let idx = gf.cell(k + off, j + off, i + off);
+                        if u == usize::MAX {
+                            assert!(!prob.active[idx] || frac[idx] < 1e-12, "activity differs");
+                            continue;
+                        }
+                        let hf = c.fine.dx;
+                        let fr = c.volume[u] / (hf * hf * hf);
+                        max_dfrac = max_dfrac.max((fr - frac[idx]).abs());
+                        let d = x[u] - pf[idx];
+                        d_all.add(d, c.volume[u]);
+                        if c.cut[u] {
+                            d_cut.add(d, c.volume[u]);
+                        }
+                        let p = c.centre[u];
+                        e_fine.add(pf[idx] - sph_u(p[0], p[1], p[2]), c.volume[u]);
+                    }
+                }
+            }
+            eprintln!(
+                "{iface:?} n {n} vs uniform {m}: max |Δfrac| {max_dfrac:.1e}; composite − fine on the patch L2 {:.3e} Linf {:.3e}; on cut cells ({}) L2 {:.3e} Linf {:.3e}; the fine grid's own error there L2 {:.3e} Linf {:.3e}",
+                d_all.l2(),
+                d_all.max,
+                d_cut.count,
+                d_cut.l2(),
+                d_cut.max,
+                e_fine.l2(),
+                e_fine.max
+            );
+            if let Some(f) = csv.as_mut() {
+                writeln!(
+                    f,
+                    "{iface:?},{n},{max_dfrac:.3e},{:.6e},{:.6e},{:.6e},{:.6e},{:.6e},{:.6e}",
+                    d_all.l2(),
+                    d_all.max,
+                    d_cut.l2(),
+                    d_cut.max,
+                    e_fine.l2(),
+                    e_fine.max
+                )
+                .ok();
+            }
+        }
+    }
+}