rtx-cfd: apply boundary conditions to u* before using its divergence

Closes the relaxation-factor dependence. Converged solutions are now
identical for velocity relaxation 0.3, 0.5, 0.7 and 0.9 -- bit for bit --
where they previously spread 18%.

The cause was ordering, not formulation. Boundary conditions were applied
only at the end of the iteration, so `copy_to_starred` snapshotted a
predicted field whose boundary faces held whatever the momentum sweep had
written there: values the wall overwrote with zero moments later. The
divergence of that field is the pressure equation's source, so those
un-constrained faces entered it as a spurious mass source, concentrated at
the two lid corners where the moving lid meets a stationary wall. The
swept value scales with the relaxation factor, so the spurious source did
too -- and so did the answer.

Diagnosis is worth recording because the symptom pointed away from the
cause. At the stalled state the interior momentum equations were satisfied
to machine precision at every relaxation factor: a fresh Gauss-Seidel
sweep moved the interior by 1e-15, the pressure correction was 1e-14, and
the momentum residual was 3.6e-16. The entire residual floor lived in the
*mass* term, and only that term varied with alpha -- 3.7e-4 at 0.3 against
1.6e-4 at 0.9. Each relaxation factor was converging honestly, to the
solution of a slightly different problem.

The correction also moves the cavity substantially closer to the reference,
because the spurious corner source had been suppressing the recirculation:

  grid    before    after     Ghia (1982)
  17^2    -0.068    -0.123    -0.2109
  33^2    -0.109    -0.154
  65^2    -0.142    -0.174
  97^2    -0.157    -0.182

Richardson extrapolation on the two finest grids now gives about -0.199
against Ghia's -0.2109, within 6%, with the remaining gap consistent with
first-order upwind's numerical viscosity. The residual floor falls roughly
linearly with mesh size (1.6e-3, 5.2e-4, 1.6e-4, 7.8e-5), which is the
signature of the corner singularity rather than of an unconverged solve --
the same one Botella & Peyret (1998) subtract analytically.

Two further fixes fell out of it:

  - The solver returned NaN rather than reporting divergence. Asked for an
    8x8 cavity at a Reynolds number of a million it now stops, says it did
    not converge, and reports the last finite residual, instead of handing
    back a field of NaN that poisons everything downstream. Previously the
    false transient's large diagonal damped that case into crawling rather
    than diverging, which hid it.

  - `apply_boundary_condition` used `start_index` where it meant
    `end_index` for the bottom wall. The other three arms are correct; with
    both indices unset the default masked it.

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

@@ -220,6 +220,22 @@ impl SimpleSolver {
         // transient term, the under-relaxation and the residual all read.
         self.momentum_prediction_step(flow_field, dt).await?;
 
+        // The predicted field must satisfy the velocity boundary conditions
+        // before its divergence is used as the pressure source.
+        //
+        // Boundary conditions were previously applied only at the end of the
+        // iteration, so `u*` carried whatever the momentum sweep happened to
+        // write on the boundary faces — values the wall then overwrote with
+        // zero anyway. Their divergence entered the pressure equation as a
+        // spurious mass source concentrated at the two lid corners, where the
+        // moving lid meets a stationary wall. Because the swept value scales
+        // with the relaxation factor, so did the spurious source, and so did
+        // the converged solution: the interior momentum equations were
+        // satisfied to machine precision at every relaxation factor, but each
+        // one satisfied them around a different corner condition.
+        flow_field.apply_boundary_conditions(boundary_conditions)?;
+        flow_field.copy_to_starred();
+
         // Step 2: Solve pressure correction equation
         let mass_residual = self.pressure_correction_step(flow_field).await?;
 
@@ -1029,6 +1045,36 @@ impl IncompressibleSolver for SimpleSolver {
 
             let total_residual =
                 (mass_residual * mass_residual + momentum_residual * momentum_residual).sqrt();
+
+            // Stop on divergence rather than returning NaN.
+            //
+            // A solver asked for something it cannot do — here an 8x8 cavity
+            // at a Reynolds number of a million — should say it did not
+            // converge, not hand back a field of NaN that silently poisons
+            // everything downstream. Reports the last finite residual so the
+            // caller can see how far it got before it blew up.
+            if !total_residual.is_finite() {
+                let solve_time = start_time.elapsed();
+                let last_finite = residual_history
+                    .iter()
+                    .rev()
+                    .copied()
+                    .find(|r: &f64| r.is_finite())
+                    .unwrap_or(f64::MAX);
+                return Ok(SimpleResult {
+                    solver_result: SolverResult {
+                        converged: false,
+                        iterations: iteration + 1,
+                        final_residual: last_finite,
+                        residual_history,
+                        solve_time,
+                    },
+                    pressure_iterations,
+                    mass_residual: last_finite,
+                    momentum_residual: last_finite,
+                });
+            }
+
             residual_history.push(total_residual);
 
             if total_residual < self.parameters.tolerance {