Merge r7-2a-composite-projection (R8/R7 phase 2 round 1; default-off, verified)

Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]>
2026-09-26 04:52:30 -05:00
parent c82f28ce21 cc68336f48
commit ac3a5e5eb1
5 changed files with 1560 additions and 9 deletions
@@ -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,
        }
    }

@@ -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,
+        }
+    }
+}
@@ -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)]
@@ -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);