rtx-cfd: ALE on a moving tensor-product grid, DGCL-exact by construction

The first brick of the Turek-Hron frontier: PISO (explicit conservative
predictor + SOR projection) generalised to a staggered grid whose x- and
y-lines move arbitrarily each step while the domain boundary stays fixed.

The discretisation choice that carries everything: time-averaged face
areas (A^n + A^{n+1})/2 in both the fluid fluxes and the face-swept
volumes. For tensor-product motion the discrete geometric conservation
law then holds as an algebraic identity, so uniform flow is a
machine-precision fixed point, not a truncation-order one:

- DGCL test: uniform (0.7, -0.4) on a 16x12 grid with interior lines
  wiggling out of phase, 400 steps: max deviation 7.9e-15 (~35 ulp).
  Negative control with end-of-step areas (per-step cell error exactly
  dw*dh/V, the cross term the identity absorbs): 1.5e-2 - a 1e12
  separation, so the test can fail.
- Degeneracy: zero motion on a uniform grid vs fixed-grid PISO over
  Taylor-Green steps: max difference 2.2e-16 - one ulp - pinning every
  geometric generalisation to the verified implementation.
- Physics under motion: Taylor-Green on the wiggling mesh, L2 error
  2.42e-2 -> 1.07e-2 (n=16 -> 32, order 1.17); moving-mesh error at
  n=32 sits below the fixed-mesh 1.1532e-2 (PISO's published value to
  four digits); energy decay unchanged by the motion.

One trap documented in the test: the projection's inner-stop floor
(0.1 * tolerance * reference_flux) at an engineering tolerance lets a
one-sweep partial p' accumulate into p, whose gradient perturbs the
velocities at ~1e-11 with the geometry blameless. The DGCL run must use
a rounding-level tolerance because machine-precision preservation is the
claim under test. Measured: 3.6e-11 at tol 1e-9, 7.9e-15 at 1e-13.

Incompressibility needs no mesh-velocity term: subtracting the GCL from
moving-cell mass conservation leaves plain div(u) = 0 on the current
geometry, so the projection is the fixed-grid one with non-uniform
coefficients.

292 rtx-cfd tests green (288 + 4).

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

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

@@ -0,0 +1,721 @@
+//! 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,
+}
+
+/// 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,
+}
+
+impl Default for AleParameters {
+    fn default() -> Self {
+        Self {
+            corrector_steps: 2,
+            tolerance: 1e-6,
+            swept_face_rule: SweptFaceRule::Trapezoidal,
+        }
+    }
+}
+
+/// 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.
+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`.
+    fn apply_boundary_normals(&self, field: &mut AleField, t: f64, x: &[f64], y: &[f64]) {
+        let (nx, ny) = (field.nx, field.ny);
+        let yc = centres(y);
+        let xc = centres(x);
+        for j in 0..ny {
+            field.u[(j, 0)] = self.boundary(x[0], yc[j], t).0;
+            field.u[(j, nx)] = self.boundary(x[nx], yc[j], t).0;
+        }
+        for i in 0..nx {
+            field.v[(0, i)] = self.boundary(xc[i], y[0], t).1;
+            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 {
+                    self.boundary(xo[i], yo[ny], t_old).0
+                };
+                let phi_s = if q_s >= 0.0 {
+                    uo[(j, i)]
+                } else if j > 0 {
+                    uo[(j - 1, i)]
+                } else {
+                    self.boundary(xo[i], yo[0], t_old).0
+                };
+                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 {
+                    let u_wall = self.boundary(xo[i], yo[ny], t_old).0;
+                    nu * (u_wall - uo[(j, i)]) / (yo[ny] - yco[j]) * w_o
+                };
+                let d_s = if j > 0 {
+                    nu * (uo[(j - 1, i)] - uo[(j, i)]) / (yco[j] - yco[j - 1]) * w_o
+                } else {
+                    let u_wall = self.boundary(xo[i], yo[0], t_old).0;
+                    nu * (u_wall - uo[(j, i)]) / (yco[j] - yo[0]) * w_o
+                };
+                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 {
+                    self.boundary(xo[nx], yo[j], t_old).1
+                };
+                let phi_w = if q_w >= 0.0 {
+                    vo[(j, i)]
+                } else if i > 0 {
+                    vo[(j, i - 1)]
+                } else {
+                    self.boundary(xo[0], yo[j], t_old).1
+                };
+                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 {
+                    let v_wall = self.boundary(xo[nx], yo[j], t_old).1;
+                    nu * (v_wall - vo[(j, i)]) / (xo[nx] - xco[i]) * h_o
+                };
+                let d_w = if i > 0 {
+                    nu * (vo[(j, i - 1)] - vo[(j, i)]) / (xco[i] - xco[i - 1]) * h_o
+                } else {
+                    let v_wall = self.boundary(xo[0], yo[j], t_old).1;
+                    nu * (v_wall - vo[(j, i)]) / (xco[i] - xo[0]) * h_o
+                };
+                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;
+            }
+        }
+
+        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 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 {
+                    if 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, where the normal velocity is data.
+                    let ae = if i + 1 == nx {
+                        0.0
+                    } else {
+                        dt * dy_j / (xcn[i + 1] - xcn[i])
+                    };
+                    let aw = if i == 0 {
+                        0.0
+                    } else {
+                        dt * dy_j / (xcn[i] - xcn[i - 1])
+                    };
+                    let an = if j + 1 == ny {
+                        0.0
+                    } else {
+                        dt * dx_i / (ycn[j + 1] - ycn[j])
+                    };
+                    let as_ = if j == 0 {
+                        0.0
+                    } else {
+                        dt * dx_i / (ycn[j] - ycn[j - 1])
+                    };
+                    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;
+            }
+        }
+        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` (which must keep the domain endpoints fixed and the
+    /// lines strictly increasing — the motion may not invert a cell).
+    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 eps_x = 1e-12 * (field.x[field.nx] - field.x[0]).abs();
+        let eps_y = 1e-12 * (field.y[field.ny] - field.y[0]).abs();
+        if (new_x[0] - field.x[0]).abs() > eps_x
+            || (new_x[field.nx] - field.x[field.nx]).abs() > eps_x
+            || (new_y[0] - field.y[0]).abs() > eps_y
+            || (new_y[field.ny] - field.y[field.ny]).abs() > eps_y
+        {
+            return Err(CfdError::invalid_parameter(
+                "domain boundary must stay fixed: only interior node lines may move",
+            ));
+        }
+
+        let t_old = self.time;
+        let t_new = t_old + dt;
+
+        // Boundary data at the start of the step, on the start-of-step
+        // geometry: this is what the explicit predictor differentiates.
+        let (x0, y0) = (field.x.clone(), field.y.clone());
+        self.apply_boundary_normals(field, t_old, &x0, &y0);
+
+        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,
+        })
+    }
+}

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

@@ -9,6 +9,8 @@ use crate::{CfdConfig, CfdError, CfdResult};
 // use nalgebra::{DMatrix, DVector};
 // use std::collections::HashMap;
 
+/// ALE solver on a moving tensor-product staggered grid
+pub mod ale;
 /// Boundary conditions
 pub mod boundary_conditions;
 /// Flow field data structures
@@ -25,6 +27,7 @@ pub mod simple;
 pub mod simple_gpu;
 
 // Re-export main types
+pub use ale::{AleField, AleParameters, AlePisoSolver, AleResult, SweptFaceRule};
 pub use boundary_conditions::{
     BoundaryCondition, BoundaryConditions, BoundaryLocation, BoundaryType,
 };