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

//! Quantitative cavity benchmark against Ghia, Ghia & Shin (1982) at Re = 400.
//!
//! The Re = 100 cavity in `tests/simple_tests.rs` asserts bands measured on
//! its own grid; this test makes the quantitative claim, on a domain that is
//! *exactly* the unit square (`n` cells of `1/n` — the older cavity test's
//! `65 x (1/64)` grid is 1.6% oversized, which is fine for banded assertions
//! and not fine for a table comparison). Velocities are sampled on the
//! staggered faces that lie exactly on the centrelines: u faces at
//! `x = 0.5`, v faces at `y = 0.5`.
//!
//! # Reference values, and where they come from
//!
//! Ghia et al. (1982), 257^2 multigrid. Sourced from two independent
//! transcriptions that agree digit for digit: the tabulation in Mramor,
//! Vertnik & Sarler, CMC vol. 36 (2013), Table 1, and the benchmark data
//! collection at gist.github.com/ivan-pi. For Re = 400:
//!
//! ```text
//!     u on x = 0.5:  min -0.32726 at y = 0.2813
//!     v on y = 0.5:  min -0.44993 at x = 0.8594,  max 0.30203 at x = 0.2266
//! ```
//!
//! plus the 15 interior u-profile stations asserted below. Worth knowing
//! when reading the tolerances: Ghia's own values carry their own
//! discretisation error — Sahin & Owens (2003) on finer grids put the u
//! minimum at -0.32838 (+0.3%) — so agreement with Ghia much tighter than a
//! few tenths of a percent would be overfitting to the reference, not
//! accuracy.
//!
//! # The stopping-tolerance trap this test walked into
//!
//! The first run used the Re = 100 test's residual tolerance of 1e-4 and
//! read u_min = -0.31987 — "converged", 2.3% shy of Ghia — and refining to
//! 192^2 made it WORSE (-0.30879, profile deviation doubled). The residual
//! had dropped below tolerance while the field was still developing, and
//! the effect grows with mesh size because SIMPLE's per-iteration
//! contraction weakens as h -> 0: at a fixed residual level, the finer mesh
//! is at an *earlier* stage of convergence. Tightening the stop, not the
//! mesh, was the fix — 3e-5 gave -0.32565 and 1.5e-5 gave -0.32667, at
//! which point the movement per halving (~0.3%, then ~0.1%) is below the
//! reference's own error and the answer is converged for this comparison.
//! This is why a fixed-tolerance grid study is not presented here, and why
//! "the residual converged" must never stand in for "the answer stopped
//! moving".

use rtx_cfd::solvers::incompressible::{
    BoundaryConditions, BoundaryLocation, BoundaryType, ConvectionScheme, FlowField,
    IncompressibleSolver, SimpleParameters, SimpleSolver,
};
use rtx_cfd::{CfdConfig, CfdResult};

/// Ghia et al. (1982), Table I, Re = 400: u on the vertical centreline at
/// the tabulated y stations (interior rows only; the 0 and 1 endpoints are
/// boundary conditions, not solution).
const GHIA_U_RE400: &[(f64, f64)] = &[
    (0.9766, 0.75837),
    (0.9688, 0.68439),
    (0.9609, 0.61756),
    (0.9531, 0.55892),
    (0.8516, 0.29093),
    (0.7344, 0.16256),
    (0.6172, 0.02135),
    (0.5000, -0.11477),
    (0.4531, -0.17119),
    (0.2813, -0.32726),
    (0.1719, -0.24299),
    (0.1016, -0.14612),
    (0.0703, -0.10338),
    (0.0625, -0.09266),
    (0.0547, -0.08186),
];

struct CavitySolution {
    field: FlowField,
    n: usize,
    iterations: usize,
}

async fn solve_cavity_re400(n: usize) -> CfdResult<CavitySolution> {
    let dx = 1.0 / n as f64;

    let config = CfdConfig::new()
        .with_density(1.0)
        .with_viscosity(1.0 / 400.0)
        .with_reference_velocity(1.0)
        .with_reference_length(1.0);

    // The corner-singularity floor on the normalised mass residual falls
    // roughly linearly with mesh size (1.6e-4 at 65^2); 1e-4 sits above the
    // ~8e-5 floor of this grid.
    let params = SimpleParameters::default()
        .with_pressure_relaxation(0.3)
        .with_velocity_relaxation(0.7)
        .with_max_iterations(40000)
        .with_convection_scheme(ConvectionScheme::TvdVanAlbada)
        // See the module docs: 1e-4 "converges" 2.3% short of Ghia, and the
        // shortfall grows with mesh size. 1.5e-5 is where the answer stops
        // moving relative to the reference's own accuracy.
        .with_tolerance(1.5e-5);

    let mut solver = SimpleSolver::new(config, params)?;

    let mut field = FlowField::new(n, n, dx, dx)?;
    let mut bcs = BoundaryConditions::new();
    for location in [
        BoundaryLocation::Left,
        BoundaryLocation::Right,
        BoundaryLocation::Bottom,
        BoundaryLocation::Top,
    ] {
        bcs.add_boundary_condition(location, BoundaryType::FreeSlipWall);
    }
    solver.set_wall_velocity(move |_x, y| if y > 0.5 { (1.0, 0.0) } else { (0.0, 0.0) });
    field.apply_boundary_conditions(&bcs)?;

    let result = solver.solve(&mut field, &bcs).await?;
    assert!(
        result.solver_result.converged,
        "cavity did not converge: residual {:.3e} after {} iterations",
        result.solver_result.final_residual, result.solver_result.iterations
    );

    Ok(CavitySolution {
        field,
        n,
        iterations: result.solver_result.iterations,
    })
}

/// u on the vertical centreline at height y, sampled from the u faces at
/// exactly x = 0.5 (face i = n/2) and linearly interpolated between the face
/// rows at y = (j + 0.5)/n.
fn u_centerline(solution: &CavitySolution, y: f64) -> f64 {
    let n = solution.n;
    let mid = n / 2;
    let position = y * n as f64 - 0.5;
    let j0 = position.floor().clamp(0.0, (n - 2) as f64) as usize;
    let t = (position - j0 as f64).clamp(0.0, 1.0);
    (1.0 - t) * solution.field.u[(j0, mid)] + t * solution.field.u[(j0 + 1, mid)]
}

#[tokio::test]
async fn cavity_re400_matches_ghia() -> CfdResult<()> {
    let solution = solve_cavity_re400(128).await?;
    let n = solution.n;

    // u extremum on the vertical centreline, from the face values directly.
    let mid = n / 2;
    let mut u_min = f64::INFINITY;
    let mut y_at_min = 0.0;
    for j in 0..n {
        let u = solution.field.u[(j, mid)];
        if u < u_min {
            u_min = u;
            y_at_min = (j as f64 + 0.5) / n as f64;
        }
    }

    // v extrema on the horizontal centreline: v faces at y = 0.5 exactly.
    let jmid = n / 2;
    let mut v_min = f64::INFINITY;
    let mut v_max = f64::NEG_INFINITY;
    let mut x_at_vmin = 0.0;
    let mut x_at_vmax = 0.0;
    for i in 0..n {
        let v = solution.field.v[(jmid, i)];
        let x = (i as f64 + 0.5) / n as f64;
        if v < v_min {
            v_min = v;
            x_at_vmin = x;
        }
        if v > v_max {
            v_max = v;
            x_at_vmax = x;
        }
    }

    println!("  {n}^2 TVD Re=400 ({} iterations)", solution.iterations);
    println!("  u_min = {u_min:.5} at y = {y_at_min:.4}   (Ghia: -0.32726 at 0.2813)");
    println!("  v_min = {v_min:.5} at x = {x_at_vmin:.4}   (Ghia: -0.44993 at 0.8594)");
    println!("  v_max = {v_max:.5} at x = {x_at_vmax:.4}   (Ghia:  0.30203 at 0.2266)");

    let mut max_profile_error: f64 = 0.0;
    for &(y, ghia_u) in GHIA_U_RE400 {
        let u = u_centerline(&solution, y);
        let error = (u - ghia_u).abs();
        max_profile_error = max_profile_error.max(error);
        println!("    y = {y:.4}   u = {u:+.5}   Ghia = {ghia_u:+.5}   diff = {error:.4}");
    }
    println!("  max |u - Ghia| over the profile = {max_profile_error:.4}");

    // Measured at 128^2, tolerance 1.5e-5: u_min = -0.32667 at y = 0.2852,
    // v_min = -0.45024 at x = 0.8633, v_max = 0.30044 at x = 0.2305, max
    // profile deviation 0.0051. Bands are set at 2-3x the measured gaps —
    // tight enough that the 1e-4-tolerance state (-0.31987, deviation
    // 0.0335) fails all of them.
    assert!(
        (u_min + 0.32726).abs() < 0.002,
        "u_min = {u_min:.5}, more than 0.6% from Ghia's -0.32726"
    );
    assert!(
        (y_at_min - 0.2813).abs() < 0.02,
        "u minimum at y = {y_at_min:.4}, Ghia puts it at 0.2813"
    );
    assert!(
        (v_min + 0.44993).abs() < 0.002,
        "v_min = {v_min:.5}, Ghia: -0.44993"
    );
    assert!(
        (v_max - 0.30203).abs() < 0.004,
        "v_max = {v_max:.5}, Ghia: 0.30203"
    );
    assert!(
        (x_at_vmin - 0.8594).abs() < 0.02,
        "v minimum at x = {x_at_vmin:.4}, Ghia puts it at 0.8594"
    );
    assert!(
        (x_at_vmax - 0.2266).abs() < 0.02,
        "v maximum at x = {x_at_vmax:.4}, Ghia puts it at 0.2266"
    );
    assert!(
        max_profile_error < 0.01,
        "max |u - Ghia| = {max_profile_error:.4} over the centreline profile, \
         above the 1% of lid speed this resolution warrants"
    );
    Ok(())
}