//! ALE (arbitrary Lagrangian–Eulerian) incompressible solver on a moving
//! tensor-product staggered grid.
//!
//! This is the PISO scheme — explicit conservative momentum predictor, then
//! pressure-correction projections — generalised to a mesh whose x-lines and
//! y-lines move arbitrarily in time while the domain boundary stays fixed.
//! Cells remain axis-aligned rectangles (tensor-product motion), so the
//! staggered MAC layout survives: `u[(j, i)]` on the x-line `x[i]` at the
//! cell-centre height, `v[(j, i)]` on the y-line `y[j]`, `p[(j, i)]` at cell
//! centres. Spacing is non-uniform in both directions and changes every step.
//!
//! # The discrete geometric conservation law
//!
//! The momentum update is the conservative ALE form
//!
//! ```text
//!   (V^{n+1} u^{n+1} - V^n u^n)/dt + sum_f q_f u_f = RHS,
//!   q_f = u_f . n A_f  -  sweptVol_f / dt
//! ```
//!
//! and its face areas are the **time-averaged** (trapezoidal) ones,
//! `A_f = (A_f^n + A_f^{n+1}) / 2`, in both the fluid flux and the swept
//! volume. For tensor-product motion that choice satisfies the geometric
//! conservation law *exactly*:
//!
//! ```text
//!   dx1 dy1 - dx0 dy0 = (dx1 - dx0)(dy0 + dy1)/2 + (dy1 - dy0)(dx0 + dx1)/2
//! ```
//!
//! is an algebraic identity, so the sum of the signed swept volumes equals
//! the cell's volume increment to rounding error and a uniform flow is an
//! exact fixed point of the discrete update on any admissible mesh motion —
//! which is what `tests/ale_dgcl.rs` asserts at 1e-12. The tempting
//! alternative — end-of-step areas, [`SweptFaceRule::EndOfStep`] — is kept
//! only as the test's negative control: it leaves a per-step relative error
//! of exactly `dw dh / V` per cell (the cross term the identity absorbs),
//! invisible to every consistency check and fatal to long FSI runs.
//!
//! # Incompressibility on a moving mesh
//!
//! Mass conservation for a moving cell is `dV/dt + sum (u - w).n A = 0`;
//! subtracting the GCL (`dV/dt = sum w.n A`) leaves `sum u.n A = 0` — plain
//! divergence-freedom in the *current* geometry, with no mesh-velocity term.
//! The projection therefore works exactly as on a fixed grid, assembled on
//! the end-of-step geometry: prescribed normal velocities on the whole
//! boundary make it pure Neumann, one cell anchors the level, and SOR at the
//! optimal Poisson factor with a true-residual stop does the inner solve
//! (both lessons inherited from the fixed-grid PISO: see its module docs).

use super::SolverResult;
use crate::{CfdConfig, CfdError, CfdResult};
use nalgebra::DMatrix;

/// Which face areas enter the fluid fluxes and swept volumes.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum SweptFaceRule {
    /// Time-averaged areas: satisfies the discrete GCL exactly for
    /// tensor-product motion. The only correct choice; the default.
    Trapezoidal,
    /// End-of-step areas: first-order consistent and GCL-violating. Exists
    /// solely as the negative control for `tests/ale_dgcl.rs`.
    EndOfStep,
}

/// What one side of the domain boundary is.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum SideBoundary {
    /// Prescribed velocity (the default): the boundary function supplies
    /// the normal component (data for the projection) and the tangential
    /// value for no-slip half-cell wall diffusion. If the boundary line
    /// moves, the prescribed normal velocity must equal the line's motion
    /// `(new - old)/dt` — a material wall — or mass bookkeeping will not
    /// close.
    #[default]
    Velocity,
    /// Impenetrable frictionless wall: normal velocity from the boundary
    /// function (usually zero), zero tangential shear.
    SlipWall,
    /// Open boundary at gauge pressure zero: the normal velocity is an
    /// unknown (zero-gradient predictor, corrected by the projection, whose
    /// `p'` takes a Dirichlet zero on the face — which also makes the
    /// Poisson system non-singular, so no cell is anchored). For a non-zero
    /// outlet pressure, shift the gauge.
    PressureOutlet,
}

/// Boundary type per domain side.
#[derive(Debug, Clone, Copy, Default)]
pub struct AleBoundaries {
    /// x = x\[0\].
    pub left: SideBoundary,
    /// x = x\[nx\].
    pub right: SideBoundary,
    /// y = y\[0\].
    pub bottom: SideBoundary,
    /// y = y\[ny\].
    pub top: SideBoundary,
}

impl AleBoundaries {
    fn any_outlet(self) -> bool {
        [self.left, self.right, self.bottom, self.top].contains(&SideBoundary::PressureOutlet)
    }
}

/// Parameters for the ALE solver.
#[derive(Debug, Clone)]
pub struct AleParameters {
    /// Projection passes per step (2 suffices with an explicit predictor;
    /// more only mop up inner-solver truncation).
    pub corrector_steps: usize,
    /// Convergence tolerance on the normalised mass imbalance after
    /// correction.
    pub tolerance: f64,
    /// Face-area rule; see [`SweptFaceRule`].
    pub swept_face_rule: SweptFaceRule,
    /// Boundary type per domain side; all [`SideBoundary::Velocity`] by
    /// default.
    pub boundaries: AleBoundaries,
}

impl Default for AleParameters {
    fn default() -> Self {
        Self {
            corrector_steps: 2,
            tolerance: 1e-6,
            swept_face_rule: SweptFaceRule::Trapezoidal,
            boundaries: AleBoundaries::default(),
        }
    }
}

/// Result of one ALE time step.
#[derive(Debug, Clone)]
pub struct AleResult {
    /// Base solver result information.
    pub solver_result: SolverResult,
    /// Number of projection passes performed.
    pub corrector_steps_performed: usize,
}

/// Staggered flow state on a moving tensor-product grid. The node lines `x`
/// (length `nx + 1`) and `y` (length `ny + 1`) are part of the state and are
/// advanced by [`AlePisoSolver::advance`]; `x_old`/`y_old` hold the previous
/// step's lines so the solver can form swept volumes.
#[derive(Debug, Clone)]
pub struct AleField {
    /// Cells in x.
    pub nx: usize,
    /// Cells in y.
    pub ny: usize,
    /// Current node lines.
    pub x: Vec<f64>,
    /// Current node lines.
    pub y: Vec<f64>,
    /// Node lines at the start of the current step.
    pub x_old: Vec<f64>,
    /// Node lines at the start of the current step.
    pub y_old: Vec<f64>,
    /// u on x-lines: `(ny, nx + 1)`.
    pub u: DMatrix<f64>,
    /// v on y-lines: `(ny + 1, nx)`.
    pub v: DMatrix<f64>,
    /// Pressure at cell centres: `(ny, nx)`.
    pub p: DMatrix<f64>,
    /// Start-of-step velocities (what the explicit predictor differentiates).
    pub u_old: DMatrix<f64>,
    /// Start-of-step velocities.
    pub v_old: DMatrix<f64>,
    /// Predicted (pre-projection) velocities.
    pub u_star: DMatrix<f64>,
    /// Predicted (pre-projection) velocities.
    pub v_star: DMatrix<f64>,
    /// Pressure correction.
    pub p_prime: DMatrix<f64>,
    /// Projection source (per-cell mass imbalance flux).
    pub sp: DMatrix<f64>,
}

fn validate_lines(lines: &[f64], name: &str) -> CfdResult<()> {
    if lines.len() < 4 {
        return Err(CfdError::invalid_parameter(format!(
            "{name}: need at least 3 cells (4 node lines), got {}",
            lines.len().saturating_sub(1)
        )));
    }
    for pair in lines.windows(2) {
        if pair[1] <= pair[0] {
            return Err(CfdError::invalid_parameter(format!(
                "{name}: node lines must be strictly increasing \
                 ({} then {} — a cell has non-positive volume)",
                pair[0], pair[1]
            )));
        }
    }
    Ok(())
}

impl AleField {
    /// Create a field on the given node lines, all values zero.
    pub fn new(x: Vec<f64>, y: Vec<f64>) -> CfdResult<Self> {
        validate_lines(&x, "x")?;
        validate_lines(&y, "y")?;
        let nx = x.len() - 1;
        let ny = y.len() - 1;
        Ok(Self {
            nx,
            ny,
            x_old: x.clone(),
            y_old: y.clone(),
            x,
            y,
            u: DMatrix::zeros(ny, nx + 1),
            v: DMatrix::zeros(ny + 1, nx),
            p: DMatrix::zeros(ny, nx),
            u_old: DMatrix::zeros(ny, nx + 1),
            v_old: DMatrix::zeros(ny + 1, nx),
            u_star: DMatrix::zeros(ny, nx + 1),
            v_star: DMatrix::zeros(ny + 1, nx),
            p_prime: DMatrix::zeros(ny, nx),
            sp: DMatrix::zeros(ny, nx),
        })
    }

    /// Uniformly spaced field on `[0, lx] x [0, ly]`.
    pub fn uniform(nx: usize, ny: usize, lx: f64, ly: f64) -> CfdResult<Self> {
        if lx <= 0.0 || ly <= 0.0 {
            return Err(CfdError::invalid_parameter(
                "domain lengths must be positive",
            ));
        }
        let x = (0..=nx).map(|i| lx * i as f64 / nx as f64).collect();
        let y = (0..=ny).map(|j| ly * j as f64 / ny as f64).collect();
        Self::new(x, y)
    }

    fn copy_to_starred(&mut self) {
        self.u_star.copy_from(&self.u);
        self.v_star.copy_from(&self.v);
    }
}

/// Cell-centre coordinates for a set of node lines.
fn centres(lines: &[f64]) -> Vec<f64> {
    lines.windows(2).map(|w| 0.5 * (w[0] + w[1])).collect()
}

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

/// The ALE PISO solver. See the module docs for the discretisation.
pub struct AlePisoSolver {
    config: CfdConfig,
    parameters: AleParameters,
    /// Prescribed velocity `(x, y, t) -> (u, v)` on the domain boundary: it
    /// supplies the normal components on boundary faces (which the
    /// projection treats as data, not unknowns) and the tangential values
    /// the near-wall half-cell diffusion needs. `None` means a closed
    /// no-slip box.
    boundary_velocity: Option<VelocityFn>,
    /// Optional volumetric momentum source `(x, y, t) -> (f_x, f_y)` per
    /// unit volume — the hook a manufactured solution enters through.
    momentum_source: Option<SourceFn>,
    /// Accumulated physical time; advances by `dt` each step.
    time: f64,
}

impl AlePisoSolver {
    /// Create a new solver.
    pub fn new(config: CfdConfig, parameters: AleParameters) -> CfdResult<Self> {
        config.validate()?;
        Ok(Self {
            config,
            parameters,
            boundary_velocity: None,
            momentum_source: None,
            time: 0.0,
        })
    }

    /// Set the boundary velocity. See [`Self::boundary_velocity`].
    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));
    }

    /// Set a volumetric momentum source. See [`Self::momentum_source`].
    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));
    }

    /// Physical time the state has been advanced to.
    pub fn time(&self) -> f64 {
        self.time
    }

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

    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))
    }

    /// Write the prescribed normal velocities onto the boundary faces of the
    /// given geometry at time `t`. Outlet faces are unknowns and are left
    /// alone.
    fn apply_boundary_normals(&self, field: &mut AleField, t: f64, x: &[f64], y: &[f64]) {
        let (nx, ny) = (field.nx, field.ny);
        let b = self.parameters.boundaries;
        let outlet = SideBoundary::PressureOutlet;
        let yc = centres(y);
        let xc = centres(x);
        for j in 0..ny {
            if b.left != outlet {
                field.u[(j, 0)] = self.boundary(x[0], yc[j], t).0;
            }
            if b.right != outlet {
                field.u[(j, nx)] = self.boundary(x[nx], yc[j], t).0;
            }
        }
        for i in 0..nx {
            if b.bottom != outlet {
                field.v[(0, i)] = self.boundary(xc[i], y[0], t).1;
            }
            if b.top != outlet {
                field.v[(ny, i)] = self.boundary(xc[i], y[ny], t).1;
            }
        }
    }

    /// Explicit conservative ALE momentum predictor. Every flux is built
    /// from `u_old`/`v_old` and the old/new node lines, so the step is
    /// genuinely explicit and independent of sweep order.
    #[allow(clippy::too_many_lines)]
    fn momentum_predictor(&self, field: &mut AleField, dt: f64, t_old: f64) -> CfdResult<()> {
        let (nx, ny) = (field.nx, field.ny);
        let rho = self.config.density;
        let nu = self.config.viscosity / rho;

        let xo = field.x_old.clone();
        let yo = field.y_old.clone();
        let xn = field.x.clone();
        let yn = field.y.clone();
        let xco = centres(&xo);
        let yco = centres(&yo);
        let xcn = centres(&xn);
        let ycn = centres(&yn);

        // Face area per the configured rule: the trapezoidal average is the
        // GCL-exact choice, end-of-step is the negative control.
        let area = |old: f64, new: f64| -> f64 {
            match self.parameters.swept_face_rule {
                SweptFaceRule::Trapezoidal => 0.5 * (old + new),
                SweptFaceRule::EndOfStep => new,
            }
        };

        // u control volumes: [xc(i-1), xc(i)] x [y_j, y_{j+1}], i = 1..nx.
        for j in 0..ny {
            for i in 1..nx {
                let uo = &field.u_old;
                let vo = &field.v_old;

                let w_o = xco[i] - xco[i - 1];
                let w_n = xcn[i] - xcn[i - 1];
                let h_o = yo[j + 1] - yo[j];
                let h_n = yn[j + 1] - yn[j];
                let v_old_cell = w_o * h_o;
                let v_new_cell = w_n * h_n;

                // Vertical faces at the cell centres east and west.
                let a_ew = area(h_o, h_n);
                let swept_e = (xcn[i] - xco[i]) * a_ew;
                let swept_w = (xcn[i - 1] - xco[i - 1]) * a_ew;

                // Horizontal faces: the CV width splits at the u-node into
                // the halves owned by the two neighbouring pressure cells,
                // which carry different v values.
                let l_half = area(xo[i] - xco[i - 1], xn[i] - xcn[i - 1]);
                let r_half = area(xco[i] - xo[i], xcn[i] - xn[i]);
                let swept_n = (yn[j + 1] - yo[j + 1]) * (l_half + r_half);
                let swept_s = (yn[j] - yo[j]) * (l_half + r_half);

                // Outward relative fluxes q = u.n A - swept/dt.
                let u_e = 0.5 * (uo[(j, i)] + uo[(j, i + 1)]);
                let u_w = 0.5 * (uo[(j, i - 1)] + uo[(j, i)]);
                let q_e = u_e * a_ew - swept_e / dt;
                let q_w = -(u_w * a_ew - swept_w / dt);
                let vn_flux = vo[(j + 1, i - 1)] * l_half + vo[(j + 1, i)] * r_half;
                let vs_flux = vo[(j, i - 1)] * l_half + vo[(j, i)] * r_half;
                let q_n = vn_flux - swept_n / dt;
                let q_s = -(vs_flux - swept_s / dt);

                // Upwinded momentum on each face; inflow across a domain
                // boundary carries the prescribed boundary value.
                let phi_e = if q_e >= 0.0 {
                    uo[(j, i)]
                } else {
                    uo[(j, i + 1)]
                };
                let phi_w = if q_w >= 0.0 {
                    uo[(j, i)]
                } else {
                    uo[(j, i - 1)]
                };
                let phi_n = if q_n >= 0.0 {
                    uo[(j, i)]
                } else if j + 1 < ny {
                    uo[(j + 1, i)]
                } else if self.parameters.boundaries.top == SideBoundary::Velocity {
                    self.boundary(xo[i], yo[ny], t_old).0
                } else {
                    // Slip wall or outlet: no prescribed tangential value;
                    // carry the interior one.
                    uo[(j, i)]
                };
                let phi_s = if q_s >= 0.0 {
                    uo[(j, i)]
                } else if j > 0 {
                    uo[(j - 1, i)]
                } else if self.parameters.boundaries.bottom == SideBoundary::Velocity {
                    self.boundary(xo[i], yo[0], t_old).0
                } else {
                    uo[(j, i)]
                };
                let conv = q_e * phi_e + q_w * phi_w + q_n * phi_n + q_s * phi_s;

                // Diffusive fluxes on the old geometry; wall-adjacent fluxes
                // act over the actual half-cell distance to the wall.
                let d_e = nu * (uo[(j, i + 1)] - uo[(j, i)]) / (xo[i + 1] - xo[i]) * h_o;
                let d_w = nu * (uo[(j, i - 1)] - uo[(j, i)]) / (xo[i] - xo[i - 1]) * h_o;
                let d_n = if j + 1 < ny {
                    nu * (uo[(j + 1, i)] - uo[(j, i)]) / (yco[j + 1] - yco[j]) * w_o
                } else if self.parameters.boundaries.top == SideBoundary::Velocity {
                    let u_wall = self.boundary(xo[i], yo[ny], t_old).0;
                    nu * (u_wall - uo[(j, i)]) / (yo[ny] - yco[j]) * w_o
                } else {
                    // Slip wall or outlet: zero tangential shear.
                    0.0
                };
                let d_s = if j > 0 {
                    nu * (uo[(j - 1, i)] - uo[(j, i)]) / (yco[j] - yco[j - 1]) * w_o
                } else if self.parameters.boundaries.bottom == SideBoundary::Velocity {
                    let u_wall = self.boundary(xo[i], yo[0], t_old).0;
                    nu * (u_wall - uo[(j, i)]) / (yco[j] - yo[0]) * w_o
                } else {
                    0.0
                };
                let diff = d_e + d_w + d_n + d_s;

                // Net pressure force on the CV; the staggered layout puts
                // the cell-centre pressures exactly on its vertical faces.
                let pres = -(field.p[(j, i)] - field.p[(j, i - 1)]) * h_o / rho;

                let src = self
                    .momentum_source
                    .as_ref()
                    .map_or(0.0, |f| f(xo[i], yco[j], t_old).0 * v_old_cell / rho);

                field.u[(j, i)] =
                    (v_old_cell * uo[(j, i)] + dt * (-conv + diff + pres + src)) / v_new_cell;
            }
        }

        // v control volumes: [x_i, x_{i+1}] x [yc(j-1), yc(j)], j = 1..ny.
        for j in 1..ny {
            for i in 0..nx {
                let uo = &field.u_old;
                let vo = &field.v_old;

                let w_o = xo[i + 1] - xo[i];
                let w_n = xn[i + 1] - xn[i];
                let h_o = yco[j] - yco[j - 1];
                let h_n = ycn[j] - ycn[j - 1];
                let v_old_cell = w_o * h_o;
                let v_new_cell = w_n * h_n;

                let a_ns = area(w_o, w_n);
                let swept_n = (ycn[j] - yco[j]) * a_ns;
                let swept_s = (ycn[j - 1] - yco[j - 1]) * a_ns;

                let b_half = area(yo[j] - yco[j - 1], yn[j] - ycn[j - 1]);
                let t_half = area(yco[j] - yo[j], ycn[j] - yn[j]);
                let swept_e = (xn[i + 1] - xo[i + 1]) * (b_half + t_half);
                let swept_w = (xn[i] - xo[i]) * (b_half + t_half);

                let v_n = 0.5 * (vo[(j, i)] + vo[(j + 1, i)]);
                let v_s = 0.5 * (vo[(j - 1, i)] + vo[(j, i)]);
                let q_n = v_n * a_ns - swept_n / dt;
                let q_s = -(v_s * a_ns - swept_s / dt);
                let ue_flux = uo[(j - 1, i + 1)] * b_half + uo[(j, i + 1)] * t_half;
                let uw_flux = uo[(j - 1, i)] * b_half + uo[(j, i)] * t_half;
                let q_e = ue_flux - swept_e / dt;
                let q_w = -(uw_flux - swept_w / dt);

                let phi_n = if q_n >= 0.0 {
                    vo[(j, i)]
                } else {
                    vo[(j + 1, i)]
                };
                let phi_s = if q_s >= 0.0 {
                    vo[(j, i)]
                } else {
                    vo[(j - 1, i)]
                };
                let phi_e = if q_e >= 0.0 {
                    vo[(j, i)]
                } else if i + 1 < nx {
                    vo[(j, i + 1)]
                } else if self.parameters.boundaries.right == SideBoundary::Velocity {
                    self.boundary(xo[nx], yo[j], t_old).1
                } else {
                    vo[(j, i)]
                };
                let phi_w = if q_w >= 0.0 {
                    vo[(j, i)]
                } else if i > 0 {
                    vo[(j, i - 1)]
                } else if self.parameters.boundaries.left == SideBoundary::Velocity {
                    self.boundary(xo[0], yo[j], t_old).1
                } else {
                    vo[(j, i)]
                };
                let conv = q_e * phi_e + q_w * phi_w + q_n * phi_n + q_s * phi_s;

                let d_n = nu * (vo[(j + 1, i)] - vo[(j, i)]) / (yo[j + 1] - yo[j]) * w_o;
                let d_s = nu * (vo[(j - 1, i)] - vo[(j, i)]) / (yo[j] - yo[j - 1]) * w_o;
                let d_e = if i + 1 < nx {
                    nu * (vo[(j, i + 1)] - vo[(j, i)]) / (xco[i + 1] - xco[i]) * h_o
                } else if self.parameters.boundaries.right == SideBoundary::Velocity {
                    let v_wall = self.boundary(xo[nx], yo[j], t_old).1;
                    nu * (v_wall - vo[(j, i)]) / (xo[nx] - xco[i]) * h_o
                } else {
                    0.0
                };
                let d_w = if i > 0 {
                    nu * (vo[(j, i - 1)] - vo[(j, i)]) / (xco[i] - xco[i - 1]) * h_o
                } else if self.parameters.boundaries.left == SideBoundary::Velocity {
                    let v_wall = self.boundary(xo[0], yo[j], t_old).1;
                    nu * (v_wall - vo[(j, i)]) / (xco[i] - xo[0]) * h_o
                } else {
                    0.0
                };
                let diff = d_e + d_w + d_n + d_s;

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

                let src = self
                    .momentum_source
                    .as_ref()
                    .map_or(0.0, |f| f(xco[i], yo[j], t_old).1 * v_old_cell / rho);

                field.v[(j, i)] =
                    (v_old_cell * vo[(j, i)] + dt * (-conv + diff + pres + src)) / v_new_cell;
            }
        }

        // Outlet faces are unknowns without a control volume of their own:
        // give them the zero-gradient (fully developed) predictor value and
        // let the projection correct them.
        let b = self.parameters.boundaries;
        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)];
            }
        }

        Ok(())
    }

    /// One projection on the end-of-step geometry: solve the
    /// pressure-correction Poisson equation and subtract
    /// `(dt/rho) grad(p')` from the predicted velocities. Structure and
    /// inner-solve safeguards are the fixed-grid PISO's (anchored Neumann,
    /// SOR at the optimal factor, true-residual stop) with the coefficients
    /// generalised to non-uniform spacing.
    fn project(&self, field: &mut AleField, dt: f64) -> CfdResult<f64> {
        let (nx, ny) = (field.nx, field.ny);
        let rho = self.config.density;
        let xn = field.x.clone();
        let yn = field.y.clone();
        let xcn = centres(&xn);
        let ycn = centres(&yn);

        field.p_prime.fill(0.0);

        let mut source_scale = 0.0;
        for j in 0..ny {
            let dy_j = yn[j + 1] - yn[j];
            for i in 0..nx {
                let dx_i = xn[i + 1] - xn[i];
                let divergence_flux = rho
                    * ((field.u_star[(j, i + 1)] - field.u_star[(j, i)]) * dy_j
                        + (field.v_star[(j + 1, i)] - field.v_star[(j, i)]) * dx_i);
                field.sp[(j, i)] = -divergence_flux;
                source_scale += divergence_flux.abs();
            }
        }

        let b = self.parameters.boundaries;
        let outlet = SideBoundary::PressureOutlet;
        let reference_flux = rho * self.config.reference_velocity * self.config.reference_length;
        let inner_stop =
            (1e-2 * source_scale).max(0.1 * self.parameters.tolerance * reference_flux) + 1e-14;
        let omega = 2.0 / (1.0 + (std::f64::consts::PI / nx.max(ny) as f64).sin());
        for _sweep in 0..2000 {
            let mut residual = 0.0;
            for j in 0..ny {
                let dy_j = yn[j + 1] - yn[j];
                for i in 0..nx {
                    // With velocity prescribed on the whole boundary the
                    // system is pure Neumann and one cell anchors the level;
                    // any outlet contributes a Dirichlet face instead, and
                    // the anchor must NOT also be imposed.
                    if !b.any_outlet() && i == 1 && j == 1 {
                        field.p_prime[(j, i)] = 0.0;
                        continue;
                    }
                    let dx_i = xn[i + 1] - xn[i];

                    // A coefficient is zero exactly when its face is a
                    // domain boundary with prescribed normal velocity; an
                    // outlet face instead carries `p' = 0` half a cell away,
                    // so its coefficient survives with no neighbour term.
                    let ae = if i + 1 < nx {
                        dt * dy_j / (xcn[i + 1] - xcn[i])
                    } else if b.right == outlet {
                        dt * dy_j / (xn[nx] - xcn[i])
                    } else {
                        0.0
                    };
                    let aw = if i > 0 {
                        dt * dy_j / (xcn[i] - xcn[i - 1])
                    } else if b.left == outlet {
                        dt * dy_j / (xcn[0] - xn[0])
                    } else {
                        0.0
                    };
                    let an = if j + 1 < ny {
                        dt * dx_i / (ycn[j + 1] - ycn[j])
                    } else if b.top == outlet {
                        dt * dx_i / (yn[ny] - ycn[j])
                    } else {
                        0.0
                    };
                    let as_ = if j > 0 {
                        dt * dx_i / (ycn[j] - ycn[j - 1])
                    } else if b.bottom == outlet {
                        dt * dx_i / (ycn[0] - yn[0])
                    } else {
                        0.0
                    };
                    let ap = ae + aw + an + as_;

                    let east = if i + 1 < nx {
                        ae * field.p_prime[(j, i + 1)]
                    } else {
                        0.0
                    };
                    let west = if i > 0 {
                        aw * field.p_prime[(j, i - 1)]
                    } else {
                        0.0
                    };
                    let north = if j + 1 < ny {
                        an * field.p_prime[(j + 1, i)]
                    } else {
                        0.0
                    };
                    let south = if j > 0 {
                        as_ * field.p_prime[(j - 1, i)]
                    } else {
                        0.0
                    };

                    let rhs = field.sp[(j, i)] + east + west + north + south;
                    let p_old = field.p_prime[(j, i)];
                    residual += (rhs - ap * p_old).abs();
                    field.p_prime[(j, i)] = (1.0 - omega) * p_old + omega * rhs / ap;
                }
            }
            if residual < inner_stop {
                break;
            }
        }

        // Correct exactly the faces the equation treated as correctable:
        // every interior face.
        for j in 0..ny {
            for i in 1..nx {
                let dp_dx =
                    (field.p_prime[(j, i)] - field.p_prime[(j, i - 1)]) / (xcn[i] - xcn[i - 1]);
                field.u[(j, i)] = field.u_star[(j, i)] - (dt / rho) * dp_dx;
            }
        }
        for j in 1..ny {
            for i in 0..nx {
                let dp_dy =
                    (field.p_prime[(j, i)] - field.p_prime[(j - 1, i)]) / (ycn[j] - ycn[j - 1]);
                field.v[(j, i)] = field.v_star[(j, i)] - (dt / rho) * dp_dy;
            }
        }
        // Outlet faces are correctable too, against the Dirichlet `p' = 0`
        // on the face itself.
        if b.right == outlet {
            for j in 0..ny {
                let dp_dx = (0.0 - field.p_prime[(j, nx - 1)]) / (xn[nx] - xcn[nx - 1]);
                field.u[(j, nx)] = field.u_star[(j, nx)] - (dt / rho) * dp_dx;
            }
        }
        if b.left == outlet {
            for j in 0..ny {
                let dp_dx = (field.p_prime[(j, 0)] - 0.0) / (xcn[0] - xn[0]);
                field.u[(j, 0)] = field.u_star[(j, 0)] - (dt / rho) * dp_dx;
            }
        }
        if b.top == outlet {
            for i in 0..nx {
                let dp_dy = (0.0 - field.p_prime[(ny - 1, i)]) / (yn[ny] - ycn[ny - 1]);
                field.v[(ny, i)] = field.v_star[(ny, i)] - (dt / rho) * dp_dy;
            }
        }
        if b.bottom == outlet {
            for i in 0..nx {
                let dp_dy = (field.p_prime[(0, i)] - 0.0) / (ycn[0] - yn[0]);
                field.v[(0, i)] = field.v_star[(0, i)] - (dt / rho) * dp_dy;
            }
        }
        for j in 0..ny {
            for i in 0..nx {
                field.p[(j, i)] += field.p_prime[(j, i)];
            }
        }

        let mut mass_imbalance = 0.0;
        for j in 0..ny {
            let dy_j = yn[j + 1] - yn[j];
            for i in 0..nx {
                let dx_i = xn[i + 1] - xn[i];
                let divergence_flux = rho
                    * ((field.u[(j, i + 1)] - field.u[(j, i)]) * dy_j
                        + (field.v[(j + 1, i)] - field.v[(j, i)]) * dx_i);
                mass_imbalance += divergence_flux.abs();
            }
        }
        Ok(if reference_flux > 0.0 {
            mass_imbalance / reference_flux
        } else {
            mass_imbalance
        })
    }

    /// Advance one time step of size `dt`, moving the mesh nodes to
    /// `new_x`/`new_y` (strictly increasing — the motion may not invert a
    /// cell). Boundary lines may move: a moving `Velocity` side is a
    /// material wall, so its prescribed normal velocity must equal the
    /// line's motion `(new - old)/dt` or discrete mass bookkeeping will
    /// not close.
    pub async fn advance(
        &mut self,
        field: &mut AleField,
        new_x: &[f64],
        new_y: &[f64],
        dt: f64,
    ) -> CfdResult<AleResult> {
        let start_time = std::time::Instant::now();
        if dt <= 0.0 {
            return Err(CfdError::invalid_parameter("dt must be positive"));
        }
        if new_x.len() != field.nx + 1 || new_y.len() != field.ny + 1 {
            return Err(CfdError::invalid_parameter(format!(
                "node-line counts must not change: expected {}+1 x-lines and {}+1 y-lines, \
                 got {} and {}",
                field.nx,
                field.ny,
                new_x.len(),
                new_y.len()
            )));
        }
        validate_lines(new_x, "new_x")?;
        validate_lines(new_y, "new_y")?;

        let t_old = self.time;
        let t_new = t_old + dt;

        // The start-of-step boundary faces are whatever the previous step's
        // end-of-step application (or the caller's initial condition) left
        // there — the fluid's actual state at t_old. Re-stamping them here
        // from the boundary function would silently substitute the *new*
        // interval's wall velocity for the old one whenever the function
        // carries per-step data (an FSI coupling does exactly that), and
        // the resulting inconsistent old state leaves an O(dt) pressure
        // artifact in the wall-adjacent cells. Found by the piston test:
        // p exact to 6e-11 everywhere except the wall cell at 4.7e-5.
        field.x_old.clone_from(&field.x);
        field.y_old.clone_from(&field.y);
        field.x.copy_from_slice(new_x);
        field.y.copy_from_slice(new_y);
        field.u_old.copy_from(&field.u);
        field.v_old.copy_from(&field.v);

        self.momentum_predictor(field, dt, t_old)?;

        // The projection enforces continuity at the end of the step, so the
        // boundary faces must already carry their end-of-step data.
        let (x1, y1) = (field.x.clone(), field.y.clone());
        self.apply_boundary_normals(field, t_new, &x1, &y1);
        field.copy_to_starred();

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

        self.time = t_new;
        Ok(AleResult {
            solver_result: SolverResult {
                converged: final_residual < self.parameters.tolerance,
                iterations: total_correctors,
                final_residual,
                residual_history,
                solve_time: start_time.elapsed(),
            },
            corrector_steps_performed: total_correctors,
        })
    }
}