rustytorch/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded/mod.rs

//! PISO on the fixed uniform staggered grid with an embedded body.
//!
//! This is the fixed-grid PISO scheme (`piso.rs`: explicit momentum
//! predictor, SOR pressure-correction projections) with three additions:
//!
//! - **per-side domain boundaries** ([`SideBoundary`]: prescribed velocity,
//!   slip wall, pressure outlet), with exactly the ALE solver's semantics
//!   — the channel of the Turek–Hron benchmark needs an outlet and the
//!   fixed-grid PISO had none;
//! - a **time-dependent boundary-velocity function** `(x, y, t)` supplying
//!   the normal data and the tangential wall values, and a source
//!   `(x, y, t)`, as on the ALE solver;
//! - an optional **embedded body** ([`EmbeddedBody`]) classified onto the
//!   grid by [`EmbeddedMask`]: the momentum predictor updates only fluid
//!   faces, the projection enforces continuity only on fluid cells with
//!   zero coefficient across every prescribed face, and after each
//!   projection the ghost faces are re-imposed from the corrected fluid
//!   field (see `embedded_body.rs` for the reconstruction and the
//!   compatibility correction).
//!
//! With no body, velocity on every side and no normal flow through the
//! sides (a closed box), `advance` is the fixed-grid PISO step to the last
//! bit — `tests/embedded_mms.rs` pins that degeneracy before it measures
//! anything else. With through-flow the two differ by design: the
//! fixed-grid PISO zeroes the transverse convective face velocity on its
//! domain sides (exact for walls), this solver takes it from the stored
//! boundary faces, which is what an inlet or outlet needs — dropping the
//! outgoing flux at an outlet let the last column accumulate momentum and
//! the Turek–Hron channel blew up at t ≈ 4 s.
//!
//! # Boundary history
//!
//! As the ALE solver learned (the twelfth defect of the campaign), the
//! start-of-step boundary faces are *not* re-stamped from the boundary
//! function: the previous step's end-of-step application is the material
//! history the explicit predictor differentiates. Only the first step
//! stamps `t = 0` data, via [`EmbeddedPisoSolver::initialize`] or lazily.

mod projection;

use super::ale::{AleBoundaries, SideBoundary};
use super::embedded_body::{EmbeddedBody, EmbeddedMask, FaceKind};
use super::poisson::{
    MgPrecision, MultigridParameters, PoissonProblem, PoissonSolverKind, solve_multigrid_pcg,
};
use super::simple::ConvectionScheme;
use super::{FlowField, SolverResult};
use crate::{CfdConfig, CfdError, CfdResult};

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

/// Parameters for the embedded-boundary PISO solver.
#[derive(Debug, Clone)]
pub struct EmbeddedParameters {
    /// Projection passes per step.
    pub corrector_steps: usize,
    /// Convergence tolerance on the normalised mass imbalance after
    /// correction.
    pub tolerance: f64,
    /// Boundary type per domain side (all prescribed velocity by default).
    /// The struct is the ALE solver's; the semantics are identical.
    pub boundaries: AleBoundaries,
    /// Inner solver of the pressure-correction system (default
    /// [`PoissonSolverKind::Sor`]). Both solve the same system to the same
    /// true-residual stop; multigrid's cost is mesh-independent.
    pub poisson_solver: PoissonSolverKind,
    /// Precision of the multigrid V-cycle (default [`MgPrecision::F64`] =
    /// bit-identical). `F32` is the M1 precision probe of
    /// `overset_metal_campaign.md` §3.2 / §5.2: the CG, its true residual
    /// and the stop stay f64; only the preconditioner runs in single
    /// precision. No effect with [`PoissonSolverKind::Sor`].
    pub poisson_precision: MgPrecision,
    /// The multigrid smoother ordering (PERF-2; default lexicographic, the
    /// recorded regime).
    pub poisson_smoother: super::poisson::MgSmoother,
    /// Convective face values in the explicit predictor (default
    /// [`ConvectionScheme::Upwind`], which is bit-identical to the fixed-grid
    /// PISO). The TVD schemes add SIMPLE's limited correction to each
    /// interior face — with an explicit predictor no deferred iteration is
    /// needed, the limited flux is just used directly. First-order upwind's
    /// numerical viscosity `|u| h / 2` exceeds the physical viscosity ten
    /// times over on the Turek–Hron CFD3 grids and suppressed the vortex
    /// shedding entirely; the limited scheme restores it. Faces whose
    /// far-upwind node lies outside the domain, and domain-side faces, fall
    /// back to pure upwind exactly as in SIMPLE; near the body the stencil
    /// reads ghost values, which encode the wall.
    pub convection_scheme: ConvectionScheme,
}

impl Default for EmbeddedParameters {
    fn default() -> Self {
        Self {
            corrector_steps: 2,
            tolerance: 1e-6,
            boundaries: AleBoundaries::default(),
            poisson_solver: PoissonSolverKind::Sor,
            poisson_precision: MgPrecision::F64,
            poisson_smoother: super::poisson::MgSmoother::Lexicographic,
            convection_scheme: ConvectionScheme::Upwind,
        }
    }
}

/// A snapshot of [`EmbeddedPisoSolver`]'s per-step state, for re-running a
/// step within a coupling subiteration. See [`EmbeddedPisoSolver::snapshot`].
#[derive(Clone)]
pub struct EmbeddedSolverState {
    mask: Option<EmbeddedMask>,
    time: f64,
    initialized: bool,
    alpha: Option<Vec<f64>>,
    fringe: Option<Vec<bool>>,
}

/// Result of one embedded PISO step.
#[derive(Debug, Clone)]
pub struct EmbeddedResult {
    /// Base solver result information.
    pub solver_result: SolverResult,
    /// Number of projection passes performed.
    pub corrector_steps_performed: usize,
    /// The per-face compatibility correction applied to the ghost faces at
    /// the end of the step (velocity units); zero without a body.
    pub ghost_correction: f64,
    /// Pressure cells that flipped solid → fluid in this step's mask
    /// rebuild (always zero for a static body).
    pub fresh_cells: usize,
}

/// The embedded-boundary PISO solver. See the module docs.
pub struct EmbeddedPisoSolver {
    config: CfdConfig,
    parameters: EmbeddedParameters,
    momentum_source: Option<SourceFn>,
    /// Reporting-only: the cells that turned fluid on the current step,
    /// kept for the `RTX_EMBEDDED_TRACE_SP` divergence trace (empty
    /// unless the env var is set).
    fresh_trace: Vec<(usize, usize)>,
    /// REFUTED on the falsifier (2026-09-03): 40× larger spikes at either
    /// sign — the wall faces already carry the swept volume; kept as the
    /// record of that measurement, never to be enabled.
    /// Swept-volume source strength (0 = off, bit-identical; ±1 = on,
    /// sign as registered by the falsifier): the fluid area fraction of
    /// every interface cell, α = clamp(½ + φ/h, 0, 1) from the body's
    /// signed distance at the cell centre, enters the continuity
    /// constraint as a source ρ (α^{n+1} − α^n) dx dy / dt, so a cell's
    /// fluid volume enters continuously as the wall sweeps instead of as
    /// a whole-cell jump at the mask flip — the fixed impulse per flip
    /// the fresh-cell falsifier measured (omni-cortex
    /// `docs/fresh_cell_gcl_campaign.md`).
    swept_volume: f64,
    /// Field extension for fresh faces (knob, default off = bit-identical;
    /// measured NO EFFECT on the falsifier 2026-09-03 — kept as the record):
    /// see `EmbeddedMask::extend_fresh_faces`.
    field_extension: bool,
    /// The previous step's fluid area fractions (moving path, knob on).
    alpha_old: Option<Vec<f64>>,
    /// This step's fractions, computed at the mask rebuild.
    alpha_new: Option<Vec<f64>>,
    boundary_velocity: Option<VelocityFn>,
    body: Option<EmbeddedBody>,
    mask: Option<EmbeddedMask>,
    /// Overset fringe (A-P2): per-cell flag, and the Dirichlet `p'` the
    /// fringe cells carry in the next projection (row-major, full size).
    fringe: Option<Vec<bool>>,
    fringe_correction: Vec<f64>,
    /// Relative part of the pressure solve's inner stop (default 1e-2 =
    /// bit-identical with the record): each solve reduces the residual to
    /// this fraction of the continuity source. The overset's Schwarz rounds
    /// exchange the solution and need it accurate below their own stop
    /// (measured: at 1e-2 the first corrector's rounds never settled below
    /// a 1e-3 relative change — 20/20 rounds every step).
    inner_stop_factor: f64,
    /// PERF-2 P0: Poisson `(setup ns, iterate ns, calls, CG iterations)`
    /// summed over every multigrid-PCG solve (`docs/perf2_campaign.md`).
    poisson_profile: std::cell::Cell<(u64, u64, u64, u64)>,
    /// PERF-2 P1.1: the prepared multigrid operator, reused across rounds
    /// and steps while the operator's coefficients and mask are unchanged.
    pcg_cache: std::cell::RefCell<super::poisson::PcgCache>,
    moving: bool,
    /// Mask hysteresis band in multiples of the min cell size (0 = off).
    mask_hysteresis: f64,
    time: f64,
    initialized: bool,
}

impl EmbeddedPisoSolver {
    /// Create the solver.
    pub fn new(config: CfdConfig, parameters: EmbeddedParameters) -> CfdResult<Self> {
        config.validate()?;
        Ok(Self {
            config,
            parameters,
            momentum_source: None,
            fresh_trace: Vec::new(),
            swept_volume: 0.0,
            field_extension: false,
            alpha_old: None,
            alpha_new: None,
            boundary_velocity: None,
            body: None,
            mask: None,
            fringe: None,
            fringe_correction: Vec::new(),
            inner_stop_factor: 1e-2,
            poisson_profile: std::cell::Cell::new((0, 0, 0, 0)),
            pcg_cache: std::cell::RefCell::new(super::poisson::PcgCache::default()),
            moving: false,
            mask_hysteresis: 0.0,
            time: 0.0,
            initialized: false,
        })
    }

    /// Precision of the multigrid V-cycle (see
    /// [`EmbeddedParameters::poisson_precision`]); the M1 probe knob.
    pub fn set_poisson_precision(&mut self, precision: MgPrecision) {
        self.parameters.poisson_precision = precision;
    }

    /// The multigrid smoother ordering (PERF-2's red-black knob); the
    /// cached operator is rebuilt on the next solve.
    pub fn set_poisson_smoother(&mut self, smoother: super::poisson::MgSmoother) {
        self.parameters.poisson_smoother = smoother;
    }

    /// Mask hysteresis for the moving-body rebuild, as a fraction of the
    /// min cell size (default 0, exactly the plain rebuild). With a band,
    /// a cell within `band * h_min` of the surface keeps the
    /// classification it has in the mask held at rebuild time — in a
    /// coupling loop that restores a [`Self::snapshot`] before each
    /// subiteration, that is the committed step-start mask, so every pass
    /// of a step classifies against ONE reference and candidate geometries
    /// within the band all see the SAME mask (the pass map stops flipping
    /// cells on sub-band candidate differences). The cost is the effective
    /// wall lagging the true surface by up to the band.
    pub fn set_mask_hysteresis(&mut self, band_in_h: f64) {
        self.mask_hysteresis = band_in_h;
    }

    /// Volumetric momentum source `(x, y, t) -> (f_x, f_y)` per unit volume.
    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));
    }

    /// Boundary velocity `(x, y, t) -> (u, v)` on the domain sides: normal
    /// component prescribed on `Velocity` and `SlipWall` sides, tangential
    /// component the no-slip value on `Velocity` sides.
    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));
    }

    /// Embed a body, treated as fixed in shape and position: the mask is
    /// built once, on the first step.
    pub fn set_body(&mut self, body: EmbeddedBody) {
        self.body = Some(body);
        self.mask = None;
        self.moving = false;
    }

    /// Embed a body whose signed distance and surface velocity depend on
    /// time. The mask is rebuilt at the end-of-step time every step; the
    /// new mask's ghost values are reconstructed from the previous
    /// corrected field, so a *stationary* body run through this path is
    /// bit-identical to [`Self::set_body`]'s. A velocity face that flips
    /// solid → fluid (a *fresh* face) enters the new interval holding
    /// exactly the ghost reconstruction the previous step left on it —
    /// a consistent near-wall value, not garbage — and a fresh pressure
    /// cell is refilled from its fluid neighbours before the predictor's
    /// gradient can read its stale value. The body must move less than a
    /// cell per step (the convective time-step limit already enforces
    /// this for a body slower than the local peak velocity).
    pub fn set_moving_body(&mut self, body: EmbeddedBody) {
        self.body = Some(body);
        self.mask = None;
        self.moving = true;
    }

    /// The overset background (A-P2): a mask from the overlap
    /// classification (active cells fluid, prescribed faces `Ghost` with no
    /// reconstruction) and the fringe flags. Fringe cells carry no
    /// continuity equation; their `p'` is Dirichlet
    /// ([`Self::set_fringe_correction`]) and their `p` and prescribed faces
    /// are stamped by the caller from the patch. No body: nothing is
    /// re-imposed at the end of a step.
    pub fn set_overlap(&mut self, mask: EmbeddedMask, fringe: Vec<bool>) {
        self.fringe_correction = vec![0.0; fringe.len()];
        self.fringe = Some(fringe);
        self.mask = Some(mask);
        self.body = None;
        self.moving = false;
    }

    /// Dirichlet `p'` of the fringe cells for the next
    /// [`Self::solve_correction`] (`cells` as `(j, i)`).
    pub fn set_fringe_correction(&mut self, cells: &[(usize, usize)], values: &[f64]) {
        let nx = self.mask.as_ref().map_or(0, |m| m.nx());
        for (&(j, i), &v) in cells.iter().zip(values) {
            self.fringe_correction[j * nx + i] = v;
        }
    }

    /// Relative part of the pressure solve's inner stop (see the field).
    /// The Poisson solves' `(setup ns, iterate ns, calls, CG iterations)` so far.
    pub fn poisson_profile(&self) -> (u64, u64, u64, u64) {
        self.poisson_profile.get()
    }

    pub fn set_inner_stop_factor(&mut self, factor: f64) {
        self.inner_stop_factor = factor;
    }
    /// The inner-stop factor.
    pub fn inner_stop_factor(&self) -> f64 {
        self.inner_stop_factor
    }

    /// Whether an overset fringe is set.
    pub fn has_fringe(&self) -> bool {
        self.fringe.is_some()
    }

    #[inline]
    pub(crate) fn is_fringe(&self, j: usize, i: usize) -> bool {
        match (&self.fringe, &self.mask) {
            (Some(f), Some(m)) => f[j * m.nx() + i],
            _ => false,
        }
    }

    /// Dirichlet `p'` of fringe cell `(j, i)`.
    #[inline]
    pub(crate) fn fringe_correction(&self, j: usize, i: usize) -> f64 {
        let nx = self.mask.as_ref().map_or(0, |m| m.nx());
        self.fringe_correction[j * nx + i]
    }

    /// Write the fringe cells' Dirichlet `p'` into `p_prime` (the face
    /// corrections across active–fringe faces read it there).
    pub(crate) fn stamp_fringe_correction(&self, p_prime: &mut nalgebra::DMatrix<f64>) {
        if let (Some(f), Some(m)) = (&self.fringe, &self.mask) {
            let nx = m.nx();
            for (idx, &is_fringe) in f.iter().enumerate() {
                if is_fringe {
                    p_prime[(idx / nx, idx % nx)] = self.fringe_correction[idx];
                }
            }
        }
    }

    /// Field extension for faces that turn fluid (moving-body path).
    pub fn set_field_extension(&mut self, on: bool) {
        self.field_extension = on;
    }

    /// Swept-volume source strength for the moving-body path (0 = off).
    pub fn set_swept_volume_source(&mut self, strength: f64) {
        self.swept_volume = strength;
    }

    /// The body, if any.
    pub fn body(&self) -> Option<&EmbeddedBody> {
        self.body.as_ref()
    }

    /// The mask, once built (after `initialize` or the first step).
    pub fn mask(&self) -> Option<&EmbeddedMask> {
        self.mask.as_ref()
    }

    /// Accumulated time.
    pub fn time(&self) -> f64 {
        self.time
    }

    /// Snapshot of the solver's own per-step state — the mask, the
    /// accumulated time and the initialization flag. A coupling
    /// subiteration re-runs one step from the same start: clone the
    /// [`FlowField`], take this snapshot, and [`Self::restore`] both
    /// before every re-run — otherwise the moving-body path's fresh-cell
    /// detection compares against the *previous subiteration's* mask
    /// instead of the committed step-start mask.
    pub fn snapshot(&self) -> EmbeddedSolverState {
        EmbeddedSolverState {
            mask: self.mask.clone(),
            time: self.time,
            initialized: self.initialized,
            alpha: self.alpha_old.clone(),
            fringe: self.fringe.clone(),
        }
    }

    /// Restore a [`Self::snapshot`]. The snapshot is cloned, so one
    /// snapshot serves any number of re-runs.
    pub fn restore(&mut self, state: &EmbeddedSolverState) {
        self.mask = state.mask.clone();
        self.time = state.time;
        self.initialized = state.initialized;
        self.alpha_old = state.alpha.clone();
        self.fringe = state.fringe.clone();
        if let Some(f) = &self.fringe {
            self.fringe_correction = vec![0.0; f.len()];
        }
    }

    /// Reset the accumulated time.
    pub fn set_time(&mut self, t: f64) {
        self.time = t;
    }

    /// Solver configuration.
    pub fn config(&self) -> &CfdConfig {
        &self.config
    }

    /// Solver parameters.
    pub fn parameters(&self) -> &EmbeddedParameters {
        &self.parameters
    }

    fn boundary(&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))
    }

    /// Build the mask (if a body is set) and stamp the `t = time` boundary
    /// data and ghost values onto the field. Called lazily by the first
    /// `advance`; call it explicitly to inspect the mask or to run
    /// diagnostics on the initial field.
    pub fn initialize(&mut self, field: &mut FlowField) -> CfdResult<()> {
        let (nx, ny, dx, dy) = field.grid_info();
        if let Some(body) = &self.body {
            if self.mask.is_none() {
                self.mask = Some(EmbeddedMask::build(body, nx, ny, dx, dy, self.time)?);
            }
        }
        let t = self.time;
        self.apply_boundary_normals(field, t);
        if let (Some(body), Some(mask)) = (&self.body, &self.mask) {
            mask.impose(body, &mut field.u, &mut field.v, t);
        }
        self.initialized = true;
        Ok(())
    }

    /// Write the prescribed normal velocities onto the boundary faces at
    /// time `t`. Outlet faces are unknowns and are left alone.
    fn apply_boundary_normals(&self, field: &mut FlowField, t: f64) {
        let (nx, ny, dx, dy) = field.grid_info();
        let b = self.parameters.boundaries;
        let outlet = SideBoundary::PressureOutlet;
        for j in 0..ny {
            let y = (j as f64 + 0.5) * dy;
            if b.left != outlet {
                field.u[(j, 0)] = self.boundary(0.0, y, t).0;
            }
            if b.right != outlet {
                field.u[(j, nx)] = self.boundary(nx as f64 * dx, y, t).0;
            }
        }
        for i in 0..nx {
            let x = (i as f64 + 0.5) * dx;
            if b.bottom != outlet {
                field.v[(0, i)] = self.boundary(x, 0.0, t).1;
            }
            if b.top != outlet {
                field.v[(ny, i)] = self.boundary(x, ny as f64 * dy, t).1;
            }
        }
    }

    #[inline]
    fn u_is_fluid(&self, j: usize, i: usize) -> bool {
        self.mask
            .as_ref()
            .is_none_or(|m| m.u_kind(j, i) == FaceKind::Fluid)
    }

    #[inline]
    fn v_is_fluid(&self, j: usize, i: usize) -> bool {
        self.mask
            .as_ref()
            .is_none_or(|m| m.v_kind(j, i) == FaceKind::Fluid)
    }

    #[inline]
    fn cell_is_fluid(&self, j: usize, i: usize) -> bool {
        self.mask.as_ref().is_none_or(|m| m.is_fluid_cell(j, i))
    }

    fn upwind(face_velocity: f64, upstream: f64, downstream: f64) -> f64 {
        if face_velocity >= 0.0 {
            upstream
        } else {
            downstream
        }
    }

    /// The predictor's right-hand side on the u face `(j, i)`, `i = 1..nx`,
    /// from `field.u_old`, `field.v_old` and `field.p`: `−conv + diff −
    /// ∇p/ρ + f/ρ`, expression for expression the fixed-grid PISO's. The
    /// predictor writes `u_old + dt · rhs` on the fluid faces; the overset
    /// momentum-residual diagnostic (P4 option B) evaluates the same
    /// operator on the prescribed faces, so "the solver's own stencil" is
    /// this function by construction.
    #[allow(clippy::too_many_lines)]
    pub(crate) fn u_rhs(&self, field: &FlowField, j: usize, i: usize, t_old: f64) -> f64 {
        let (nx, ny, dx, dy) = field.grid_info();
        let rho = self.config.density;
        let nu = self.config.viscosity / rho;
        let b = self.parameters.boundaries;
        let velocity = SideBoundary::Velocity;
        let uo = &field.u_old;
        let vo = &field.v_old;
        let u_p = uo[(j, i)];

        let ue_face = 0.5 * (uo[(j, i)] + uo[(j, i + 1)]);
        let uw_face = 0.5 * (uo[(j, i - 1)] + uo[(j, i)]);

        let south_is_wall = j == 0;
        let north_is_wall = j + 1 == ny;

        // Transverse face velocities from the stored v faces — on a
        // domain side these are the prescribed boundary normals
        // (zero on a wall, the outflow on an outlet). The fixed-grid
        // PISO zeroes them on its walls, which is the same number on
        // a wall and wrong on an outlet: the outgoing mass flux must
        // carry momentum out, or the last row accumulates it.
        let vn_face = 0.5 * (vo[(j + 1, i - 1)] + vo[(j + 1, i)]);
        let vs_face = 0.5 * (vo[(j, i - 1)] + vo[(j, i)]);

        // Upwind value across a domain side: the boundary function's
        // tangential value on a Velocity side, the interior value
        // otherwise (zero-gradient).
        let beyond_north = if b.top == velocity {
            self.boundary(i as f64 * dx, ny as f64 * dy, t_old).0
        } else {
            u_p
        };
        let beyond_south = if b.bottom == velocity {
            self.boundary(i as f64 * dx, 0.0, t_old).0
        } else {
            u_p
        };

        let conv_x = (ue_face * Self::upwind(ue_face, uo[(j, i)], uo[(j, i + 1)])
            - uw_face * Self::upwind(uw_face, uo[(j, i - 1)], uo[(j, i)]))
            / dx;
        let conv_y = (vn_face
            * if north_is_wall {
                Self::upwind(vn_face, u_p, beyond_north)
            } else {
                Self::upwind(vn_face, uo[(j, i)], uo[(j + 1, i)])
            }
            - vs_face
                * if south_is_wall {
                    Self::upwind(vs_face, beyond_south, u_p)
                } else {
                    Self::upwind(vs_face, uo[(j - 1, i)], uo[(j, i)])
                })
            / dy;

        // Limited (TVD) corrections to the four convective face
        // values; exactly zero-cost on the default upwind scheme.
        let scheme = self.parameters.convection_scheme;
        let mut conv_x = conv_x;
        let mut conv_y = conv_y;
        if scheme != ConvectionScheme::Upwind {
            let delta_e = if ue_face >= 0.0 {
                scheme.face_correction(Some(uo[(j, i - 1)]), uo[(j, i)], uo[(j, i + 1)])
            } else {
                let far = (i + 2 <= nx).then(|| uo[(j, i + 2)]);
                scheme.face_correction(far, uo[(j, i + 1)], uo[(j, i)])
            };
            let delta_w = if uw_face >= 0.0 {
                let far = (i >= 2).then(|| uo[(j, i - 2)]);
                scheme.face_correction(far, uo[(j, i - 1)], uo[(j, i)])
            } else {
                scheme.face_correction(Some(uo[(j, i + 1)]), uo[(j, i)], uo[(j, i - 1)])
            };
            let delta_n = if north_is_wall {
                0.0
            } else if vn_face >= 0.0 {
                let far = (j >= 1).then(|| uo[(j - 1, i)]);
                scheme.face_correction(far, uo[(j, i)], uo[(j + 1, i)])
            } else {
                let far = (j + 2 < ny).then(|| uo[(j + 2, i)]);
                scheme.face_correction(far, uo[(j + 1, i)], uo[(j, i)])
            };
            let delta_s = if south_is_wall {
                0.0
            } else if vs_face >= 0.0 {
                let far = (j >= 2).then(|| uo[(j - 2, i)]);
                scheme.face_correction(far, uo[(j - 1, i)], uo[(j, i)])
            } else {
                let far = (j + 1 < ny).then(|| uo[(j + 1, i)]);
                scheme.face_correction(far, uo[(j, i)], uo[(j - 1, i)])
            };
            conv_x += (ue_face * delta_e - uw_face * delta_w) / dx;
            conv_y += (vn_face * delta_n - vs_face * delta_s) / dy;
        }

        let diff_x = nu * (uo[(j, i + 1)] - 2.0 * u_p + uo[(j, i - 1)]) / (dx * dx);

        // Wall-adjacent diffusive fluxes act over half a cell on a
        // Velocity side; a slip wall or outlet carries no shear.
        let flux_north = if north_is_wall {
            if b.top == velocity {
                let u_wall = self.boundary(i as f64 * dx, ny as f64 * dy, t_old).0;
                nu * (u_wall - u_p) / (0.5 * dy)
            } else {
                0.0
            }
        } else {
            nu * (uo[(j + 1, i)] - u_p) / dy
        };
        let flux_south = if south_is_wall {
            if b.bottom == velocity {
                let u_wall = self.boundary(i as f64 * dx, 0.0, t_old).0;
                nu * (u_p - u_wall) / (0.5 * dy)
            } else {
                0.0
            }
        } else {
            nu * (u_p - uo[(j - 1, i)]) / dy
        };
        let diff_y = (flux_north - flux_south) / dy;

        let pressure_gradient = -(field.p[(j, i)] - field.p[(j, i - 1)]) / (rho * dx);

        let body_force = self.momentum_source.as_ref().map_or(0.0, |f| {
            f(i as f64 * dx, (j as f64 + 0.5) * dy, t_old).0 / rho
        });

        -conv_x - conv_y + diff_x + diff_y + pressure_gradient + body_force
    }

    /// See [`Self::u_rhs`]: the v face `(j, i)`, `j = 1..ny`.
    #[allow(clippy::too_many_lines)]
    pub(crate) fn v_rhs(&self, field: &FlowField, j: usize, i: usize, t_old: f64) -> f64 {
        let (nx, ny, dx, dy) = field.grid_info();
        let rho = self.config.density;
        let nu = self.config.viscosity / rho;
        let b = self.parameters.boundaries;
        let velocity = SideBoundary::Velocity;
        let uo = &field.u_old;
        let vo = &field.v_old;
        let v_p = vo[(j, i)];

        let vn_face = 0.5 * (vo[(j, i)] + vo[(j + 1, i)]);
        let vs_face = 0.5 * (vo[(j - 1, i)] + vo[(j, i)]);

        let west_is_wall = i == 0;
        let east_is_wall = i + 1 == nx;

        let ue_face = 0.5 * (uo[(j - 1, i + 1)] + uo[(j, i + 1)]);
        let uw_face = 0.5 * (uo[(j - 1, i)] + uo[(j, i)]);
        let beyond_east = if b.right == velocity {
            self.boundary(nx as f64 * dx, j as f64 * dy, t_old).1
        } else {
            v_p
        };
        let beyond_west = if b.left == velocity {
            self.boundary(0.0, j as f64 * dy, t_old).1
        } else {
            v_p
        };

        let conv_y = (vn_face * Self::upwind(vn_face, vo[(j, i)], vo[(j + 1, i)])
            - vs_face * Self::upwind(vs_face, vo[(j - 1, i)], vo[(j, i)]))
            / dy;
        let conv_x = (ue_face
            * if east_is_wall {
                Self::upwind(ue_face, v_p, beyond_east)
            } else {
                Self::upwind(ue_face, vo[(j, i)], vo[(j, i + 1)])
            }
            - uw_face
                * if west_is_wall {
                    Self::upwind(uw_face, beyond_west, v_p)
                } else {
                    Self::upwind(uw_face, vo[(j, i - 1)], vo[(j, i)])
                })
            / dx;

        let scheme = self.parameters.convection_scheme;
        let mut conv_x = conv_x;
        let mut conv_y = conv_y;
        if scheme != ConvectionScheme::Upwind {
            let delta_n = if vn_face >= 0.0 {
                scheme.face_correction(Some(vo[(j - 1, i)]), vo[(j, i)], vo[(j + 1, i)])
            } else {
                let far = (j + 2 <= ny).then(|| vo[(j + 2, i)]);
                scheme.face_correction(far, vo[(j + 1, i)], vo[(j, i)])
            };
            let delta_s = if vs_face >= 0.0 {
                let far = (j >= 2).then(|| vo[(j - 2, i)]);
                scheme.face_correction(far, vo[(j - 1, i)], vo[(j, i)])
            } else {
                scheme.face_correction(Some(vo[(j + 1, i)]), vo[(j, i)], vo[(j - 1, i)])
            };
            let delta_e = if east_is_wall {
                0.0
            } else if ue_face >= 0.0 {
                let far = (i >= 1).then(|| vo[(j, i - 1)]);
                scheme.face_correction(far, vo[(j, i)], vo[(j, i + 1)])
            } else {
                let far = (i + 2 < nx).then(|| vo[(j, i + 2)]);
                scheme.face_correction(far, vo[(j, i + 1)], vo[(j, i)])
            };
            let delta_w = if west_is_wall {
                0.0
            } else if uw_face >= 0.0 {
                let far = (i >= 2).then(|| vo[(j, i - 2)]);
                scheme.face_correction(far, vo[(j, i - 1)], vo[(j, i)])
            } else {
                let far = (i + 1 < nx).then(|| vo[(j, i + 1)]);
                scheme.face_correction(far, vo[(j, i)], vo[(j, i - 1)])
            };
            conv_y += (vn_face * delta_n - vs_face * delta_s) / dy;
            conv_x += (ue_face * delta_e - uw_face * delta_w) / dx;
        }

        let diff_y = nu * (vo[(j + 1, i)] - 2.0 * v_p + vo[(j - 1, i)]) / (dy * dy);

        let flux_east = if east_is_wall {
            if b.right == velocity {
                let v_wall = self.boundary(nx as f64 * dx, j as f64 * dy, t_old).1;
                nu * (v_wall - v_p) / (0.5 * dx)
            } else {
                0.0
            }
        } else {
            nu * (vo[(j, i + 1)] - v_p) / dx
        };
        let flux_west = if west_is_wall {
            if b.left == velocity {
                let v_wall = self.boundary(0.0, j as f64 * dy, t_old).1;
                nu * (v_p - v_wall) / (0.5 * dx)
            } else {
                0.0
            }
        } else {
            nu * (v_p - vo[(j, i - 1)]) / dx
        };
        let diff_x = (flux_east - flux_west) / dx;

        let pressure_gradient = -(field.p[(j, i)] - field.p[(j - 1, i)]) / (rho * dy);

        let body_force = self.momentum_source.as_ref().map_or(0.0, |f| {
            f((i as f64 + 0.5) * dx, j as f64 * dy, t_old).1 / rho
        });

        -conv_x - conv_y + diff_x + diff_y + pressure_gradient + body_force
    }

    /// Explicit momentum predictor on the fluid faces, expression for
    /// expression the fixed-grid PISO's (so the no-body case is identical
    /// to the bit), plus the slip-wall / outlet arms of the ALE solver on
    /// the domain sides. Non-fluid faces keep their prescribed values.
    #[allow(clippy::too_many_lines)]
    fn momentum_predictor(&self, field: &mut FlowField, dt: f64, t_old: f64) -> CfdResult<()> {
        let (nx, ny, _, _) = field.grid_info();
        let b = self.parameters.boundaries;

        for j in 0..ny {
            for i in 1..nx {
                if !self.u_is_fluid(j, i) {
                    continue;
                }
                let rhs = self.u_rhs(field, j, i, t_old);
                field.u[(j, i)] = field.u_old[(j, i)] + dt * rhs;
            }
        }

        for j in 1..ny {
            for i in 0..nx {
                if !self.v_is_fluid(j, i) {
                    continue;
                }
                let rhs = self.v_rhs(field, j, i, t_old);
                field.v[(j, i)] = field.v_old[(j, i)] + dt * rhs;
            }
        }

        // Outlet faces: zero-gradient predictor value, corrected by the
        // projection.
        if b.left == SideBoundary::PressureOutlet {
            for j in 0..ny {
                field.u[(j, 0)] = field.u[(j, 1)];
            }
        }
        if b.right == SideBoundary::PressureOutlet {
            for j in 0..ny {
                field.u[(j, nx)] = field.u[(j, nx - 1)];
            }
        }
        if b.bottom == SideBoundary::PressureOutlet {
            for i in 0..nx {
                field.v[(0, i)] = field.v[(1, i)];
            }
        }
        if b.top == SideBoundary::PressureOutlet {
            for i in 0..nx {
                field.v[(ny, i)] = field.v[(ny - 1, i)];
            }
        }

        field.copy_to_starred();
        Ok(())
    }

    /// Advance one step of `dt`: [`Self::begin_step`], the correctors,
    /// [`Self::end_step`].
    pub async fn advance(&mut self, field: &mut FlowField, dt: f64) -> CfdResult<EmbeddedResult> {
        let start = self.begin_step(field, dt)?;

        let mut residual_history = Vec::new();
        let mut total_corrector_steps = 0;
        let mut final_residual = f64::INFINITY;
        for corrector in 0..self.parameters.corrector_steps.max(1) {
            let mass_residual = self.project(field, dt, corrector == 0)?;
            residual_history.push(mass_residual);
            final_residual = mass_residual;
            total_corrector_steps += 1;
            if mass_residual < self.parameters.tolerance {
                break;
            }
            field.copy_to_starred();
        }

        Ok(self.end_step(
            field,
            &start,
            residual_history,
            total_corrector_steps,
            final_residual,
        ))
    }

    /// Everything before the correctors: validation, lazy initialisation,
    /// the history shift, the explicit predictor, the new interval's
    /// boundary data, the moving-body mask rebuild (fresh-cell refill,
    /// ghost re-imposition), and `u* = u`. The overset coupling runs this
    /// on the background, then drives the correctors itself.
    pub(crate) fn begin_step(&mut self, field: &mut FlowField, dt: f64) -> CfdResult<StepStart> {
        if dt <= 0.0 || !dt.is_finite() {
            return Err(CfdError::invalid_parameter(format!(
                "time step must be positive and finite, got {dt}"
            )));
        }
        if !self.initialized {
            self.initialize(field)?;
        }
        let start_time = std::time::Instant::now();
        let t_old = self.time;
        let t_new = t_old + dt;

        // The start-of-step state is whatever the previous step left on the
        // boundary and ghost faces — no re-stamping (see module docs).
        field.update_old_values();
        self.momentum_predictor(field, dt, t_old)?;
        // Boundary data for the new interval; the predictor's `u` holds the
        // old boundary values until now.
        self.apply_boundary_normals(field, t_new);

        // A moving body: rebuild the mask at the end-of-step geometry,
        // refill the pressure of cells that just became fluid (their stored
        // p is stale by their time inside the body — the next predictor
        // would read its gradient), and impose the new mask's ghost values
        // from the previous corrected field.
        let mut fresh_cells = 0usize;
        let trace_sp = std::env::var("RTX_EMBEDDED_TRACE_SP").is_ok();
        self.fresh_trace.clear();
        if self.moving {
            if let Some(body) = &self.body {
                let (nx, ny, dx, dy) = field.grid_info();
                let new_mask = EmbeddedMask::build_with_reference(
                    body,
                    nx,
                    ny,
                    dx,
                    dy,
                    t_new,
                    self.mask.as_ref(),
                    self.mask_hysteresis * dx.min(dy),
                )?;
                if let Some(old_mask) = &self.mask {
                    for j in 0..ny {
                        for i in 0..nx {
                            if new_mask.is_fluid_cell(j, i) && !old_mask.is_fluid_cell(j, i) {
                                fresh_cells += 1;
                                if trace_sp {
                                    self.fresh_trace.push((j, i));
                                }
                                let mut sum = 0.0;
                                let mut count = 0usize;
                                let mut visit = |jj: usize, ii: usize| {
                                    if new_mask.is_fluid_cell(jj, ii)
                                        && old_mask.is_fluid_cell(jj, ii)
                                    {
                                        sum += field.p[(jj, ii)];
                                        count += 1;
                                    }
                                };
                                if i + 1 < nx {
                                    visit(j, i + 1);
                                }
                                if i > 0 {
                                    visit(j, i - 1);
                                }
                                if j + 1 < ny {
                                    visit(j + 1, i);
                                }
                                if j > 0 {
                                    visit(j - 1, i);
                                }
                                if count > 0 {
                                    field.p[(j, i)] = sum / count as f64;
                                }
                            }
                        }
                    }
                }
                if self.swept_volume != 0.0 {
                    let h = dx.min(dy);
                    let mut alpha = vec![1.0f64; nx * ny];
                    for j in 0..ny {
                        for i in 0..nx {
                            let phi = body.phi((i as f64 + 0.5) * dx, (j as f64 + 0.5) * dy, t_new);
                            alpha[j * nx + i] = (0.5 + phi / h).clamp(0.0, 1.0);
                        }
                    }
                    if self.alpha_old.is_none() {
                        self.alpha_old = Some(alpha.clone());
                    }
                    self.alpha_new = Some(alpha);
                }
                let u_history = field.u_old.clone();
                let v_history = field.v_old.clone();
                new_mask.impose_from(
                    body,
                    &u_history,
                    &v_history,
                    &mut field.u,
                    &mut field.v,
                    t_new,
                );
                if self.field_extension {
                    if let Some(old_mask) = &self.mask {
                        old_mask.extend_fresh_faces(
                            &new_mask,
                            body,
                            t_new,
                            &u_history,
                            &v_history,
                            &mut field.u,
                            &mut field.v,
                            &mut field.u_old,
                            &mut field.v_old,
                        );
                    }
                }
                self.mask = Some(new_mask);
            }
        }
        field.copy_to_starred();
        Ok(StepStart {
            start_time,
            t_new,
            fresh_cells,
        })
    }

    /// Everything after the correctors: ghost re-imposition from the
    /// corrected field, the clock, the swept-volume bookkeeping, the
    /// result.
    pub(crate) fn end_step(
        &mut self,
        field: &mut FlowField,
        start: &StepStart,
        residual_history: Vec<f64>,
        total_corrector_steps: usize,
        final_residual: f64,
    ) -> EmbeddedResult {
        let t_new = start.t_new;
        // Ghost faces follow the corrected fluid field; they are the
        // stencil and flux data of the next step.
        let ghost_correction = match (&self.body, &self.mask) {
            (Some(body), Some(mask)) => mask.impose(body, &mut field.u, &mut field.v, t_new),
            _ => 0.0,
        };

        self.time = t_new;
        if let Some(a) = self.alpha_new.take() {
            self.alpha_old = Some(a);
        }
        EmbeddedResult {
            solver_result: SolverResult {
                converged: final_residual < self.parameters.tolerance,
                iterations: total_corrector_steps,
                final_residual,
                residual_history,
                solve_time: start.start_time.elapsed(),
            },
            corrector_steps_performed: total_corrector_steps,
            ghost_correction,
            fresh_cells: start.fresh_cells,
        }
    }
}

/// What [`EmbeddedPisoSolver::begin_step`] hands to
/// [`EmbeddedPisoSolver::end_step`].
#[derive(Debug, Clone, Copy)]
pub(crate) struct StepStart {
    start_time: std::time::Instant,
    /// End-of-step time.
    pub(crate) t_new: f64,
    /// Pressure cells that flipped solid → fluid in this step's rebuild.
    pub(crate) fresh_cells: usize,
}