rtx-cfd: manufactured solution finds the diffusion conductances were 1/h too

large

Applies MMS to the SIMPLE solver. It found a major discretisation error on
the first run, which is the point of the method.

The diffusion conductances read `mu / dx` and `mu / dy`. Finite volume
requires `Gamma * A / delta` — the face area over the distance between the
nodes it separates — so they should be `mu * dy / dx` and `mu * dx / dy`.
The face area was missing entirely, making viscosity too large by a factor
of `1/h`: sixty-five times on a 65x65 mesh. Every other term in the
equation was already a force (`dp * dy` for pressure, `rho u dy` for the
convective flux), so the mismatch was confined to diffusion.

The consequence was that the solver ran at an effective Reynolds number
far below the one requested. Before the fix the manufactured-solution
error did not reduce under refinement at all — observed order about -0.05,
because the spurious viscosity grows with the mesh. After it, the error
falls monotonically.

This also explains an apparent regression that is really a correction.
The cavity vortex position moved from y = 0.484 to y = 0.391 against
Ghia's 0.4531, which reads as worse agreement. It is not: a strongly
over-diffusive cavity approaches Stokes flow, whose vortex sits near
mid-height, so the old number was closer to the reference than the scheme
deserved. Correcting the viscosity exposed the discretisation's own error.
The test now states that disagreement plainly rather than asserting a band
around the reference.

What MMS reports now, and it is not yet good enough:

    n = 16   L2 velocity error = 2.586104e-1   order    -
    n = 32   L2 velocity error = 1.797373e-1   order 0.52
    n = 64   L2 velocity error = 1.277188e-1   order 0.49

First-order upwind should give 1. It gives about 0.5, and the u component
is markedly further from exact than v on the same mesh. Both say there is
at least one more defect in the discretisation or its boundary treatment,
and the asymmetry between the two momentum equations is the clue. The test
asserts only monotone error reduction — what is established — and records
the shortfall, because asserting a rate the solver does not achieve would
either redden the suite or invite someone to weaken it later.

This changes the plan: raising the observed order to 1 is now a
precondition for the second-order convection work rather than a
consequence of it. There is no value in adding a higher-order scheme to a
discretisation that has not demonstrated first order.

Supporting changes:

  - `SimpleSolver::set_momentum_source` applies a volumetric body force,
    which is what lets a manufactured solution be imposed at all.
  - Divergence is now detected by growth, not only by NaN. The 8x8 case at
    Reynolds 10^6 reached 1e149 before anything caught it, because
    `is_finite` stays true right up until it does not.
  - `test_simple_solver_workflow` specified water properties on a unit
    domain, which is Reynolds 10^6 on ten cells: no steady laminar
    solution exists and the solver diverges on it, correctly. It passed
    only while the excess diffusion stabilised it. Now set to Reynolds 100.

561 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

@@ -151,6 +151,15 @@ pub struct SimpleSolver {
     workspace: LinearAlgebraWorkspace,
     /// Turbulence model (optional)
     turbulence_model: Option<KEpsilonModel>,
+    /// Optional volumetric momentum source `f(x, y) -> (f_x, f_y)`, per unit
+    /// volume.
+    ///
+    /// Exists so a manufactured solution can be imposed: given any velocity
+    /// and pressure field, the residual of the momentum equations *is* the
+    /// body force that makes that field exact, and applying it turns the
+    /// solver into something whose exact answer is known in closed form.
+    #[allow(clippy::type_complexity)]
+    momentum_source: Option<Box<dyn Fn(f64, f64) -> (f64, f64) + Send + Sync>>,
 }
 
 /// Workspace for linear algebra operations
@@ -183,6 +192,10 @@ struct MomentumCoefficients {
 }
 
 impl SimpleSolver {
+    /// How far the residual may rise above its best value before the solve is
+    /// declared divergent.
+    const DIVERGENCE_GROWTH: f64 = 1e6;
+
     /// Create new SIMPLE solver
     pub fn new(config: CfdConfig, parameters: SimpleParameters) -> CfdResult<Self> {
         config.validate()?;
@@ -205,9 +218,37 @@ impl SimpleSolver {
                 momentum_coefficients: None,
             },
             turbulence_model,
+            momentum_source: None,
         })
     }
 
+    /// Set a volumetric momentum source. See [`Self::momentum_source`].
+    pub fn set_momentum_source<F>(&mut self, source: F)
+    where
+        F: Fn(f64, f64) -> (f64, f64) + Send + Sync + 'static,
+    {
+        self.momentum_source = Some(Box::new(source));
+    }
+
+    /// Momentum source contribution for a u-face, already multiplied by the
+    /// control volume so it is a force, matching the pressure-gradient term.
+    ///
+    /// On this staggered layout u-face `i` sits at `x = i dx`, mid-height of
+    /// row `j`, i.e. `y = (j + 0.5) dy`.
+    fn u_source_term(&self, i: usize, j: usize, dx: f64, dy: f64) -> f64 {
+        self.momentum_source
+            .as_ref()
+            .map_or(0.0, |f| f(i as f64 * dx, (j as f64 + 0.5) * dy).0 * dx * dy)
+    }
+
+    /// Momentum source contribution for a v-face, at `x = (i + 0.5) dx`,
+    /// `y = j dy`.
+    fn v_source_term(&self, i: usize, j: usize, dx: f64, dy: f64) -> f64 {
+        self.momentum_source
+            .as_ref()
+            .map_or(0.0, |f| f((i as f64 + 0.5) * dx, j as f64 * dy).1 * dx * dy)
+    }
+
     /// Solve one SIMPLE iteration
     pub async fn solve_simple_iteration(
         &mut self,
@@ -736,10 +777,20 @@ impl SimpleSolver {
         let mu_eff = self.compute_effective_viscosity(flow_field, i, j, mu);
 
         // Diffusion coefficients using effective viscosity
-        let gamma_e = mu_eff / dx;
+        // Diffusion conductances: `Gamma * A / delta`, the face area over the
-        let gamma_w = mu_eff / dx;
+        // distance between the nodes it separates.
-        let gamma_n = mu_eff / dy;
+        //
-        let gamma_s = mu_eff / dy;
+        // These previously read `mu / dx` and `mu / dy`, omitting the face
+        // area entirely. On a square grid that makes them a factor `1/h` too
+        // large — 65 times too much diffusion on a 65x65 mesh — so the solver
+        // ran at an effective Reynolds number far below the one requested.
+        // Every other term is already a force: the pressure term is
+        // `dp * dy`, the convective flux is `rho u dy`, so the mismatch was
+        // confined to diffusion.
+        let gamma_e = mu_eff * dy / dx;
+        let gamma_w = mu_eff * dy / dx;
+        let gamma_n = mu_eff * dx / dy;
+        let gamma_s = mu_eff * dx / dy;
 
         // 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
@@ -827,10 +878,12 @@ impl SimpleSolver {
         let mu_eff = self.compute_effective_viscosity(flow_field, i, j, mu);
 
         // Similar to u-momentum but for v-component
-        let gamma_e = mu_eff / dx;
+        // See the note in the u-momentum routine: `Gamma * A / delta`, not
-        let gamma_w = mu_eff / dx;
+        // `Gamma / delta`.
-        let gamma_n = mu_eff / dy;
+        let gamma_e = mu_eff * dy / dx;
-        let gamma_s = mu_eff / dy;
+        let gamma_w = mu_eff * dy / dx;
+        let gamma_n = mu_eff * dx / dy;
+        let gamma_s = mu_eff * dx / dy;
 
         // 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
@@ -864,6 +917,7 @@ impl SimpleSolver {
         let ap = ap_unrelaxed / alpha;
         let source = pressure_gradient
             + time_term
+            + self.v_source_term(i, j, dx, dy)
             + (1.0 - alpha) / alpha * ap_unrelaxed * flow_field.v_old[(j, i)];
 
         Ok(MomentumEquationCoeffs {
@@ -1037,6 +1091,7 @@ impl IncompressibleSolver for SimpleSolver {
         let start_time = Instant::now();
         let mut residual_history = Vec::new();
         let pressure_iterations = Vec::new();
+        let mut best_residual = f64::INFINITY;
 
         for iteration in 0..self.parameters.max_iterations {
             let (mass_residual, momentum_residual) = self
@@ -1046,6 +1101,33 @@ impl IncompressibleSolver for SimpleSolver {
             let total_residual =
                 (mass_residual * mass_residual + momentum_residual * momentum_residual).sqrt();
 
+            // Stop on divergence rather than running on to overflow.
+            //
+            // A residual that has grown by orders of magnitude above its best
+            // value is diverging, and continuing only turns a large number
+            // into an enormous one — the 8x8 cavity at a Reynolds number of a
+            // million reached 1e149 before anything caught it, because
+            // `is_finite` stays true right up to the moment it does not.
+            if total_residual.is_finite()
+                && best_residual.is_finite()
+                && total_residual > best_residual * Self::DIVERGENCE_GROWTH
+            {
+                let solve_time = start_time.elapsed();
+                return Ok(SimpleResult {
+                    solver_result: SolverResult {
+                        converged: false,
+                        iterations: iteration + 1,
+                        final_residual: total_residual,
+                        residual_history,
+                        solve_time,
+                    },
+                    pressure_iterations,
+                    mass_residual,
+                    momentum_residual,
+                });
+            }
+            best_residual = best_residual.min(total_residual);
+
             // Stop on divergence rather than returning NaN.
             //
             // A solver asked for something it cannot do — here an 8x8 cavity

crates/specialized/rtx-cfd/tests/cfd_workflow_tests.rs

@@ -16,9 +16,16 @@ use std::time::Duration;
 #[tokio::test]
 async fn test_simple_solver_workflow() -> CfdResult<()> {
     // Create CFD configuration
+    // Reynolds number 100, not 10^6.
+    //
+    // The nominal water properties (rho = 1000, mu = 1e-3) give
+    // Re = rho U L / mu = 10^6 on a unit domain, which has no steady laminar
+    // solution and could not be resolved by ten cells if it did. The solver
+    // diverges on it, correctly. It only appeared to work while the diffusion
+    // conductances were a factor 1/h too large, which quietly stabilised it.
     let config = CfdConfig::new()
-        .with_density(1000.0) // Water density
+        .with_density(1000.0)
-        .with_viscosity(1e-3) // Water viscosity
+        .with_viscosity(10.0)
         .with_reference_velocity(1.0)
         .with_reference_length(1.0);
 

crates/specialized/rtx-cfd/tests/mms_navier_stokes.rs

@@ -0,0 +1,234 @@
+//! Code verification of the SIMPLE solver by manufactured solution.
+//!
+//! Cavity benchmarks tell you whether the answer looks like the picture in the
+//! paper. This tells you the rate at which the discretisation converges to an
+//! exact solution, which is a statement about the code rather than about the
+//! flow, and which no amount of tuning can fake.
+//!
+//! # The manufactured solution
+//!
+//! ```text
+//!     u(x, y) =  sin(pi x) cos(pi y)
+//!     v(x, y) = -cos(pi x) sin(pi y)
+//!     p(x, y) =  sin(pi x) sin(pi y)
+//! ```
+//!
+//! The velocity field is divergence-free by construction —
+//! `du/dx = pi cos cos` and `dv/dy = -pi cos cos` — which it must be, or the
+//! pressure-correction equation is being asked to solve an inconsistent
+//! problem and the measurement means nothing.
+//!
+//! The pressure is deliberately *not* the one that goes with this velocity in
+//! a force-free flow. With `p = (cos 2pi x + cos 2pi y)/4` the convective and
+//! pressure terms cancel identically at unit density, which would leave the
+//! convection discretisation untested. This choice keeps all three terms
+//! present in the source.
+//!
+//! # Grid convention
+//!
+//! On this staggered layout, with `nx` by `ny` cells over the unit square:
+//!
+//! ```text
+//!     u[(j, i)]  at  ( i dx,          (j + 0.5) dy )     i = 0..=nx
+//!     v[(j, i)]  at  ( (i + 0.5) dx,  j dy         )     j = 0..=ny
+//!     p[(j, i)]  at  ( (i + 0.5) dx,  (j + 0.5) dy )
+//! ```
+//!
+//! Cell `i` is bounded by u-faces `i` and `i + 1`, which is the convention
+//! `compute_mass_source` uses to form the divergence and therefore the one
+//! that defines the grid.
+//!
+//! # What this currently measures — and it is not yet good enough
+//!
+//! First-order upwind should give an observed order of 1. **It measures about
+//! 0.5**, and the `u` component is markedly further from the exact solution
+//! than `v` on the same mesh. Both facts say there is at least one more defect
+//! in the discretisation or its boundary treatment, and the asymmetry between
+//! the two momentum equations is the strongest clue as to where.
+//!
+//! This test therefore asserts what is established — that the error falls
+//! monotonically under refinement, which it did *not* do before the diffusion
+//! conductances were corrected — and records the shortfall rather than
+//! asserting a rate the solver does not achieve. Raising it to 1 is the
+//! precondition for the second-order convection work, not a consequence of it:
+//! there is no point adding a higher-order scheme to a discretisation that has
+//! not yet demonstrated first order.
+//!
+//! Found by this test already: the diffusion conductances omitted the face
+//! area, reading `mu / dx` where finite volume requires `mu * A / delta`, so
+//! viscosity was too large by a factor of `1/h` — 65 times on a 65x65 mesh.
+//! Before that fix the error did not reduce under refinement at all.
+
+use rtx_cfd::solvers::incompressible::{
+    BoundaryConditions, FlowField, SimpleParameters, SimpleSolver,
+};
+use rtx_cfd::{CfdConfig, CfdResult};
+use std::f64::consts::PI;
+
+const RHO: f64 = 1.0;
+const MU: f64 = 0.05;
+
+fn u_exact(x: f64, y: f64) -> f64 {
+    (PI * x).sin() * (PI * y).cos()
+}
+
+fn v_exact(x: f64, y: f64) -> f64 {
+    -(PI * x).cos() * (PI * y).sin()
+}
+
+/// Momentum source `f = rho (u.grad)u - mu lap(u) + grad p`, derived by hand.
+///
+/// The convective terms collapse neatly:
+/// `(u.grad)u = pi sin(pi x) cos(pi x) = (pi/2) sin(2 pi x)` and likewise
+/// `(u.grad)v = (pi/2) sin(2 pi y)`, because `sin^2 + cos^2` factors out.
+/// The Laplacians are `lap(u) = -2 pi^2 u` and `lap(v) = -2 pi^2 v`.
+fn source(x: f64, y: f64) -> (f64, f64) {
+    let fx = RHO * 0.5 * PI * (2.0 * PI * x).sin()
+        + 2.0 * PI * PI * MU * u_exact(x, y)
+        + PI * (PI * x).cos() * (PI * y).sin();
+    let fy = RHO * 0.5 * PI * (2.0 * PI * y).sin()
+        + 2.0 * PI * PI * MU * v_exact(x, y)
+        + PI * (PI * x).sin() * (PI * y).cos();
+    (fx, fy)
+}
+
+/// Solve the manufactured problem on an `n` by `n` grid, returning the L2
+/// error of the velocity field over the interior faces.
+async fn l2_error(n: usize) -> CfdResult<f64> {
+    let dx = 1.0 / n as f64;
+    let dy = dx;
+
+    let config = CfdConfig::new()
+        .with_density(RHO)
+        .with_viscosity(MU)
+        .with_reference_velocity(1.0)
+        .with_reference_length(1.0);
+    let params = SimpleParameters::default()
+        .with_max_iterations(40000)
+        .with_tolerance(1e-9);
+    let mut solver = SimpleSolver::new(config, params)?;
+    solver.set_momentum_source(source);
+
+    let mut field = FlowField::new(n, n, dx, dy)?;
+
+    // Impose the exact solution on the outermost layer of faces.
+    //
+    // The momentum sweeps run over interior faces only — `1..nx` for u and
+    // `1..ny` for v — so this layer is never overwritten and an empty
+    // boundary-condition set leaves it untouched for the whole solve. That is
+    // exactly the Dirichlet problem the manufactured solution defines.
+    let set_boundary = |field: &mut FlowField| {
+        for j in 0..n {
+            let y = (j as f64 + 0.5) * dy;
+            field.u[(j, 0)] = u_exact(0.0, y);
+            field.u[(j, n)] = u_exact(1.0, y);
+        }
+        for i in 0..n {
+            let x = (i as f64 + 0.5) * dx;
+            field.v[(0, i)] = v_exact(x, 0.0);
+            field.v[(n, i)] = v_exact(x, 1.0);
+        }
+        // u on the top and bottom rows, and v on the left and right columns,
+        // are also outside the swept range.
+        for i in 0..=n {
+            let x = i as f64 * dx;
+            field.u[(0, i)] = u_exact(x, 0.5 * dy);
+            field.u[(n - 1, i)] = u_exact(x, (n as f64 - 0.5) * dy);
+        }
+        for j in 0..=n {
+            let y = j as f64 * dy;
+            field.v[(j, 0)] = v_exact(0.5 * dx, y);
+            field.v[(j, n - 1)] = v_exact((n as f64 - 0.5) * dx, y);
+        }
+    };
+    set_boundary(&mut field);
+
+    let empty = BoundaryConditions::new();
+    let runtime_iterations = 40000;
+    for _ in 0..runtime_iterations {
+        let (mass, momentum) = solver
+            .solve_simple_iteration(&mut field, &empty, 0.01)
+            .await?;
+        if (mass * mass + momentum * momentum).sqrt() < 1e-9 {
+            break;
+        }
+    }
+
+    // L2 error over interior faces, weighted by cell volume.
+    let mut squared = 0.0;
+    let mut volume = 0.0;
+    for j in 1..n - 1 {
+        for i in 1..n {
+            let e = field.u[(j, i)] - u_exact(i as f64 * dx, (j as f64 + 0.5) * dy);
+            squared += e * e * dx * dy;
+            volume += dx * dy;
+        }
+    }
+    for j in 1..n {
+        for i in 1..n - 1 {
+            let e = field.v[(j, i)] - v_exact((i as f64 + 0.5) * dx, j as f64 * dy);
+            squared += e * e * dx * dy;
+            volume += dx * dy;
+        }
+    }
+
+    Ok((squared / volume).sqrt())
+}
+
+/// The velocity error must fall under refinement.
+///
+/// Deliberately weaker than the order-1 assertion this solver ought to
+/// satisfy. See the module documentation: the observed order is about 0.5, and
+/// asserting 1 here would either fail the suite or invite someone to weaken it
+/// later. Monotone reduction still has teeth — it is exactly what failed
+/// before the diffusion conductances were fixed, when the error grew with
+/// refinement.
+#[tokio::test]
+async fn observed_order_matches_the_convection_scheme() -> CfdResult<()> {
+    let resolutions = [16usize, 32, 64];
+    let mut errors: Vec<f64> = Vec::new();
+    for &n in &resolutions {
+        errors.push(l2_error(n).await?);
+    }
+    let rates: Vec<f64> = errors
+        .windows(2)
+        .map(|pair| (pair[0] / pair[1]).log2())
+        .collect();
+
+    for (i, &n) in resolutions.iter().enumerate() {
+        let rate = if i == 0 {
+            String::from("    -")
+        } else {
+            format!("{:5.2}", rates[i - 1])
+        };
+        println!(
+            "  n = {n:3}   L2 velocity error = {:.6e}   observed order = {rate}",
+            errors[i]
+        );
+    }
+
+    assert!(
+        errors.windows(2).all(|pair| pair[1] < pair[0]),
+        "the error must fall under refinement; got {errors:?}"
+    );
+
+    for (i, &rate) in rates.iter().enumerate() {
+        assert!(
+            rate > 0.3,
+            "refinement {} -> {}: observed order {rate:.3}. The discretisation \
+             has essentially stopped converging. Errors: {errors:?}",
+            resolutions[i],
+            resolutions[i + 1]
+        );
+        assert!(
+            rate < 2.3,
+            "refinement {} -> {}: observed order {rate:.3}, above what a \
+             first-order convection scheme can deliver — suspect the error \
+             measure rather than celebrating. Errors: {errors:?}",
+            resolutions[i],
+            resolutions[i + 1]
+        );
+    }
+
+    Ok(())
+}

crates/specialized/rtx-cfd/tests/simple_tests.rs

@@ -142,12 +142,27 @@ mod simple_tests {
              the solution is a shear layer, not a cavity"
         );
 
-        // The vortex *position* is not scheme-limited in the way its strength
+        // Ghia puts the centreline minimum at y = 0.4531. This solver reads
-        // is, so it can be asserted tightly: Ghia put it at y = 0.4531.
+        // about 0.39 on a 65^2 grid, and the band below is set accordingly
+        // rather than around the reference — an honest record of a real
+        // disagreement, not a claim of agreement.
+        //
+        // It used to read 0.484, which looked better. That was partly luck:
+        // the diffusion conductances omitted the face area and were a factor
+        // 1/h too large, so the solver was running at a Reynolds number far
+        // below 100. A strongly over-diffusive cavity approaches Stokes flow,
+        // whose vortex sits near mid-height, which happened to land closer to
+        // Ghia than the scheme deserved. Correcting the viscosity exposed the
+        // discretisation's own error.
+        //
+        // `tests/mms_navier_stokes.rs` measures that error directly and finds
+        // the observed order of accuracy is about 0.5 where first-order upwind
+        // should give 1. Quantitative agreement with Ghia is not expected
+        // until that is resolved.
         assert!(
-            (0.40..0.52).contains(&y_at_min),
+            (0.34..0.55).contains(&y_at_min),
-            "the primary vortex should sit near Ghia's y = 0.4531; \
+            "the primary vortex is at y = {y_at_min:.4}, outside even the wide \
-             found the minimum at y = {y_at_min:.4}"
+             band this solver currently warrants (Ghia: 0.4531)"
         );
 
         // The *strength* is limited by first-order upwind's numerical