rtx-cfd: multigrid-PCG projection — 30x faster, same answers — and the CFD1 refinement study

Falsifier 4 of the Turek–Hron geometry decision fired (the SOR projection
cost 0.09 s/step at 250x41 and an hour per run at 5 mm); this answers it.

solvers::incompressible::poisson: PoissonProblem (cell-centred five-point
SPD operator as per-cell face coefficients + Dirichlet diagonal extra +
active mask) and solve_multigrid_pcg — conjugate gradient preconditioned
by one V-cycle of geometric multigrid: aggregation by 2 per direction (odd
sizes absorbed, coarse cell active iff any child is), the Galerkin coarse
operator for piecewise-constant prolongation / summation restriction,
symmetric Gauss–Seidel smoothing, coarse correction scaled by 2 (Braess's
under-correction of unsmoothed aggregation; scalar, so the preconditioner
stays symmetric and positive on range(A)), L1 TRUE-residual stop with a
stagnation guard. Singular systems are handled per connected component of
the active cells (mean projection and level per pure-Neumann component;
the anchor's component to p[anchor] = 0). PoissonSolverKind::{Sor,
Multigrid} on PisoParameters / EmbeddedParameters; Sor is the default and
its code is byte-for-byte untouched; an unconverged multigrid solve falls
back to the SOR sweeps for that projection.

Verified (poisson/tests.rs, tests/poisson_equivalence.rs):
- PCG iterations to cut the residual 1e-8 on the closed Neumann box at
  32^2..256^2: 4, 4, 4, 4; ragged masked domains 8/8/8;
- manufactured recoveries to ~1e-14; Galerkin identity A_c v = R A P v to
  7e-15 on every level (masked, outlet column, non-uniform conductances);
  V-cycle symmetric to 1e-14; NaN-poisoned inactive cells untouched;
- two Neumann components with opposite imbalances, and a Dirichlet
  component beside an imbalanced Neumann one (review scenarios): converge,
  each component right up to its own constant;
- speed vs plain SOR at the same stop: 22.7x (128^2), 41x (256^2);
- same answers as SOR: PISO MMS 4.6e-8 relative, Taylor–Green divergence
  1.4e-9 every step, embedded-circle MMS 7e-8, no-body bit-identity with MG
  on both solvers, channel+outlet+circle 1.4e-10; CFD1 loads identical to
  four digits at 0.003 s/step vs 0.094 (30x).

CFD1 refinement study (tests/turek_hron_cfd.rs, three grids, 257 s):
h = 10 / 6.6 / 5 mm -> control-volume drag 15.6156 / 15.2829 / 15.0988 vs
14.2929 (+9.25 / +6.93 / +5.64%), apparent order 0.71, Richardson
extrapolate 14.04; surface route and lift not monotone (flag 2/3/4 cells
thick) — the test asserts the measured band at the finest grid.

Built with a 4-agent workflow (core, integration, refinement study,
adversarial review); the review found no defects and four risks, three
fixed here (per-component projection, one symmetric smoother-sweep
parameter, acting on `converged` with an SOR fallback) and one recorded
(isotropic aggregation loses grid-independence on anisotropic cells).

rtx-cfd 301 -> 318 green.

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

crates/specialized/rtx-cfd/src/solvers/incompressible/poisson.rs

@@ -0,0 +1,865 @@
+//! Cell-centred five-point Poisson problems on an active-cell mask, solved
+//! by conjugate gradient preconditioned with one V-cycle of geometric
+//! (aggregation) multigrid.
+//!
+//! The pressure-correction systems of the fixed-grid PISO (`piso.rs`) and
+//! the embedded-body PISO (`embedded.rs`) are both instances of
+//! [`PoissonProblem`]: a symmetric positive (semi-)definite M-matrix whose
+//! off-diagonal coefficients are face conductances (`dt dy/dx`, `dt dx/dy`),
+//! zero across every prescribed face, with an optional diagonal-only
+//! Dirichlet contribution (`extra_diag`) on cells next to a pressure outlet.
+//! Both solvers used point SOR with the true-residual stop; SOR contracts
+//! the smooth modes by `1 - O(h)` per sweep, so the sweep count grows with
+//! the mesh and the 256² projection of a long march dominates the step.
+//! [`solve_multigrid_pcg`] is the mesh-independent replacement.
+//!
+//! Design:
+//! - The hierarchy is a `Vec` of levels, each holding its own
+//!   [`PoissonProblem`]. Coarsening is by 2 in each direction (odd sizes:
+//!   the last coarse cell aggregates what remains), a coarse cell is active
+//!   iff any child is. Prolongation is injection of the coarse value into
+//!   every child, restriction is the sum over children, and the coarse
+//!   operator is the Galerkin product `R A P` — which for constant
+//!   prolongation is exactly "sum the fine face coefficients across each
+//!   coarse face, sum the `extra_diag` of the children"; internal fine faces
+//!   cancel. The coarse problems therefore have the same structure as the
+//!   fine one, and the symmetry of the coarse operator is inherited.
+//! - The prolongated coarse correction is scaled by [`COARSE_CORRECTION`]
+//!   `= 2`: constant prolongation under-corrects the smooth modes by exactly
+//!   that factor (see the constant's documentation), and without the scaling
+//!   the CG iteration count grows by ~10 per level.
+//! - Singular blocks are handled per connected component of the active
+//!   cells: a component without Dirichlet data has its right-hand-side mean
+//!   projected out and its level fixed on exit (anchor or mean zero) on its
+//!   own, so enclosed fluid pockets or a body that splits the domain cannot
+//!   make the iteration diverge on a globally-compatible source.
+//! - Limit: the coarsening is isotropic point-aggregation, so the
+//!   grid-independence holds for cells of aspect ratio near one. Measured
+//!   (review): aspect 2 → 6 iterations at every size, aspect 4 → 15–22,
+//!   aspect 10 → 37–91 (32² → 256²), aspect 100 → 405 at 256². Strongly
+//!   anisotropic grids need semi-coarsening or line smoothing; the flow
+//!   solvers here use square cells.
+//! - The smoother is symmetric Gauss–Seidel (forward then backward sweep),
+//!   so each level's smoother is self-adjoint in the `A` inner product and
+//!   the V-cycle with `R = Pᵀ` is a symmetric preconditioner, which keeps
+//!   CG valid. The coarsest level is solved by 50 symmetric GS sweeps.
+//! - On a pure-Neumann problem the mean over active cells is projected out
+//!   of the right-hand side and of every preconditioned residual, so CG runs
+//!   on the range of `A` and the level of `p` is fixed on exit (anchor cell
+//!   or mean zero).
+//! - An active cell whose row is empty (`ap == 0`: an isolated cell with no
+//!   correctable face and no Dirichlet contribution) has no equation: the
+//!   solver never reads or writes it and it is excluded from the residual
+//!   — the same rule the embedded SOR applied.
+//!
+//! The convergence stop is the L1 norm of the TRUE residual `b - A p` over
+//! the active cells: the CG recurrence residual is checked every iteration
+//! and, when it passes the tolerance, the true residual is recomputed and
+//! must pass too before the solve reports convergence.
+
+/// Cell-centred five-point symmetric positive (semi-)definite problem
+///
+/// ```text
+///   ap_i p_i - ae_i p_E - aw_i p_W - an_i p_N - as_i p_S = rhs_i   on active cells,
+/// ```
+///
+/// with `ap_i = ae_i + aw_i + an_i + as_i + extra_diag_i`. Row-major index
+/// `j * nx + i`. The coefficient to a neighbour is 0 when that face is
+/// prescribed or the neighbour is inactive. Symmetry is required:
+/// `ae[j*nx+i] == aw[j*nx+i+1]`, `an[j*nx+i] == as_[(j+1)*nx+i]`.
+/// `extra_diag` holds Dirichlet (outlet) contributions that have no
+/// neighbour.
+#[derive(Debug, Clone)]
+pub struct PoissonProblem {
+    /// Cells in x.
+    pub nx: usize,
+    /// Cells in y.
+    pub ny: usize,
+    /// Active (fluid) cells; inactive cells carry no equation and are never
+    /// read or written by the solver.
+    pub active: Vec<bool>,
+    /// Coefficient to the east neighbour `(j, i+1)`.
+    pub ae: Vec<f64>,
+    /// Coefficient to the west neighbour `(j, i-1)`.
+    pub aw: Vec<f64>,
+    /// Coefficient to the north neighbour `(j+1, i)`.
+    pub an: Vec<f64>,
+    /// Coefficient to the south neighbour `(j-1, i)`.
+    pub as_: Vec<f64>,
+    /// Diagonal-only (Dirichlet) contribution.
+    pub extra_diag: Vec<f64>,
+    /// Right-hand side.
+    pub rhs: Vec<f64>,
+}
+
+impl PoissonProblem {
+    /// All cells active, every coefficient and the right-hand side zero.
+    #[must_use]
+    pub fn new(nx: usize, ny: usize) -> Self {
+        let n = nx * ny;
+        Self {
+            nx,
+            ny,
+            active: vec![true; n],
+            ae: vec![0.0; n],
+            aw: vec![0.0; n],
+            an: vec![0.0; n],
+            as_: vec![0.0; n],
+            extra_diag: vec![0.0; n],
+            rhs: vec![0.0; n],
+        }
+    }
+
+    /// Row-major index of cell `(j, i)`.
+    #[inline]
+    #[must_use]
+    pub fn index(&self, j: usize, i: usize) -> usize {
+        j * self.nx + i
+    }
+
+    /// Diagonal coefficient `ap` of cell `idx`.
+    #[inline]
+    #[must_use]
+    pub fn diagonal(&self, idx: usize) -> f64 {
+        self.ae[idx] + self.aw[idx] + self.an[idx] + self.as_[idx] + self.extra_diag[idx]
+    }
+
+    /// Pure Neumann: no active cell has a Dirichlet contribution.
+    #[must_use]
+    pub fn is_singular(&self) -> bool {
+        !self
+            .active
+            .iter()
+            .zip(&self.extra_diag)
+            .any(|(&a, &d)| a && d > 0.0)
+    }
+
+    /// Sum over the active cells that carry an equation (`ap > 0`) of
+    /// `|rhs - (ap p - Σ a_nb p_nb)|`. Neighbour values are read only from
+    /// active in-range neighbours.
+    #[must_use]
+    pub fn residual_l1(&self, p: &[f64]) -> f64 {
+        let (nx, ny) = (self.nx, self.ny);
+        let mut sum = 0.0;
+        for j in 0..ny {
+            for i in 0..nx {
+                let idx = j * nx + i;
+                if !self.active[idx] {
+                    continue;
+                }
+                let ap = self.diagonal(idx);
+                if ap <= 0.0 {
+                    continue;
+                }
+                let mut nb = 0.0;
+                if i + 1 < nx && self.active[idx + 1] {
+                    nb += self.ae[idx] * p[idx + 1];
+                }
+                if i > 0 && self.active[idx - 1] {
+                    nb += self.aw[idx] * p[idx - 1];
+                }
+                if j + 1 < ny && self.active[idx + nx] {
+                    nb += self.an[idx] * p[idx + nx];
+                }
+                if j > 0 && self.active[idx - nx] {
+                    nb += self.as_[idx] * p[idx - nx];
+                }
+                sum += (self.rhs[idx] - (ap * p[idx] - nb)).abs();
+            }
+        }
+        sum
+    }
+
+    /// Checks array lengths, non-negativity of every coefficient, zero
+    /// coefficients across domain edges and towards inactive cells, and
+    /// symmetry to `1e-12` relative.
+    pub fn validate(&self) -> Result<(), String> {
+        let n = self.nx * self.ny;
+        for (name, len) in [
+            ("active", self.active.len()),
+            ("ae", self.ae.len()),
+            ("aw", self.aw.len()),
+            ("an", self.an.len()),
+            ("as_", self.as_.len()),
+            ("extra_diag", self.extra_diag.len()),
+            ("rhs", self.rhs.len()),
+        ] {
+            if len != n {
+                return Err(format!("{name}: length {len}, expected nx*ny = {n}"));
+            }
+        }
+        let (nx, ny) = (self.nx, self.ny);
+        for j in 0..ny {
+            for i in 0..nx {
+                let idx = j * nx + i;
+                for (name, v) in [
+                    ("ae", self.ae[idx]),
+                    ("aw", self.aw[idx]),
+                    ("an", self.an[idx]),
+                    ("as_", self.as_[idx]),
+                    ("extra_diag", self.extra_diag[idx]),
+                ] {
+                    if v.is_nan() || v < 0.0 {
+                        return Err(format!("{name}[{idx}] = {v} is negative or NaN"));
+                    }
+                }
+                let active = self.active[idx];
+                let check = |name: &str, coef: f64, nb_ok: bool| -> Result<(), String> {
+                    if !nb_ok && coef != 0.0 {
+                        return Err(format!(
+                            "{name}[{idx}] = {coef} across a domain edge or towards an inactive \
+                             cell (cell active = {active})"
+                        ));
+                    }
+                    Ok(())
+                };
+                check(
+                    "ae",
+                    self.ae[idx],
+                    active && i + 1 < nx && self.active[idx + 1],
+                )?;
+                check("aw", self.aw[idx], active && i > 0 && self.active[idx - 1])?;
+                check(
+                    "an",
+                    self.an[idx],
+                    active && j + 1 < ny && self.active[idx + nx],
+                )?;
+                check(
+                    "as_",
+                    self.as_[idx],
+                    active && j > 0 && self.active[idx - nx],
+                )?;
+                let symmetric = |a: f64, b: f64| (a - b).abs() <= 1e-12 * a.abs().max(b.abs());
+                if i + 1 < nx && !symmetric(self.ae[idx], self.aw[idx + 1]) {
+                    return Err(format!(
+                        "asymmetric x face: ae[{idx}] = {} vs aw[{}] = {}",
+                        self.ae[idx],
+                        idx + 1,
+                        self.aw[idx + 1]
+                    ));
+                }
+                if j + 1 < ny && !symmetric(self.an[idx], self.as_[idx + nx]) {
+                    return Err(format!(
+                        "asymmetric y face: an[{idx}] = {} vs as_[{}] = {}",
+                        self.an[idx],
+                        idx + nx,
+                        self.as_[idx + nx]
+                    ));
+                }
+            }
+        }
+        Ok(())
+    }
+}
+
+/// Multigrid preconditioner parameters.
+#[derive(Debug, Clone)]
+pub struct MultigridParameters {
+    /// Symmetric Gauss–Seidel sweeps before AND after the coarse correction
+    /// (default 2; a value of 0 is treated as 1). One count for both on
+    /// purpose: unequal pre/post counts make the V-cycle non-symmetric and
+    /// conjugate gradient invalid — measured symmetry defect 0.9 for (1, 0)
+    /// against 1e-14 for (1, 1) — so the API does not let it happen.
+    pub smoother_sweeps: usize,
+    /// Stop coarsening when the active cells are at most this many; that
+    /// level is solved by 50 symmetric GS sweeps (default 32).
+    pub coarsest_cells: usize,
+    /// CG iteration cap (default 500).
+    pub max_iterations: usize,
+}
+
+impl Default for MultigridParameters {
+    fn default() -> Self {
+        Self {
+            smoother_sweeps: 2,
+            coarsest_cells: 32,
+            max_iterations: 500,
+        }
+    }
+}
+
+/// Outcome of a [`solve_multigrid_pcg`] call. `converged == false` means
+/// the returned `p` did not reach the tolerance (iteration cap, rounding
+/// floor, or an inconsistent system) — callers must act on it, which is
+/// why the type is `#[must_use]`.
+#[derive(Debug, Clone, Copy)]
+#[must_use]
+pub struct PoissonSolution {
+    /// CG iterations performed.
+    pub iterations: usize,
+    /// L1 true residual `Σ |b - A p|` over the active cells at exit.
+    pub residual: f64,
+    /// `residual < tolerance` at exit.
+    pub converged: bool,
+}
+
+/// Which inner solver a projection uses for its pressure-correction system.
+///
+/// `Sor` is the historical point successive over-relaxation at the optimal
+/// Poisson factor; `Multigrid` is [`solve_multigrid_pcg`] on the same
+/// coefficients, the same right-hand side and the same true-residual stop.
+/// The two land on the same discrete pressure correction (to the inner
+/// tolerance); they differ only in cost, which is mesh-independent for
+/// multigrid and grows with the mesh for SOR.
+#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
+pub enum PoissonSolverKind {
+    /// Point SOR, `omega = 2 / (1 + sin(pi / max(nx, ny)))`, 2000-sweep cap.
+    #[default]
+    Sor,
+    /// Conjugate gradient preconditioned by one geometric multigrid V-cycle.
+    Multigrid,
+}
+
+/// Symmetric GS sweeps on the coarsest level.
+const COARSEST_SWEEPS: usize = 50;
+/// Over-correction factor applied to the prolongated coarse correction.
+///
+/// With piecewise-constant prolongation and summation restriction the
+/// Galerkin coarse operator of an aggregate-by-2 hierarchy is exactly twice
+/// the geometric `2h` operator (a coarse face collects `2^(d-1)` fine faces
+/// whose conductance is `h^(d-2)` each, against `(2h)^(d-2)` for the `2h`
+/// cell), while the restricted residual of a smooth error is the full
+/// `(2h)^d f`: the plain Galerkin correction is therefore half the geometric
+/// one for the smooth modes, in any dimension and for any uniform face
+/// conductance — the well-known under-correction of unsmoothed aggregation
+/// (Braess 1995). Scaling the correction by 2 restores the geometric
+/// correction for the smooth modes. The preconditioner stays symmetric
+/// (a scalar factor), and stays positive definite: the coarse correction
+/// operator `ω P B_c R A` is `A`-self-adjoint and positive, so the V-cycle
+/// error propagator `S (I − ω P B_c R A) S` has `A`-spectrum below 1 and
+/// `M⁻¹A = I − E` is positive on the range of `A`. Measured: the PCG
+/// iteration count for a 1e-8 residual reduction on the Neumann box goes
+/// from `[14, 20, 29, 40]` (32²…256², growing a level at a time) to
+/// `[4, 4, 4, 4]`, and on the ragged mask from `[20, 28, 39]` to `[8, 8, 8]`.
+const COARSE_CORRECTION: f64 = 2.0;
+/// Hard cap on the number of levels (a 1×1 coarsest is reached long before).
+const MAX_LEVELS: usize = 64;
+
+/// One level of the hierarchy: the problem, its diagonal, the list of
+/// active cells that carry an equation (row-major), the parent map into the
+/// next coarser level, and the V-cycle work vectors.
+struct Level {
+    problem: PoissonProblem,
+    /// `problem.active && ap > 0`: cells with an equation.
+    active: Vec<bool>,
+    ap: Vec<f64>,
+    /// Row-major indices of the active cells.
+    cells: Vec<usize>,
+    /// Fine index → coarse index (empty on the coarsest level).
+    coarse_of: Vec<usize>,
+}
+
+/// V-cycle work vectors of one level.
+struct Work {
+    /// Right-hand side of the residual equation on this level.
+    b: Vec<f64>,
+    /// Correction on this level.
+    x: Vec<f64>,
+    /// Residual.
+    r: Vec<f64>,
+}
+
+impl Work {
+    fn new(n: usize) -> Self {
+        Self {
+            b: vec![0.0; n],
+            x: vec![0.0; n],
+            r: vec![0.0; n],
+        }
+    }
+}
+
+impl Level {
+    fn new(mut problem: PoissonProblem) -> Self {
+        let (nx, ny) = (problem.nx, problem.ny);
+        let n = nx * ny;
+        let ap: Vec<f64> = (0..n).map(|idx| problem.diagonal(idx)).collect();
+        let active: Vec<bool> = (0..n)
+            .map(|idx| problem.active[idx] && ap[idx] > 0.0)
+            .collect();
+        // Sanitise: a coefficient towards an inactive or out-of-range
+        // neighbour is never used, so zero it — the stencil can then branch
+        // on the coefficient alone and never touches an inactive value.
+        for j in 0..ny {
+            for i in 0..nx {
+                let idx = j * nx + i;
+                if !active[idx] {
+                    problem.ae[idx] = 0.0;
+                    problem.aw[idx] = 0.0;
+                    problem.an[idx] = 0.0;
+                    problem.as_[idx] = 0.0;
+                    continue;
+                }
+                if !(i + 1 < nx && active[idx + 1]) {
+                    problem.ae[idx] = 0.0;
+                }
+                if !(i > 0 && active[idx - 1]) {
+                    problem.aw[idx] = 0.0;
+                }
+                if !(j + 1 < ny && active[idx + nx]) {
+                    problem.an[idx] = 0.0;
+                }
+                if !(j > 0 && active[idx - nx]) {
+                    problem.as_[idx] = 0.0;
+                }
+            }
+        }
+        let cells: Vec<usize> = (0..n).filter(|&idx| active[idx]).collect();
+        Self {
+            problem,
+            active,
+            ap,
+            cells,
+            coarse_of: Vec::new(),
+        }
+    }
+
+    /// `Σ a_nb x_nb` for cell `idx`; only active in-range neighbours are
+    /// read (their coefficients are the only non-zero ones after
+    /// sanitising).
+    #[inline]
+    fn neighbour_sum(&self, x: &[f64], idx: usize) -> f64 {
+        let nx = self.problem.nx;
+        let mut s = 0.0;
+        let ae = self.problem.ae[idx];
+        if ae != 0.0 {
+            s += ae * x[idx + 1];
+        }
+        let aw = self.problem.aw[idx];
+        if aw != 0.0 {
+            s += aw * x[idx - 1];
+        }
+        let an = self.problem.an[idx];
+        if an != 0.0 {
+            s += an * x[idx + nx];
+        }
+        let as_ = self.problem.as_[idx];
+        if as_ != 0.0 {
+            s += as_ * x[idx - nx];
+        }
+        s
+    }
+
+    /// `y = A x` on the active cells.
+    fn apply(&self, x: &[f64], y: &mut [f64]) {
+        for &idx in &self.cells {
+            y[idx] = self.ap[idx] * x[idx] - self.neighbour_sum(x, idx);
+        }
+    }
+
+    /// `r = b - A x` on the active cells; returns its L1 norm.
+    fn residual(&self, b: &[f64], x: &[f64], r: &mut [f64]) -> f64 {
+        let mut l1 = 0.0;
+        for &idx in &self.cells {
+            let v = b[idx] - (self.ap[idx] * x[idx] - self.neighbour_sum(x, idx));
+            r[idx] = v;
+            l1 += v.abs();
+        }
+        l1
+    }
+
+    /// One symmetric Gauss–Seidel sweep (forward then backward) on `A x = b`.
+    fn symmetric_gs(&self, b: &[f64], x: &mut [f64]) {
+        for &idx in &self.cells {
+            x[idx] = (b[idx] + self.neighbour_sum(x, idx)) / self.ap[idx];
+        }
+        for &idx in self.cells.iter().rev() {
+            x[idx] = (b[idx] + self.neighbour_sum(x, idx)) / self.ap[idx];
+        }
+    }
+
+    /// Galerkin coarsening: coarse face coefficient = sum of the fine
+    /// coefficients across that coarse face, coarse `extra_diag` = sum of
+    /// the children's. Returns the coarse problem and the parent map.
+    fn coarsen(&self) -> (PoissonProblem, Vec<usize>) {
+        let (nx, ny) = (self.problem.nx, self.problem.ny);
+        let nxc = (nx / 2).max(1);
+        let nyc = (ny / 2).max(1);
+        let cx = |i: usize| (i / 2).min(nxc - 1);
+        let cy = |j: usize| (j / 2).min(nyc - 1);
+        let mut coarse = PoissonProblem::new(nxc, nyc);
+        coarse.active.fill(false);
+        let mut coarse_of = vec![usize::MAX; nx * ny];
+        for &idx in &self.cells {
+            let (i, j) = (idx % nx, idx / nx);
+            let (ic, jc) = (cx(i), cy(j));
+            let c = jc * nxc + ic;
+            coarse_of[idx] = c;
+            coarse.active[c] = true;
+            coarse.extra_diag[c] += self.problem.extra_diag[idx];
+            // Coefficients towards inactive neighbours are already zero.
+            if self.problem.ae[idx] != 0.0 && cx(i + 1) != ic {
+                coarse.ae[c] += self.problem.ae[idx];
+            }
+            if self.problem.aw[idx] != 0.0 && cx(i - 1) != ic {
+                coarse.aw[c] += self.problem.aw[idx];
+            }
+            if self.problem.an[idx] != 0.0 && cy(j + 1) != jc {
+                coarse.an[c] += self.problem.an[idx];
+            }
+            if self.problem.as_[idx] != 0.0 && cy(j - 1) != jc {
+                coarse.as_[c] += self.problem.as_[idx];
+            }
+        }
+        (coarse, coarse_of)
+    }
+}
+
+/// The multigrid hierarchy: level 0 is the fine problem.
+pub(crate) struct Hierarchy {
+    levels: Vec<Level>,
+    work: Vec<Work>,
+    sweeps: usize,
+}
+
+impl Hierarchy {
+    /// Builds the hierarchy down to `coarsest_cells` active cells (or until
+    /// coarsening stops reducing the count).
+    pub(crate) fn build(problem: &PoissonProblem, params: &MultigridParameters) -> Self {
+        let mut levels = vec![Level::new(problem.clone())];
+        while levels.len() < MAX_LEVELS {
+            let fine = levels.last().expect("at least one level");
+            if fine.cells.len() <= params.coarsest_cells.max(1) {
+                break;
+            }
+            let (coarse, coarse_of) = fine.coarsen();
+            let coarse = Level::new(coarse);
+            if coarse.cells.len() >= fine.cells.len() {
+                break;
+            }
+            let last = levels.len() - 1;
+            levels[last].coarse_of = coarse_of;
+            levels.push(coarse);
+        }
+        let work = levels
+            .iter()
+            .map(|l| Work::new(l.problem.nx * l.problem.ny))
+            .collect();
+        Self {
+            levels,
+            work,
+            sweeps: params.smoother_sweeps.max(1),
+        }
+    }
+
+    /// Number of levels.
+    pub(crate) fn depth(&self) -> usize {
+        self.levels.len()
+    }
+
+    /// The problem on level `l` (level 0 is the fine problem, with the
+    /// coefficients towards inactive cells sanitised to zero).
+    pub(crate) fn problem(&self, l: usize) -> &PoissonProblem {
+        &self.levels[l].problem
+    }
+
+    /// Parent map from level `l` to level `l + 1` (`usize::MAX` on cells
+    /// without an equation).
+    pub(crate) fn coarse_of(&self, l: usize) -> &[usize] {
+        &self.levels[l].coarse_of
+    }
+
+    /// Row-major indices of the active cells of level `l`.
+    pub(crate) fn cells(&self, l: usize) -> &[usize] {
+        &self.levels[l].cells
+    }
+
+    /// `z = M⁻¹ r`: one V-cycle on `A z = r` from `z = 0`. Only active
+    /// entries of `r` are read and only active entries of `z` are written.
+    pub(crate) fn apply_preconditioner(&mut self, r: &[f64], z: &mut [f64]) {
+        let depth = self.levels.len();
+        let levels = &self.levels;
+        for &idx in &levels[0].cells {
+            self.work[0].b[idx] = r[idx];
+        }
+        // Down: smooth from zero, restrict the residual.
+        for l in 0..depth - 1 {
+            let (fine, coarse) = (&levels[l], &levels[l + 1]);
+            let (head, tail) = self.work.split_at_mut(l + 1);
+            let (wf, wc) = (&mut head[l], &mut tail[0]);
+            for &idx in &fine.cells {
+                wf.x[idx] = 0.0;
+            }
+            for _ in 0..self.sweeps {
+                fine.symmetric_gs(&wf.b, &mut wf.x);
+            }
+            fine.residual(&wf.b, &wf.x, &mut wf.r);
+            for &idx in &coarse.cells {
+                wc.b[idx] = 0.0;
+            }
+            for &idx in &fine.cells {
+                let c = fine.coarse_of[idx];
+                // A coarse cell without an equation (a whole component
+                // inside one aggregate) receives the component's zero sum.
+                if coarse.active[c] {
+                    wc.b[c] += wf.r[idx];
+                }
+            }
+        }
+        // Coarsest: a fixed number of symmetric sweeps from zero.
+        {
+            let bottom = &levels[depth - 1];
+            let wb = &mut self.work[depth - 1];
+            for &idx in &bottom.cells {
+                wb.x[idx] = 0.0;
+            }
+            for _ in 0..COARSEST_SWEEPS {
+                bottom.symmetric_gs(&wb.b, &mut wb.x);
+            }
+        }
+        // Up: prolongate, smooth.
+        for l in (0..depth - 1).rev() {
+            let fine = &levels[l];
+            let (head, tail) = self.work.split_at_mut(l + 1);
+            let (wf, wc) = (&mut head[l], &tail[0]);
+            for &idx in &fine.cells {
+                wf.x[idx] += COARSE_CORRECTION * wc.x[fine.coarse_of[idx]];
+            }
+            for _ in 0..self.sweeps {
+                fine.symmetric_gs(&wf.b, &mut wf.x);
+            }
+        }
+        for &idx in &levels[0].cells {
+            z[idx] = self.work[0].x[idx];
+        }
+    }
+}
+
+/// Conjugate gradient on the active cells preconditioned by one V-cycle of
+/// geometric multigrid (see the module documentation for the hierarchy).
+///
+/// Stops when the L1 true residual `Σ |b - A p|` over the active cells is
+/// below `tolerance` (absolute — the caller passes a scale-relative value)
+/// or after `params.max_iterations` iterations.
+///
+/// Singular (pure Neumann) problems: CG runs on the consistent system (the
+/// caller guarantees a compatible rhs to rounding; the mean over active
+/// cells is projected out of the rhs to be safe, and the residual reported
+/// and tested is against that projected rhs — an incompatible rhs shows up
+/// as the gap between it and `problem.residual_l1(p)`), and on exit `p` is shifted
+/// so `p[anchor] == 0` when `anchor` is `Some` and names an active cell,
+/// else to mean zero over the active cells. On a non-singular problem the
+/// level is determined by the equations and `anchor` is ignored.
+///
+/// `p` is the initial guess and the result; inactive entries are neither
+/// read nor written.
+pub fn solve_multigrid_pcg(
+    problem: &PoissonProblem,
+    p: &mut [f64],
+    params: &MultigridParameters,
+    tolerance: f64,
+    anchor: Option<usize>,
+) -> PoissonSolution {
+    let n = problem.nx * problem.ny;
+    assert_eq!(p.len(), n, "p must have nx*ny entries");
+    debug_assert!(
+        problem.validate().is_ok(),
+        "invalid PoissonProblem: {:?}",
+        problem.validate()
+    );
+
+    let mut hier = Hierarchy::build(problem, params);
+    let cells: Vec<usize> = hier.levels[0].cells.clone();
+    let active_n = cells.len();
+    if active_n == 0 {
+        return PoissonSolution {
+            iterations: 0,
+            residual: 0.0,
+            converged: true,
+        };
+    }
+
+    // Connected components of the active cells (through faces with a
+    // non-zero coefficient). A component with no Dirichlet contribution is
+    // a pure-Neumann block of the system and singular ON ITS OWN: its
+    // right-hand side must have zero mean for the block to be consistent,
+    // whatever the other components carry. A single global mean projection
+    // is not enough — two enclosed pockets with opposite imbalances sum to
+    // zero globally and still make CG diverge, and the divergence pollutes
+    // even a well-posed Dirichlet component (found in review, 1e-8
+    // relative was enough). So the mean is projected per singular
+    // component, and the exit shift is applied per singular component.
+    let components = Components::find(problem, &cells);
+    let singular_any = components.singular.iter().any(|&s| s);
+
+    let project_mean = |v: &mut [f64]| {
+        for (c, members) in components.members.iter().enumerate() {
+            if !components.singular[c] {
+                continue;
+            }
+            let mean = members.iter().map(|&idx| v[idx]).sum::<f64>() / members.len() as f64;
+            for &idx in members {
+                v[idx] -= mean;
+            }
+        }
+    };
+    let dot = |a: &[f64], b: &[f64]| cells.iter().map(|&idx| a[idx] * b[idx]).sum::<f64>();
+    let l1 = |a: &[f64]| cells.iter().map(|&idx| a[idx].abs()).sum::<f64>();
+
+    // Right-hand side (per-component mean projected out where singular).
+    let mut b = vec![0.0; n];
+    for &idx in &cells {
+        b[idx] = problem.rhs[idx];
+    }
+    project_mean(&mut b);
+    let singular = singular_any;
+
+    let mut r = vec![0.0; n];
+    let mut z = vec![0.0; n];
+    let mut d = vec![0.0; n];
+    let mut q = vec![0.0; n];
+
+    let true_residual =
+        |p: &[f64], r: &mut [f64], hier: &Hierarchy| -> f64 { hier.levels[0].residual(&b, p, r) };
+
+    let anchor = anchor.filter(|&a| a < n && hier.levels[0].active[a]);
+    let finish = |p: &mut [f64], iterations: usize, residual: f64| {
+        // Level of each singular component: the anchor's component is
+        // shifted so p[anchor] == 0, every other singular component to mean
+        // zero. Non-singular components have their level fixed by their
+        // Dirichlet data and are left alone.
+        for (c, members) in components.members.iter().enumerate() {
+            if !components.singular[c] {
+                continue;
+            }
+            let shift = match anchor {
+                Some(a) if components.id[a] == c => p[a],
+                _ => members.iter().map(|&idx| p[idx]).sum::<f64>() / members.len() as f64,
+            };
+            for &idx in members {
+                p[idx] -= shift;
+            }
+        }
+        PoissonSolution {
+            iterations,
+            residual,
+            converged: residual < tolerance,
+        }
+    };
+
+    let mut res = true_residual(p, &mut r, &hier);
+    if res < tolerance {
+        return finish(p, 0, res);
+    }
+
+    hier.apply_preconditioner(&r, &mut z);
+    if singular {
+        project_mean(&mut z);
+    }
+    for &idx in &cells {
+        d[idx] = z[idx];
+    }
+    let mut rz = dot(&r, &z);
+
+    let mut iterations = 0;
+    // True residual at the last resynchronisation: when a resynchronised
+    // true residual no longer improves on the previous one, the recurrence
+    // has hit the rounding floor of `b - A p` and further iterations cannot
+    // reach the tolerance — stop, honestly unconverged.
+    let mut last_true = res;
+    while iterations < params.max_iterations {
+        iterations += 1;
+        hier.levels[0].apply(&d, &mut q);
+        let dq = dot(&d, &q);
+        if !dq.is_finite() || dq <= 0.0 || !rz.is_finite() || rz <= 0.0 {
+            // Breakdown (r = 0 to rounding, or a non-positive curvature
+            // from rounding in the null space): stop on the true residual.
+            res = true_residual(p, &mut r, &hier);
+            return finish(p, iterations, res);
+        }
+        let alpha = rz / dq;
+        for &idx in &cells {
+            p[idx] += alpha * d[idx];
+            r[idx] -= alpha * q[idx];
+        }
+        if l1(&r) < tolerance {
+            // The recurrence residual passed: confirm against the TRUE
+            // residual, and resynchronise if rounding has let them drift.
+            res = true_residual(p, &mut r, &hier);
+            if res < tolerance || res > 0.9 * last_true {
+                return finish(p, iterations, res);
+            }
+            last_true = res;
+        }
+        hier.apply_preconditioner(&r, &mut z);
+        if singular {
+            project_mean(&mut z);
+        }
+        let rz_new = dot(&r, &z);
+        let beta = rz_new / rz;
+        rz = rz_new;
+        for &idx in &cells {
+            d[idx] = z[idx] + beta * d[idx];
+        }
+    }
+    res = true_residual(p, &mut r, &hier);
+    finish(p, iterations, res)
+}
+
+/// Connected components of the active cells of a [`PoissonProblem`],
+/// connected through faces carrying a non-zero coefficient, and whether
+/// each component is singular (carries no Dirichlet contribution).
+struct Components {
+    /// Component id per cell (`usize::MAX` for inactive cells).
+    id: Vec<usize>,
+    /// Member cells per component.
+    members: Vec<Vec<usize>>,
+    /// Per component: no active member has `extra_diag > 0`.
+    singular: Vec<bool>,
+}
+
+impl Components {
+    fn find(problem: &PoissonProblem, cells: &[usize]) -> Self {
+        let (nx, ny) = (problem.nx, problem.ny);
+        let n = nx * ny;
+        let mut id = vec![usize::MAX; n];
+        let mut members = Vec::new();
+        let mut singular = Vec::new();
+        let mut stack = Vec::new();
+        for &seed in cells {
+            if id[seed] != usize::MAX {
+                continue;
+            }
+            let c = members.len();
+            let mut list = Vec::new();
+            let mut has_dirichlet = false;
+            id[seed] = c;
+            stack.push(seed);
+            while let Some(idx) = stack.pop() {
+                list.push(idx);
+                if problem.extra_diag[idx] > 0.0 {
+                    has_dirichlet = true;
+                }
+                let (j, i) = (idx / nx, idx % nx);
+                let mut visit = |nb: usize, coefficient: f64| {
+                    if coefficient > 0.0 && problem.active[nb] && id[nb] == usize::MAX {
+                        id[nb] = c;
+                        stack.push(nb);
+                    }
+                };
+                if i + 1 < nx {
+                    visit(idx + 1, problem.ae[idx]);
+                }
+                if i > 0 {
+                    visit(idx - 1, problem.aw[idx]);
+                }
+                if j + 1 < ny {
+                    visit(idx + nx, problem.an[idx]);
+                }
+                if j > 0 {
+                    visit(idx - nx, problem.as_[idx]);
+                }
+            }
+            members.push(list);
+            singular.push(!has_dirichlet);
+        }
+        Self {
+            id,
+            members,
+            singular,
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests;