rtx-cfd: overset A-P1 GATED — the curvilinear patch moves and deforms under an exact 2-D DGCL

StepGeometry (motion.rs): time-averaged face vectors S̄_f = ½(S^n + S^{n+1}) and
swept volumes δV_f = S̄_f·δc_f — exact for linear node motion on any quad, so
Σ sign δV_f = V^{n+1} − V^n is algebra (1.8e-14 measured; the EndOfStep control
2.1e-3). CurvilinearPisoSolver::set_mesh(next) names the end-of-step geometry;
advance swaps it in, rebuilds operators + pressure matrix on it (L_f, LSQ
gradients, no mesh-velocity term in the projection), keeps the old mesh for
V^n and the explicit boundary data; predictor in the conservative ALE form
V^{n+1} û = V^n u^n + dt(−Σ sign (F − δV/dt) u_f + ν D + f V^n), written as
u·(V^n/V^{n+1}) + … so a stationary mesh is bitwise the static path; fluxes on
S̄_f; snapshot carries both meshes; swept_face_rule knob (Trapezoidal default,
EndOfStep = negative control). Stokes limit keeps the mesh flux (was dropped
with convection) and centres it (upwinding it cost an order: 1.06/1.00).

Gates (tests/curvilinear_ale.rs, 11 tests, 83 s): uniform flow on a wiggling
AND bending annulus 4.44e-15 over 400 steps, p exactly 0, 0 pressure
iterations; control deviates 4.9e-5; stationary mesh through the moving path
bit-identical (both diffusion variants); snapshot/restore on the moving mesh
bit-identical; Taylor–Green orders unchanged — upwind 0.995/0.976 vs fixed
0.987/0.975 at 1.05× error, Stokes 1.92/1.97 vs 1.94/1.99 at 2.4×, moving
annulus 2.11/2.02; linear-field falsifier 1.95/1.92 (annulus), 1.91/1.43
(square, sliding wall nodes). P0 ladders re-run identical to every digit.

Rule from the diagnosis: start a moving run ON the t = 0 mesh and sweep less
than a cell per step — a first step that jumped 2–4 cells imprinted an
O(displacement) error no refinement removed (dt- and motion-independent).

Co-Authored-By: Claude Fable 5.1 <[email protected]>
Claude-Session: https://claude.ai/code/session_01X2GmJXeQ2njUecEKiJZ1G2

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

@@ -19,16 +19,29 @@
 //! re-project the STORED fluxes. `L_f` is the orthogonal difference plus a
 //! node-based tangential correction (a 9-point stencil), assembled and
 //! applied by the same code.
+//!
+//! Phase P1: the patch MOVES. [`CurvilinearPisoSolver::set_mesh`] names the
+//! geometry at the end of the coming step; the step then runs in the
+//! conservative ALE form `(V^{n+1} û − V^n u^n)/dt = …` with the
+//! convecting flux relative to the mesh, `F_f − δV_f/dt`, on the
+//! time-averaged face vectors `S̄_f` (`motion.rs`: the 2-D discrete GCL is
+//! exact algebra under that rule). `L_f`, the cell gradients and the
+//! pressure matrix live on the end-of-step geometry; the projection is the
+//! static one (no mesh-velocity term). A stationary mesh through this path
+//! is bit-identical to the static path.
 
+mod motion;
 mod operators;
 mod predictor;
 mod projection;
 
+pub use motion::StepGeometry;
 pub use operators::Operators;
 
 use crate::mesh::{PatchMesh, PatchSide};
+use crate::solvers::incompressible::ale::SweptFaceRule;
 use crate::solvers::incompressible::sparse_bicgstab::CsrMatrix;
-use crate::{CfdConfig, CfdResult};
+use crate::{CfdConfig, CfdError, CfdResult};
 
 /// What one side of the patch is.
 #[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
@@ -110,6 +123,10 @@ pub struct CurvilinearParameters {
     pub boundaries: PatchBoundaries,
     /// BiCGSTAB iteration cap.
     pub max_poisson_iterations: usize,
+    /// Face-area rule on a moving mesh: `Trapezoidal` (the DGCL-exact
+    /// choice, the default) or `EndOfStep` (the negative control, GCL-
+    /// violating). Irrelevant when the mesh does not move.
+    pub swept_face_rule: SweptFaceRule,
 }
 
 impl Default for CurvilinearParameters {
@@ -121,6 +138,7 @@ impl Default for CurvilinearParameters {
             normal_diffusion: NormalDiffusion::Explicit,
             boundaries: PatchBoundaries::default(),
             max_poisson_iterations: 5000,
+            swept_face_rule: SweptFaceRule::Trapezoidal,
         }
     }
 }
@@ -166,10 +184,15 @@ pub struct CurvilinearResult {
     pub boundary_flux_adjustment: f64,
 }
 
-/// Restorable solver state (the coupling re-runs a step).
+/// Restorable solver state (the coupling re-runs a step): the time, the
+/// current mesh and any mesh already named for the next step. Restoring
+/// rebuilds the operators and drops the cached matrix, so a re-run from
+/// the snapshot is bit-identical to the first run.
 #[derive(Debug, Clone)]
 pub struct CurvilinearSolverState {
     time: f64,
+    mesh: PatchMesh,
+    pending: Option<PatchMesh>,
 }
 
 type VelocityFn = Box<dyn Fn(f64, f64, f64) -> (f64, f64) + Send + Sync>;
@@ -179,6 +202,8 @@ pub struct CurvilinearPisoSolver {
     config: CfdConfig,
     params: CurvilinearParameters,
     mesh: PatchMesh,
+    /// The mesh the next `advance` ends on (`set_mesh`), if it moves.
+    pending: Option<PatchMesh>,
     ops: Operators,
     boundary_velocity: Option<VelocityFn>,
     momentum_source: Option<VelocityFn>,
@@ -198,6 +223,7 @@ impl CurvilinearPisoSolver {
             config,
             params,
             mesh,
+            pending: None,
             ops,
             boundary_velocity: None,
             momentum_source: None,
@@ -220,10 +246,37 @@ impl CurvilinearPisoSolver {
     {
         self.momentum_source = Some(Box::new(f));
     }
-    /// The mesh.
+    /// The mesh the field currently lives on (the end of the last step).
     pub fn mesh(&self) -> &PatchMesh {
         &self.mesh
     }
+    /// The mesh named for the end of the next step, if any.
+    pub fn next_mesh(&self) -> Option<&PatchMesh> {
+        self.pending.as_ref()
+    }
+    /// Name the geometry the NEXT `advance` ends on. Same topology as the
+    /// current mesh (`ns`, `nn`, periodicity); the node motion between the
+    /// two is taken as linear in time. Calling it again before `advance`
+    /// replaces the earlier choice (the coupling loop re-tries a step);
+    /// not calling it leaves the mesh where it is.
+    pub fn set_mesh(&mut self, next: PatchMesh) -> CfdResult<()> {
+        if next.ns() != self.mesh.ns()
+            || next.nn() != self.mesh.nn()
+            || next.periodic().is_some() != self.mesh.periodic().is_some()
+        {
+            return Err(CfdError::mesh(format!(
+                "set_mesh: topology changed ({}x{}, periodic {}) -> ({}x{}, periodic {})",
+                self.mesh.ns(),
+                self.mesh.nn(),
+                self.mesh.periodic().is_some(),
+                next.ns(),
+                next.nn(),
+                next.periodic().is_some()
+            )));
+        }
+        self.pending = Some(next);
+        Ok(())
+    }
     /// The operators.
     pub fn operators(&self) -> &Operators {
         &self.ops
@@ -244,13 +297,21 @@ impl CurvilinearPisoSolver {
     pub fn set_time(&mut self, t: f64) {
         self.time = t;
     }
-    /// Capture the state.
+    /// Capture the state (time and both meshes).
     pub fn snapshot(&self) -> CurvilinearSolverState {
-        CurvilinearSolverState { time: self.time }
+        CurvilinearSolverState {
+            time: self.time,
+            mesh: self.mesh.clone(),
+            pending: self.pending.clone(),
+        }
     }
-    /// Restore a captured state.
+    /// Restore a captured state: operators rebuilt, matrix cache dropped.
     pub fn restore(&mut self, state: &CurvilinearSolverState) {
         self.time = state.time;
+        self.mesh = state.mesh.clone();
+        self.pending = state.pending.clone();
+        self.ops = Operators::new(&self.mesh, &self.params.boundaries);
+        self.matrix = None;
     }
 
     pub(crate) fn boundary_velocity(&self, x: f64, y: f64, t: f64) -> (f64, f64) {
@@ -278,8 +339,15 @@ impl CurvilinearPisoSolver {
             field.v[c] = v;
         }
         let zero = vec![0.0; mesh.cell_count()];
-        field.flux =
-            self.predicted_fluxes(&field.u.clone(), &field.v.clone(), &zero, 0.0, self.time);
+        let geo = StepGeometry::stationary(mesh);
+        field.flux = self.predicted_fluxes(
+            &field.u.clone(),
+            &field.v.clone(),
+            &zero,
+            0.0,
+            self.time,
+            &geo,
+        );
     }
 
     /// Advance one step of `dt`.
@@ -291,10 +359,28 @@ impl CurvilinearPisoSolver {
         let t_old = self.time;
         let t_new = t_old + dt;
         let rho = self.config.density;
+
+        // The mesh moves: `self.mesh` becomes the end-of-step geometry
+        // (operators and matrix follow it); the start-of-step geometry is
+        // kept for this step's volumes and swept faces only.
+        let old = match self.pending.take() {
+            Some(next) => {
+                let old = std::mem::replace(&mut self.mesh, next);
+                self.ops = Operators::new(&self.mesh, &self.params.boundaries);
+                self.matrix = None;
+                Some(old)
+            }
+            None => None,
+        };
+        let geo = match &old {
+            Some(o) => StepGeometry::new(o, &self.mesh, self.params.swept_face_rule),
+            None => StepGeometry::stationary(&self.mesh),
+        };
+        let old_mesh = old.as_ref().unwrap_or(&self.mesh);
         let mesh = &self.mesh;
 
-        let (uh, vh) = self.predict(field, dt, t_old);
-        let mut flux = self.predicted_fluxes(&uh, &vh, &field.p, dt, t_new);
+        let (uh, vh) = self.predict(field, dt, t_old, old_mesh, &geo);
+        let mut flux = self.predicted_fluxes(&uh, &vh, &field.p, dt, t_new, &geo);
         let adjustment = self.adjust_boundary_flux(&mut flux);
         for c in 0..mesh.cell_count() {
             let g = self.pressure_gradient(&field.p, c);

crates/specialized/rtx-cfd/src/solvers/incompressible/curvilinear/motion.rs

@@ -0,0 +1,106 @@
+//! The moving patch (overset A-P1): the geometry one step consumes when the
+//! mesh at the end of the step differs from the mesh at its start.
+//!
+//! In two dimensions, for a face whose end nodes move linearly in time from
+//! `x^n` to `x^{n+1}`, the area it sweeps is EXACTLY
+//!
+//! ```text
+//!   δV_f = S̄_f · (c_f^{n+1} − c_f^n),   S̄_f = ½ (S_f^n + S_f^{n+1}),
+//! ```
+//!
+//! (`S_f(τ) · ċ_f` is linear in `τ`, so the midpoint rule integrates it
+//! exactly, and `S_f` is linear in the node coordinates, so the midpoint
+//! configuration's area vector is the average). Summed with the outward
+//! signs around a cell the same integrand is `dV/dτ`, hence
+//!
+//! ```text
+//!   Σ_f sign_f δV_f = V^{n+1} − V^n
+//! ```
+//!
+//! is an algebraic identity for every quadrilateral — the discrete
+//! geometric conservation law of Thomas & Lombard (1979), in the form
+//! Farhat, Geuzaine & Grandmont (2001) prove needs only one configuration
+//! in 2-D. The fluid flux uses the same `S̄_f`, so a uniform flow `U` gives
+//! `Σ_f sign_f F_f = U · Σ_f sign_f S̄_f = 0` (a closed polygon) and its
+//! relative flux `F_f − δV_f/dt` convects exactly `−U (V^{n+1} − V^n)/dt`,
+//! which the conservative update `(V^{n+1} u − V^n u)/dt` absorbs: uniform
+//! flow is a fixed point to rounding on any admissible motion
+//! (`tests/curvilinear_ale.rs`). Mass conservation for the moving cell,
+//! `ΔV/dt + Σ (u − w)·S = 0`, minus that identity leaves `Σ_f sign_f F_f =
+//! 0`: the projection has no mesh-velocity term (`ale.rs` module docs).
+//!
+//! `SweptFaceRule::EndOfStep` — `S_f^{n+1}` in both places — is kept as the
+//! negative control: it is first-order consistent and violates the
+//! identity by the cross term of the node displacements, visibly.
+
+use crate::mesh::PatchMesh;
+use crate::solvers::incompressible::ale::SweptFaceRule;
+
+/// Face geometry for one step from `old` to `new`.
+#[derive(Debug, Clone)]
+pub struct StepGeometry {
+    /// Face area vector the fluxes see (`S̄_f`, or `S_f^{n+1}` under the
+    /// control rule); `S_f` itself when the mesh is stationary.
+    pub s_bar: Vec<[f64; 2]>,
+    /// Volume swept by each face, oriented +s / +n; exactly `0.0` when the
+    /// mesh is stationary.
+    pub swept: Vec<f64>,
+    /// Whether `old` and `new` are the same geometry (the static path).
+    pub stationary: bool,
+}
+
+impl StepGeometry {
+    /// The stationary geometry of `mesh`.
+    pub fn stationary(mesh: &PatchMesh) -> Self {
+        Self {
+            s_bar: mesh.faces().iter().map(|f| f.s).collect(),
+            swept: vec![0.0; mesh.faces().len()],
+            stationary: true,
+        }
+    }
+
+    /// The geometry of a step from `old` to `new` (same topology). With
+    /// `old` and `new` identical the result is bitwise the stationary one:
+    /// `½ (S + S) = S` and `S̄ · 0 = 0` exactly.
+    pub fn new(old: &PatchMesh, new: &PatchMesh, rule: SweptFaceRule) -> Self {
+        assert_eq!(
+            old.faces().len(),
+            new.faces().len(),
+            "meshes differ in topology"
+        );
+        let mut s_bar = Vec::with_capacity(new.faces().len());
+        let mut swept = Vec::with_capacity(new.faces().len());
+        for (fo, fn_) in old.faces().iter().zip(new.faces()) {
+            let s = match rule {
+                SweptFaceRule::Trapezoidal => {
+                    [0.5 * (fo.s[0] + fn_.s[0]), 0.5 * (fo.s[1] + fn_.s[1])]
+                }
+                SweptFaceRule::EndOfStep => fn_.s,
+            };
+            let dc = [fn_.centre[0] - fo.centre[0], fn_.centre[1] - fo.centre[1]];
+            s_bar.push(s);
+            swept.push(s[0] * dc[0] + s[1] * dc[1]);
+        }
+        Self {
+            s_bar,
+            swept,
+            stationary: false,
+        }
+    }
+
+    /// The largest cell violation of the discrete GCL, `|Σ_f sign_f δV_f −
+    /// (V^{n+1} − V^n)|`, relative to the cell volume. Rounding-level under
+    /// the trapezoidal rule; the displacement cross term under the control.
+    pub fn gcl_defect(&self, old: &PatchMesh, new: &PatchMesh) -> f64 {
+        (0..new.cell_count())
+            .map(|c| {
+                let sum: f64 = new
+                    .cell_faces(c)
+                    .iter()
+                    .map(|&(f, sign)| sign * self.swept[f])
+                    .sum();
+                (sum - (new.area(c) - old.area(c))).abs() / new.area(c)
+            })
+            .fold(0.0, f64::max)
+    }
+}

crates/specialized/rtx-cfd/src/solvers/incompressible/curvilinear/predictor.rs

@@ -1,13 +1,35 @@
 //! The momentum predictor: explicit forward Euler from the old field
 //! (matching the background PISO), or line-implicit in `n` for the
 //! orthogonal part of the across-patch diffusion.
+//!
+//! On a moving mesh the update is the conservative ALE form
+//!
+//! ```text
+//!   V^{n+1} û = V^n u^n + dt (−Σ_f sign_f (F_f − δV_f/dt) u_f + ν D(u^n) + f V^n),
+//! ```
+//!
+//! with the relative flux deciding the upwind side, `D` and (when line-
+//! implicit) the implicit operator on the end-of-step geometry, and the
+//! boundary data of the explicit terms at the start-of-step face centres.
+//! Written as `u^n (V^n/V^{n+1}) + …/V^{n+1}` so a stationary mesh (ratio
+//! exactly `1.0`, `δV_f` exactly `0.0`) reproduces the static path bitwise.
 
-use super::{CurvilinearPisoSolver, NormalDiffusion, PatchConvection, PatchField, SideBc};
-use crate::mesh::PatchSide;
+use super::{
+    CurvilinearPisoSolver, NormalDiffusion, PatchConvection, PatchField, SideBc, StepGeometry,
+};
+use crate::mesh::{PatchMesh, PatchSide};
 
 impl CurvilinearPisoSolver {
-    /// `û = u^n + dt (−C(F^n, u^n) + ν D(u^n) + f/ρ)`, no pressure.
-    pub(super) fn predict(&self, field: &PatchField, dt: f64, t_old: f64) -> (Vec<f64>, Vec<f64>) {
+    /// `û` without the pressure; `old` is the start-of-step mesh (`self.mesh`
+    /// itself when stationary), `geo` the step's face geometry.
+    pub(super) fn predict(
+        &self,
+        field: &PatchField,
+        dt: f64,
+        t_old: f64,
+        old: &PatchMesh,
+        geo: &StepGeometry,
+    ) -> (Vec<f64>, Vec<f64>) {
         let mesh = &self.mesh;
         let n = mesh.cell_count();
         let rho = self.config.density;
@@ -31,21 +53,43 @@ impl CurvilinearPisoSolver {
             for (f, sign) in mesh.cell_faces(c) {
                 let face = &mesh.faces()[f];
                 let side = mesh.side(f);
-                let bval = side.and_then(|s| bvel(s, face.centre));
-                // Convection: outward flux times the upwind face value.
-                if matches!(self.params.convection, PatchConvection::Upwind) {
-                    let out = sign * field.flux[f];
+                let bval = side.and_then(|s| bvel(s, old.faces()[f].centre));
+                // Convection: outward flux RELATIVE to the moving face times
+                // the upwind face value. In the Stokes limit the fluid flux
+                // is dropped but the mesh flux stays: the conservative
+                // update needs `−Σ sign δV_f u_f` whenever the mesh moves.
+                let fluid = match self.params.convection {
+                    PatchConvection::Upwind => field.flux[f],
+                    PatchConvection::None => 0.0,
+                };
+                let out = sign * (fluid - geo.swept[f] / dt);
+                if out != 0.0 {
                     let (uf, vf) = match (face.owner, face.neigh) {
-                        (Some(p), Some(q)) => {
-                            let up = if out >= 0.0 {
-                                c
-                            } else if p == c {
-                                q
-                            } else {
-                                p
-                            };
-                            (field.u[up], field.v[up])
-                        }
+                        (Some(p), Some(q)) => match self.params.convection {
+                            PatchConvection::Upwind => {
+                                let up = if out >= 0.0 {
+                                    c
+                                } else if p == c {
+                                    q
+                                } else {
+                                    p
+                                };
+                                (field.u[up], field.v[up])
+                            }
+                            // The Stokes limit has no upwind scheme to be
+                            // consistent with; the mesh flux takes the linear
+                            // face value and keeps its second order (upwinding
+                            // it measured first order on moving Taylor–Green:
+                            // 6.59e-3 / 3.17e-3 / 1.58e-3 against a
+                            // second-order fixed ladder).
+                            PatchConvection::None => {
+                                let w = face.w;
+                                (
+                                    w * field.u[p] + (1.0 - w) * field.u[q],
+                                    w * field.v[p] + (1.0 - w) * field.v[q],
+                                )
+                            }
+                        },
                         _ => match bval {
                             Some(b) => b,
                             None => (field.u[c], field.v[c]),
@@ -92,19 +136,28 @@ impl CurvilinearPisoSolver {
                 dv += lv;
             }
             let a = mesh.area(c);
-            let (fx, fy) = self.source_at(mesh.centre(c), t_old);
-            uh[c] = field.u[c] + dt * ((-cu + nu * du) / a + fx / rho);
-            vh[c] = field.v[c] + dt * ((-cv + nu * dv) / a + fy / rho);
+            let ratio = old.area(c) / a;
+            let (fx, fy) = self.source_at(old.centre(c), t_old);
+            uh[c] = field.u[c] * ratio + dt * ((-cu + nu * du) / a + fx / rho * ratio);
+            vh[c] = field.v[c] * ratio + dt * ((-cv + nu * dv) / a + fy / rho * ratio);
         }
 
         if implicit_n {
-            self.solve_lines(&mut uh, &mut vh, dt, nu, t_old);
+            self.solve_lines(&mut uh, &mut vh, dt, nu, t_old, old);
         }
         (uh, vh)
     }
 
-    /// `(I − dt ν L_n/A) û = rhs` along every s-line, Thomas algorithm.
-    fn solve_lines(&self, uh: &mut [f64], vh: &mut [f64], dt: f64, nu: f64, t_old: f64) {
+    /// `(I − dt ν L_n/V^{n+1}) û = rhs` along every s-line, Thomas algorithm.
+    fn solve_lines(
+        &self,
+        uh: &mut [f64],
+        vh: &mut [f64],
+        dt: f64,
+        nu: f64,
+        t_old: f64,
+        old: &PatchMesh,
+    ) {
         let mesh = &self.mesh;
         let (ns, nn) = (mesh.ns(), mesh.nn());
         let mut lower = vec![0.0; nn];
@@ -134,7 +187,7 @@ impl CurvilinearPisoSolver {
                         }
                         Some(s) => match self.params.boundaries.get(s) {
                             SideBc::Velocity => {
-                                let face = &mesh.faces()[f];
+                                let face = &old.faces()[f];
                                 let b =
                                     self.boundary_velocity(face.centre[0], face.centre[1], t_old);
                                 d += coef;

crates/specialized/rtx-cfd/src/solvers/incompressible/curvilinear/projection.rs

@@ -3,7 +3,7 @@
 //! pressure-correction equation on the 9-point operator, and the flux and
 //! velocity corrections.
 
-use super::{CurvilinearPisoSolver, PatchField, SideBc};
+use super::{CurvilinearPisoSolver, PatchField, SideBc, StepGeometry};
 use crate::mesh::PatchSide;
 use crate::solvers::incompressible::sparse_bicgstab::{
     BicgstabResult, CsrMatrix, bicgstab_jacobi, project_mean,
@@ -48,8 +48,9 @@ impl CurvilinearPisoSolver {
         })
     }
 
-    /// `F* = interp(û)·S − (dt/ρ) L_f(p^n)` on interior and outlet faces,
-    /// the prescribed flux on velocity faces (at `t_new`).
+    /// `F* = interp(û)·S̄ − (dt/ρ) L_f(p^n)` on interior and outlet faces,
+    /// the prescribed flux `u_b · S̄` on velocity faces (at `t_new`, on the
+    /// end-of-step face centres). `S̄` is the step's face vector (`geo`).
     pub(super) fn predicted_fluxes(
         &self,
         uh: &[f64],
@@ -57,6 +58,7 @@ impl CurvilinearPisoSolver {
         p: &[f64],
         dt: f64,
         t_new: f64,
+        geo: &StepGeometry,
     ) -> Vec<f64> {
         let mesh = &self.mesh;
         let rho = self.config.density;
@@ -64,22 +66,25 @@ impl CurvilinearPisoSolver {
         mesh.faces()
             .iter()
             .enumerate()
-            .map(|(f, face)| match (face.owner, face.neigh) {
-                (Some(o), Some(n)) => {
-                    let w = face.w;
-                    let uf = w * uh[o] + (1.0 - w) * uh[n];
-                    let vf = w * vh[o] + (1.0 - w) * vh[n];
-                    uf * face.s[0] + vf * face.s[1] - dt / rho * lp[f]
-                }
-                _ => {
-                    let c = mesh.boundary_cell(f);
-                    match self.params.boundaries.get(mesh.side(f).expect("boundary")) {
-                        SideBc::Velocity => {
-                            let (ub, vb) =
-                                self.boundary_velocity(face.centre[0], face.centre[1], t_new);
-                            ub * face.s[0] + vb * face.s[1]
+            .map(|(f, face)| {
+                let s = geo.s_bar[f];
+                match (face.owner, face.neigh) {
+                    (Some(o), Some(n)) => {
+                        let w = face.w;
+                        let uf = w * uh[o] + (1.0 - w) * uh[n];
+                        let vf = w * vh[o] + (1.0 - w) * vh[n];
+                        uf * s[0] + vf * s[1] - dt / rho * lp[f]
+                    }
+                    _ => {
+                        let c = mesh.boundary_cell(f);
+                        match self.params.boundaries.get(mesh.side(f).expect("boundary")) {
+                            SideBc::Velocity => {
+                                let (ub, vb) =
+                                    self.boundary_velocity(face.centre[0], face.centre[1], t_new);
+                                ub * s[0] + vb * s[1]
+                            }
+                            SideBc::Outlet => uh[c] * s[0] + vh[c] * s[1] - dt / rho * lp[f],
                         }
-                        SideBc::Outlet => uh[c] * face.s[0] + vh[c] * face.s[1] - dt / rho * lp[f],
                     }
                 }
             })

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

@@ -46,7 +46,7 @@ pub use boundary_conditions::{
 };
 pub use curvilinear::{
     CurvilinearParameters, CurvilinearPisoSolver, CurvilinearResult, CurvilinearSolverState,
-    NormalDiffusion, Operators, PatchBoundaries, PatchConvection, PatchField, SideBc,
+    NormalDiffusion, Operators, PatchBoundaries, PatchConvection, PatchField, SideBc, StepGeometry,
 };
 pub use embedded::{EmbeddedParameters, EmbeddedPisoSolver, EmbeddedResult, EmbeddedSolverState};
 pub use embedded_body::{