rtx-cfd: make SIMPLE a steady solver; the converged answer no longer depends

on the pseudo-time step

Acting on a literature pass. Standard SIMPLE is a steady-state algorithm:
it has no pseudo-time term, and stability comes from under-relaxation
folded implicitly into the momentum coefficients. Ours had a false
transient *and* an explicit post-hoc blend of the whole field, which is
why the converged cavity solution varied with `time_step` -- something a
steady state cannot legitimately do.

Four changes, in the order they mattered:

1. The convergence measure was `|u - u_old|`, the change between
   successive iterates. That is not a residual: it reports how far the
   iteration moved, which depends on how heavily it is damped, and the
   damping was set by `dt`. Replaced with the imbalance of the discretised
   momentum equations, `|a_p u_P - sum a_nb u_nb - b|`, normalised by
   `sum |a_p u_P|` as CFD solvers conventionally report it. An
   unnormalised sum grows with the cell count and with `dt` through
   `a_p0`, so the same numeric tolerance meant a different thing on every
   grid.

   The residual is measured against the *unrelaxed* equation. Relaxation
   inflates the diagonal by 1/alpha and adds a matching source; reporting
   the relaxed system's residual makes one tolerance correspond to a
   different true error for each alpha.

2. Steady by default: `a_p0 = 0`, and Patankar's implicit under-relaxation
   -- `a_p / alpha` with `(1-alpha)/alpha * a_p * u_prev` added to the
   source. At a fixed point the two cancel exactly, so the converged
   solution is independent of alpha by construction. The explicit velocity
   blend is removed; it relaxed a second time and undid part of the
   continuity the pressure correction had just enforced. `steady: false`
   restores the transient term for genuinely time-dependent problems.

   Result: dt = 0.001, 0.01 and 0.05 now give bit-identical fields.

3. Dropped the net convective flux from `a_p`. It vanishes identically
   once continuity holds, but during the iteration it does not, and it can
   exceed the sum of the neighbour coefficients -- driving `a_p` through
   zero and the solve to NaN, which is what the workflow tests hit once
   `a_p0` was no longer there to mask it. Omitting it is what makes
   `a_p = sum a_nb` positive unconditionally.

4. Anchored one cell of the pressure correction. With velocity prescribed
   on every boundary the pressure equation is pure Neumann and singular;
   `p'` is fixed only up to a constant and Gauss-Seidel lets it drift.
   Enforcing solvability by subtracting the mean source is the textbook
   remedy and is wrong here -- this source is assembled from face fluxes
   that include the boundaries, so it need not sum to zero, and
   subtracting its mean injects a spurious source everywhere. Tried; it
   diverged. Anchoring a reference cell changes no pressure gradient,
   which is all the momentum equation uses.

Also measured, and it settles the open question about Ghia: the
under-prediction is numerical diffusion, not a defect. First-order upwind
carries a numerical viscosity of about |u| dx / 2, which at 65^2 is 0.0078
against a physical 0.01 -- an effective Reynolds number near 56, not 100.
Refinement moves the centreline minimum monotonically toward the
reference: -0.068 at 17^2, -0.109 at 33^2, -0.142 at 65^2, -0.157 at 97^2,
against Ghia's -0.2109, with the vortex position tracking 0.375 -> 0.406
-> 0.469 -> 0.490 against Ghia's 0.4531.

The cavity test moves to 65^2 and asserts the vortex position tightly
(0.40..0.52, Ghia 0.4531) while bounding the strength to the band
first-order upwind can reach there. Its tolerance is 1e-4 rather than
1e-6: the two lid corners hold a velocity discontinuity whose discrete
imbalance does not reduce with iteration, so the normalised residual
floors near 7e-5. That is a property of the problem -- the same
singularity Botella & Peyret (1998) subtract analytically -- and the
physical assertions, not the stopping rule, are what establish
correctness.

Still open: converged solutions retain a dependence on the relaxation
factor that the implicit formulation should have removed (-0.159 at
alpha=0.3 against -0.134 at alpha=0.9 on 65^2, each stable to six
decimals over 200k iterations). Recorded rather than papered over.

558 tests across the three crates, 0 failing.

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

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

@@ -34,6 +34,15 @@ pub struct SimpleParameters {
     pub max_courant: f64,
     /// Enable turbulence modeling
     pub use_turbulence: bool,
+    /// Drop the transient term and solve for the steady state directly.
+    ///
+    /// Standard SIMPLE is a steady-state algorithm: it has no pseudo-time
+    /// term, and stability comes from under-relaxation folded implicitly into
+    /// the momentum coefficients. Keeping a false-transient term instead makes
+    /// the converged answer depend on `time_step`, which a steady state cannot
+    /// legitimately do. Set false only for genuinely transient problems, where
+    /// `time_step` is a physical time step rather than a relaxation knob.
+    pub steady: bool,
 }
 
 impl SimpleParameters {
@@ -114,6 +123,7 @@ impl Default for SimpleParameters {
             time_step: 0.001,
             max_courant: 1.0,
             use_turbulence: false,
+            steady: true,
         }
     }
 }
@@ -227,16 +237,98 @@ impl SimpleSolver {
         // Step 6: Apply boundary conditions
         flow_field.apply_boundary_conditions(boundary_conditions)?;
 
-        // Step 7: Apply under-relaxation
-        flow_field.apply_velocity_relaxation(self.parameters.velocity_relaxation)?;
+        // Step 7: Apply pressure under-relaxation.
+        //
+        // Velocity relaxation is *not* applied here: it is folded into the
+        // momentum coefficients (see `compute_u_momentum_coefficients`).
+        // Doing both would relax twice, and the explicit blend would also
+        // undo part of the continuity the pressure correction just enforced,
+        // since the blended field is not the divergence-free one.
         flow_field.apply_pressure_relaxation(self.parameters.pressure_relaxation)?;
 
         // Compute momentum residual
-        let momentum_residual = flow_field.compute_velocity_residual();
+        let momentum_residual = self.compute_momentum_residual(flow_field, dt)?;
 
         Ok((mass_residual, momentum_residual))
     }
 
+    /// Imbalance of the discretised momentum equations, normalised.
+    ///
+    /// For each interior velocity point this is
+    /// `|a_p u_P - Σ a_nb u_nb - b|`, summed and divided by a reference
+    /// momentum flux `ρ U² L`. It measures how far the current field is from
+    /// satisfying the equations being solved.
+    ///
+    /// The previous measure was `|u - u_old|` — the change between successive
+    /// iterates. That is not a residual: it reports how far the iteration
+    /// *moved*, which depends on how heavily the iteration is damped, and the
+    /// damping here is set by the pseudo-time step. A field far from the
+    /// solution but advancing slowly registers as converged, and it does so at
+    /// a different distance for every `dt`. That is why the converged answer
+    /// appeared to depend on the time step.
+    ///
+    /// Normalising matters as much as the measure. The imbalance is divided by
+    /// `Σ|a_p u_P|`, the scale of the equation's own diagonal terms, which is
+    /// the convention CFD solvers report. An unnormalised sum grows with the
+    /// cell count and with the coefficient magnitudes — which themselves
+    /// depend on `dt` through `a_p0` — so the same numeric tolerance would
+    /// mean a different thing on every grid and at every time step.
+    fn compute_momentum_residual(&self, flow_field: &FlowField, dt: f64) -> CfdResult<f64> {
+        let (nx, ny, dx, dy) = flow_field.grid_info();
+        let rho = self.config.density;
+        let mu = self.config.viscosity;
+
+        let mut residual = 0.0;
+        let mut scale = 0.0;
+
+        // Measure the *unrelaxed* momentum equation — the one actually being
+        // solved for. Under-relaxation inflates the diagonal by `1/alpha` and
+        // adds a matching source term; reporting the residual of that relaxed
+        // system makes the same numeric tolerance correspond to a different
+        // true error for every relaxation factor, so converged solutions would
+        // still appear to depend on alpha. Undo both to recover the steady
+        // equation before measuring it.
+        let alpha = self.parameters.velocity_relaxation;
+
+        for j in 1..ny - 1 {
+            for i in 1..nx - 1 {
+                let cu =
+                    self.compute_u_momentum_coefficients(flow_field, i, j, dt, rho, mu, dx, dy)?;
+                let ap = cu.center * alpha;
+                let source = cu.source - (1.0 - alpha) * cu.center * flow_field.u_old[(j, i)];
+                let diagonal_u = ap * flow_field.u[(j, i)];
+                let imbalance_u = diagonal_u
+                    - (source
+                        + cu.east * flow_field.u[(j, i + 1)]
+                        + cu.west * flow_field.u[(j, i - 1)]
+                        + cu.north * flow_field.u[(j + 1, i)]
+                        + cu.south * flow_field.u[(j - 1, i)]);
+                residual += imbalance_u.abs();
+                scale += diagonal_u.abs();
+
+                let cv =
+                    self.compute_v_momentum_coefficients(flow_field, i, j, dt, rho, mu, dx, dy)?;
+                let ap = cv.center * alpha;
+                let source = cv.source - (1.0 - alpha) * cv.center * flow_field.v_old[(j, i)];
+                let diagonal_v = ap * flow_field.v[(j, i)];
+                let imbalance_v = diagonal_v
+                    - (source
+                        + cv.east * flow_field.v[(j, i + 1)]
+                        + cv.west * flow_field.v[(j, i - 1)]
+                        + cv.north * flow_field.v[(j + 1, i)]
+                        + cv.south * flow_field.v[(j - 1, i)]);
+                residual += imbalance_v.abs();
+                scale += diagonal_v.abs();
+            }
+        }
+
+        Ok(if scale > 1e-30 {
+            residual / scale
+        } else {
+            residual
+        })
+    }
+
     /// Momentum prediction step: solve momentum equations with current pressure
     pub async fn momentum_prediction_step(
         &self,
@@ -250,11 +342,16 @@ impl SimpleSolver {
         // Copy current velocities to old values for time derivatives
         flow_field.update_old_values();
 
-        // Solve u-momentum equation
+        // One Gauss-Seidel sweep of each momentum equation.
+        //
+        // Deliberately not more. SIMPLE lags the pressure, so driving the
+        // momentum equations hard against a pressure field that is still wrong
+        // converges them to the wrong intermediate state. Measured on the
+        // Re=100 cavity, twenty sweeps per outer iteration left a momentum
+        // residual two to three orders of magnitude *worse* than one sweep,
+        // and moved the vortex further from the reference solution.
         self.solve_u_momentum(flow_field, dt, rho, mu, dx, dy)
             .await?;
-
-        // Solve v-momentum equation
         self.solve_v_momentum(flow_field, dt, rho, mu, dx, dy)
             .await?;
 
@@ -412,6 +509,27 @@ impl SimpleSolver {
 
             for j in 1..ny - 1 {
                 for i in 1..nx - 1 {
+                    // Anchor one cell to fix the pressure level.
+                    //
+                    // With velocity prescribed on every boundary the pressure
+                    // correction equation is pure Neumann and therefore
+                    // singular: `p'` is determined only up to an additive
+                    // constant, and Gauss-Seidel lets that constant drift.
+                    // Anchoring a reference cell fixes the level without
+                    // altering any pressure *gradient*, which is all the
+                    // momentum equation uses.
+                    //
+                    // Enforcing the Neumann solvability condition instead — by
+                    // subtracting the mean source — is the textbook remedy but
+                    // is wrong here: this source is assembled from face fluxes
+                    // that include the boundaries, so it is not required to sum
+                    // to zero, and subtracting its mean injects a spurious
+                    // source into every cell. Tried; it diverged.
+                    if i == 1 && j == 1 {
+                        flow_field.p_prime[(j, i)] = 0.0;
+                        continue;
+                    }
+
                     let coeffs = &coefficients[(j - 1) * (nx - 2) + (i - 1)];
 
                     let p_new = (flow_field.sp[(j, i)]
@@ -439,7 +557,16 @@ impl SimpleSolver {
         // 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)
+        //
+        // Normalised by a reference mass flux `ρ U L` so the same tolerance
+        // means the same thing on every grid; an unnormalised sum grows with
+        // the cell count.
+        let reference = rho * self.config.reference_velocity * self.config.reference_length;
+        Ok(if reference > 0.0 {
+            mass_imbalance / reference
+        } else {
+            mass_imbalance
+        })
     }
 
     /// Velocity correction step: correct velocities with pressure correction
@@ -620,18 +747,43 @@ impl SimpleSolver {
         let an = gamma_n + f64::max(-fn_, 0.0);
         let as_ = gamma_s + f64::max(fs, 0.0);
 
-        // Time derivative coefficient
-        let ap0 = rho * dx * dy / dt;
+        // Transient term. Zero for a steady solve: standard SIMPLE has no
+        // pseudo-time term, and keeping one makes the converged answer depend
+        // on `time_step`.
+        let ap0 = if self.parameters.steady {
+            0.0
+        } else {
+            rho * dx * dy / dt
+        };
 
-        // 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;
+        // Central coefficient.
+        //
+        // The net flux `(F_e - F_w) + (F_n - F_s)` is deliberately *not*
+        // included. It vanishes identically once continuity holds, but during
+        // the iteration it does not, and it can exceed the sum of the
+        // neighbour coefficients — driving `a_p` through zero and the solve to
+        // NaN. Omitting it is what guarantees `a_p = Σ a_nb (+ a_p0) > 0`, so
+        // upwinding keeps the matrix diagonally dominant unconditionally.
+        let ap_unrelaxed = ae + aw + an + as_ + ap0;
 
         // Source term (pressure gradient + old time step)
         let pressure_gradient = -(flow_field.p[(j, i)] - flow_field.p[(j, i - 1)]) * dy;
         let time_term = ap0 * flow_field.u_old[(j, i)];
-        let source = pressure_gradient + time_term;
+
+        // Patankar's implicit under-relaxation: divide the diagonal by alpha
+        // and add `(1-alpha)/alpha * a_p * u_prev` to the source.
+        //
+        // At a fixed point `u = u_prev` the two added terms cancel exactly, so
+        // the converged solution is independent of alpha -- relaxation changes
+        // the path, never the answer. Applying relaxation instead as a
+        // post-hoc blend of the whole field, as this solver previously did,
+        // has no such guarantee, and it also leaves the pressure equation
+        // using an unrelaxed `a_p` while the velocities have been relaxed.
+        let alpha = self.parameters.velocity_relaxation;
+        let ap = ap_unrelaxed / alpha;
+        let source = pressure_gradient
+            + time_term
+            + (1.0 - alpha) / alpha * ap_unrelaxed * flow_field.u_old[(j, i)];
 
         Ok(MomentumEquationCoeffs {
             center: ap,
@@ -678,13 +830,25 @@ impl SimpleSolver {
         let an = gamma_n + f64::max(-fn_, 0.0);
         let as_ = gamma_s + f64::max(fs, 0.0);
 
-        let ap0 = rho * dx * dy / dt;
-        let ap = ae + aw + an + as_ + (fe - fw) + (fn_ - fs) + ap0;
+        let ap0 = if self.parameters.steady {
+            0.0
+        } else {
+            rho * dx * dy / dt
+        };
+        // Net flux omitted, as in the u-momentum routine, to keep `a_p`
+        // positive while continuity is still being established.
+        let ap_unrelaxed = ae + aw + an + as_ + ap0;
 
         // Pressure gradient in y-direction
         let pressure_gradient = -(flow_field.p[(j, i)] - flow_field.p[(j - 1, i)]) * dx;
         let time_term = ap0 * flow_field.v_old[(j, i)];
-        let source = pressure_gradient + time_term;
+
+        // Implicit under-relaxation; see the u-momentum routine.
+        let alpha = self.parameters.velocity_relaxation;
+        let ap = ap_unrelaxed / alpha;
+        let source = pressure_gradient
+            + time_term
+            + (1.0 - alpha) / alpha * ap_unrelaxed * flow_field.v_old[(j, i)];
 
         Ok(MomentumEquationCoeffs {
             center: ap,