rtx-cfd: second-order convection by deferred-correction TVD; MMS order 1.84, cavity closes on Ghia

First-order upwind's O(h) numerical viscosity was the measured limit on the
whole discretisation: MMS order ~0.9 at Re = 20 against 2.05 in the Stokes
limit. This adds a ConvectionScheme parameter to SimPLE — Upwind (default,
behaviour unchanged), TvdVanAlbada, TvdVanLeer — implemented by deferred
correction: the upwind operator stays implicit, so a_p = sum(a_nb) and
diagonal dominance survive unconditionally, and the limited
high-order-minus-upwind flux difference enters the source explicitly at the
current iterate. At a fixed point the two agree, so the converged answer is
the TVD discretisation. Faces whose far-upwind node lies outside the domain
fall back to pure upwind; wall faces pass no mass, so no correction enters.

Measured by the manufactured solution (van Albada, 16 -> 32 -> 64):

    L2 velocity   1.325e-3   4.406e-4   1.232e-4    orders 1.59, 1.84
    (upwind)      3.516e-2   1.954e-2   1.038e-2    orders 0.85, 0.91

The error is 27x to 84x below upwind's at equal resolution, the order climbs
toward 2 (the shortfall is limiter clipping plus the boundary fallback, both
of which shrink with h), the pressure error falls at the same rate, and
continuity still holds to solver tolerance in every cell.

On the Re = 100 lid-driven cavity at 65^2 the centreline minimum moves from
-0.1932 (upwind) to -0.2036 against Ghia's -0.2109 — 59% of the remaining
gap closed at equal resolution, converged in 790 iterations — and the vortex
position moves from 0.5000 to 0.4844 toward Ghia's 0.4531. Both new cavity
bounds exclude the upwind values, so falling back to first order fails them.

284 tests, 0 failing.

Co-Authored-By: Claude Fable 5 <[email protected]>

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

@@ -32,7 +32,7 @@ pub use flow_field::FlowField;
 pub use piso::{PisoParameters, PisoResult, PisoSolver};
 #[cfg(feature = "cuda")]
 pub use piso_gpu::PisoGpuSolver;
-pub use simple::{SimpleParameters, SimpleResult, SimpleSolver};
+pub use simple::{ConvectionScheme, SimpleParameters, SimpleResult, SimpleSolver};
 #[cfg(feature = "cuda")]
 pub use simple_gpu::SimpleGpuSolver;
 

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

@@ -17,6 +17,65 @@ use async_trait::async_trait;
 use nalgebra::{DMatrix, DVector, Vector3};
 use std::time::Instant;
 
+/// Discretisation of the convective term in the momentum equations.
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub enum ConvectionScheme {
+    /// First-order upwind. Unconditionally bounded, but carries a numerical
+    /// viscosity of about `|u| dx / 2`, which caps the observed order of the
+    /// whole discretisation at 1 whenever convection matters.
+    Upwind,
+    /// Deferred-correction TVD with the van Albada limiter
+    /// `psi(r) = (r^2 + r) / (r^2 + 1)` (0 for `r <= 0`).
+    ///
+    /// The upwind operator stays implicit, so `a_p = sum(a_nb)` and diagonal
+    /// dominance survive unconditionally; the limited high-order-minus-upwind
+    /// flux difference is added explicitly to the source, evaluated at the
+    /// current iterate. At a converged state the two agree, so the fixed
+    /// point is the TVD discretisation — relaxation changes the path, never
+    /// the answer. Faces whose far-upwind node lies outside the domain fall
+    /// back to pure upwind, the standard TVD boundary treatment.
+    TvdVanAlbada,
+    /// Deferred-correction TVD with the van Leer limiter
+    /// `psi(r) = (r + |r|) / (1 + |r|)`. Same construction as
+    /// [`ConvectionScheme::TvdVanAlbada`].
+    TvdVanLeer,
+}
+
+impl ConvectionScheme {
+    /// Flux limiter `psi(r)`. Zero recovers pure upwind, one recovers central
+    /// differencing; both TVD limiters satisfy `psi(1) = 1`, which is what
+    /// makes them second order in smooth regions.
+    fn limiter(self, r: f64) -> f64 {
+        match self {
+            Self::Upwind => 0.0,
+            Self::TvdVanAlbada => {
+                if r > 0.0 {
+                    (r * r + r) / (r * r + 1.0)
+                } else {
+                    0.0
+                }
+            }
+            Self::TvdVanLeer => (r + r.abs()) / (1.0 + r.abs()),
+        }
+    }
+
+    /// The limited correction `u_face_HO - u_face_upwind` for one face, given
+    /// the far-upwind, upwind and downwind values along the flow direction.
+    /// `None` for the far-upwind value means it lies outside the domain, and
+    /// the face falls back to pure upwind.
+    fn face_correction(self, far_upwind: Option<f64>, upwind: f64, downwind: f64) -> f64 {
+        let Some(far) = far_upwind else {
+            return 0.0;
+        };
+        let denominator = downwind - upwind;
+        if denominator.abs() < 1e-300 {
+            return 0.0;
+        }
+        let r = (upwind - far) / denominator;
+        0.5 * self.limiter(r) * denominator
+    }
+}
+
 /// Parameters for SIMPLE algorithm
 #[derive(Debug, Clone)]
 pub struct SimpleParameters {
@@ -24,6 +83,8 @@ pub struct SimpleParameters {
     pub pressure_relaxation: f64,
     /// Under-relaxation factor for velocity (typically 0.5-0.8)
     pub velocity_relaxation: f64,
+    /// Convection discretisation. Defaults to first-order upwind.
+    pub convection_scheme: ConvectionScheme,
     /// Maximum number of iterations
     pub max_iterations: usize,
     /// Convergence tolerance for residuals
@@ -66,6 +127,13 @@ impl SimpleParameters {
         self
     }
 
+    /// Set the convection scheme
+    #[must_use]
+    pub fn with_convection_scheme(mut self, scheme: ConvectionScheme) -> Self {
+        self.convection_scheme = scheme;
+        self
+    }
+
     /// Set maximum iterations
     #[must_use]
     pub fn with_max_iterations(mut self, max_iter: usize) -> Self {
@@ -118,6 +186,7 @@ impl Default for SimpleParameters {
         Self {
             pressure_relaxation: 0.3,
             velocity_relaxation: 0.7,
+            convection_scheme: ConvectionScheme::Upwind,
             max_iterations: 1000,
             tolerance: 1e-6,
             time_step: 0.001,
@@ -929,6 +998,140 @@ impl SimpleSolver {
         }
     }
 
+    /// Deferred-correction source for the u-momentum equation: the limited
+    /// high-order convective fluxes minus their upwind counterparts, moved to
+    /// the right-hand side with the sign that puts convection on the left.
+    ///
+    /// Face stencils run along the flow direction: for each face the upwind
+    /// node `C`, downwind node `D` and far-upwind node `U` define
+    /// `r = (C - U) / (D - C)`, and the correction is
+    /// `psi(r) (D - C) / 2`. A face whose far-upwind node lies outside the
+    /// domain falls back to pure upwind, and a wall face has zero mass flux,
+    /// so its correction never enters.
+    #[allow(clippy::too_many_arguments)]
+    fn u_deferred_correction(
+        &self,
+        flow_field: &FlowField,
+        i: usize,
+        j: usize,
+        nx: usize,
+        ny: usize,
+        fe: f64,
+        fw: f64,
+        fn_: f64,
+        fs: f64,
+    ) -> f64 {
+        let scheme = self.parameters.convection_scheme;
+        if scheme == ConvectionScheme::Upwind {
+            return 0.0;
+        }
+        let u = &flow_field.u;
+
+        // East face of the u control volume, between u faces `i` and `i + 1`.
+        let delta_e = if fe >= 0.0 {
+            // `i >= 1` for every unknown, so the far-upwind node exists.
+            scheme.face_correction(Some(u[(j, i - 1)]), u[(j, i)], u[(j, i + 1)])
+        } else {
+            let far = (i + 2 <= nx).then(|| u[(j, i + 2)]);
+            scheme.face_correction(far, u[(j, i + 1)], u[(j, i)])
+        };
+
+        // West face, between u faces `i - 1` and `i`.
+        let delta_w = if fw >= 0.0 {
+            let far = (i >= 2).then(|| u[(j, i - 2)]);
+            scheme.face_correction(far, u[(j, i - 1)], u[(j, i)])
+        } else {
+            scheme.face_correction(Some(u[(j, i + 1)]), u[(j, i)], u[(j, i - 1)])
+        };
+
+        // North face, between rows `j` and `j + 1`; a wall face passes no mass.
+        let delta_n = if j + 1 >= ny {
+            0.0
+        } else if fn_ >= 0.0 {
+            let far = (j >= 1).then(|| u[(j - 1, i)]);
+            scheme.face_correction(far, u[(j, i)], u[(j + 1, i)])
+        } else {
+            let far = (j + 2 < ny).then(|| u[(j + 2, i)]);
+            scheme.face_correction(far, u[(j + 1, i)], u[(j, i)])
+        };
+
+        // South face, between rows `j - 1` and `j`.
+        let delta_s = if j == 0 {
+            0.0
+        } else if fs >= 0.0 {
+            let far = (j >= 2).then(|| u[(j - 2, i)]);
+            scheme.face_correction(far, u[(j - 1, i)], u[(j, i)])
+        } else {
+            let far = (j + 1 < ny).then(|| u[(j + 1, i)]);
+            scheme.face_correction(far, u[(j, i)], u[(j - 1, i)])
+        };
+
+        -(fe * delta_e - fw * delta_w + fn_ * delta_n - fs * delta_s)
+    }
+
+    /// Deferred-correction source for the v-momentum equation; mirrors
+    /// [`Self::u_deferred_correction`] with the roles of the axes swapped.
+    #[allow(clippy::too_many_arguments)]
+    fn v_deferred_correction(
+        &self,
+        flow_field: &FlowField,
+        i: usize,
+        j: usize,
+        nx: usize,
+        ny: usize,
+        fe: f64,
+        fw: f64,
+        fn_: f64,
+        fs: f64,
+    ) -> f64 {
+        let scheme = self.parameters.convection_scheme;
+        if scheme == ConvectionScheme::Upwind {
+            return 0.0;
+        }
+        let v = &flow_field.v;
+
+        // North face of the v control volume, between v faces `j` and `j + 1`.
+        let delta_n = if fn_ >= 0.0 {
+            scheme.face_correction(Some(v[(j - 1, i)]), v[(j, i)], v[(j + 1, i)])
+        } else {
+            let far = (j + 2 <= ny).then(|| v[(j + 2, i)]);
+            scheme.face_correction(far, v[(j + 1, i)], v[(j, i)])
+        };
+
+        // South face, between v faces `j - 1` and `j`.
+        let delta_s = if fs >= 0.0 {
+            let far = (j >= 2).then(|| v[(j - 2, i)]);
+            scheme.face_correction(far, v[(j - 1, i)], v[(j, i)])
+        } else {
+            scheme.face_correction(Some(v[(j + 1, i)]), v[(j, i)], v[(j - 1, i)])
+        };
+
+        // East face, between columns `i` and `i + 1`; a wall face passes no
+        // mass.
+        let delta_e = if i + 1 >= nx {
+            0.0
+        } else if fe >= 0.0 {
+            let far = (i >= 1).then(|| v[(j, i - 1)]);
+            scheme.face_correction(far, v[(j, i)], v[(j, i + 1)])
+        } else {
+            let far = (i + 2 < nx).then(|| v[(j, i + 2)]);
+            scheme.face_correction(far, v[(j, i + 1)], v[(j, i)])
+        };
+
+        // West face, between columns `i - 1` and `i`.
+        let delta_w = if i == 0 {
+            0.0
+        } else if fw >= 0.0 {
+            let far = (i >= 2).then(|| v[(j, i - 2)]);
+            scheme.face_correction(far, v[(j, i - 1)], v[(j, i)])
+        } else {
+            let far = (i + 1 < nx).then(|| v[(j, i + 1)]);
+            scheme.face_correction(far, v[(j, i)], v[(j, i - 1)])
+        };
+
+        -(fe * delta_e - fw * delta_w + fn_ * delta_n - fs * delta_s)
+    }
+
     /// Compute coefficients for u-momentum equation
     fn compute_u_momentum_coefficients(
         &self,
@@ -941,7 +1144,7 @@ impl SimpleSolver {
         dx: f64,
         dy: f64,
     ) -> CfdResult<MomentumEquationCoeffs> {
-        let (_nx, ny, _, _) = flow_field.grid_info();
+        let (nx, ny, _, _) = flow_field.grid_info();
 
         // Compute effective viscosity (molecular + turbulent)
         let mu_eff = self.compute_effective_viscosity(flow_field, i, j, mu);
@@ -1074,6 +1277,7 @@ impl SimpleSolver {
             + time_term
             + wall_source
             + self.u_source_term(i, j, dx, dy)
+            + self.u_deferred_correction(flow_field, i, j, nx, ny, fe, fw, fn_, fs)
             + (1.0 - alpha) / alpha * ap_unrelaxed * flow_field.u_old[(j, i)];
 
         Ok(MomentumEquationCoeffs {
@@ -1098,7 +1302,7 @@ impl SimpleSolver {
         dx: f64,
         dy: f64,
     ) -> CfdResult<MomentumEquationCoeffs> {
-        let (nx, _ny, _, _) = flow_field.grid_info();
+        let (nx, ny, _, _) = flow_field.grid_info();
 
         // Compute effective viscosity (molecular + turbulent)
         let mu_eff = self.compute_effective_viscosity(flow_field, i, j, mu);
@@ -1177,6 +1381,7 @@ impl SimpleSolver {
             + time_term
             + wall_source
             + self.v_source_term(i, j, dx, dy)
+            + self.v_deferred_correction(flow_field, i, j, nx, ny, fe, fw, fn_, fs)
             + (1.0 - alpha) / alpha * ap_unrelaxed * flow_field.v_old[(j, i)];
 
         Ok(MomentumEquationCoeffs {