rtx-cfd: TVD convection in the embedded predictor — the wake sheds — and Turek–Hron CFD2/CFD3

First-order upwind's numerical viscosity |u| h / 2 is ~10x the physical
viscosity on the Turek–Hron grids: the effective Reynolds number lands
near 20 and CFD3 (Re 200) produced NO vortex shedding at all — one lift
zero-crossing in three seconds at h = 10 mm. The physics, not a bug.

EmbeddedParameters gains `convection_scheme` (default Upwind, bit-
identical — the no-body degeneracy test still reads 0.0): the TVD branch
adds SIMPLE's limited face corrections (van Albada / van Leer,
`face_correction` now pub(crate)) directly in the explicit predictor —
no deferred iteration needed in an explicit step. Domain-side faces and
faces whose far-upwind node is outside fall back to upwind exactly as in
SIMPLE; near the body the stencil reads ghost values, which encode the
wall. Verified: the embedded-circle MMS error drops 10–16x below upwind
(8.16e-4 vs 8.49e-3 at n = 32) at observed order 1.56 (SIMPLE's TVD
measured 1.59–1.84).

tests/turek_hron_cfd23.rs — CFD2 (Re 100, steady) and CFD3 (Re 200,
periodic), both with the benchmark's inflow ramp, both measured as time
statistics over a window (never a snapshot), surface route primary and
the control volume printed as the diagnostic (its central-difference
evaluation truncation grows with the convective flux: the routes agree
to 0.6% at Re 20 and differ 15–25% at Re 100–200 on these grids).

Measured across h = 10 / 6.6 / 5 mm:
- CFD3 shedding frequency 4.2746 / 4.3400 / 4.3939 Hz vs the reference
  4.3956 — converging −2.8% -> −1.3% -> −0.04%;
- CFD3 lift mean −184 / +160 / −2.6 vs −11.9 — lands on the reference;
  lift amplitude ±438 / ±556 / ±557 vs ±437.8 — +27% at the finer grids,
  unconverged (the flag is 2/3/4 cells thick);
- CFD2 control-volume drag 152.4 / 143.3 / 139.4 vs 136.700 — +2.0% at
  5 mm; CFD2 surface drag sits ~−10% (the boundary layer is ~one cell);
  CFD2 lift −3.4 / +30.2 / +8.4 vs 10.53.
Suite defaults run CFD2 at ny = 62 and CFD3 at ny = 41 (cost); the
asserted bands are the measured ones (frequency 10%, mean drag 15%,
amplitude 35%), not accuracy claims; RTX_CFD2_NY / RTX_CFD3_NY run the
studies.

Also recorded: the CFD1 refinement study extended to h = 3.3 mm
(RTX_CFD1_NY): control-volume drag 14.8996 (+4.25%), apparent order
~0.70 sustained over four grids, control-volume lift 1.1332 vs 1.11905
(+1.3%).

rtx-cfd 318 -> 321 green (full suite 321 passed / 0 failed).

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

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

@@ -39,6 +39,7 @@
 use super::ale::{AleBoundaries, SideBoundary};
 use super::embedded_body::{EmbeddedBody, EmbeddedMask, FaceKind};
 use super::poisson::{MultigridParameters, PoissonProblem, PoissonSolverKind, solve_multigrid_pcg};
+use super::simple::ConvectionScheme;
 use super::{FlowField, SolverResult};
 use crate::{CfdConfig, CfdError, CfdResult};
 
@@ -60,6 +61,18 @@ pub struct EmbeddedParameters {
     /// [`PoissonSolverKind::Sor`]). Both solve the same system to the same
     /// true-residual stop; multigrid's cost is mesh-independent.
     pub poisson_solver: PoissonSolverKind,
+    /// Convective face values in the explicit predictor (default
+    /// [`ConvectionScheme::Upwind`], which is bit-identical to the fixed-grid
+    /// PISO). The TVD schemes add SIMPLE's limited correction to each
+    /// interior face — with an explicit predictor no deferred iteration is
+    /// needed, the limited flux is just used directly. First-order upwind's
+    /// numerical viscosity `|u| h / 2` exceeds the physical viscosity ten
+    /// times over on the Turek–Hron CFD3 grids and suppressed the vortex
+    /// shedding entirely; the limited scheme restores it. Faces whose
+    /// far-upwind node lies outside the domain, and domain-side faces, fall
+    /// back to pure upwind exactly as in SIMPLE; near the body the stencil
+    /// reads ghost values, which encode the wall.
+    pub convection_scheme: ConvectionScheme,
 }
 
 impl Default for EmbeddedParameters {
@@ -69,6 +82,7 @@ impl Default for EmbeddedParameters {
             tolerance: 1e-6,
             boundaries: AleBoundaries::default(),
             poisson_solver: PoissonSolverKind::Sor,
+            convection_scheme: ConvectionScheme::Upwind,
         }
     }
 }
@@ -315,6 +329,46 @@ impl EmbeddedPisoSolver {
                         })
                     / dy;
 
+                // Limited (TVD) corrections to the four convective face
+                // values; exactly zero-cost on the default upwind scheme.
+                let scheme = self.parameters.convection_scheme;
+                let mut conv_x = conv_x;
+                let mut conv_y = conv_y;
+                if scheme != ConvectionScheme::Upwind {
+                    let delta_e = if ue_face >= 0.0 {
+                        scheme.face_correction(Some(uo[(j, i - 1)]), uo[(j, i)], uo[(j, i + 1)])
+                    } else {
+                        let far = (i + 2 <= nx).then(|| uo[(j, i + 2)]);
+                        scheme.face_correction(far, uo[(j, i + 1)], uo[(j, i)])
+                    };
+                    let delta_w = if uw_face >= 0.0 {
+                        let far = (i >= 2).then(|| uo[(j, i - 2)]);
+                        scheme.face_correction(far, uo[(j, i - 1)], uo[(j, i)])
+                    } else {
+                        scheme.face_correction(Some(uo[(j, i + 1)]), uo[(j, i)], uo[(j, i - 1)])
+                    };
+                    let delta_n = if north_is_wall {
+                        0.0
+                    } else if vn_face >= 0.0 {
+                        let far = (j >= 1).then(|| uo[(j - 1, i)]);
+                        scheme.face_correction(far, uo[(j, i)], uo[(j + 1, i)])
+                    } else {
+                        let far = (j + 2 < ny).then(|| uo[(j + 2, i)]);
+                        scheme.face_correction(far, uo[(j + 1, i)], uo[(j, i)])
+                    };
+                    let delta_s = if south_is_wall {
+                        0.0
+                    } else if vs_face >= 0.0 {
+                        let far = (j >= 2).then(|| uo[(j - 2, i)]);
+                        scheme.face_correction(far, uo[(j - 1, i)], uo[(j, i)])
+                    } else {
+                        let far = (j + 1 < ny).then(|| uo[(j + 1, i)]);
+                        scheme.face_correction(far, uo[(j, i)], uo[(j - 1, i)])
+                    };
+                    conv_x += (ue_face * delta_e - uw_face * delta_w) / dx;
+                    conv_y += (vn_face * delta_n - vs_face * delta_s) / dy;
+                }
+
                 let diff_x = nu * (uo[(j, i + 1)] - 2.0 * u_p + uo[(j, i - 1)]) / (dx * dx);
 
                 // Wall-adjacent diffusive fluxes act over half a cell on a
@@ -397,6 +451,44 @@ impl EmbeddedPisoSolver {
                         })
                     / dx;
 
+                let scheme = self.parameters.convection_scheme;
+                let mut conv_x = conv_x;
+                let mut conv_y = conv_y;
+                if scheme != ConvectionScheme::Upwind {
+                    let delta_n = if vn_face >= 0.0 {
+                        scheme.face_correction(Some(vo[(j - 1, i)]), vo[(j, i)], vo[(j + 1, i)])
+                    } else {
+                        let far = (j + 2 <= ny).then(|| vo[(j + 2, i)]);
+                        scheme.face_correction(far, vo[(j + 1, i)], vo[(j, i)])
+                    };
+                    let delta_s = if vs_face >= 0.0 {
+                        let far = (j >= 2).then(|| vo[(j - 2, i)]);
+                        scheme.face_correction(far, vo[(j - 1, i)], vo[(j, i)])
+                    } else {
+                        scheme.face_correction(Some(vo[(j + 1, i)]), vo[(j, i)], vo[(j - 1, i)])
+                    };
+                    let delta_e = if east_is_wall {
+                        0.0
+                    } else if ue_face >= 0.0 {
+                        let far = (i >= 1).then(|| vo[(j, i - 1)]);
+                        scheme.face_correction(far, vo[(j, i)], vo[(j, i + 1)])
+                    } else {
+                        let far = (i + 2 < nx).then(|| vo[(j, i + 2)]);
+                        scheme.face_correction(far, vo[(j, i + 1)], vo[(j, i)])
+                    };
+                    let delta_w = if west_is_wall {
+                        0.0
+                    } else if uw_face >= 0.0 {
+                        let far = (i >= 2).then(|| vo[(j, i - 2)]);
+                        scheme.face_correction(far, vo[(j, i - 1)], vo[(j, i)])
+                    } else {
+                        let far = (i + 1 < nx).then(|| vo[(j, i + 1)]);
+                        scheme.face_correction(far, vo[(j, i)], vo[(j, i - 1)])
+                    };
+                    conv_y += (vn_face * delta_n - vs_face * delta_s) / dy;
+                    conv_x += (ue_face * delta_e - uw_face * delta_w) / dx;
+                }
+
                 let diff_y = nu * (vo[(j + 1, i)] - 2.0 * v_p + vo[(j - 1, i)]) / (dy * dy);
 
                 let flux_east = if east_is_wall {

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

@@ -62,8 +62,14 @@ impl ConvectionScheme {
     /// 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 {
+    /// the face falls back to pure upwind. `pub(crate)` so the embedded
+    /// solver's explicit predictor can use the same limited fluxes.
+    pub(crate) fn face_correction(
+        self,
+        far_upwind: Option<f64>,
+        upwind: f64,
+        downwind: f64,
+    ) -> f64 {
         let Some(far) = far_upwind else {
             return 0.0;
         };