rtx-cfd: curvilinear collocated PISO on a structured patch (overset A-P0, WIP) — PatchMesh (right-handed s,n; periodic seam with shift; face metrics), patch generators (TFI, skewed annulus, sheared/varying-skew channels), CSR + Jacobi-BiCGSTAB, the Zang–Street–Koseff incremental step with the node-based 9-point L_f, LSQ gradients, explicit and line-implicit-n predictors, adjustPhi; tests: mesh metrics (5 green), operators exact on linear fields incl. the seam (green), sparse (2 green), MMS ladder (Cartesian 16/32: 1.37–1.39x the staggered error, order 0.83; n=64 stalls at a |du/dt| floor 2e-4 — open, tolerance-scaling hypothesis), annulus/Poiseuille not yet run

Co-Authored-By: Claude Fable 5.1 <[email protected]> Claude-Session: https://claude.ai/code/session_01X2GmJXeQ2njUecEKiJZ1G2
2026-09-04 05:00:08 -07:00
parent 1347bc6772
commit 52da75a3a9
13 changed files with 2659 additions and 0 deletions
@@ -0,0 +1,333 @@
+//! Collocated finite-volume PISO on a structured curvilinear patch — the
+//! body-fitted half of the overset hybrid (`docs/overset_metal_campaign.md`
+//! §2.1, §5.3). Phase P0: a STATIC patch, manufactured-solution gated.
+//!
+//! The step is the Zang–Street–Koseff (1994) fractional step in the
+//! incremental form the background PISO uses:
+//!
+//! ```text
+//! û  = u^n + dt (−C(F^n, u^n) + ν D(u^n) + f/ρ)
+//! u* = û − (dt/ρ) G_c p^n                 cell gradient (least squares)
+//! F* = interp(û)·S_f − (dt/ρ) L_f(p^n)    compact face operator
+//! Σ_f (dt/ρ) L_f(p') = Σ_f F*             pressure correction
+//! F  = F* − (dt/ρ) L_f(p'),  u = u* − (dt/ρ) G_c p',  p += p'
+//! ```
+//!
+//! so the face flux — the primary variable for convection, divergence and
+//! the next step — sees the whole pressure through one compact operator
+//! (no checkerboard), and the cell velocity is a slave. Further correctors
+//! re-project the STORED fluxes. `L_f` is the orthogonal difference plus a
+//! node-based tangential correction (a 9-point stencil), assembled and
+//! applied by the same code.
+
+mod operators;
+mod predictor;
+mod projection;
+
+pub use operators::Operators;
+
+use crate::mesh::{PatchMesh, PatchSide};
+use crate::solvers::incompressible::sparse_bicgstab::CsrMatrix;
+use crate::{CfdConfig, CfdResult};
+
+/// What one side of the patch is.
+#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
+pub enum SideBc {
+    /// Prescribed velocity (a wall or an inflow): the boundary function
+    /// gives the velocity; mass flux `u_b · S_f`; pressure Neumann.
+    #[default]
+    Velocity,
+    /// Open boundary at gauge pressure zero: zero-gradient velocity,
+    /// flux from the predictor, corrected by the projection.
+    Outlet,
+}
+
+/// The four sides. `s_start`/`s_end` are ignored on a periodic patch.
+#[derive(Debug, Clone, Copy, Default)]
+pub struct PatchBoundaries {
+    /// `k = 0` (the wall of an O-grid).
+    pub inner: SideBc,
+    /// `k = nn`.
+    pub outer: SideBc,
+    /// `i = 0`.
+    pub s_start: SideBc,
+    /// `i = ns`.
+    pub s_end: SideBc,
+}
+
+impl PatchBoundaries {
+    /// The condition on a side.
+    pub fn get(&self, side: PatchSide) -> SideBc {
+        match side {
+            PatchSide::Inner => self.inner,
+            PatchSide::Outer => self.outer,
+            PatchSide::SStart => self.s_start,
+            PatchSide::SEnd => self.s_end,
+        }
+    }
+}
+
+/// Convection treatment.
+#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
+pub enum PatchConvection {
+    /// First-order upwind on the face fluxes (the background's scheme).
+    #[default]
+    Upwind,
+    /// No convection: the Stokes limit, for the second-order MMS gate.
+    None,
+}
+
+/// How the across-patch diffusion is time-stepped.
+#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
+pub enum NormalDiffusion {
+    /// Forward Euler, like the background (limit `~Δn²/(4ν)`).
+    #[default]
+    Explicit,
+    /// The orthogonal part of the n-face diffusion implicit along each
+    /// s-line (tridiagonal); the tangential part and everything else
+    /// explicit.
+    LineImplicit,
+}
+
+/// Solver parameters.
+#[derive(Debug, Clone)]
+pub struct CurvilinearParameters {
+    /// Projections per step (the first removes the divergence; the rest
+    /// mop up inner-solver truncation).
+    pub corrector_steps: usize,
+    /// Pressure-correction stop, relative to the step's flux scale.
+    pub tolerance: f64,
+    /// Convection scheme.
+    pub convection: PatchConvection,
+    /// Across-patch diffusion treatment.
+    pub normal_diffusion: NormalDiffusion,
+    /// Side conditions.
+    pub boundaries: PatchBoundaries,
+    /// BiCGSTAB iteration cap.
+    pub max_poisson_iterations: usize,
+}
+
+impl Default for CurvilinearParameters {
+    fn default() -> Self {
+        Self {
+            corrector_steps: 2,
+            tolerance: 1e-8,
+            convection: PatchConvection::Upwind,
+            normal_diffusion: NormalDiffusion::Explicit,
+            boundaries: PatchBoundaries::default(),
+            max_poisson_iterations: 5000,
+        }
+    }
+}
+
+/// The patch state: cell-centred velocity and pressure, face mass fluxes.
+#[derive(Debug, Clone)]
+pub struct PatchField {
+    /// Cell x-velocity.
+    pub u: Vec<f64>,
+    /// Cell y-velocity.
+    pub v: Vec<f64>,
+    /// Cell pressure.
+    pub p: Vec<f64>,
+    /// Volume flux through every face, oriented +s / +n.
+    pub flux: Vec<f64>,
+}
+
+impl PatchField {
+    /// Zero field on `mesh`.
+    pub fn new(mesh: &PatchMesh) -> Self {
+        let n = mesh.cell_count();
+        Self {
+            u: vec![0.0; n],
+            v: vec![0.0; n],
+            p: vec![0.0; n],
+            flux: vec![0.0; mesh.faces().len()],
+        }
+    }
+}
+
+/// What one step reports.
+#[derive(Debug, Clone, Copy)]
+pub struct CurvilinearResult {
+    /// Correctors performed.
+    pub corrector_steps_performed: usize,
+    /// Largest cell mass imbalance after the last corrector.
+    pub max_divergence: f64,
+    /// Pressure-solver iterations, summed over the correctors.
+    pub poisson_iterations: usize,
+    /// Whether every pressure solve reached its tolerance.
+    pub poisson_converged: bool,
+    /// Boundary-flux defect removed on a closed patch (zero if an outlet exists).
+    pub boundary_flux_adjustment: f64,
+}
+
+/// Restorable solver state (the coupling re-runs a step).
+#[derive(Debug, Clone)]
+pub struct CurvilinearSolverState {
+    time: f64,
+}
+
+type VelocityFn = Box<dyn Fn(f64, f64, f64) -> (f64, f64) + Send + Sync>;
+
+/// PISO on a curvilinear patch.
+pub struct CurvilinearPisoSolver {
+    config: CfdConfig,
+    params: CurvilinearParameters,
+    mesh: PatchMesh,
+    ops: Operators,
+    boundary_velocity: Option<VelocityFn>,
+    momentum_source: Option<VelocityFn>,
+    time: f64,
+    matrix: Option<(f64, CsrMatrix, Option<usize>)>,
+}
+
+impl CurvilinearPisoSolver {
+    /// Build on `mesh`.
+    pub fn new(
+        config: CfdConfig,
+        params: CurvilinearParameters,
+        mesh: PatchMesh,
+    ) -> CfdResult<Self> {
+        let ops = Operators::new(&mesh, &params.boundaries);
+        Ok(Self {
+            config,
+            params,
+            mesh,
+            ops,
+            boundary_velocity: None,
+            momentum_source: None,
+            time: 0.0,
+            matrix: None,
+        })
+    }
+
+    /// Velocity on every `Velocity` side, `(x, y, t) -> (u, v)`.
+    pub fn set_boundary_velocity<F>(&mut self, f: F)
+    where
+        F: Fn(f64, f64, f64) -> (f64, f64) + Send + Sync + 'static,
+    {
+        self.boundary_velocity = Some(Box::new(f));
+    }
+    /// Body force per unit volume, `(x, y, t) -> (fx, fy)`.
+    pub fn set_momentum_source<F>(&mut self, f: F)
+    where
+        F: Fn(f64, f64, f64) -> (f64, f64) + Send + Sync + 'static,
+    {
+        self.momentum_source = Some(Box::new(f));
+    }
+    /// The mesh.
+    pub fn mesh(&self) -> &PatchMesh {
+        &self.mesh
+    }
+    /// The operators.
+    pub fn operators(&self) -> &Operators {
+        &self.ops
+    }
+    /// Parameters.
+    pub fn parameters(&self) -> &CurvilinearParameters {
+        &self.params
+    }
+    /// Configuration.
+    pub fn config(&self) -> &CfdConfig {
+        &self.config
+    }
+    /// Current time.
+    pub fn time(&self) -> f64 {
+        self.time
+    }
+    /// Set the time.
+    pub fn set_time(&mut self, t: f64) {
+        self.time = t;
+    }
+    /// Capture the state.
+    pub fn snapshot(&self) -> CurvilinearSolverState {
+        CurvilinearSolverState { time: self.time }
+    }
+    /// Restore a captured state.
+    pub fn restore(&mut self, state: &CurvilinearSolverState) {
+        self.time = state.time;
+    }
+
+    pub(crate) fn boundary_velocity(&self, x: f64, y: f64, t: f64) -> (f64, f64) {
+        self.boundary_velocity
+            .as_ref()
+            .map_or((0.0, 0.0), |f| f(x, y, t))
+    }
+    pub(crate) fn source_at(&self, xy: [f64; 2], t: f64) -> (f64, f64) {
+        self.momentum_source
+            .as_ref()
+            .map_or((0.0, 0.0), |f| f(xy[0], xy[1], t))
+    }
+
+    /// Set the cell velocities from a function and make the face fluxes
+    /// consistent (interpolated; prescribed on velocity sides).
+    pub fn initialize<F>(&self, field: &mut PatchField, velocity: F)
+    where
+        F: Fn(f64, f64) -> (f64, f64),
+    {
+        let mesh = &self.mesh;
+        for c in 0..mesh.cell_count() {
+            let xy = mesh.centre(c);
+            let (u, v) = velocity(xy[0], xy[1]);
+            field.u[c] = u;
+            field.v[c] = v;
+        }
+        let zero = vec![0.0; mesh.cell_count()];
+        field.flux =
+            self.predicted_fluxes(&field.u.clone(), &field.v.clone(), &zero, 0.0, self.time);
+    }
+
+    /// Advance one step of `dt`.
+    pub async fn advance(
+        &mut self,
+        field: &mut PatchField,
+        dt: f64,
+    ) -> CfdResult<CurvilinearResult> {
+        let t_old = self.time;
+        let t_new = t_old + dt;
+        let rho = self.config.density;
+        let mesh = &self.mesh;
+
+        let (uh, vh) = self.predict(field, dt, t_old);
+        let mut flux = self.predicted_fluxes(&uh, &vh, &field.p, dt, t_new);
+        let adjustment = self.adjust_boundary_flux(&mut flux);
+        for c in 0..mesh.cell_count() {
+            let g = self.pressure_gradient(&field.p, c);
+            field.u[c] = uh[c] - dt / rho * g[0];
+            field.v[c] = vh[c] - dt / rho * g[1];
+        }
+        field.flux = flux;
+
+        if self.matrix.as_ref().is_none_or(|(d, _, _)| *d != dt) {
+            let (m, anchor) = self.assemble_pressure_matrix(dt);
+            self.matrix = Some((dt, m, anchor));
+        }
+        let (_, matrix, anchor) = self.matrix.as_ref().expect("assembled");
+        let flux_scale: f64 = field.flux.iter().map(|f| f.abs()).sum::<f64>().max(1e-300);
+        let tolerance = self.params.tolerance * flux_scale;
+
+        let mut iterations = 0;
+        let mut converged = true;
+        let mut performed = 0;
+        let mut max_div = self.max_divergence(&field.flux);
+        for _ in 0..self.params.corrector_steps {
+            let (pc, out) = self.solve_pressure_correction(matrix, *anchor, &field.flux, tolerance);
+            iterations += out.iterations;
+            converged &= out.converged;
+            self.apply_correction(field, &pc, dt);
+            performed += 1;
+            max_div = self.max_divergence(&field.flux);
+            if max_div <= tolerance {
+                break;
+            }
+        }
+        self.time = t_new;
+        Ok(CurvilinearResult {
+            corrector_steps_performed: performed,
+            max_divergence: max_div,
+            poisson_iterations: iterations,
+            poisson_converged: converged,
+            boundary_flux_adjustment: adjustment,
+        })
+    }
+}
@@ -0,0 +1,306 @@
+//! Discrete operators on the patch: linear-exact nodal reconstruction,
+//! the compact face-gradient operator `L_f` (orthogonal part plus the
+//! node-based tangential correction), weighted least-squares cell
+//! gradients, and face interpolation. One `L_f` serves velocity diffusion,
+//! the pressure term in the face flux, and the pressure equation's matrix
+//! — the consistency the collocated scheme depends on.
+
+use super::{PatchBoundaries, SideBc};
+use crate::mesh::{PatchMesh, PatchSide};
+
+/// How a node's value is reconstructed from the cells around it.
+#[derive(Debug, Clone)]
+pub struct NodeStencil {
+    /// Cells touching the node.
+    pub cells: Vec<usize>,
+    /// Their weights (sum to one; exact on linear fields for 4-cell nodes).
+    pub weights: Vec<f64>,
+    /// The boundary the node lies on, if any (`Inner`/`Outer` win at corners).
+    pub side: Option<PatchSide>,
+}
+
+/// Per-face constants of `L_f`.
+#[derive(Debug, Clone, Copy)]
+pub struct FaceCoefs {
+    /// `|S_f|² / (d_f · S_f)`: the orthogonal (implicit) coefficient.
+    pub alpha: f64,
+    /// `(d_f · t̂) / |f|`: the tangential correction is `−alpha · tang · (φ_n1 − φ_n0)`.
+    pub tang: f64,
+}
+
+/// Precomputed operator data for one patch geometry.
+#[derive(Debug, Clone)]
+pub struct Operators {
+    nodes: Vec<NodeStencil>,
+    faces: Vec<FaceCoefs>,
+}
+
+impl Operators {
+    /// Build for `mesh`; `boundaries` decides which nodes are Dirichlet
+    /// for pressure (outlet sides).
+    pub fn new(mesh: &PatchMesh, _boundaries: &PatchBoundaries) -> Self {
+        let (ns, nn) = (mesh.ns(), mesh.nn());
+        let mut nodes = Vec::with_capacity((nn + 1) * (ns + 1));
+        for k in 0..=nn {
+            for i in 0..=ns {
+                nodes.push(Self::node_stencil(mesh, k, i));
+            }
+        }
+        let faces = mesh
+            .faces()
+            .iter()
+            .map(|f| {
+                let s2 = f.s[0] * f.s[0] + f.s[1] * f.s[1];
+                let ds = f.d[0] * f.s[0] + f.d[1] * f.s[1];
+                let alpha = if ds != 0.0 { s2 / ds } else { 0.0 };
+                let t = [
+                    mesh.node_xy(f.n1)[0] - mesh.node_xy(f.n0)[0],
+                    mesh.node_xy(f.n1)[1] - mesh.node_xy(f.n0)[1],
+                ];
+                let len = (t[0] * t[0] + t[1] * t[1]).sqrt();
+                let tang = if len > 0.0 {
+                    (f.d[0] * t[0] + f.d[1] * t[1]) / (len * len)
+                } else {
+                    0.0
+                };
+                FaceCoefs { alpha, tang }
+            })
+            .collect();
+        Self { nodes, faces }
+    }
+
+    fn node_stencil(mesh: &PatchMesh, k: usize, i: usize) -> NodeStencil {
+        let (ns, nn) = (mesh.ns(), mesh.nn());
+        let side = if k == 0 {
+            Some(PatchSide::Inner)
+        } else if k == nn {
+            Some(PatchSide::Outer)
+        } else if mesh.periodic().is_none() && i == 0 {
+            Some(PatchSide::SStart)
+        } else if mesh.periodic().is_none() && i == ns {
+            Some(PatchSide::SEnd)
+        } else {
+            None
+        };
+        let cells = mesh.node_cells(k, i);
+        let xn = mesh.node_xy(mesh.node(k, i));
+        // Cell centres relative to the node, with the periodic shift for
+        // cells reached across the seam (a cell whose column is ns-1 seen
+        // from node column 0, or column 0 seen from node column ns).
+        let rel: Vec<[f64; 2]> = cells
+            .iter()
+            .map(|&c| {
+                let (_, ci) = mesh.cell_ki(c);
+                let mut xc = mesh.centre(c);
+                if let Some(shift) = mesh.periodic() {
+                    if i == 0 && ci == ns - 1 {
+                        xc = [xc[0] - shift[0], xc[1] - shift[1]];
+                    } else if i == ns && ci == 0 {
+                        xc = [xc[0] + shift[0], xc[1] + shift[1]];
+                    }
+                }
+                [xc[0] - xn[0], xc[1] - xn[1]]
+            })
+            .collect();
+        let weights = linear_exact_weights(&rel);
+        NodeStencil {
+            cells,
+            weights,
+            side,
+        }
+    }
+
+    /// Node stencils.
+    pub fn nodes(&self) -> &[NodeStencil] {
+        &self.nodes
+    }
+    /// Face coefficients.
+    pub fn face(&self, f: usize) -> FaceCoefs {
+        self.faces[f]
+    }
+
+    /// Node values of a cell field. `boundary(side, xy)` supplies a
+    /// prescribed value at nodes on a boundary side (`None` = extrapolate
+    /// from the cells like an interior node).
+    pub fn node_values(
+        &self,
+        mesh: &PatchMesh,
+        cell_vals: &[f64],
+        boundary: &dyn Fn(PatchSide, [f64; 2]) -> Option<f64>,
+    ) -> Vec<f64> {
+        self.nodes
+            .iter()
+            .enumerate()
+            .map(|(n, st)| {
+                if let Some(side) = st.side {
+                    if let Some(v) = boundary(side, mesh.node_xy(n)) {
+                        return v;
+                    }
+                }
+                st.cells
+                    .iter()
+                    .zip(&st.weights)
+                    .map(|(&c, &w)| w * cell_vals[c])
+                    .sum()
+            })
+            .collect()
+    }
+
+    /// `L_f(φ) ≈ (∇φ)_f · S_f` for one face, given cell and node values;
+    /// `boundary_value` is the face's Dirichlet value on a boundary face
+    /// (`None` = zero-gradient: the flux is zero).
+    pub fn face_gradient_flux(
+        &self,
+        mesh: &PatchMesh,
+        f: usize,
+        cell_vals: &[f64],
+        node_vals: &[f64],
+        boundary_value: Option<f64>,
+    ) -> f64 {
+        let face = &mesh.faces()[f];
+        let FaceCoefs { alpha, tang } = self.faces[f];
+        let (phi_p, phi_n) = match (face.owner, face.neigh) {
+            (Some(p), Some(n)) => (cell_vals[p], cell_vals[n]),
+            (Some(p), None) => match boundary_value {
+                Some(b) => (cell_vals[p], b),
+                None => return 0.0,
+            },
+            (None, Some(n)) => match boundary_value {
+                Some(b) => (b, cell_vals[n]),
+                None => return 0.0,
+            },
+            (None, None) => unreachable!(),
+        };
+        alpha * ((phi_n - phi_p) - tang * (node_vals[face.n1] - node_vals[face.n0]))
+    }
+
+    /// The coefficient list of `L_f` in terms of cell unknowns for the
+    /// pressure-correction equation: interior faces expand the nodal term
+    /// through the node stencils (skipping nodes on Dirichlet sides, whose
+    /// correction is zero); Dirichlet boundary faces contribute `∓alpha`
+    /// on the interior cell; Neumann faces contribute nothing.
+    pub fn face_gradient_coeffs(
+        &self,
+        mesh: &PatchMesh,
+        boundaries: &PatchBoundaries,
+        f: usize,
+        out: &mut Vec<(usize, f64)>,
+    ) {
+        out.clear();
+        let face = &mesh.faces()[f];
+        let FaceCoefs { alpha, tang } = self.faces[f];
+        match (face.owner, face.neigh) {
+            (Some(p), Some(n)) => {
+                out.push((p, -alpha));
+                out.push((n, alpha));
+                for (node, sign) in [(face.n0, alpha * tang), (face.n1, -alpha * tang)] {
+                    let st = &self.nodes[node];
+                    if st.side.is_some_and(|s| boundaries.get(s) == SideBc::Outlet) {
+                        continue;
+                    }
+                    for (&c, &w) in st.cells.iter().zip(&st.weights) {
+                        out.push((c, sign * w));
+                    }
+                }
+            }
+            (Some(p), None) => {
+                if boundaries.get(mesh.side(f).expect("boundary")) == SideBc::Outlet {
+                    out.push((p, -alpha));
+                }
+            }
+            (None, Some(n)) => {
+                if boundaries.get(mesh.side(f).expect("boundary")) == SideBc::Outlet {
+                    out.push((n, alpha));
+                }
+            }
+            (None, None) => unreachable!(),
+        }
+    }
+
+    /// Weighted least-squares gradient of a cell field at cell `c`.
+    /// `boundary_value(face)` gives the value at a boundary face centre
+    /// (`None` = the face is left out of the fit).
+    pub fn gradient(
+        &self,
+        mesh: &PatchMesh,
+        c: usize,
+        cell_vals: &[f64],
+        boundary_value: &dyn Fn(usize) -> Option<f64>,
+    ) -> [f64; 2] {
+        let (mut sxx, mut sxy, mut syy, mut bx, mut by) = (0.0, 0.0, 0.0, 0.0, 0.0);
+        let phi_p = cell_vals[c];
+        for (f, sign) in mesh.cell_faces(c) {
+            let face = &mesh.faces()[f];
+            let (r, dphi) = match (face.owner, face.neigh) {
+                (Some(p), Some(n)) => {
+                    let other = if p == c { n } else { p };
+                    (
+                        [sign * face.d[0], sign * face.d[1]],
+                        cell_vals[other] - phi_p,
+                    )
+                }
+                _ => match boundary_value(f) {
+                    Some(b) => ([sign * face.d[0], sign * face.d[1]], b - phi_p),
+                    None => continue,
+                },
+            };
+            let w = 1.0 / (r[0] * r[0] + r[1] * r[1]);
+            sxx += w * r[0] * r[0];
+            sxy += w * r[0] * r[1];
+            syy += w * r[1] * r[1];
+            bx += w * r[0] * dphi;
+            by += w * r[1] * dphi;
+        }
+        let det = sxx * syy - sxy * sxy;
+        if det.abs() <= 1e-300 {
+            return [0.0, 0.0];
+        }
+        [(syy * bx - sxy * by) / det, (sxx * by - sxy * bx) / det]
+    }
+}
+
+/// Weights over points `rel` (relative to the reconstruction point) that
+/// sum to one and reproduce linear fields when the points allow it
+/// (Holmes–Connell pseudo-Laplacian); falls back to inverse distance if
+/// the fit is singular or produces a negative weight, and to the 1-D
+/// linear fit for two points.
+fn linear_exact_weights(rel: &[[f64; 2]]) -> Vec<f64> {
+    let n = rel.len();
+    if n == 1 {
+        return vec![1.0];
+    }
+    if n == 2 {
+        let d = [rel[1][0] - rel[0][0], rel[1][1] - rel[0][1]];
+        let dd = d[0] * d[0] + d[1] * d[1];
+        let t = if dd > 0.0 {
+            -(rel[0][0] * d[0] + rel[0][1] * d[1]) / dd
+        } else {
+            0.5
+        };
+        return vec![1.0 - t, t];
+    }
+    let (mut rx, mut ry, mut ixx, mut ixy, mut iyy) = (0.0, 0.0, 0.0, 0.0, 0.0);
+    for r in rel {
+        rx += r[0];
+        ry += r[1];
+        ixx += r[0] * r[0];
+        ixy += r[0] * r[1];
+        iyy += r[1] * r[1];
+    }
+    let det = ixx * iyy - ixy * ixy;
+    let mut weights: Vec<f64> = if det.abs() > 1e-300 {
+        let lx = -(iyy * rx - ixy * ry) / det;
+        let ly = -(ixx * ry - ixy * rx) / det;
+        rel.iter().map(|r| 1.0 + lx * r[0] + ly * r[1]).collect()
+    } else {
+        vec![-1.0]
+    };
+    if weights.iter().any(|&w| w < 0.0) {
+        weights = rel
+            .iter()
+            .map(|r| 1.0 / (r[0] * r[0] + r[1] * r[1]).sqrt().max(1e-300))
+            .collect();
+    }
+    let sum: f64 = weights.iter().sum();
+    weights.iter().map(|w| w / sum).collect()
+}
@@ -0,0 +1,171 @@
+//! The momentum predictor: explicit forward Euler from the old field
+//! (matching the background PISO), or line-implicit in `n` for the
+//! orthogonal part of the across-patch diffusion.
+
+use super::{CurvilinearPisoSolver, NormalDiffusion, PatchConvection, PatchField, SideBc};
+use crate::mesh::PatchSide;
+
+impl CurvilinearPisoSolver {
+    /// `û = u^n + dt (−C(F^n, u^n) + ν D(u^n) + f/ρ)`, no pressure.
+    pub(super) fn predict(&self, field: &PatchField, dt: f64, t_old: f64) -> (Vec<f64>, Vec<f64>) {
+        let mesh = &self.mesh;
+        let n = mesh.cell_count();
+        let rho = self.config.density;
+        let nu = self.config.viscosity / rho;
+        let ops = &self.ops;
+
+        let bvel = |side: PatchSide, xy: [f64; 2]| -> Option<(f64, f64)> {
+            match self.params.boundaries.get(side) {
+                SideBc::Velocity => Some(self.boundary_velocity(xy[0], xy[1], t_old)),
+                SideBc::Outlet => None,
+            }
+        };
+        let un = ops.node_values(mesh, &field.u, &|s, xy| bvel(s, xy).map(|v| v.0));
+        let vn = ops.node_values(mesh, &field.v, &|s, xy| bvel(s, xy).map(|v| v.1));
+
+        let implicit_n = matches!(self.params.normal_diffusion, NormalDiffusion::LineImplicit);
+        let mut uh = vec![0.0; n];
+        let mut vh = vec![0.0; n];
+        for c in 0..n {
+            let (mut cu, mut cv, mut du, mut dv) = (0.0, 0.0, 0.0, 0.0);
+            for (f, sign) in mesh.cell_faces(c) {
+                let face = &mesh.faces()[f];
+                let side = mesh.side(f);
+                let bval = side.and_then(|s| bvel(s, face.centre));
+                // Convection: outward flux times the upwind face value.
+                if matches!(self.params.convection, PatchConvection::Upwind) {
+                    let out = sign * field.flux[f];
+                    let (uf, vf) = match (face.owner, face.neigh) {
+                        (Some(p), Some(q)) => {
+                            let up = if out >= 0.0 {
+                                c
+                            } else if p == c {
+                                q
+                            } else {
+                                p
+                            };
+                            (field.u[up], field.v[up])
+                        }
+                        _ => match bval {
+                            Some(b) => b,
+                            None => (field.u[c], field.v[c]),
+                        },
+                    };
+                    cu += out * uf;
+                    cv += out * vf;
+                }
+                // Diffusion: ν L_f(u), the orthogonal n-face part deferred
+                // to the tridiagonal solve when line-implicit.
+                let is_n = !mesh.is_sface(f);
+                let (bu, bv) = match bval {
+                    Some(b) => (Some(b.0), Some(b.1)),
+                    None => (None, None),
+                };
+                let mut lu = ops.face_gradient_flux(mesh, f, &field.u, &un, bu);
+                let mut lv = ops.face_gradient_flux(mesh, f, &field.v, &vn, bv);
+                if implicit_n && is_n && (side.is_none() || bval.is_some()) {
+                    let alpha = ops.face(f).alpha;
+                    let (pu, pv) = match (face.owner, face.neigh) {
+                        (Some(p), Some(q)) => {
+                            let other = if p == c { q } else { p };
+                            (
+                                sign * alpha * (field.u[other] - field.u[c]),
+                                sign * alpha * (field.v[other] - field.v[c]),
+                            )
+                        }
+                        _ => {
+                            let b = bval.expect("dirichlet");
+                            (
+                                sign * alpha * (b.0 - field.u[c]),
+                                sign * alpha * (b.1 - field.v[c]),
+                            )
+                        }
+                    };
+                    lu = sign * lu - pu;
+                    lv = sign * lv - pv;
+                } else {
+                    lu *= sign;
+                    lv *= sign;
+                }
+                du += lu;
+                dv += lv;
+            }
+            let a = mesh.area(c);
+            let (fx, fy) = self.source_at(mesh.centre(c), t_old);
+            uh[c] = field.u[c] + dt * ((-cu + nu * du) / a + fx / rho);
+            vh[c] = field.v[c] + dt * ((-cv + nu * dv) / a + fy / rho);
+        }
+
+        if implicit_n {
+            self.solve_lines(&mut uh, &mut vh, dt, nu, t_old);
+        }
+        (uh, vh)
+    }
+
+    /// `(I − dt ν L_n/A) û = rhs` along every s-line, Thomas algorithm.
+    fn solve_lines(&self, uh: &mut [f64], vh: &mut [f64], dt: f64, nu: f64, t_old: f64) {
+        let mesh = &self.mesh;
+        let (ns, nn) = (mesh.ns(), mesh.nn());
+        let mut lower = vec![0.0; nn];
+        let mut diag = vec![0.0; nn];
+        let mut upper = vec![0.0; nn];
+        let mut ru = vec![0.0; nn];
+        let mut rv = vec![0.0; nn];
+        for i in 0..ns {
+            for k in 0..nn {
+                let c = mesh.cell(k, i);
+                let a = mesh.area(c);
+                let mut d = 1.0;
+                let (mut lo, mut up) = (0.0, 0.0);
+                ru[k] = uh[c];
+                rv[k] = vh[c];
+                for (f, sign) in [(mesh.nface(k, i), -1.0), (mesh.nface(k + 1, i), 1.0)] {
+                    let alpha = self.ops.face(f).alpha;
+                    let coef = dt * nu * alpha / a;
+                    match mesh.side(f) {
+                        None => {
+                            d += coef;
+                            if sign < 0.0 {
+                                lo = -coef;
+                            } else {
+                                up = -coef;
+                            }
+                        }
+                        Some(s) => match self.params.boundaries.get(s) {
+                            SideBc::Velocity => {
+                                let face = &mesh.faces()[f];
+                                let b =
+                                    self.boundary_velocity(face.centre[0], face.centre[1], t_old);
+                                d += coef;
+                                ru[k] += coef * b.0;
+                                rv[k] += coef * b.1;
+                            }
+                            SideBc::Outlet => {}
+                        },
+                    }
+                }
+                lower[k] = lo;
+                diag[k] = d;
+                upper[k] = up;
+            }
+            // Thomas.
+            for k in 1..nn {
+                let m = lower[k] / diag[k - 1];
+                diag[k] -= m * upper[k - 1];
+                ru[k] -= m * ru[k - 1];
+                rv[k] -= m * rv[k - 1];
+            }
+            ru[nn - 1] /= diag[nn - 1];
+            rv[nn - 1] /= diag[nn - 1];
+            for k in (0..nn - 1).rev() {
+                ru[k] = (ru[k] - upper[k] * ru[k + 1]) / diag[k];
+                rv[k] = (rv[k] - upper[k] * rv[k + 1]) / diag[k];
+            }
+            for k in 0..nn {
+                let c = mesh.cell(k, i);
+                uh[c] = ru[k];
+                vh[c] = rv[k];
+            }
+        }
+    }
+}
@@ -0,0 +1,223 @@
+//! The pressure step: face fluxes from the predictor with the compact
+//! pressure term, boundary-flux adjustment for closed patches, the
+//! pressure-correction equation on the 9-point operator, and the flux and
+//! velocity corrections.
+
+use super::{CurvilinearPisoSolver, PatchField, SideBc};
+use crate::mesh::PatchSide;
+use crate::solvers::incompressible::sparse_bicgstab::{
+    BicgstabResult, CsrMatrix, bicgstab_jacobi, project_mean,
+};
+
+impl CurvilinearPisoSolver {
+    /// Node values of a pressure-like cell field: zero on outlet sides,
+    /// extrapolated elsewhere.
+    pub(super) fn pressure_nodes(&self, p: &[f64]) -> Vec<f64> {
+        let b = &self.params.boundaries;
+        self.ops
+            .node_values(&self.mesh, p, &|side: PatchSide, _| match b.get(side) {
+                SideBc::Outlet => Some(0.0),
+                SideBc::Velocity => None,
+            })
+    }
+
+    /// `L_f(p)` on every face (zero on Neumann faces, outlet value zero).
+    pub(super) fn pressure_face_gradients(&self, p: &[f64]) -> Vec<f64> {
+        let mesh = &self.mesh;
+        let pn = self.pressure_nodes(p);
+        (0..mesh.faces().len())
+            .map(|f| {
+                let bval = mesh
+                    .side(f)
+                    .and_then(|s| match self.params.boundaries.get(s) {
+                        SideBc::Outlet => Some(0.0),
+                        SideBc::Velocity => None,
+                    });
+                self.ops.face_gradient_flux(mesh, f, p, &pn, bval)
+            })
+            .collect()
+    }
+
+    /// Least-squares cell gradient of a pressure-like field (outlet faces
+    /// at zero, Neumann faces left out).
+    pub(super) fn pressure_gradient(&self, p: &[f64], c: usize) -> [f64; 2] {
+        let mesh = &self.mesh;
+        self.ops.gradient(mesh, c, p, &|f| match mesh.side(f) {
+            Some(s) if self.params.boundaries.get(s) == SideBc::Outlet => Some(0.0),
+            _ => None,
+        })
+    }
+
+    /// `F* = interp(û)·S − (dt/ρ) L_f(p^n)` on interior and outlet faces,
+    /// the prescribed flux on velocity faces (at `t_new`).
+    pub(super) fn predicted_fluxes(
+        &self,
+        uh: &[f64],
+        vh: &[f64],
+        p: &[f64],
+        dt: f64,
+        t_new: f64,
+    ) -> Vec<f64> {
+        let mesh = &self.mesh;
+        let rho = self.config.density;
+        let lp = self.pressure_face_gradients(p);
+        mesh.faces()
+            .iter()
+            .enumerate()
+            .map(|(f, face)| match (face.owner, face.neigh) {
+                (Some(o), Some(n)) => {
+                    let w = face.w;
+                    let uf = w * uh[o] + (1.0 - w) * uh[n];
+                    let vf = w * vh[o] + (1.0 - w) * vh[n];
+                    uf * face.s[0] + vf * face.s[1] - dt / rho * lp[f]
+                }
+                _ => {
+                    let c = mesh.boundary_cell(f);
+                    match self.params.boundaries.get(mesh.side(f).expect("boundary")) {
+                        SideBc::Velocity => {
+                            let (ub, vb) =
+                                self.boundary_velocity(face.centre[0], face.centre[1], t_new);
+                            ub * face.s[0] + vb * face.s[1]
+                        }
+                        SideBc::Outlet => uh[c] * face.s[0] + vh[c] * face.s[1] - dt / rho * lp[f],
+                    }
+                }
+            })
+            .collect()
+    }
+
+    /// On a patch with no outlet the prescribed boundary fluxes must sum
+    /// to zero for the projection to be solvable; the O(h²) defect of
+    /// face-centre sampling is spread over the velocity faces by area
+    /// (OpenFOAM's `adjustPhi`). Returns the defect removed.
+    pub(super) fn adjust_boundary_flux(&self, flux: &mut [f64]) -> f64 {
+        let mesh = &self.mesh;
+        let has_outlet = [
+            PatchSide::Inner,
+            PatchSide::Outer,
+            PatchSide::SStart,
+            PatchSide::SEnd,
+        ]
+        .iter()
+        .any(|&s| self.params.boundaries.get(s) == SideBc::Outlet);
+        if has_outlet {
+            return 0.0;
+        }
+        let (mut net, mut total_len) = (0.0, 0.0);
+        for (f, face) in mesh.faces().iter().enumerate() {
+            if mesh.side(f).is_some() {
+                let out_sign = if face.owner.is_some() { 1.0 } else { -1.0 };
+                net += out_sign * flux[f];
+                total_len += (face.s[0] * face.s[0] + face.s[1] * face.s[1]).sqrt();
+            }
+        }
+        if total_len == 0.0 {
+            return net;
+        }
+        for (f, face) in mesh.faces().iter().enumerate() {
+            if mesh.side(f).is_some() {
+                let out_sign = if face.owner.is_some() { 1.0 } else { -1.0 };
+                let len = (face.s[0] * face.s[0] + face.s[1] * face.s[1]).sqrt();
+                flux[f] -= out_sign * net * len / total_len;
+            }
+        }
+        net
+    }
+
+    /// Assemble `−Σ_f sign (dt/ρ) L_f` (positive diagonal) and pick the
+    /// anchor for the pure-Neumann case.
+    pub(super) fn assemble_pressure_matrix(&self, dt: f64) -> (CsrMatrix, Option<usize>) {
+        let mesh = &self.mesh;
+        let rho = self.config.density;
+        let n = mesh.cell_count();
+        let mut tri = Vec::with_capacity(n * 12);
+        let mut coefs = Vec::new();
+        let mut any_dirichlet = false;
+        for c in 0..n {
+            for (f, sign) in mesh.cell_faces(c) {
+                self.ops
+                    .face_gradient_coeffs(mesh, &self.params.boundaries, f, &mut coefs);
+                if mesh.side(f).is_some() && !coefs.is_empty() {
+                    any_dirichlet = true;
+                }
+                for &(col, v) in &coefs {
+                    tri.push((c, col, -sign * dt / rho * v));
+                }
+            }
+            tri.push((c, c, 0.0)); // guarantee a diagonal entry
+        }
+        let mut a = CsrMatrix::from_triplets(n, &tri);
+        let anchor = if any_dirichlet {
+            None
+        } else {
+            // An interior cell away from the seam: (1, 1).
+            let a_cell = mesh.cell(1.min(mesh.nn() - 1), 1.min(mesh.ns() - 1));
+            a.set_row_identity(a_cell);
+            Some(a_cell)
+        };
+        (a, anchor)
+    }
+
+    /// Solve `−Σ sign (dt/ρ) L_f(p') = −Σ sign F` for `p'` (zero start).
+    pub(super) fn solve_pressure_correction(
+        &self,
+        matrix: &CsrMatrix,
+        anchor: Option<usize>,
+        flux: &[f64],
+        tolerance: f64,
+    ) -> (Vec<f64>, BicgstabResult) {
+        let mesh = &self.mesh;
+        let n = mesh.cell_count();
+        let mut rhs = vec![0.0; n];
+        for c in 0..n {
+            let mut div = 0.0;
+            for (f, sign) in mesh.cell_faces(c) {
+                div += sign * flux[f];
+            }
+            rhs[c] = -div;
+        }
+        if let Some(a) = anchor {
+            project_mean(&mut rhs);
+            rhs[a] = 0.0;
+        }
+        let mut pc = vec![0.0; n];
+        let out = bicgstab_jacobi(
+            matrix,
+            &rhs,
+            &mut pc,
+            tolerance,
+            self.params.max_poisson_iterations,
+        );
+        (pc, out)
+    }
+
+    /// `F −= (dt/ρ) L_f(p')`, `u −= (dt/ρ) ∇p'`, `p += p'`.
+    pub(super) fn apply_correction(&self, field: &mut PatchField, pc: &[f64], dt: f64) {
+        let mesh = &self.mesh;
+        let rho = self.config.density;
+        let lp = self.pressure_face_gradients(pc);
+        for f in 0..mesh.faces().len() {
+            field.flux[f] -= dt / rho * lp[f];
+        }
+        for c in 0..mesh.cell_count() {
+            let g = self.pressure_gradient(pc, c);
+            field.u[c] -= dt / rho * g[0];
+            field.v[c] -= dt / rho * g[1];
+            field.p[c] += pc[c];
+        }
+    }
+
+    /// Largest cell mass imbalance `|Σ sign F_f|`.
+    pub(super) fn max_divergence(&self, flux: &[f64]) -> f64 {
+        let mesh = &self.mesh;
+        (0..mesh.cell_count())
+            .map(|c| {
+                mesh.cell_faces(c)
+                    .iter()
+                    .map(|&(f, sign)| sign * flux[f])
+                    .sum::<f64>()
+                    .abs()
+            })
+            .fold(0.0, f64::max)
+    }
+}
@@ -13,6 +13,8 @@ use crate::{CfdConfig, CfdError, CfdResult};
 pub mod ale;
 /// Boundary conditions
 pub mod boundary_conditions;
+/// Collocated PISO on a structured curvilinear patch (the overset patch)
+pub mod curvilinear;
 /// PISO on the fixed grid with an embedded body
 pub mod embedded;
 /// Embedded-body geometry, classification and loads
@@ -32,6 +34,8 @@ pub mod simple;
 /// GPU-accelerated SIMPLE algorithm implementation
 #[cfg(feature = "cuda")]
 pub mod simple_gpu;
+/// CSR matrix + Jacobi-BiCGSTAB for the curvilinear pressure equation
+pub mod sparse_bicgstab;

 // Re-export main types
 pub use ale::{
@@ -40,6 +44,10 @@ pub use ale::{
 pub use boundary_conditions::{
    BoundaryCondition, BoundaryConditions, BoundaryLocation, BoundaryType,
 };
+pub use curvilinear::{
+    CurvilinearParameters, CurvilinearPisoSolver, CurvilinearResult, CurvilinearSolverState,
+    NormalDiffusion, Operators, PatchBoundaries, PatchConvection, PatchField, SideBc,
+};
 pub use embedded::{EmbeddedParameters, EmbeddedPisoSolver, EmbeddedResult, EmbeddedSolverState};
 pub use embedded_body::{
    EmbeddedBody, EmbeddedMask, FaceKind, SurfaceForce, SurfaceSample, polygon_interface_velocity,
@@ -0,0 +1,315 @@
+//! A small CSR matrix and a Jacobi-preconditioned BiCGSTAB for the
+//! curvilinear patch's pressure equation (`docs/overset_metal_campaign.md`
+//! §5.3): the non-orthogonal operator is not symmetric, the patch is a few
+//! thousand cells, and `solve_multigrid_pcg` is hard-wired to the
+//! five-point Cartesian stencil.
+
+/// Compressed sparse rows.
+#[derive(Debug, Clone)]
+pub struct CsrMatrix {
+    n: usize,
+    row_ptr: Vec<usize>,
+    col: Vec<usize>,
+    val: Vec<f64>,
+}
+
+impl CsrMatrix {
+    /// Build from `(row, col, value)` triplets; duplicates accumulate,
+    /// columns are sorted within each row.
+    pub fn from_triplets(n: usize, triplets: &[(usize, usize, f64)]) -> Self {
+        let mut rows: Vec<Vec<(usize, f64)>> = vec![Vec::new(); n];
+        for &(r, c, v) in triplets {
+            rows[r].push((c, v));
+        }
+        let mut row_ptr = Vec::with_capacity(n + 1);
+        let mut col = Vec::with_capacity(triplets.len());
+        let mut val = Vec::with_capacity(triplets.len());
+        row_ptr.push(0);
+        for row in rows.iter_mut() {
+            row.sort_by_key(|e| e.0);
+            let mut last: Option<usize> = None;
+            for &(c, v) in row.iter() {
+                if last == Some(c) {
+                    *val.last_mut().expect("entry") += v;
+                } else {
+                    col.push(c);
+                    val.push(v);
+                    last = Some(c);
+                }
+            }
+            row_ptr.push(col.len());
+        }
+        Self {
+            n,
+            row_ptr,
+            col,
+            val,
+        }
+    }
+
+    /// Dimension.
+    pub fn n(&self) -> usize {
+        self.n
+    }
+    /// Stored entries.
+    pub fn nnz(&self) -> usize {
+        self.val.len()
+    }
+    /// `y = A x`.
+    pub fn matvec(&self, x: &[f64], y: &mut [f64]) {
+        for r in 0..self.n {
+            let mut acc = 0.0;
+            for k in self.row_ptr[r]..self.row_ptr[r + 1] {
+                acc += self.val[k] * x[self.col[k]];
+            }
+            y[r] = acc;
+        }
+    }
+    /// The diagonal (zero where absent).
+    pub fn diagonal(&self) -> Vec<f64> {
+        let mut d = vec![0.0; self.n];
+        for r in 0..self.n {
+            for k in self.row_ptr[r]..self.row_ptr[r + 1] {
+                if self.col[k] == r {
+                    d[r] = self.val[k];
+                }
+            }
+        }
+        d
+    }
+    /// Replace row `r` by the identity row (`x_r = b_r`): the anchor of a
+    /// pure-Neumann problem.
+    pub fn set_row_identity(&mut self, r: usize) {
+        let mut has_diag = false;
+        for k in self.row_ptr[r]..self.row_ptr[r + 1] {
+            self.val[k] = if self.col[k] == r {
+                has_diag = true;
+                1.0
+            } else {
+                0.0
+            };
+        }
+        assert!(has_diag, "row {r} has no diagonal entry to anchor");
+    }
+    /// Row `r` as `(columns, values)`.
+    pub fn row(&self, r: usize) -> (&[usize], &[f64]) {
+        let (a, b) = (self.row_ptr[r], self.row_ptr[r + 1]);
+        (&self.col[a..b], &self.val[a..b])
+    }
+}
+
+/// Outcome of a BiCGSTAB solve.
+#[derive(Debug, Clone, Copy)]
+#[must_use]
+pub struct BicgstabResult {
+    /// Iterations taken.
+    pub iterations: usize,
+    /// L1 norm of the true residual `b − A x` at exit.
+    pub residual: f64,
+    /// Whether the residual reached the tolerance.
+    pub converged: bool,
+}
+
+/// Subtract the mean from `v` (the consistency projection for a singular
+/// right-hand side).
+pub fn project_mean(v: &mut [f64]) {
+    let mean = v.iter().sum::<f64>() / v.len() as f64;
+    for x in v.iter_mut() {
+        *x -= mean;
+    }
+}
+
+/// Jacobi-preconditioned BiCGSTAB (van der Vorst 1992) with an L1
+/// true-residual stop; `x` is the initial guess and the result.
+pub fn bicgstab_jacobi(
+    a: &CsrMatrix,
+    b: &[f64],
+    x: &mut [f64],
+    tolerance: f64,
+    max_iterations: usize,
+) -> BicgstabResult {
+    let n = a.n();
+    let diag = a.diagonal();
+    let inv_diag: Vec<f64> = diag
+        .iter()
+        .map(|&d| if d != 0.0 { 1.0 / d } else { 1.0 })
+        .collect();
+    let l1 = |v: &[f64]| v.iter().map(|t| t.abs()).sum::<f64>();
+    let dot = |u: &[f64], v: &[f64]| u.iter().zip(v).map(|(p, q)| p * q).sum::<f64>();
+
+    let mut r = vec![0.0; n];
+    a.matvec(x, &mut r);
+    for i in 0..n {
+        r[i] = b[i] - r[i];
+    }
+    let mut res = l1(&r);
+    if res <= tolerance {
+        return BicgstabResult {
+            iterations: 0,
+            residual: res,
+            converged: true,
+        };
+    }
+    let r0 = r.clone();
+    let mut p = vec![0.0; n];
+    let mut v = vec![0.0; n];
+    let mut s = vec![0.0; n];
+    let mut t = vec![0.0; n];
+    let mut y = vec![0.0; n];
+    let mut z = vec![0.0; n];
+    let (mut rho_old, mut alpha, mut omega) = (1.0, 1.0, 1.0);
+
+    for it in 1..=max_iterations {
+        let rho = dot(&r0, &r);
+        if rho == 0.0 || !rho.is_finite() {
+            break;
+        }
+        let beta = (rho / rho_old) * (alpha / omega);
+        for i in 0..n {
+            p[i] = r[i] + beta * (p[i] - omega * v[i]);
+        }
+        for i in 0..n {
+            y[i] = inv_diag[i] * p[i];
+        }
+        a.matvec(&y, &mut v);
+        let r0v = dot(&r0, &v);
+        if r0v == 0.0 || !r0v.is_finite() {
+            break;
+        }
+        alpha = rho / r0v;
+        for i in 0..n {
+            s[i] = r[i] - alpha * v[i];
+        }
+        if l1(&s) <= tolerance {
+            for i in 0..n {
+                x[i] += alpha * y[i];
+            }
+            a.matvec(x, &mut r);
+            for i in 0..n {
+                r[i] = b[i] - r[i];
+            }
+            res = l1(&r);
+            return BicgstabResult {
+                iterations: it,
+                residual: res,
+                converged: res <= tolerance,
+            };
+        }
+        for i in 0..n {
+            z[i] = inv_diag[i] * s[i];
+        }
+        a.matvec(&z, &mut t);
+        let tt = dot(&t, &t);
+        omega = if tt > 0.0 { dot(&t, &s) / tt } else { 0.0 };
+        for i in 0..n {
+            x[i] += alpha * y[i] + omega * z[i];
+            r[i] = s[i] - omega * t[i];
+        }
+        rho_old = rho;
+        // The recurrence residual drifts from the true one; check the true
+        // residual whenever the recurrence claims convergence.
+        if l1(&r) <= tolerance {
+            a.matvec(x, &mut r);
+            for i in 0..n {
+                r[i] = b[i] - r[i];
+            }
+            res = l1(&r);
+            if res <= tolerance {
+                return BicgstabResult {
+                    iterations: it,
+                    residual: res,
+                    converged: true,
+                };
+            }
+        }
+        if omega == 0.0 {
+            break;
+        }
+    }
+    a.matvec(x, &mut r);
+    for i in 0..n {
+        r[i] = b[i] - r[i];
+    }
+    res = l1(&r);
+    BicgstabResult {
+        iterations: max_iterations,
+        residual: res,
+        converged: res <= tolerance,
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn a_small_nonsymmetric_system_is_solved_to_rounding() {
+        // Diagonally dominant, non-symmetric.
+        let dense = [
+            [4.0, -1.0, 0.0, 0.5, 0.0],
+            [-0.5, 5.0, -1.0, 0.0, 0.2],
+            [0.0, -1.5, 6.0, -1.0, 0.0],
+            [0.1, 0.0, -1.0, 4.0, -1.0],
+            [0.0, 0.3, 0.0, -0.5, 3.0],
+        ];
+        let mut tri = Vec::new();
+        for (r, row) in dense.iter().enumerate() {
+            for (c, &v) in row.iter().enumerate() {
+                if v != 0.0 {
+                    tri.push((r, c, v));
+                }
+            }
+        }
+        // Duplicate entry: must accumulate.
+        tri.push((0, 1, -0.5));
+        tri.push((0, 1, 0.5));
+        let a = CsrMatrix::from_triplets(5, &tri);
+        assert_eq!(a.nnz(), 17);
+        let x_true = [1.0, -2.0, 3.0, 0.5, -1.5];
+        let mut b = vec![0.0; 5];
+        a.matvec(&x_true, &mut b);
+        let mut x = vec![0.0; 5];
+        let out = bicgstab_jacobi(&a, &b, &mut x, 1e-13, 100);
+        assert!(out.converged, "{out:?}");
+        for i in 0..5 {
+            assert!(
+                (x[i] - x_true[i]).abs() < 1e-11,
+                "x[{i}] = {} vs {}",
+                x[i],
+                x_true[i]
+            );
+        }
+    }
+
+    #[test]
+    fn an_anchored_periodic_laplacian_recovers_a_periodic_field_up_to_a_constant() {
+        // 1-D periodic second difference (singular): anchor one row, project
+        // the mean out of the right-hand side.
+        let n = 64;
+        let h = 1.0 / n as f64;
+        let mut tri = Vec::new();
+        for i in 0..n {
+            tri.push((i, i, 2.0 / (h * h)));
+            tri.push((i, (i + 1) % n, -1.0 / (h * h)));
+            tri.push((i, (i + n - 1) % n, -1.0 / (h * h)));
+        }
+        let mut a = CsrMatrix::from_triplets(n, &tri);
+        let phi = |x: f64| (2.0 * std::f64::consts::PI * x).sin();
+        let exact: Vec<f64> = (0..n).map(|i| phi((i as f64 + 0.5) * h)).collect();
+        let mut b = vec![0.0; n];
+        a.matvec(&exact, &mut b);
+        b[3] += 1e-3; // an inconsistent perturbation the projection must remove
+        project_mean(&mut b);
+        a.set_row_identity(0);
+        b[0] = 0.0;
+        let mut x = vec![0.0; n];
+        // 1e-9 absolute: |A| ~ 1/h² = 4e3 and |x| ~ 1 put the rounding floor near 1e-10.
+        let out = bicgstab_jacobi(&a, &b, &mut x, 1e-9, 2000);
+        assert!(out.converged, "{out:?}");
+        let shift = exact[0] - x[0];
+        let worst = (0..n)
+            .map(|i| (x[i] + shift - exact[i]).abs())
+            .fold(0.0, f64::max);
+        assert!(worst < 2e-3, "worst {worst:.3e}"); // the 1e-3 perturbation's response bounds it
+    }
+}