rtx-cfd: repair the pressure-velocity coupling, LBM walls and mesh quality

Clears the rest of the quarantine. All three crates now run 558 tests
with 0 failures and no `#[ignore]` markers.

SIMPLE could not converge, and the reason was not slow convergence but
wrong physics.

The pressure correction equation used a bare Laplacian, 1/dx^2 and
1/dy^2, while the velocity correction divided by a_p = rho dx dy / dt.
SIMPLE requires these to be each other's inverse: substituting the
corrected velocities into continuity must reproduce the pressure
equation, which fixes a_E = rho d dy/dx with d = dV/a_p. The two
disagreed by roughly 1/(h^2 dt) -- about 2e4 on a 16x16 cavity -- so the
pressure correction was that many times too weak to enforce continuity.

The consequence was visible and specific. A lid-driven cavity at Re=100
produced a monotonic profile rising from 0 at the floor to 1 at the lid:
Couette flow, with no recirculation anywhere, and a peak pressure of
1.6e-4 against the rho U^2 scale of 1. The return flow in a cavity is
driven entirely by the pressure gradient, so with the pressure pinned
near zero there was nothing to turn the flow around. With the
coefficients made consistent the profile recirculates, the peak pressure
is 2.9, and the solver converges.

Also in SIMPLE:
  - `p'` was never reset between outer iterations. It is a correction
    that `pressure_update_step` folds into `p`, so carrying it forward
    applied the same correction twice.
  - The convergence measure was the inner Gauss-Seidel residual, which
    goes to zero whether or not the flow satisfies continuity. Now the
    mass imbalance.
  - The velocity correction used only the transient part of a_p,
    `rho dV/dt`, rather than the diagonal the momentum equation was
    actually solved with.
  - All four convective face fluxes were computed from a single
    cell-centred velocity, so `fe` and `fw` were the same number, as were
    `fn` and `fs`. Upwinding then picked the same direction on opposite
    faces of the control volume. Now interpolated per face on the
    staggered grid.

Not claimed: agreement with Ghia, Ghia & Shin (1982). The vortex centre
moves toward their y = 0.4531 under refinement (0.400 at 16^2, 0.419 at
32^2, 0.460 at 64^2) but the minimum centreline velocity reaches only
-0.130 against their -0.2109, and the converged field still depends
slightly on the pseudo-time step, which a true steady state cannot. The
cavity test therefore asserts what is established -- convergence,
recirculation, vortex position, and an O(1) pressure field -- and the
remaining gap is recorded in omni-cortex/docs/solver_status.md rather
than papered over with a loose tolerance.

LBM bounce-back was doing neither of the things its name claims. It was
written as assignment (`f[2] = f[4]`) rather than a swap, discarding the
population being reflected -- bounce-back is a permutation and conserves
mass exactly, so the domain leaked 0.013% of its mass every 100 steps and
would have kept draining. And the pairs used were 5<->8 and 6<->7, which
reverse only the wall-normal component: that is specular reflection, a
free-slip wall, so the no-slip condition the walls were supposed to
impose never held.

Mesh quality:
  - Quadrilateral aspect ratio included the diagonals in the maximum but
    not the minimum, so it could never return 1: a unit square reported
    sqrt(2) and a 2:1 rectangle sqrt(5).
  - Triangle aspect ratio used longest-over-shortest edge, which does not
    detect the failure mode that matters. A sliver with vertices (0,0),
    (10,0), (5,0.1) scores 2.0 -- indistinguishable from a healthy 2:1
    triangle -- while its area is a twentieth of what its edges suggest.
    Now the radius ratio R/2r, which is 1 for equilateral and 1250 for
    that sliver, and which also fixes the quality histogram.
  - StructuredMesh aspect ratio took bounding-box extents and guarded the
    z-extent with `.max(1e-10)`. On a 2-D mesh the depth is exactly zero,
    so the guard became the minimum and a unit square reported 2e10.

Mesh refinement produced meshes that failed their own validation.
`subdivide_triangle` reserved midpoint ids as `next_node_id + k`, then
advanced the counter by 3, after which `refine_cells` called `add_node`
and advanced it three more -- so every refined cell referenced vertices
three ids away from the ones actually created. Separately, the position
lookup selected by slot rather than by id ("This is simplified, should
look up correct midpoint"), so three of four sub-triangles had their
areas computed from the wrong points; the quadrilateral version mapped
every new id to the cell centre.

Fixtures corrected rather than tolerances loosened: a structured mesh
test asserted 0.16 for the average cell volume while the comment beside
it computed 0.25 from the node-count convention the code actually uses;
the Zou-He pressure test built a *velocity* boundary at u = 1.2, far
above the lattice speed of sound, making the density negative; and the
cavity-setup test required the lid to influence the domain centre 16
rows away in 10 steps, which exceeds the lattice propagation speed.

Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>

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

@@ -205,7 +205,9 @@ impl SimpleSolver {
         boundary_conditions: &BoundaryConditions,
         dt: f64,
     ) -> CfdResult<(f64, f64)> {
-        // Step 1: Solve momentum equations with current pressure field
+        // Step 1: Solve momentum equations with current pressure field.
+        // This begins by storing the current iterate in `u_old`, which the
+        // transient term, the under-relaxation and the residual all read.
         self.momentum_prediction_step(flow_field, dt).await?;
 
         // Step 2: Solve pressure correction equation
@@ -337,24 +339,80 @@ impl SimpleSolver {
         let (nx, ny, dx, dy) = flow_field.grid_info();
         let rho = self.config.density;
 
+        // The pressure correction starts from zero every outer iteration.
+        //
+        // `p'` is a correction to the current pressure field, not a field in
+        // its own right: `pressure_update_step` folds it into `p` at the end
+        // of the iteration, so carrying it into the next one applies the same
+        // correction twice.
+        flow_field.p_prime.fill(0.0);
+
         // Setup pressure correction equation: ∇²p' = ρ/Δt * ∇·u*
+        let mut mass_imbalance: f64 = 0.0;
         for j in 1..ny - 1 {
             for i in 1..nx - 1 {
                 // Compute mass source (continuity equation residual)
                 let mass_source = self.compute_mass_source(flow_field, i, j, dx, dy, rho)?;
                 flow_field.sp[(j, i)] = mass_source;
+                mass_imbalance += mass_source.abs();
+            }
+        }
+
+        // Neighbour coefficients, evaluated per face from the momentum
+        // equation's own diagonal at that face — the same `a_p` the velocity
+        // correction divides by, so the two remain each other's inverse.
+        // Cached because `a_p` depends on the velocity field, which does not
+        // change during the inner sweeps.
+        let mut coefficients = Vec::with_capacity((nx - 2) * (ny - 2));
+        for j in 1..ny - 1 {
+            for i in 1..nx - 1 {
+                let volume = dx * dy;
+                let d_east = volume
+                    / self.compute_u_momentum_center_coefficient(
+                        flow_field,
+                        i + 1,
+                        j,
+                        dx,
+                        dy,
+                        rho,
+                    )?;
+                let d_west = volume
+                    / self.compute_u_momentum_center_coefficient(flow_field, i, j, dx, dy, rho)?;
+                let d_north = volume
+                    / self.compute_v_momentum_center_coefficient(
+                        flow_field,
+                        i,
+                        j + 1,
+                        dx,
+                        dy,
+                        rho,
+                    )?;
+                let d_south = volume
+                    / self.compute_v_momentum_center_coefficient(flow_field, i, j, dx, dy, rho)?;
+
+                let ae = rho * d_east * dy / dx;
+                let aw = rho * d_west * dy / dx;
+                let an = rho * d_north * dx / dy;
+                let as_ = rho * d_south * dx / dy;
+                coefficients.push(MomentumEquationCoeffs {
+                    center: ae + aw + an + as_,
+                    east: ae,
+                    west: aw,
+                    north: an,
+                    south: as_,
+                    source: 0.0,
+                });
             }
         }
 
         // Solve pressure correction equation using Gauss-Seidel
-        let mut max_residual = 0.0;
-        for _iteration in 0..100 {
+        for _iteration in 0..200 {
             // Inner pressure correction iterations
             let mut residual = 0.0;
 
             for j in 1..ny - 1 {
                 for i in 1..nx - 1 {
-                    let coeffs = self.compute_pressure_coefficients(dx, dy)?;
+                    let coeffs = &coefficients[(j - 1) * (nx - 2) + (i - 1)];
 
                     let p_new = (flow_field.sp[(j, i)]
                         + coeffs.east * flow_field.p_prime[(j, i + 1)]
@@ -369,15 +427,19 @@ impl SimpleSolver {
                 }
             }
 
-            residual = residual.sqrt();
-            max_residual = residual;
-
-            if residual < 1e-8 {
+            if residual.sqrt() < 1e-10 {
                 break;
             }
         }
 
-        Ok(max_residual)
+        // Report the mass imbalance, not the inner Gauss-Seidel residual.
+        //
+        // The outer loop treats this as its convergence measure, and the
+        // inner residual only says how well the pressure-correction equation
+        // was solved — it goes to zero whether or not the flow satisfies
+        // continuity, so the solver could report convergence while the field
+        // was still divergent.
+        Ok(mass_imbalance)
     }
 
     /// Velocity correction step: correct velocities with pressure correction
@@ -536,12 +598,21 @@ impl SimpleSolver {
         let gamma_n = mu_eff / dy;
         let gamma_s = mu_eff / dy;
 
-        // Convection coefficients (using upwind)
-        let (u_center, v_center) = flow_field.get_velocity_at(i, j)?;
-        let fe = rho * u_center * dy; // East face mass flux
-        let fw = rho * u_center * dy; // West face mass flux
-        let fn_ = rho * v_center * dx; // North face mass flux
-        let fs = rho * v_center * dx; // South face mass flux
+        // Convective mass fluxes through the four faces of the u control
+        // volume, which on a staggered grid is centred on the u face `i` and
+        // spans from cell centre `i-1` to cell centre `i`. Its east and west
+        // faces therefore sit at those cell centres, where the velocity is the
+        // average of the two neighbouring u values.
+        //
+        // All four fluxes were previously taken from a single cell-centred
+        // velocity, so `fe` and `fw` were literally the same number, as were
+        // `fn_` and `fs`. Upwinding then chose the same direction on opposite
+        // faces of the volume, which cannot represent transport across it: the
+        // scheme reduced to diffusion plus a spurious diagonal term.
+        let fe = rho * 0.5 * (flow_field.u[(j, i)] + flow_field.u[(j, i + 1)]) * dy;
+        let fw = rho * 0.5 * (flow_field.u[(j, i - 1)] + flow_field.u[(j, i)]) * dy;
+        let fn_ = rho * 0.5 * (flow_field.v[(j + 1, i - 1)] + flow_field.v[(j + 1, i)]) * dx;
+        let fs = rho * 0.5 * (flow_field.v[(j, i - 1)] + flow_field.v[(j, i)]) * dx;
 
         // Compute coefficients with upwind scheme
         let ae = gamma_e + f64::max(-fe, 0.0);
@@ -552,8 +623,10 @@ impl SimpleSolver {
         // Time derivative coefficient
         let ap0 = rho * dx * dy / dt;
 
-        // Central coefficient
-        let ap = ae + aw + an + as_ + ap0;
+        // Central coefficient. The net flux term vanishes for a
+        // divergence-free field but is retained so the equation stays
+        // conservative while continuity is still being enforced.
+        let ap = ae + aw + an + as_ + (fe - fw) + (fn_ - fs) + ap0;
 
         // Source term (pressure gradient + old time step)
         let pressure_gradient = -(flow_field.p[(j, i)] - flow_field.p[(j, i - 1)]) * dy;
@@ -591,11 +664,14 @@ impl SimpleSolver {
         let gamma_n = mu_eff / dy;
         let gamma_s = mu_eff / dy;
 
-        let (u_center, v_center) = flow_field.get_velocity_at(i, j)?;
-        let fe = rho * u_center * dy;
-        let fw = rho * u_center * dy;
-        let fn_ = rho * v_center * dx;
-        let fs = rho * v_center * dx;
+        // Face fluxes for the v control volume, centred on the v face `j` and
+        // spanning cell centre `j-1` to cell centre `j`. Mirrors the u case
+        // above; see the note there on why a single cell-centred velocity for
+        // all four faces cannot represent transport.
+        let fn_ = rho * 0.5 * (flow_field.v[(j, i)] + flow_field.v[(j + 1, i)]) * dx;
+        let fs = rho * 0.5 * (flow_field.v[(j - 1, i)] + flow_field.v[(j, i)]) * dx;
+        let fe = rho * 0.5 * (flow_field.u[(j - 1, i + 1)] + flow_field.u[(j, i + 1)]) * dy;
+        let fw = rho * 0.5 * (flow_field.u[(j - 1, i)] + flow_field.u[(j, i)]) * dy;
 
         let ae = gamma_e + f64::max(-fe, 0.0);
         let aw = gamma_w + f64::max(fw, 0.0);
@@ -603,7 +679,7 @@ impl SimpleSolver {
         let as_ = gamma_s + f64::max(fs, 0.0);
 
         let ap0 = rho * dx * dy / dt;
-        let ap = ae + aw + an + as_ + ap0;
+        let ap = ae + aw + an + as_ + (fe - fw) + (fn_ - fs) + ap0;
 
         // Pressure gradient in y-direction
         let pressure_gradient = -(flow_field.p[(j, i)] - flow_field.p[(j - 1, i)]) * dx;
@@ -641,14 +717,46 @@ impl SimpleSolver {
         Ok(-mass_flux_imbalance) // Negative because we want ∇²p' = -∇·u*
     }
 
-    /// Compute coefficients for pressure correction equation
-    fn compute_pressure_coefficients(&self, dx: f64, dy: f64) -> CfdResult<MomentumEquationCoeffs> {
-        // Pressure correction equation: ∇²p' = S
-        // Standard 5-point stencil with unit coefficients
-        let ae = 1.0 / (dx * dx);
-        let aw = 1.0 / (dx * dx);
-        let an = 1.0 / (dy * dy);
-        let as_ = 1.0 / (dy * dy);
+    /// Compute coefficients for the pressure correction equation.
+    ///
+    /// The coefficients are not free: SIMPLE requires that substituting the
+    /// corrected velocities back into the continuity equation *reproduces*
+    /// this equation. With the velocity correction
+    /// `u_e = u*_e - d (p'_E - p'_P) / dx` and `d = ΔV / a_p`, continuity
+    /// gives
+    ///
+    /// ```text
+    ///     a_E = a_W = rho d dy/dx,   a_N = a_S = rho d dx/dy
+    /// ```
+    ///
+    /// so the neighbour coefficients carry `d` — the momentum equation's own
+    /// diagonal — and the mass imbalance in `compute_mass_source` is the
+    /// source. Any other scaling breaks the link between the pressure the
+    /// equation produces and the velocity correction it is supposed to drive.
+    ///
+    /// This previously used a bare Laplacian, `1/dx²` and `1/dy²`, while
+    /// `velocity_correction_step` divided by `a_p = rho dx dy / dt`. The two
+    /// therefore disagreed by a factor of roughly `1 / (h² dt)` — about
+    /// 2 x 10^4 on a 16 x 16 cavity — so the pressure correction was that many
+    /// times too small to enforce continuity. The consequence was not slow
+    /// convergence but the wrong physics: with the pressure field pinned near
+    /// zero, a lid-driven cavity produced a monotonic Couette profile with no
+    /// recirculation at all, since the return flow in a cavity is created
+    /// entirely by the pressure gradient.
+    fn compute_pressure_coefficients(
+        &self,
+        dx: f64,
+        dy: f64,
+        rho: f64,
+    ) -> CfdResult<MomentumEquationCoeffs> {
+        // d = ΔV / a_p, with a_p = rho dx dy / dt as used by the velocity
+        // correction, so d = dt / rho and `rho * d` is just the time step.
+        let rho_d = rho * (dx * dy) / (rho * dx * dy / self.parameters.time_step);
+
+        let ae = rho_d * dy / dx;
+        let aw = ae;
+        let an = rho_d * dx / dy;
+        let as_ = an;
         let ap = ae + aw + an + as_;
 
         Ok(MomentumEquationCoeffs {
@@ -661,32 +769,62 @@ impl SimpleSolver {
         })
     }
 
-    /// Compute center coefficient for u-momentum equation
+    /// Diagonal coefficient of the u-momentum equation, as used by the
+    /// velocity correction and the pressure equation.
+    ///
+    /// This must be the *same* `a_p` the momentum equation was solved with —
+    /// convection and diffusion included, not only the transient term — or the
+    /// correction `u = u* - (ΔV / a_p) ∂p'/∂x` does not undo the momentum
+    /// imbalance it is meant to.
+    ///
+    /// It previously returned `rho dx dy / dt`, which is only the transient
+    /// contribution `a_p0`. On a Re = 100 cavity the convective and diffusive
+    /// terms are of the same order as `a_p0`, so the correction was roughly
+    /// twice as large as it should have been.
     fn compute_u_momentum_center_coefficient(
         &self,
-        _flow_field: &FlowField,
-        _i: usize,
-        _j: usize,
+        flow_field: &FlowField,
+        i: usize,
+        j: usize,
         dx: f64,
         dy: f64,
         rho: f64,
     ) -> CfdResult<f64> {
-        // Simplified calculation for velocity correction
-        // This would be the diagonal coefficient from momentum discretization
-        Ok(rho * dx * dy / self.parameters.time_step)
+        let coeffs = self.compute_u_momentum_coefficients(
+            flow_field,
+            i,
+            j,
+            self.parameters.time_step,
+            rho,
+            self.config.viscosity,
+            dx,
+            dy,
+        )?;
+        Ok(coeffs.center)
     }
 
-    /// Compute center coefficient for v-momentum equation
+    /// Diagonal coefficient of the v-momentum equation. See
+    /// [`Self::compute_u_momentum_center_coefficient`].
     fn compute_v_momentum_center_coefficient(
         &self,
-        _flow_field: &FlowField,
-        _i: usize,
-        _j: usize,
+        flow_field: &FlowField,
+        i: usize,
+        j: usize,
         dx: f64,
         dy: f64,
         rho: f64,
     ) -> CfdResult<f64> {
-        Ok(rho * dx * dy / self.parameters.time_step)
+        let coeffs = self.compute_v_momentum_coefficients(
+            flow_field,
+            i,
+            j,
+            self.parameters.time_step,
+            rho,
+            self.config.viscosity,
+            dx,
+            dy,
+        )?;
+        Ok(coeffs.center)
     }
 }