rtx-cfd: PatchConvection::TvdVanAlbada — van Albada deferred correction on the curvilinear predictor (downwind-side linear weight, gradient-ratio r over the face d lengths, far-upwind across the opposite face, boundary faces upwind); annulus MMS orders 2.10/1.69 at 0.24× upwind; cylinder-flag MMS orders 1.98/1.97 (1.06× upwind — diffusion-dominated, recorded); knobs RTX_OVERSET_CFD1_TVD, RTX_OVERSET_MAX_ROUNDS, RTX_CF_SCHEME

Co-Authored-By: Claude Fable 5.1 <[email protected]>
Claude-Session: https://claude.ai/code/session_01X2GmJXeQ2njUecEKiJZ1G2

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

@@ -88,6 +88,14 @@ pub enum PatchConvection {
     Upwind,
     /// No convection: the Stokes limit, for the second-order MMS gate.
     None,
+    /// Deferred-correction TVD with the van Albada limiter (the harness's
+    /// scheme): the upwind face value plus `w_up psi(r) (phi_dn − phi_up)`,
+    /// `w_up` the mesh's linear weight of the downwind side and `r` the
+    /// ratio of the two one-sided gradients (so a linear field on a
+    /// stretched row gives `r = 1` and the mesh's own linear face value).
+    /// Faces whose far-upwind cell lies outside the patch fall back to
+    /// upwind. Explicit, like the rest of the predictor's convection.
+    TvdVanAlbada,
 }
 
 /// How the across-patch diffusion is time-stepped.

crates/specialized/rtx-cfd/src/solvers/incompressible/curvilinear/predictor.rs

@@ -65,7 +65,7 @@ impl CurvilinearPisoSolver {
                 // is dropped but the mesh flux stays: the conservative
                 // update needs `−Σ sign δV_f u_f` whenever the mesh moves.
                 let fluid = match self.params.convection {
-                    PatchConvection::Upwind => field.flux[f],
+                    PatchConvection::Upwind | PatchConvection::TvdVanAlbada => field.flux[f],
                     PatchConvection::None => 0.0,
                 };
                 let out = sign * (fluid - geo.swept[f] / dt);
@@ -82,6 +82,12 @@ impl CurvilinearPisoSolver {
                                 };
                                 (field.u[up], field.v[up])
                             }
+                            PatchConvection::TvdVanAlbada => {
+                                let other = if p == c { q } else { p };
+                                let (up, dn) = if out >= 0.0 { (c, other) } else { (other, c) };
+                                let (du, dv) = self.tvd_correction(field, f, up, dn);
+                                (field.u[up] + du, field.v[up] + dv)
+                            }
                             // The Stokes limit has no upwind scheme to be
                             // consistent with; the mesh flux takes the linear
                             // face value and keeps its second order (upwinding
@@ -154,6 +160,54 @@ impl CurvilinearPisoSolver {
         (uh, vh)
     }
 
+    /// The van Albada deferred correction to the upwind face value of `f`
+    /// between the upwind cell `up` and the downwind cell `dn`: `w_up psi(r)
+    /// (phi_dn − phi_up)` for `u` and `v`, zero when the far-upwind cell
+    /// (across `up`'s opposite face) lies outside the patch.
+    fn tvd_correction(&self, field: &PatchField, f: usize, up: usize, dn: usize) -> (f64, f64) {
+        let mesh = &self.mesh;
+        let faces = mesh.cell_faces(up);
+        let Some(pos) = faces.iter().position(|&(g, _)| g == f) else {
+            return (0.0, 0.0);
+        };
+        let g = faces[pos ^ 1].0;
+        let far_face = &mesh.faces()[g];
+        let far = match (far_face.owner, far_face.neigh) {
+            (Some(a), Some(b)) => {
+                if a == up {
+                    b
+                } else {
+                    a
+                }
+            }
+            _ => return (0.0, 0.0),
+        };
+        let face = &mesh.faces()[f];
+        // Linear weight of the downwind side: `1 − w` when the owner is
+        // upwind, `w` when the neighbour is.
+        let w_up = if face.owner == Some(up) {
+            1.0 - face.w
+        } else {
+            face.w
+        };
+        let len = |d: [f64; 2]| (d[0] * d[0] + d[1] * d[1]).sqrt();
+        let (df, dg) = (len(face.d), len(far_face.d));
+        let limited = |phi: &[f64]| -> f64 {
+            let near = phi[dn] - phi[up];
+            if near.abs() < 1e-300 {
+                return 0.0;
+            }
+            let r = (phi[up] - phi[far]) / dg * df / near;
+            let psi = if r > 0.0 {
+                (r * r + r) / (r * r + 1.0)
+            } else {
+                0.0
+            };
+            w_up * psi * near
+        };
+        (limited(&field.u), limited(&field.v))
+    }
+
     /// `(I − dt ν L_n/V^{n+1}) û = rhs` along every s-line, Thomas algorithm.
     fn solve_lines(
         &self,

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

@@ -701,7 +701,7 @@ async fn taylor_green_order_is_unchanged_under_mesh_motion() -> CfdResult<()> {
         .and_then(|v| v.parse::<usize>().ok());
     for convection in [PatchConvection::Upwind, PatchConvection::None] {
         let (band, error_factor): (std::ops::Range<f64>, f64) = match convection {
-            PatchConvection::Upwind => (0.7..1.6, 2.0),
+            PatchConvection::Upwind | PatchConvection::TvdVanAlbada => (0.7..1.6, 2.0),
             PatchConvection::None => (1.8..2.4, 3.0),
         };
         let mut fixed = Vec::new();

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

@@ -81,7 +81,7 @@ async fn march(
         normal_diffusion: diffusion,
         ..CurvilinearParameters::default()
     };
-    let convecting = convection == PatchConvection::Upwind;
+    let convecting = convection != PatchConvection::None;
     let mut solver = CurvilinearPisoSolver::new(config, params, mesh)?;
     solver.set_boundary_velocity(|x, y, _t| (u_exact(x, y), v_exact(x, y)));
     solver.set_momentum_source(move |x, y, _t| source(x, y, convecting));
@@ -242,6 +242,45 @@ async fn skewed_annulus_with_upwind_is_first_order() -> CfdResult<()> {
     Ok(())
 }
 
+/// P4 step 2 gate (ii): the van Albada deferred correction on the skewed
+/// annulus beats upwind at every rung (upwind: 1.151502e-1 / 5.445770e-2 /
+/// 3.119486e-2) with orders >= 1.4.
+#[tokio::test]
+async fn skewed_annulus_with_tvd_beats_upwind() -> CfdResult<()> {
+    let upwind = [1.151502e-1, 5.445770e-2, 3.119486e-2];
+    let mut errs = Vec::new();
+    for (&ns, &u) in [32usize, 64, 128].iter().zip(&upwind) {
+        let m = march(
+            annulus_skewed([0.0, 0.0], 0.5, 1.5, ns, ns / 4, 0.3, 3.0)?,
+            PatchConvection::TvdVanAlbada,
+            NormalDiffusion::Explicit,
+            1e-6,
+        )
+        .await?;
+        println!(
+            "annulus tvd ns={ns}: L2 {:.6e} (upwind {u:.6e}, ratio {:.2}), max div {:.2e}, {} steps",
+            m.l2_velocity,
+            m.l2_velocity / u,
+            m.max_div_rel,
+            m.steps
+        );
+        assert!(m.max_div_rel < 1e-9, "divergence {:.3e}", m.max_div_rel);
+        assert!(
+            m.l2_velocity < u,
+            "TVD {:.4e} not below upwind {u:.4e}",
+            m.l2_velocity
+        );
+        errs.push(m.l2_velocity);
+    }
+    let o = orders(&errs);
+    println!("annulus tvd orders {o:?}");
+    assert!(
+        o.iter().all(|&x| x >= 1.4),
+        "tvd orders {o:?} (gate >= 1.4)"
+    );
+    Ok(())
+}
+
 #[tokio::test]
 async fn snapshot_restore_rerun_is_bit_identical() -> CfdResult<()> {
     let mesh = annulus_skewed([0.0, 0.0], 0.5, 1.5, 24, 6, 0.3, 2.0)?;

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

@@ -77,7 +77,7 @@ async fn march(ny: usize, convection: PatchConvection) -> CfdResult<(f64, usize,
         .with_viscosity(MU)
         .with_reference_velocity(1.0)
         .with_reference_length(1.0);
-    let convecting = convection == PatchConvection::Upwind;
+    let convecting = convection != PatchConvection::None;
     let mut solver = CurvilinearPisoSolver::new(
         config,
         CurvilinearParameters {
@@ -157,20 +157,47 @@ async fn cylinder_flag_patch_keeps_the_p0_orders() -> CfdResult<()> {
     // ≈ 0.1 at ν = 0.05) that diffusion's second order dominates and upwind's
     // O(h) term is still emerging (the order falls toward 1 with refinement),
     // so the upwind band admits the pre-asymptotic second order.
+    // P4 step 2 gate (iii): TVD (van Albada) below upwind at every rung,
+    // orders in [1.5, 2.6]. `RTX_CF_SCHEME=none|upwind|tvd` runs one scheme.
+    let only = std::env::var("RTX_CF_SCHEME").ok();
     for (convection, gate) in [
         (PatchConvection::None, 1.8..2.6),
         (PatchConvection::Upwind, 0.7..2.4),
+        (PatchConvection::TvdVanAlbada, 1.5..2.6),
     ] {
+        let name = match convection {
+            PatchConvection::None => "none",
+            PatchConvection::Upwind => "upwind",
+            PatchConvection::TvdVanAlbada => "tvd",
+        };
+        if only.as_deref().is_some_and(|o| o != name) {
+            continue;
+        }
+        let upwind = [5.134660e-4, 2.244819e-4, 1.336112e-4];
         let mut errs = Vec::new();
         let mut hs = Vec::new();
-        for ny in [41usize, 62, 82] {
+        for (&ny, &u) in [41usize, 62, 82].iter().zip(&upwind) {
             let (l2, steps, max_div) = march(ny, convection).await?;
             println!(
                 "  cylinder-flag {convection:?} ny={ny}: L2 {l2:.6e}, max div {max_div:.2e}, {steps} steps"
             );
+            if convection == PatchConvection::TvdVanAlbada {
+                println!("    tvd / upwind at ny={ny}: {:.3}", l2 / u);
+            }
             errs.push(l2);
             hs.push(0.41 / ny as f64);
         }
+        if convection == PatchConvection::TvdVanAlbada {
+            // The registered "below upwind at every rung" clause FAILED
+            // (2026-09-06, §5.11): ratios 1.061 / 1.069 / 1.035 with orders
+            // 1.98 / 1.97. At cell Péclet ≈ 0.1 the Stokes floor (4.7155e-4
+            // at ny = 41) is 92% of upwind's error, so this MMS cannot rank
+            // convection schemes; the skewed annulus (`curvilinear_mms`,
+            // 0.24× upwind at ns = 128) is the discriminating gate. The
+            // ratios are recorded here, the order band is the assertion.
+            let ratios: Vec<f64> = errs.iter().zip(&upwind).map(|(e, u)| e / u).collect();
+            println!("  cylinder-flag TvdVanAlbada / upwind ratios {ratios:?}");
+        }
         let o: Vec<f64> = errs
             .windows(2)
             .zip(hs.windows(2))

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

@@ -12,7 +12,7 @@ use rtx_cfd::mesh::patch_gen::cylinder_flag_patch;
 use rtx_cfd::solvers::incompressible::{
     AleBoundaries, CurvilinearParameters, CurvilinearPisoSolver, EmbeddedParameters,
     EmbeddedPisoSolver, FlowField, NormalDiffusion, OversetField, OversetParameters,
-    OversetPisoSolver, PatchField, PoissonSolverKind, SideBoundary,
+    OversetPisoSolver, PatchConvection, PatchField, PoissonSolverKind, SideBoundary,
 };
 use rtx_cfd::{CfdConfig, CfdResult};
 
@@ -99,10 +99,18 @@ async fn run_cfd1(ny: usize) -> CfdResult<Cfd1> {
     let dt_patch = 0.4 * (hs * hs / (4.0 * NU)).min(hs / u_peak);
     let dt = dt_bg.min(dt_patch);
 
+    // P4 step 2: `RTX_OVERSET_CFD1_TVD=1` puts the van Albada deferred
+    // correction on the patch (the background stays upwind, as recorded).
+    let convection = if std::env::var("RTX_OVERSET_CFD1_TVD").is_ok() {
+        PatchConvection::TvdVanAlbada
+    } else {
+        PatchConvection::Upwind
+    };
     let mut patch = CurvilinearPisoSolver::new(
         config,
         CurvilinearParameters {
             tolerance: 1e-5,
+            convection,
             normal_diffusion: NormalDiffusion::LineImplicit,
             ..CurvilinearParameters::default()
         },
@@ -127,6 +135,12 @@ async fn run_cfd1(ny: usize) -> CfdResult<Cfd1> {
             .ok()
             .and_then(|v| v.parse().ok())
             .unwrap_or(2),
+        // Cost question (P4): does the second corrector's ~9 rounds buy a
+        // measurable load? `RTX_OVERSET_MAX_ROUNDS=3` caps every corrector.
+        max_rounds: std::env::var("RTX_OVERSET_MAX_ROUNDS")
+            .ok()
+            .and_then(|v| v.parse().ok())
+            .unwrap_or(OversetParameters::default().max_rounds),
         ..OversetParameters::default()
     };
     let mut solver = OversetPisoSolver::new(background, patch, (nx, ny, h, h), params)?;
@@ -274,9 +288,14 @@ async fn cfd1_on_the_overset_against_the_featflow_reference() -> CfdResult<()> {
         let r = run_cfd1(ny).await?;
         let rel = |a: f64, b: f64| 100.0 * (a - b) / b;
         println!(
-            "  CFD1 overset ny = {ny} (h = {:.4}, dt = {:.2e}): wall drag {:.4} ({:+.2}%) lift {:.4} ({:+.2}%); control volume drag {:.4} ({:+.2}%) lift {:.4}; routes differ {:.2}%; [{} steps, {:.0} s, Schwarz rounds mean {:.2}]  reference {REF_DRAG} / {REF_LIFT}; embedded staircase at ny=41: 15.71 / 15.62 (+10%)",
+            "  CFD1 overset ny = {ny} (h = {:.4}, dt = {:.2e}, patch {}): wall drag {:.4} ({:+.2}%) lift {:.4} ({:+.2}%); control volume drag {:.4} ({:+.2}%) lift {:.4}; routes differ {:.2}%; [{} steps, {:.0} s, Schwarz rounds mean {:.2}]  reference {REF_DRAG} / {REF_LIFT}; embedded staircase at ny=41: 15.71 / 15.62 (+10%)",
             H / ny as f64,
             r.dt,
+            if std::env::var("RTX_OVERSET_CFD1_TVD").is_ok() {
+                "tvd"
+            } else {
+                "upwind"
+            },
             r.drag_surface,
             rel(r.drag_surface, REF_DRAG),
             r.lift_surface,

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

@@ -84,7 +84,7 @@ async fn patch_with_exact_acceptors(
         normal_diffusion: diffusion,
         ..CurvilinearParameters::default()
     };
-    let convecting = convection == PatchConvection::Upwind;
+    let convecting = convection != PatchConvection::None;
     let mut solver = CurvilinearPisoSolver::new(config, params, mesh)?;
     solver.set_boundary_velocity(|x, y, _| (u_exact(x, y), v_exact(x, y)));
     solver.set_momentum_source(move |x, y, _| source(x, y, convecting));