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

//! 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(())
}