rtx-cfd: overset P4 — CFD1 on the composite (tests/overset_cfd1.rs): ny=41 wall drag 15.2156 (+6.46%), lift 1.0879 (−2.78%), CV drag 15.528 (+8.64%), routes 2.0% apart (staircase +10%); 28550 steps, 2099 s at dt 4.57e-4 with 5.0 Schwarz rounds mean (stall rule off) — steady-march stall rule on for the next rungs

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

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

@@ -0,0 +1,257 @@
+//! A-P4 (`docs/overset_metal_campaign.md` §2.2 P4, §5.11): Turek–Hron CFD1
+//! (Re = 20, steady) on the OVERSET — the rigid harness's background
+//! (`turek_hron_cfd.rs`: parabolic inflow, outlet, multigrid, upwind) with
+//! the cylinder–flag O-grid as a static patch (no-slip wall, line-implicit
+//! across). Loads by the patch's wall stress (`surface_force`) and by the
+//! background's control-volume momentum balance (the two-route rule).
+//! Reference (FEATFLOW level 6): drag 14.2929, lift 1.11905. The embedded
+//! staircase measured drag 15.71 (surface) / 15.62 (CV) at ny = 41 (+10%).
+
+use rtx_cfd::mesh::PatchSide;
+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,
+};
+use rtx_cfd::{CfdConfig, CfdResult};
+
+const L: f64 = 2.5;
+const H: f64 = 0.41;
+const RHO: f64 = 1000.0;
+const NU: f64 = 1e-3;
+const U_MEAN: f64 = 0.2;
+const REF_DRAG: f64 = 14.2929;
+const REF_LIFT: f64 = 1.11905;
+
+fn inflow(y: f64) -> f64 {
+    1.5 * U_MEAN * y * (H - y) / (0.5 * H).powi(2)
+}
+
+struct Cfd1 {
+    drag_surface: f64,
+    lift_surface: f64,
+    drag_cv: f64,
+    lift_cv: f64,
+    steps: usize,
+    seconds: f64,
+    rounds_mean: f64,
+    dt: f64,
+}
+
+async fn run_cfd1(ny: usize) -> CfdResult<Cfd1> {
+    let h = H / ny as f64;
+    let nx = (L / h).round() as usize;
+    let mu = RHO * NU;
+    let config = CfdConfig::new()
+        .with_density(RHO)
+        .with_viscosity(mu)
+        .with_reference_velocity(U_MEAN)
+        .with_reference_length(0.1);
+    let mut background = EmbeddedPisoSolver::new(
+        config.clone(),
+        EmbeddedParameters {
+            corrector_steps: 2,
+            tolerance: 1e-7,
+            boundaries: AleBoundaries {
+                left: SideBoundary::Velocity,
+                right: SideBoundary::PressureOutlet,
+                bottom: SideBoundary::Velocity,
+                top: SideBoundary::Velocity,
+            },
+            poisson_solver: PoissonSolverKind::Multigrid,
+            ..EmbeddedParameters::default()
+        },
+    )?;
+    background.set_boundary_velocity(|x, y, _| {
+        if x <= 0.0 {
+            (inflow(y), 0.0)
+        } else {
+            (0.0, 0.0)
+        }
+    });
+
+    let (mesh, _) = cylinder_flag_patch(
+        [0.2, 0.2],
+        0.05,
+        0.01,
+        0.6,
+        h,
+        0.5 * 0.41 / 41.0,
+        6.0 * h,
+        12,
+        4.0,
+        500,
+    )?;
+    // The patch's explicit along-body diffusion limit (its wall row is
+    // line-implicit); the harness's combined criterion for the background.
+    let mut hs = f64::INFINITY;
+    for c in 0..mesh.cell_count() {
+        for (f, _) in mesh.cell_faces(c) {
+            if mesh.is_sface(f) {
+                let d = mesh.faces()[f].d;
+                hs = hs.min((d[0] * d[0] + d[1] * d[1]).sqrt());
+            }
+        }
+    }
+    let u_peak = 1.5 * 1.5 * U_MEAN;
+    let dt_bg = 0.25 / (2.0 * u_peak / h + 4.0 * NU / (h * h));
+    let dt_patch = 0.4 * (hs * hs / (4.0 * NU)).min(hs / u_peak);
+    let dt = dt_bg.min(dt_patch);
+
+    let mut patch = CurvilinearPisoSolver::new(
+        config,
+        CurvilinearParameters {
+            tolerance: 1e-5,
+            normal_diffusion: NormalDiffusion::LineImplicit,
+            ..CurvilinearParameters::default()
+        },
+        mesh,
+    )?;
+    patch.set_side_velocity(PatchSide::Inner, |_, _, _| (0.0, 0.0));
+    let mut patch_field = PatchField::new(patch.mesh());
+    patch.initialize(&mut patch_field, |_, _| (0.0, 0.0));
+
+    let mut bg_field = FlowField::new(nx, ny, h, h)?;
+    for j in 0..ny {
+        let u0 = inflow((j as f64 + 0.5) * h);
+        for i in 0..=nx {
+            bg_field.u[(j, i)] = u0;
+        }
+    }
+    let mut solver = OversetPisoSolver::new(
+        background,
+        patch,
+        (nx, ny, h, h),
+        OversetParameters::default(),
+    )?;
+    let mut field = OversetField {
+        background: bg_field,
+        patch: patch_field,
+    };
+    solver.initialize(&mut field)?;
+
+    let cv = (
+        (0.10 / h).round() as usize,
+        (0.75 / h).round() as usize,
+        (0.05 / h).round() as usize,
+        (0.36 / h).round() as usize,
+    );
+    let cv_force = |field: &OversetField, solver: &OversetPisoSolver| {
+        solver
+            .background()
+            .mask()
+            .expect("mask")
+            .control_volume_force(
+                &field.background.u,
+                &field.background.v,
+                &field.background.p,
+                &field.background.u_old,
+                &field.background.v_old,
+                dt,
+                RHO,
+                mu,
+                None,
+                cv,
+            )
+    };
+
+    let start = std::time::Instant::now();
+    let flow_through = L / U_MEAN;
+    let min_steps = (flow_through / dt).ceil() as usize;
+    let mut history: Vec<f64> = Vec::new();
+    let mut steps = 0;
+    let mut rounds_total = 0usize;
+    let mut correctors_total = 0usize;
+    loop {
+        let r = solver.advance(&mut field, dt).await?;
+        steps += 1;
+        rounds_total += r.rounds.iter().sum::<usize>();
+        correctors_total += r.rounds.len();
+        if steps % 50 == 0 {
+            let (fx, _) = cv_force(&field, &solver);
+            history.push(fx);
+            let load = solver
+                .patch()
+                .surface_force(&field.patch, PatchSide::Inner, solver.time());
+            let umax = field
+                .background
+                .u
+                .iter()
+                .fold(0.0_f64, |m, v| m.max(v.abs()));
+            if steps % 500 == 0 || !umax.is_finite() {
+                println!(
+                    "    ny = {ny}: step {steps} t = {:.2} s  drag_cv {fx:.4}  drag_wall {:.4} lift_wall {:.4}  max|u| {umax:.3}  rounds {:?} bg res {:.1e} patch div {:.1e}  [{:.0} s]",
+                    solver.time(),
+                    load.total()[0],
+                    load.total()[1],
+                    r.rounds,
+                    r.background_residual,
+                    r.patch_max_divergence,
+                    start.elapsed().as_secs_f64()
+                );
+            }
+            assert!(
+                umax.is_finite(),
+                "velocity became non-finite at step {steps}"
+            );
+            if steps >= min_steps && history.len() > 4 {
+                let now = history[history.len() - 1];
+                let then = history[history.len() - 5];
+                if ((now - then) / now).abs() < 1e-4 {
+                    break;
+                }
+            }
+        }
+        assert!(steps < 2_000_000, "CFD1 at ny = {ny} did not settle");
+    }
+    let seconds = start.elapsed().as_secs_f64();
+    let load = solver
+        .patch()
+        .surface_force(&field.patch, PatchSide::Inner, solver.time());
+    let (drag_cv, lift_cv) = cv_force(&field, &solver);
+    Ok(Cfd1 {
+        drag_surface: load.total()[0],
+        lift_surface: load.total()[1],
+        drag_cv,
+        lift_cv,
+        steps,
+        seconds,
+        rounds_mean: rounds_total as f64 / correctors_total.max(1) as f64,
+        dt,
+    })
+}
+
+#[tokio::test]
+async fn cfd1_on_the_overset_against_the_featflow_reference() -> CfdResult<()> {
+    let resolutions: Vec<usize> = std::env::var("RTX_OVERSET_CFD1_NY").ok().map_or_else(
+        || vec![41usize],
+        |list| {
+            list.split(',')
+                .map(|t| t.trim().parse().expect("ny list"))
+                .collect()
+        },
+    );
+    for &ny in &resolutions {
+        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%)",
+            H / ny as f64,
+            r.dt,
+            r.drag_surface,
+            rel(r.drag_surface, REF_DRAG),
+            r.lift_surface,
+            rel(r.lift_surface, REF_LIFT),
+            r.drag_cv,
+            rel(r.drag_cv, REF_DRAG),
+            r.lift_cv,
+            100.0 * ((r.drag_surface - r.drag_cv) / r.drag_cv).abs(),
+            r.steps,
+            r.seconds,
+            r.rounds_mean
+        );
+        assert!(r.drag_surface.is_finite() && r.drag_cv.is_finite());
+    }
+    Ok(())
+}