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

//! 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, CellClass, CurvilinearParameters, CurvilinearPisoSolver, EmbeddedParameters,
    EmbeddedPisoSolver, FlowField, 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 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);

    // 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()
        },
        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;
        }
    }
    // A steady march: stop a corrector's rounds when two rounds make no
    // progress (the noise floor); measured at ny = 41 without it: 5.0 rounds
    // per corrector on average (the second corrector 8 every step), 2099 s.
    let params = OversetParameters {
        stall_rounds: std::env::var("RTX_OVERSET_STALL")
            .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)?;
    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;
    let max_steps: usize = std::env::var("RTX_OVERSET_CFD1_MAX_STEPS")
        .ok()
        .and_then(|v| v.parse().ok())
        .unwrap_or(2_000_000);
    let trace_first = std::env::var("RTX_OVERSET_CFD1_TRACE").is_ok();
    loop {
        let r = solver.advance(&mut field, dt).await?;
        steps += 1;
        let every: usize = std::env::var("RTX_OVERSET_CFD1_TRACE_EVERY")
            .ok()
            .and_then(|v| v.parse().ok())
            .unwrap_or(0);
        if (trace_first && steps <= 6) || (every > 0 && steps % every == 0) {
            let pmax = field
                .background
                .p
                .iter()
                .fold(0.0_f64, |m, v| m.max(v.abs()));
            let ppmax = field.patch.p.iter().fold(0.0_f64, |m, v| m.max(v.abs()));
            let upmax = field.patch.u.iter().fold(0.0_f64, |m, v| m.max(v.abs()));
            println!(
                "      step {steps}: rounds {:?} converged {} stalled {} bg res {:.2e} patch div {:.2e} patch iters {} conv {} | max|p| bg {pmax:.3e} patch {ppmax:.3e} max|u| patch {upmax:.3e} | defect bg {:.2e} patch {:.2e}",
                r.rounds,
                r.schwarz_converged,
                r.schwarz_stalled,
                r.background_residual,
                r.patch_max_divergence,
                r.patch_poisson_iterations,
                r.patch_converged,
                r.background_mass_defect,
                r.patch_mass_defect
            );
        }
        if steps >= max_steps {
            break;
        }
        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);
    // P4 momentum-defect measurement: the force the background transmits
    // into the ring (fringe + hole), into the hole alone, and the patch's
    // wall force — consecutive differences are the active region's
    // residual, the fringe ring's momentum defect, and the patch region's.
    let ring = solver
        .overlap()
        .region_force(&field.background, RHO, mu, |c| c != CellClass::Active);
    let hole = solver
        .overlap()
        .region_force(&field.background, RHO, mu, |c| c == CellClass::Hole);
    let wall = load.total();
    println!(
        "    momentum routes ny = {ny}: CV box ({drag_cv:.4}, {lift_cv:.4}) | ring outer ({:.4}, {:.4}) | hole boundary ({:.4}, {:.4}) | wall ({:.4}, {:.4}); defects [% of wall drag]: active {:+.2} fringe ring {:+.2} patch region {:+.2}",
        ring.0,
        ring.1,
        hole.0,
        hole.1,
        wall[0],
        wall[1],
        100.0 * (drag_cv - ring.0) / wall[0],
        100.0 * (ring.0 - hole.0) / wall[0],
        100.0 * (hole.0 - wall[0]) / wall[0],
    );
    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}, 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,
            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(())
}