rtx-cfd: overset_cfd23 — Turek–Hron CFD2/CFD3 on the overset (the harness's ramped inflow from rest, TVD background and patch, 3-round cap), loads by the patch wall stress with the solver-flux-form box and the CV formula beside it as window statistics, CFD3 frequency from lift crossings, the solver-metric momentum chain at the final state (solved-face residual is the pin), save/load of the fields by case; momentum_residual checks the far-upwind neighbours explicitly under a limited background (the van Albada limiter swallows a NaN silently)

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

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

@@ -0,0 +1,572 @@
+//! P4 close-out (`docs/overset_metal_campaign.md` §5.11): Turek–Hron CFD2
+//! (steady, Re = 100) and CFD3 (periodic shedding, Re = 200) on the
+//! OVERSET — the rigid harness's background (`turek_hron_cfd23.rs`: the
+//! benchmark's inflow ramp from rest, outlet, multigrid, van Albada TVD)
+//! with the cylinder–flag O-grid as a static patch under TVD, every
+//! corrector capped at 3 Schwarz rounds (the P5 budget setting). Loads by
+//! the patch's wall stress (the route the momentum audit of §5.11 settled
+//! on) with the background's box in the solver's own flux form and the CV
+//! formula recorded beside it, all as time statistics over the
+//! benchmark's window; the solver-metric momentum chain is printed at the
+//! final state as the audit (its solved-face residual is the pin).
+//!
+//! References: CFD2 drag 136.700, lift 10.5343; CFD3 drag 439.45 ± 5.62,
+//! lift −11.893 ± 437.81, f = 4.3956 Hz (level 4, dt = 0.005). The
+//! embedded staircase read CFD2 −12.3 / −11.2 / −10.3 % and CFD3 −6.9 /
+//! −6.0 / −10.3 % (drag) at ny = 41 / 62 / 82. Numbers are recorded, not
+//! asserted, until the ladder is seen; `RTX_OVERSET_CFD2_NY` /
+//! `RTX_OVERSET_CFD3_NY` (default 41), `RTX_OVERSET_CFD23_T_END` (smoke),
+//! `RTX_OVERSET_MAX_ROUNDS` (3), `RTX_OVERSET_ROWS` (4),
+//! `RTX_OVERSET_CFD1_SAVE` / `_LOAD` (fields, tagged by case).
+
+use rtx_cfd::mesh::patch_gen::cylinder_flag_patch;
+use rtx_cfd::mesh::PatchSide;
+use rtx_cfd::solvers::incompressible::{
+    AleBoundaries, ConvectionScheme, CurvilinearParameters, CurvilinearPisoSolver,
+    EmbeddedParameters, EmbeddedPisoSolver, FlowField, MgPrecision, NormalDiffusion, OversetField,
+    OversetParameters, OversetPisoSolver, PatchConvection, 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 CFD2_U: f64 = 1.0;
+const CFD2_REF_DRAG: f64 = 136.700;
+const CFD2_REF_LIFT: f64 = 10.5343;
+
+const CFD3_U: f64 = 2.0;
+const CFD3_REF_DRAG_MEAN: f64 = 439.45;
+const CFD3_REF_DRAG_AMP: f64 = 5.6183;
+const CFD3_REF_LIFT_MEAN: f64 = -11.893;
+const CFD3_REF_LIFT_AMP: f64 = 437.81;
+const CFD3_REF_FREQUENCY: f64 = 4.3956;
+
+/// The ramped parabolic inflow of the benchmark definition.
+fn inflow(u_mean: f64, y: f64, t: f64) -> f64 {
+    let ramp = if t < 2.0 {
+        0.5 * (1.0 - (std::f64::consts::PI * t / 2.0).cos())
+    } else {
+        1.0
+    };
+    ramp * 1.5 * u_mean * y * (H - y) / (0.5 * H).powi(2)
+}
+
+fn env_usize(var: &str, default: usize) -> usize {
+    std::env::var(var)
+        .ok()
+        .and_then(|v| v.parse().ok())
+        .unwrap_or(default)
+}
+
+fn ny_list(var: &str, default: &[usize]) -> Vec<usize> {
+    std::env::var(var)
+        .ok()
+        .map(|s| {
+            s.split(',')
+                .map(|t| t.trim().parse().expect("integer ny"))
+                .collect()
+        })
+        .unwrap_or_else(|| default.to_vec())
+}
+
+fn overlap_rows() -> usize {
+    env_usize(
+        "RTX_OVERSET_ROWS",
+        OversetParameters::default().overlap_rows,
+    )
+}
+
+fn field_tag(case: &str, ny: usize) -> String {
+    let rows = overlap_rows();
+    format!(
+        "{case}_ny{ny}_tvd{}",
+        if rows == OversetParameters::default().overlap_rows {
+            String::new()
+        } else {
+            format!("_rows{rows}")
+        }
+    )
+}
+
+struct Composite {
+    solver: OversetPisoSolver,
+    field: OversetField,
+    ny: usize,
+    h: f64,
+    dt: f64,
+    mu: f64,
+    cv: (usize, usize, usize, usize),
+}
+
+impl Composite {
+    fn new(u_mean: f64, ny: usize) -> CfdResult<Self> {
+        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,
+                poisson_precision: MgPrecision::F64,
+                // The harness's finding: upwind's numerical viscosity
+                // suppressed CFD3's shedding entirely.
+                convection_scheme: ConvectionScheme::TvdVanAlbada,
+            },
+        )?;
+        background.set_boundary_velocity(move |x, y, t| {
+            if x <= 0.0 {
+                (inflow(u_mean, y, t), 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,
+        )?;
+        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,
+                convection: PatchConvection::TvdVanAlbada,
+                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));
+
+        // At rest: the ramp brings the inflow up from zero.
+        let bg_field = FlowField::new(nx, ny, h, h)?;
+        let params = OversetParameters {
+            stall_rounds: env_usize("RTX_OVERSET_STALL", 2),
+            max_rounds: env_usize("RTX_OVERSET_MAX_ROUNDS", 3),
+            overlap_rows: overlap_rows(),
+            ..OversetParameters::default()
+        };
+        let mut solver = OversetPisoSolver::new(background, patch, (nx, ny, h, h), params)?;
+        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,
+        );
+        Ok(Self {
+            solver,
+            field,
+            ny,
+            h,
+            dt,
+            mu,
+            cv,
+        })
+    }
+
+    fn wall(&self) -> [f64; 2] {
+        self.solver
+            .patch()
+            .surface_force(&self.field.patch, PatchSide::Inner, self.solver.time())
+            .total()
+    }
+
+    fn cv_force(&self) -> (f64, f64) {
+        self.solver
+            .background()
+            .mask()
+            .expect("mask")
+            .control_volume_force(
+                &self.field.background.u,
+                &self.field.background.v,
+                &self.field.background.p,
+                &self.field.background.u_old,
+                &self.field.background.v_old,
+                self.dt,
+                RHO,
+                self.mu,
+                None,
+                self.cv,
+            )
+    }
+
+    /// The solver-metric momentum chain at the current state (§5.11):
+    /// residual buckets, the box in the solver's flux form, the patch's
+    /// balance. Returns the solved-far Σ|r| (the pin).
+    fn chain(&self, case: &str) -> f64 {
+        let (ny, h, dt) = (self.ny, self.h, self.dt);
+        let wall = self.wall();
+        let mr = self.solver.momentum_residual(&self.field, dt);
+        let ring = mr.fringe_fringe.fx + mr.fringe_hole.fx;
+        let mut bands = [
+            (0.0_f64, 0.20, 0.0_f64),
+            (0.20, 0.30, 0.0),
+            (0.30, 0.55, 0.0),
+            (0.55, 1.0, 0.0),
+        ];
+        for f in mr.prescribed.iter().filter(|f| f.is_u && f.r.is_finite()) {
+            let x = f.i as f64 * h;
+            if let Some(b) = bands.iter_mut().find(|b| x >= b.0 && x < b.1) {
+                b.2 += f.r;
+            }
+        }
+        let boxes = [
+            (0.10, 0.75, 0.05, 0.36),
+            (0.09, 0.70, 0.07, 0.34),
+            (0.08, 1.00, 0.03, 0.38),
+        ];
+        let forces: Vec<(f64, f64)> = boxes
+            .iter()
+            .map(|&(x0, x1, y0, y1)| {
+                self.solver.solver_metric_force(
+                    &self.field,
+                    dt,
+                    (
+                        (x0 / h).round() as usize,
+                        (x1 / h).round() as usize,
+                        (y0 / h).round() as usize,
+                        (y1 / h).round() as usize,
+                    ),
+                )
+            })
+            .collect();
+        let sf = forces[0];
+        let spread = forces
+            .iter()
+            .fold(0.0_f64, |m, f| m.max((f.0 - sf.0).abs()));
+        let (cvx, _) = self.cv_force();
+        let pb = self
+            .solver
+            .patch()
+            .momentum_balance(&self.field.patch, self.solver.time());
+        let ff = pb.flux_force();
+        let fw = pb.wall_force();
+        let bal = pb.balance();
+        let hole = sf.0 + ring;
+        println!(
+            "    {case} chain ny = {ny} at t = {:.3} [N/m]: solved far Σ|r| ({:.2e}, {:.2e}) {}/{} | near ring Σr ({:+.1e}, {:+.1e}) δ {:.2e} Pa | box (solver flux form) {:.3} [3 boxes spread {:.1e}; CV formula {:.3}] → ring Σr {:+.3} ({} + {} faces; x-bands {}) → hole flux {:.3} → band {:+.3} → patch interface {:.3} → interior {:+.3} (δP {:+.3}, balance residual {:+.3} — CFD3's unsteady term is not stored) → wall, scheme fluxes {:.3} → wall formula {:+.3} → wall {:.3}; total wall − box {:+.3} ({:+.2} %)",
+            self.solver.time(),
+            mr.solved_far.abs_x,
+            mr.solved_far.abs_y,
+            mr.solved_far.evaluated,
+            mr.solved_far.total,
+            mr.solved_near.fx,
+            mr.solved_near.fy,
+            mr.level_offset(h),
+            sf.0,
+            spread,
+            cvx,
+            ring,
+            mr.fringe_fringe.evaluated,
+            mr.fringe_hole.evaluated,
+            bands
+                .iter()
+                .map(|b| format!("{:.2}–{:.2}: {:+.3}", b.0, b.1, b.2))
+                .collect::<Vec<_>>()
+                .join(", "),
+            hole,
+            ff[0] - hole,
+            ff[0],
+            fw[0] - ff[0],
+            pb.pressure_defect()[0],
+            bal[0],
+            fw[0],
+            wall[0] - fw[0],
+            wall[0],
+            wall[0] - sf.0,
+            100.0 * (wall[0] - sf.0) / wall[0],
+        );
+        mr.solved_far.abs_x.max(mr.solved_far.abs_y)
+    }
+
+    fn save_or_load(&mut self, case: &str) -> CfdResult<bool> {
+        let tag = field_tag(case, self.ny);
+        if let Ok(dir) = std::env::var("RTX_OVERSET_CFD1_LOAD") {
+            let dir = std::path::Path::new(&dir);
+            self.field.background = FlowField::load(&dir.join(format!("bg_{tag}.bin")))?;
+            let read = |name: &str| -> Vec<f64> {
+                let bytes = std::fs::read(dir.join(format!("patch_{tag}_{name}.bin")))
+                    .unwrap_or_else(|e| panic!("load patch {name}: {e}"));
+                bytes
+                    .chunks_exact(8)
+                    .map(|c| f64::from_le_bytes(c.try_into().expect("8 bytes")))
+                    .collect()
+            };
+            self.field.patch.u = read("u");
+            self.field.patch.v = read("v");
+            self.field.patch.p = read("p");
+            self.field.patch.flux = read("flux");
+            println!("    loaded {tag} from {}", dir.display());
+            return Ok(true);
+        }
+        Ok(false)
+    }
+
+    fn save(&self, case: &str) -> CfdResult<()> {
+        if let Ok(dir) = std::env::var("RTX_OVERSET_CFD1_SAVE") {
+            let tag = field_tag(case, self.ny);
+            let dir = std::path::Path::new(&dir);
+            std::fs::create_dir_all(dir).expect("save dir");
+            self.field
+                .background
+                .save(&dir.join(format!("bg_{tag}.bin")))?;
+            for (name, vals) in [
+                ("u", &self.field.patch.u),
+                ("v", &self.field.patch.v),
+                ("p", &self.field.patch.p),
+                ("flux", &self.field.patch.flux),
+            ] {
+                let bytes: Vec<u8> = vals.iter().flat_map(|x| x.to_le_bytes()).collect();
+                std::fs::write(dir.join(format!("patch_{tag}_{name}.bin")), bytes)
+                    .expect("save patch");
+            }
+            println!("    saved {tag} to {}", dir.display());
+        }
+        Ok(())
+    }
+}
+
+/// One sampled series of the three load routes.
+struct Series {
+    times: Vec<f64>,
+    wall_drag: Vec<f64>,
+    wall_lift: Vec<f64>,
+    box_drag: Vec<f64>,
+    box_lift: Vec<f64>,
+    cv_drag: Vec<f64>,
+    steps: usize,
+    seconds: f64,
+    rounds_mean: f64,
+    dt: f64,
+    pin: f64,
+}
+
+/// March to `t_end`, sampling every 25 steps once `t >= t_start`; the
+/// chain at the end (`RTX_OVERSET_CFD23_T_END` overrides the end for a
+/// smoke run; a loaded field set skips the march).
+async fn run_sampled(
+    case: &str,
+    u_mean: f64,
+    ny: usize,
+    t_start: f64,
+    t_end: f64,
+) -> CfdResult<Series> {
+    let mut c = Composite::new(u_mean, ny)?;
+    let t_end = std::env::var("RTX_OVERSET_CFD23_T_END")
+        .ok()
+        .and_then(|v| v.parse().ok())
+        .unwrap_or(t_end);
+    let loaded = c.save_or_load(case)?;
+    let start = std::time::Instant::now();
+    let mut s = Series {
+        times: Vec::new(),
+        wall_drag: Vec::new(),
+        wall_lift: Vec::new(),
+        box_drag: Vec::new(),
+        box_lift: Vec::new(),
+        cv_drag: Vec::new(),
+        steps: 0,
+        seconds: 0.0,
+        rounds_mean: 0.0,
+        dt: c.dt,
+        pin: 0.0,
+    };
+    let (mut rounds_total, mut correctors_total) = (0usize, 0usize);
+    while !loaded && c.solver.time() < t_end {
+        let r = c.solver.advance(&mut c.field, c.dt).await?;
+        s.steps += 1;
+        rounds_total += r.rounds.iter().sum::<usize>();
+        correctors_total += r.rounds.len();
+        let umax = c
+            .field
+            .background
+            .u
+            .iter()
+            .fold(0.0_f64, |m, v| m.max(v.abs()));
+        assert!(
+            umax.is_finite(),
+            "{case} ny = {ny}: velocity became non-finite at t = {:.3}",
+            c.solver.time()
+        );
+        if s.steps % 25 == 0 && c.solver.time() >= t_start {
+            let w = c.wall();
+            let b = c.solver.solver_metric_force(&c.field, c.dt, c.cv);
+            let (cvx, _) = c.cv_force();
+            s.times.push(c.solver.time());
+            s.wall_drag.push(w[0]);
+            s.wall_lift.push(w[1]);
+            s.box_drag.push(b.0);
+            s.box_lift.push(b.1);
+            s.cv_drag.push(cvx);
+        }
+        if s.steps % 2000 == 0 {
+            let w = c.wall();
+            println!(
+                "    {case} ny = {ny}: step {} t = {:.3} s  wall drag {:.3} lift {:.3}  max|u| {umax:.3}  rounds {:?}  [{:.0} s]",
+                s.steps,
+                c.solver.time(),
+                w[0],
+                w[1],
+                r.rounds,
+                start.elapsed().as_secs_f64()
+            );
+        }
+    }
+    s.seconds = start.elapsed().as_secs_f64();
+    s.rounds_mean = rounds_total as f64 / correctors_total.max(1) as f64;
+    s.pin = c.chain(case);
+    c.save(case)?;
+    Ok(s)
+}
+
+/// Mid-range mean and half-range amplitude of a series.
+fn mid_amp(series: &[f64]) -> (f64, f64) {
+    if series.is_empty() {
+        return (f64::NAN, f64::NAN);
+    }
+    let max = series.iter().copied().fold(f64::MIN, f64::max);
+    let min = series.iter().copied().fold(f64::MAX, f64::min);
+    (0.5 * (max + min), 0.5 * (max - min))
+}
+
+/// Frequency from linearly-interpolated upward zero crossings about the
+/// mid-range; `None` with fewer than four crossings.
+fn crossing_frequency(times: &[f64], series: &[f64]) -> Option<f64> {
+    let (mean, _) = mid_amp(series);
+    let mut crossings: Vec<f64> = Vec::new();
+    for k in 1..series.len() {
+        let (a, b) = (series[k - 1] - mean, series[k] - mean);
+        if a < 0.0 && b >= 0.0 {
+            let frac = a / (a - b);
+            crossings.push(times[k - 1] + frac * (times[k] - times[k - 1]));
+        }
+    }
+    (crossings.len() >= 4).then(|| {
+        (crossings.len() - 1) as f64 / (crossings.last().unwrap() - crossings.first().unwrap())
+    })
+}
+
+fn pct(a: f64, b: f64) -> f64 {
+    100.0 * (a - b) / b
+}
+
+#[tokio::test]
+async fn cfd2_on_the_overset() -> CfdResult<()> {
+    for &ny in &ny_list("RTX_OVERSET_CFD2_NY", &[41]) {
+        let r = run_sampled("cfd2", CFD2_U, ny, 8.0, 10.0).await?;
+        let (drag, drag_amp) = mid_amp(&r.wall_drag);
+        let (lift, lift_amp) = mid_amp(&r.wall_lift);
+        let (bdrag, _) = mid_amp(&r.box_drag);
+        let (blift, _) = mid_amp(&r.box_lift);
+        let (cdrag, _) = mid_amp(&r.cv_drag);
+        println!(
+            "  CFD2 overset ny = {ny} (h = {:.4}, dt = {:.2e}, tvd/tvd, rows {}): wall drag {drag:.3} ± {drag_amp:.3} ({:+.2} %) lift {lift:.3} ± {lift_amp:.3} ({:+.2} %); box (solver flux form) drag {bdrag:.3} ({:+.2} %) lift {blift:.3}; CV formula drag {cdrag:.3}; {} samples [{} steps, {:.0} s, rounds mean {:.2}]  reference {CFD2_REF_DRAG} / {CFD2_REF_LIFT}; embedded staircase −12.3 / −11.2 / −10.3 % at ny 41 / 62 / 82",
+            H / ny as f64,
+            r.dt,
+            overlap_rows(),
+            pct(drag, CFD2_REF_DRAG),
+            pct(lift, CFD2_REF_LIFT),
+            pct(bdrag, CFD2_REF_DRAG),
+            r.times.len(),
+            r.steps,
+            r.seconds,
+            r.rounds_mean
+        );
+        assert!(
+            r.pin <= 1e-9 * drag.abs().max(1.0),
+            "solved-face residual {:.3e}",
+            r.pin
+        );
+        if !r.times.is_empty() {
+            assert!(drag.is_finite() && lift.is_finite());
+        }
+    }
+    Ok(())
+}
+
+#[tokio::test]
+async fn cfd3_on_the_overset() -> CfdResult<()> {
+    for &ny in &ny_list("RTX_OVERSET_CFD3_NY", &[41]) {
+        let r = run_sampled("cfd3", CFD3_U, ny, 6.0, 9.0).await?;
+        let (drag, drag_amp) = mid_amp(&r.wall_drag);
+        let (lift, lift_amp) = mid_amp(&r.wall_lift);
+        let f = crossing_frequency(&r.times, &r.wall_lift);
+        let half = r.wall_lift.len() / 2;
+        let (_, amp_first) = mid_amp(&r.wall_lift[..half]);
+        let (_, amp_second) = mid_amp(&r.wall_lift[half..]);
+        let (bdrag, bdrag_amp) = mid_amp(&r.box_drag);
+        let (blift, blift_amp) = mid_amp(&r.box_lift);
+        let fb = crossing_frequency(&r.times, &r.box_lift);
+        println!(
+            "  CFD3 overset ny = {ny} (h = {:.4}, dt = {:.2e}, tvd/tvd, rows {}): wall drag {drag:.2} ± {drag_amp:.2} ({:+.2} % / amp {:+.1} %), lift {lift:.2} ± {lift_amp:.2} (amp {:+.2} %), f = {f:?} Hz ({:+.2} %); half-window lift amps {amp_first:.2} / {amp_second:.2}; box (solver flux form) drag {bdrag:.2} ± {bdrag_amp:.2} lift {blift:.2} ± {blift_amp:.2} f {fb:?}; {} samples [{} steps, {:.0} s, rounds mean {:.2}]  reference drag {CFD3_REF_DRAG_MEAN} ± {CFD3_REF_DRAG_AMP}, lift {CFD3_REF_LIFT_MEAN} ± {CFD3_REF_LIFT_AMP}, f {CFD3_REF_FREQUENCY}; embedded staircase drag −6.9 / −6.0 / −10.3 %, f −2.8 / −1.3 / −0.04 % at ny 41 / 62 / 82",
+            H / ny as f64,
+            r.dt,
+            overlap_rows(),
+            pct(drag, CFD3_REF_DRAG_MEAN),
+            pct(drag_amp, CFD3_REF_DRAG_AMP),
+            pct(lift_amp, CFD3_REF_LIFT_AMP),
+            f.map_or(f64::NAN, |f| pct(f, CFD3_REF_FREQUENCY)),
+            r.times.len(),
+            r.steps,
+            r.seconds,
+            r.rounds_mean
+        );
+        assert!(
+            r.pin <= 1e-9 * drag.abs().max(1.0),
+            "solved-face residual {:.3e}",
+            r.pin
+        );
+        if !r.times.is_empty() && std::env::var("RTX_OVERSET_CFD23_T_END").is_err() {
+            assert!(
+                f.is_some(),
+                "the wake must shed: fewer than four lift zero-crossings"
+            );
+        }
+    }
+    Ok(())
+}