//! 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.
//!
//! Phase P1: the patch MOVES. [`CurvilinearPisoSolver::set_mesh`] names the
//! geometry at the end of the coming step; the step then runs in the
//! conservative ALE form `(V^{n+1} û − V^n u^n)/dt = …` with the
//! convecting flux relative to the mesh, `F_f − δV_f/dt`, on the
//! time-averaged face vectors `S̄_f` (`motion.rs`: the 2-D discrete GCL is
//! exact algebra under that rule). `L_f`, the cell gradients and the
//! pressure matrix live on the end-of-step geometry; the projection is the
//! static one (no mesh-velocity term). A stationary mesh through this path
//! is bit-identical to the static path.

mod balance;
mod motion;
mod operators;
mod predictor;
mod projection;

pub use balance::PatchBalance;
pub use motion::StepGeometry;
pub use operators::Operators;

use crate::mesh::{PatchMesh, PatchSide};
use crate::solvers::incompressible::ale::SweptFaceRule;
use crate::solvers::incompressible::sparse_bicgstab::CsrMatrix;
use crate::{CfdConfig, CfdError, 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,
    /// Deferred-correction TVD with the van Albada limiter (the harness's
    /// scheme): the upwind face value plus `w_up psi(r) (phi_dn − phi_up)`,
    /// `w_up` the mesh's linear weight of the downwind side and `r` the
    /// ratio of the two one-sided gradients (so a linear field on a
    /// stretched row gives `r = 1` and the mesh's own linear face value).
    /// Faces whose far-upwind cell lies outside the patch fall back to
    /// upwind. Explicit, like the rest of the predictor's convection.
    TvdVanAlbada,
}

/// 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: each projection reduces the L1 cell mass
    /// imbalance by this factor (plus a rounding floor of 1e-15 × Σ|F|).
    /// Relative to the incoming divergence, not to the flux scale, so the
    /// stop tightens as the flow settles — an absolute stop leaves a
    /// velocity-noise floor that grows with the grid (measured: |du/dt|
    /// floored at 2e-4 on the 64² Cartesian MMS with `1e-10 × Σ|F|`).
    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,
    /// Face-area rule on a moving mesh: `Trapezoidal` (the DGCL-exact
    /// choice, the default) or `EndOfStep` (the negative control, GCL-
    /// violating). Irrelevant when the mesh does not move.
    pub swept_face_rule: SweptFaceRule,
}

impl Default for CurvilinearParameters {
    fn default() -> Self {
        Self {
            corrector_steps: 2,
            tolerance: 1e-4,
            convection: PatchConvection::Upwind,
            normal_diffusion: NormalDiffusion::Explicit,
            boundaries: PatchBoundaries::default(),
            max_poisson_iterations: 5000,
            swept_face_rule: SweptFaceRule::Trapezoidal,
        }
    }
}

/// 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): the time, the
/// current mesh and any mesh already named for the next step. Restoring
/// rebuilds the operators and drops the cached matrix, so a re-run from
/// the snapshot is bit-identical to the first run.
#[derive(Debug, Clone)]
pub struct CurvilinearSolverState {
    time: f64,
    mesh: PatchMesh,
    pending: Option<PatchMesh>,
}

type VelocityFn = Box<dyn Fn(f64, f64, f64) -> (f64, f64) + Send + Sync>;

/// The acceptor ring of an overset patch (A-P2): the outer row of cells
/// (`k = nn − 1`) carries values stamped from the background — `u, v, p`
/// at the end of every step, the pressure correction `p'` during each
/// projection — and no momentum or continuity equation of its own.
#[derive(Debug, Clone)]
struct AcceptorRing {
    /// Dirichlet `p'` per acceptor cell (column order), for the next solve.
    correction: Vec<f64>,
}

/// PISO on a curvilinear patch.
pub struct CurvilinearPisoSolver {
    config: CfdConfig,
    params: CurvilinearParameters,
    mesh: PatchMesh,
    /// The mesh the next `advance` ends on (`set_mesh`), if it moves.
    pending: Option<PatchMesh>,
    ops: Operators,
    boundary_velocity: Option<VelocityFn>,
    /// Per-side overrides of `boundary_velocity` (Inner, Outer, SStart, SEnd).
    side_velocity: [Option<VelocityFn>; 4],
    momentum_source: Option<VelocityFn>,
    acceptors: Option<AcceptorRing>,
    time: f64,
    matrix: Option<PressureSystem>,
    /// The Robin wall on the Inner side, if any.
    robin: Option<RobinWall>,
    /// Global face index of every Inner face (Inner-face order).
    robin_faces: Vec<usize>,
    /// `(t_f − datum) / alpha` per Inner face for the current step.
    robin_offset: Vec<[f64; 2]>,
}

/// A Robin wall on the `Inner` side (P6-b, the coupler with the added
/// mass built in — `docs/overset_metal_campaign.md` §5.19 in omni-cortex):
/// the wall velocity is the prescribed one plus `(t_f − datum) / alpha`,
/// with `t_f` the fluid's traction on the body (the [`Self::wall_tractions`]
/// convention) and `datum` the traction the structure was loaded with —
/// a wall of impedance `alpha` (Pa·s/m) that recedes when the fluid
/// pushes harder than the structure expects. Explicit in the predictor
/// (the start-of-step traction); IMPLICIT in the pressure: the wall flux
/// answers the pressure correction with `|S| p' / alpha` (a compliant
/// wall), which is the term that carries the added-mass operator into
/// the fluid's own response per subiterate. At the coupled fixed point
/// `t_f = datum` and the wall is the Dirichlet one.
#[derive(Debug, Clone)]
pub struct RobinWall {
    /// Impedance, Pa·s/m (`ρ_s h_s / Δt` for a plate of thickness `h_s`).
    pub alpha: f64,
    /// The structure's traction per Inner face, in Inner-face order (the
    /// order [`CurvilinearPisoSolver::wall_tractions`] returns).
    pub datum: Vec<[f64; 2]>,
}

/// The assembled pressure-correction system for one `dt` and geometry.
pub(crate) struct PressureSystem {
    /// The `dt` it was assembled for.
    pub(crate) dt: f64,
    /// `−Σ_f sign (dt/ρ) L_f` with identity rows on acceptor cells.
    pub(crate) matrix: CsrMatrix,
    /// The anchored cell of a pure-Neumann patch.
    pub(crate) anchor: Option<usize>,
    /// `(row, acceptor cell, coefficient)`: the interior rows' couplings to
    /// acceptor cells, eliminated to the right-hand side at solve time so
    /// the residual stays in flux units.
    pub(crate) links: Vec<(usize, usize, f64)>,
}

fn side_index(side: PatchSide) -> usize {
    match side {
        PatchSide::Inner => 0,
        PatchSide::Outer => 1,
        PatchSide::SStart => 2,
        PatchSide::SEnd => 3,
    }
}

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,
            pending: None,
            ops,
            boundary_velocity: None,
            side_velocity: [None, None, None, None],
            momentum_source: None,
            acceptors: None,
            time: 0.0,
            matrix: None,
            robin: None,
            robin_faces: Vec::new(),
            robin_offset: Vec::new(),
        })
    }

    /// Put a [`RobinWall`] on the Inner side (`datum` per Inner face, in
    /// [`Self::wall_tractions`] order) or replace its datum; the pressure
    /// matrix is rebuilt. `None` restores the Dirichlet wall bit for bit.
    pub fn set_robin_wall(&mut self, wall: Option<RobinWall>) {
        let faces: Vec<usize> = (0..self.mesh.faces().len())
            .filter(|&f| self.mesh.side(f) == Some(PatchSide::Inner))
            .collect();
        if let Some(w) = &wall {
            assert_eq!(
                w.datum.len(),
                faces.len(),
                "Robin datum must have one entry per Inner face"
            );
        }
        if self.robin_offset.len() != faces.len() {
            self.robin_offset = vec![[0.0, 0.0]; faces.len()];
        }
        self.robin_faces = faces;
        self.robin = wall;
        self.matrix = None;
    }
    /// The Robin wall's current velocity offsets per Inner face.
    pub fn robin_offsets(&self) -> &[[f64; 2]] {
        &self.robin_offset
    }
    /// `(t_f − datum) / alpha` on every Inner face from the field's current
    /// tractions (the explicit part of the Robin wall), called at the start
    /// of a step; with no Robin wall the offsets stay zero.
    fn refresh_robin_offsets(&mut self, field: &PatchField, t: f64) {
        let Some(w) = &self.robin else { return };
        let alpha = w.alpha;
        let datum = w.datum.clone();
        let mu = self.config.viscosity;
        let tractions = self.wall_tractions(field, PatchSide::Inner, t);
        let faces: Vec<usize> = self.robin_faces.clone();
        for (j, (_, _, _, tf)) in tractions.iter().enumerate() {
            // The wall's own viscous stress answers the wall velocity as
            // μ/d; taken implicitly in the update (a plain explicit
            // (t_f − datum)/α has gain (μ/d)/α and blew up at α = μ/h), the
            // fixed point unchanged: offset = (t_f − datum)/α.
            let f = faces[j];
            let c = self.mesh.boundary_cell(f);
            let xc = self.mesh.centre(c);
            let xf = self.mesh.faces()[f].centre;
            let d = ((xf[0] - xc[0]).powi(2) + (xf[1] - xc[1]).powi(2))
                .sqrt()
                .max(1e-300);
            let g = mu / (alpha * d);
            let old = self.robin_offset[j];
            self.robin_offset[j] = [
                ((tf[0] - datum[j][0]) / alpha + g * old[0]) / (1.0 + g),
                ((tf[1] - datum[j][1]) / alpha + g * old[1]) / (1.0 + g),
            ];
        }
    }
    /// The Robin offset at a point of the Inner side (the nearest face).
    fn robin_offset_at(&self, x: f64, y: f64) -> (f64, f64) {
        let mut best = (f64::INFINITY, [0.0, 0.0]);
        for (j, &f) in self.robin_faces.iter().enumerate() {
            let c = self.mesh.faces()[f].centre;
            let d = (c[0] - x).powi(2) + (c[1] - y).powi(2);
            if d < best.0 {
                best = (d, self.robin_offset[j]);
            }
        }
        (best.1[0], best.1[1])
    }

    /// 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));
    }
    /// Velocity on ONE `Velocity` side, overriding [`Self::set_boundary_velocity`]
    /// there (the overset patch: the wall on `Inner`, the background on
    /// `Outer`).
    pub fn set_side_velocity<F>(&mut self, side: PatchSide, f: F)
    where
        F: Fn(f64, f64, f64) -> (f64, f64) + Send + Sync + 'static,
    {
        self.side_velocity[side_index(side)] = Some(Box::new(f));
    }
    /// Turn the outer row of cells into acceptors (see `AcceptorRing`) or
    /// back into ordinary cells. The pressure matrix is rebuilt.
    pub fn set_acceptor_ring(&mut self, on: bool) {
        self.acceptors = on.then(|| AcceptorRing {
            correction: vec![0.0; self.mesh.ns()],
        });
        self.matrix = None;
    }
    /// Whether the outer row is an acceptor ring.
    pub fn has_acceptor_ring(&self) -> bool {
        self.acceptors.is_some()
    }
    /// Is cell `c` an acceptor (no equation of its own)?
    pub fn is_acceptor(&self, c: usize) -> bool {
        self.acceptors.is_some() && self.mesh.cell_ki(c).0 == self.mesh.nn() - 1
    }
    /// Stamp `(u, v, p)` onto the acceptor cells, column order `i = 0..ns`.
    pub fn stamp_acceptors(&self, field: &mut PatchField, values: &[(f64, f64, f64)]) {
        let nn = self.mesh.nn();
        for (i, &(u, v, p)) in values.iter().enumerate() {
            let c = self.mesh.cell(nn - 1, i);
            field.u[c] = u;
            field.v[c] = v;
            field.p[c] = p;
        }
    }
    /// The Dirichlet `p'` of the acceptor cells for the next
    /// [`Self::solve_correction`] (column order).
    pub fn set_acceptor_correction(&mut self, values: &[f64]) {
        if let Some(ring) = &mut self.acceptors {
            ring.correction.clear();
            ring.correction.extend_from_slice(values);
        }
    }
    /// One pressure-correction SOLVE (no application): `p'` on every cell
    /// (acceptor rows hold their Dirichlet values), with the BiCGSTAB
    /// report. `None` when the incoming divergence is already at the
    /// rounding floor. For the overset's Schwarz rounds.
    pub(crate) fn solve_correction(
        &self,
        field: &PatchField,
        dt: f64,
    ) -> Option<(
        Vec<f64>,
        crate::solvers::incompressible::sparse_bicgstab::BicgstabResult,
    )> {
        let system = self
            .matrix
            .as_ref()
            .expect("begin_step assembled the matrix");
        debug_assert_eq!(system.dt, dt);
        let flux_scale: f64 = field.flux.iter().map(|f| f.abs()).sum::<f64>().max(1e-300);
        let floor = 1e-15 * flux_scale;
        let (rhs, incoming) = self.pressure_rhs(system, &field.flux);
        if incoming <= floor {
            return None;
        }
        let tolerance = self.params.tolerance * incoming + floor;
        Some(self.solve_pressure_correction(system, rhs, tolerance))
    }
    /// Apply a correction `pc` (from [`Self::solve_correction`]).
    pub(crate) fn apply_correction_pub(&self, field: &mut PatchField, pc: &[f64], dt: f64) {
        self.apply_correction(field, pc, dt);
    }
    /// Largest cell mass imbalance over the equation-carrying cells.
    pub(crate) fn max_divergence_pub(&self, flux: &[f64]) -> f64 {
        self.max_divergence(flux)
    }
    /// Overlap mass defect on the patch side: `Σ_acceptors |Σ_f sign F_f|`
    /// (the acceptors carry no continuity) with the OUTER face flux taken
    /// from `outer_velocity` (the background's velocity at that face
    /// centre, column order), and the flux scale `Σ |F_f|` over the faces
    /// between the acceptor ring and the interior.
    pub fn acceptor_mass_defect(
        &self,
        field: &PatchField,
        outer_velocity: &[(f64, f64)],
    ) -> (f64, f64) {
        let mesh = &self.mesh;
        let nn = mesh.nn();
        if self.acceptors.is_none() || nn < 2 {
            return (0.0, 0.0);
        }
        let mut defect = 0.0;
        let mut scale = 0.0;
        for i in 0..mesh.ns() {
            let c = mesh.cell(nn - 1, i);
            let outer = mesh.nface(nn, i);
            let mut div = 0.0;
            for (f, sign) in mesh.cell_faces(c) {
                if f == outer {
                    let s = mesh.faces()[f].s;
                    let (uo, vo) = outer_velocity.get(i).copied().unwrap_or((0.0, 0.0));
                    div += sign * (uo * s[0] + vo * s[1]);
                } else {
                    div += sign * field.flux[f];
                }
            }
            defect += div.abs();
            scale += field.flux[mesh.nface(nn - 1, i)].abs();
        }
        (defect, scale)
    }
    /// 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 the field currently lives on (the end of the last step).
    pub fn mesh(&self) -> &PatchMesh {
        &self.mesh
    }
    /// The mesh named for the end of the next step, if any.
    pub fn next_mesh(&self) -> Option<&PatchMesh> {
        self.pending.as_ref()
    }
    /// Name the geometry the NEXT `advance` ends on. Same topology as the
    /// current mesh (`ns`, `nn`, periodicity); the node motion between the
    /// two is taken as linear in time. Calling it again before `advance`
    /// replaces the earlier choice (the coupling loop re-tries a step);
    /// not calling it leaves the mesh where it is.
    pub fn set_mesh(&mut self, next: PatchMesh) -> CfdResult<()> {
        if next.ns() != self.mesh.ns()
            || next.nn() != self.mesh.nn()
            || next.periodic().is_some() != self.mesh.periodic().is_some()
        {
            return Err(CfdError::mesh(format!(
                "set_mesh: topology changed ({}x{}, periodic {}) -> ({}x{}, periodic {})",
                self.mesh.ns(),
                self.mesh.nn(),
                self.mesh.periodic().is_some(),
                next.ns(),
                next.nn(),
                next.periodic().is_some()
            )));
        }
        self.pending = Some(next);
        Ok(())
    }
    /// 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 (time and both meshes).
    pub fn snapshot(&self) -> CurvilinearSolverState {
        CurvilinearSolverState {
            time: self.time,
            mesh: self.mesh.clone(),
            pending: self.pending.clone(),
        }
    }
    /// Restore a captured state: operators rebuilt, matrix cache dropped.
    pub fn restore(&mut self, state: &CurvilinearSolverState) {
        self.time = state.time;
        self.mesh = state.mesh.clone();
        self.pending = state.pending.clone();
        self.ops = Operators::new(&self.mesh, &self.params.boundaries);
        self.matrix = None;
    }

    pub(crate) fn boundary_velocity(&self, side: PatchSide, x: f64, y: f64, t: f64) -> (f64, f64) {
        let (u, v) = self.side_velocity[side_index(side)]
            .as_ref()
            .or(self.boundary_velocity.as_ref())
            .map_or((0.0, 0.0), |f| f(x, y, t));
        if side == PatchSide::Inner && self.robin.is_some() {
            let (ou, ov) = self.robin_offset_at(x, y);
            (u + ou, v + ov)
        } else {
            (u, v)
        }
    }
    /// Dirichlet `p'` of acceptor cell `c`, if it is one.
    pub(crate) fn acceptor_correction(&self, c: usize) -> Option<f64> {
        let ring = self.acceptors.as_ref()?;
        let (k, i) = self.mesh.cell_ki(c);
        (k == self.mesh.nn() - 1).then(|| ring.correction.get(i).copied().unwrap_or(0.0))
    }
    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()];
        let geo = StepGeometry::stationary(mesh);
        field.flux = self.predicted_fluxes(
            &field.u.clone(),
            &field.v.clone(),
            &zero,
            0.0,
            self.time,
            &geo,
        );
    }

    /// Advance one step of `dt`: [`Self::begin_step`], the correctors,
    /// [`Self::end_step`].
    pub async fn advance(
        &mut self,
        field: &mut PatchField,
        dt: f64,
    ) -> CfdResult<CurvilinearResult> {
        let start = self.begin_step(field, dt)?;
        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 Some((pc, out)) = self.solve_correction(field, dt) else {
                break;
            };
            iterations += out.iterations;
            converged &= out.converged;
            self.apply_correction(field, &pc, dt);
            performed += 1;
            max_div = self.max_divergence(&field.flux);
        }
        Ok(self.end_step(&start, performed, max_div, iterations, converged))
    }

    /// Everything before the correctors: the mesh swap (if `set_mesh` named
    /// one), the step geometry, the predictor, the predicted fluxes with
    /// the closed-patch adjustment, `u*`, and the pressure matrix on the
    /// end-of-step geometry. The overset coupling runs this on the patch,
    /// then drives the correctors itself.
    pub(crate) fn begin_step(&mut self, field: &mut PatchField, dt: f64) -> CfdResult<StepStart> {
        let t_old = self.time;
        let t_new = t_old + dt;
        let rho = self.config.density;

        // The mesh moves: `self.mesh` becomes the end-of-step geometry
        // (operators and matrix follow it); the start-of-step geometry is
        // kept for this step's volumes and swept faces only.
        let old = match self.pending.take() {
            Some(next) => {
                let old = std::mem::replace(&mut self.mesh, next);
                self.ops = Operators::new(&self.mesh, &self.params.boundaries);
                self.matrix = None;
                Some(old)
            }
            None => None,
        };
        let geo = match &old {
            Some(o) => StepGeometry::new(o, &self.mesh, self.params.swept_face_rule),
            None => StepGeometry::stationary(&self.mesh),
        };
        if self.robin.is_some() {
            self.refresh_robin_offsets(field, t_old);
        }
        let old_mesh = old.as_ref().unwrap_or(&self.mesh);
        let mesh = &self.mesh;

        let (uh, vh) = self.predict(field, dt, t_old, old_mesh, &geo);
        let mut flux = self.predicted_fluxes(&uh, &vh, &field.p, dt, t_new, &geo);
        let adjustment = self.adjust_boundary_flux(&mut flux);
        for c in 0..mesh.cell_count() {
            if self.is_acceptor(c) {
                continue;
            }
            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(|s| s.dt != dt) {
            self.matrix = Some(self.assemble_pressure_matrix(dt));
        }
        Ok(StepStart { t_new, adjustment })
    }

    /// Everything after the correctors: the clock and the result.
    pub(crate) fn end_step(
        &mut self,
        start: &StepStart,
        performed: usize,
        max_div: f64,
        iterations: usize,
        converged: bool,
    ) -> CurvilinearResult {
        self.time = start.t_new;
        CurvilinearResult {
            corrector_steps_performed: performed,
            max_divergence: max_div,
            poisson_iterations: iterations,
            poisson_converged: converged,
            boundary_flux_adjustment: start.adjustment,
        }
    }
}

/// What [`CurvilinearPisoSolver::begin_step`] hands to
/// [`CurvilinearPisoSolver::end_step`].
#[derive(Debug, Clone, Copy)]
pub(crate) struct StepStart {
    /// End-of-step time.
    pub(crate) t_new: f64,
    /// Boundary-flux defect removed on a closed patch.
    pub(crate) adjustment: f64,
}

/// Fluid force on a boundary side of the patch (per unit depth).
#[derive(Debug, Clone, Copy)]
pub struct PatchLoad {
    /// Pressure part.
    pub pressure: [f64; 2],
    /// Viscous part.
    pub viscous: [f64; 2],
    /// Faces integrated.
    pub faces: usize,
}

impl PatchLoad {
    /// Total force.
    pub fn total(&self) -> [f64; 2] {
        [
            self.pressure[0] + self.viscous[0],
            self.pressure[1] + self.viscous[1],
        ]
    }
}

impl CurvilinearPisoSolver {
    /// The fluid force on the body along `side` (the wall of an O-grid:
    /// `Inner`): `F = Σ_f (−p_f S_f + μ (∇u + ∇uᵀ)_f · S_f)` with `S_f`
    /// pointing from the body into the fluid. The face pressure is the
    /// wall cell's pressure extrapolated linearly with its least-squares
    /// gradient; the face velocity gradient is the wall cell's least-squares
    /// gradient with the face's Dirichlet value in the fit (the wall shear
    /// enters through the wall point itself). No traction reconstruction
    /// through a staircase — the second thing the overset buys (§2.1).
    pub fn surface_force(&self, field: &PatchField, side: PatchSide, t: f64) -> PatchLoad {
        let (mut fp, mut fv) = ([0.0; 2], [0.0; 2]);
        let mut faces = 0usize;
        for (_, (p_f, s, tau)) in self.wall_face_terms(field, side, t) {
            fp[0] -= p_f * s[0];
            fp[1] -= p_f * s[1];
            fv[0] += tau[0];
            fv[1] += tau[1];
            faces += 1;
        }
        PatchLoad {
            pressure: fp,
            viscous: fv,
            faces,
        }
    }

    /// Per wall face of `side`: the face pressure, the area vector INTO the
    /// fluid, and the viscous force `μ (∇u + ∇uᵀ) · S` — the terms
    /// [`Self::surface_force`] sums. Iterates in face order.
    fn wall_face_terms<'a>(
        &'a self,
        field: &'a PatchField,
        side: PatchSide,
        t: f64,
    ) -> impl Iterator<Item = (usize, (f64, [f64; 2], [f64; 2]))> + 'a {
        let mesh = &self.mesh;
        let mu = self.config.viscosity;
        mesh.faces()
            .iter()
            .enumerate()
            .filter(move |(f, _)| mesh.side(*f) == Some(side))
            .map(move |(f, face)| {
                let c = mesh.boundary_cell(f);
                // Outward from the body = into the fluid: for the Inner side
                // (k = 0, the cell is the neighbour) that is +S; for the Outer
                // side (the cell is the owner) it is −S.
                let sign = if face.neigh.is_some() { 1.0 } else { -1.0 };
                let s = [sign * face.s[0], sign * face.s[1]];
                let xc = mesh.centre(c);
                let dxf = [face.centre[0] - xc[0], face.centre[1] - xc[1]];
                let gp = self.pressure_gradient(&field.p, c);
                let p_f = field.p[c] + gp[0] * dxf[0] + gp[1] * dxf[1];
                let wall = |ff: usize| -> Option<(f64, f64)> {
                    let fc = &mesh.faces()[ff];
                    let sd = mesh.side(ff)?;
                    match self.params.boundaries.get(sd) {
                        SideBc::Velocity => {
                            Some(self.boundary_velocity(sd, fc.centre[0], fc.centre[1], t))
                        }
                        SideBc::Outlet => None,
                    }
                };
                let gu = self
                    .ops
                    .gradient(mesh, c, &field.u, &|ff| wall(ff).map(|w| w.0));
                let gv = self
                    .ops
                    .gradient(mesh, c, &field.v, &|ff| wall(ff).map(|w| w.1));
                // τ = μ (∇u + ∇uᵀ): τxx = 2 u_x, τxy = u_y + v_x, τyy = 2 v_y.
                let tau = [
                    mu * (2.0 * gu[0] * s[0] + (gu[1] + gv[0]) * s[1]),
                    mu * ((gu[1] + gv[0]) * s[0] + 2.0 * gv[1] * s[1]),
                ];
                (f, (p_f, s, tau))
            })
    }

    /// The traction on every wall face of `side` (P5, the load transfer):
    /// `(face centre, unit normal into the fluid, face length, traction
    /// per unit length = (−p_f S + μ (∇u + ∇uᵀ) · S) / |S|)`, in face
    /// order — the same terms as [`Self::surface_force`].
    pub fn wall_tractions(
        &self,
        field: &PatchField,
        side: PatchSide,
        t: f64,
    ) -> Vec<([f64; 2], [f64; 2], f64, [f64; 2])> {
        let mesh = &self.mesh;
        self.wall_face_terms(field, side, t)
            .map(|(f, (p_f, s, tau))| {
                let len = (s[0] * s[0] + s[1] * s[1]).sqrt().max(1e-300);
                (
                    mesh.faces()[f].centre,
                    [s[0] / len, s[1] / len],
                    len,
                    [(-p_f * s[0] + tau[0]) / len, (-p_f * s[1] + tau[1]) / len],
                )
            })
            .collect()
    }
}