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

//! A-P4-0 gate 3 (`docs/overset_metal_campaign.md` §5.11): the P0
//! manufactured solution on the cylinder–flag O-grid (exact acceptors, the
//! S3 harness) — do the Stokes-limit and upwind orders survive the junction
//! fillets' skew? Rungs at the benchmark's h = 0.41 / 41, 62, 82.

use rtx_cfd::mesh::PatchMesh;
use rtx_cfd::mesh::patch_gen::{FlagEdges, cylinder_flag_patch, cylinder_flag_patch_deformed};
use rtx_cfd::solvers::incompressible::{
    CurvilinearParameters, CurvilinearPisoSolver, NormalDiffusion, PatchConvection, PatchField,
};
use rtx_cfd::{CfdConfig, CfdResult};
use std::f64::consts::PI;

const RHO: f64 = 1.0;
const MU: f64 = 0.05;

fn u_exact(x: f64, y: f64) -> f64 {
    (PI * x).sin() * (PI * y).cos()
}
fn v_exact(x: f64, y: f64) -> f64 {
    -(PI * x).cos() * (PI * y).sin()
}
fn p_exact(x: f64, y: f64) -> f64 {
    (PI * x).sin() * (PI * y).sin()
}
fn source(x: f64, y: f64, convecting: bool) -> (f64, f64) {
    let conv = if convecting { RHO * 0.5 * PI } else { 0.0 };
    (
        conv * (2.0 * PI * x).sin()
            + 2.0 * PI * PI * MU * u_exact(x, y)
            + PI * (PI * x).cos() * (PI * y).sin(),
        conv * (2.0 * PI * y).sin()
            + 2.0 * PI * PI * MU * v_exact(x, y)
            + PI * (PI * x).sin() * (PI * y).cos(),
    )
}

fn patch(ny: usize) -> CfdResult<PatchMesh> {
    let h = 0.41 / ny as f64;
    let sweeps: usize = std::env::var("RTX_CF_SWEEPS")
        .ok()
        .and_then(|v| v.parse().ok())
        .unwrap_or(500);
    // Fixed geometry across the ladder: the fillet of the coarsest rung.
    let fillet = 0.5 * 0.41 / 41.0;
    // P5-0 gate (iv), `RTX_CF_BEND=a`: the same MMS on the patch around
    // the flag bent to tip deflection `a` (the cantilever shape).
    if let Some(a) = std::env::var("RTX_CF_BEND")
        .ok()
        .and_then(|v| v.parse::<f64>().ok())
    {
        return Ok(cylinder_flag_patch_deformed(
            [0.2, 0.2],
            0.05,
            0.01,
            &bent_edges(a, 35, 2),
            0.6,
            h,
            fillet,
            6.0 * h,
            12,
            4.0,
            sweeps,
        )?
        .0);
    }
    Ok(cylinder_flag_patch(
        [0.2, 0.2],
        0.05,
        0.01,
        0.6,
        h,
        fillet,
        6.0 * h,
        12,
        4.0,
        sweeps,
    )?
    .0)
}

/// The FEA flag's wetted edges under the cantilever end-load shape with
/// tip deflection `a` (as `patch_cylinder_flag_deformed.rs`).
fn bent_edges(a: f64, nx: usize, ny: usize) -> FlagEdges {
    let (x0, x1, t, cy) = (0.25, 0.6, 0.01, 0.2);
    let len = x1 - x0;
    let centre = |x: f64| {
        let xi = (x - x0) / len;
        (
            a * xi * xi * (3.0 - xi) / 2.0,
            a * (6.0 * xi - 3.0 * xi * xi) / 2.0 / len,
        )
    };
    let edge = |x: f64, side: f64| -> [f64; 2] {
        let (y, dy) = centre(x);
        let n = (1.0 + dy * dy).sqrt();
        [x - side * t * dy / n, cy + y + side * t / n]
    };
    let m = 2 * nx;
    let bottom: Vec<[f64; 2]> = (0..=m)
        .map(|i| edge(x0 + len * i as f64 / m as f64, -1.0))
        .collect();
    let top: Vec<[f64; 2]> = (0..=m)
        .rev()
        .map(|i| edge(x0 + len * i as f64 / m as f64, 1.0))
        .collect();
    let (b, tp) = (bottom[m], top[0]);
    let k = 2 * ny;
    let tip: Vec<[f64; 2]> = (0..=k)
        .map(|j| {
            let f = j as f64 / k as f64;
            [b[0] + f * (tp[0] - b[0]), b[1] + f * (tp[1] - b[1])]
        })
        .collect();
    FlagEdges { bottom, tip, top }
}

async fn march(ny: usize, convection: PatchConvection) -> CfdResult<(f64, usize, f64)> {
    let mesh = patch(ny)?;
    let h = 0.41 / ny as f64;
    let nu = MU / RHO;
    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 dt = 0.4 * (hs * hs / (4.0 * nu)).min(h);
    let config = CfdConfig::new()
        .with_density(RHO)
        .with_viscosity(MU)
        .with_reference_velocity(1.0)
        .with_reference_length(1.0);
    let convecting = convection != PatchConvection::None;
    let mut solver = CurvilinearPisoSolver::new(
        config,
        CurvilinearParameters {
            tolerance: 1e-5,
            convection,
            normal_diffusion: NormalDiffusion::LineImplicit,
            ..CurvilinearParameters::default()
        },
        mesh,
    )?;
    solver.set_boundary_velocity(|x, y, _| (u_exact(x, y), v_exact(x, y)));
    solver.set_momentum_source(move |x, y, _| source(x, y, convecting));
    solver.set_acceptor_ring(true);
    let (ns, nn) = (solver.mesh().ns(), solver.mesh().nn());
    let acc: Vec<(f64, f64, f64)> = (0..ns)
        .map(|i| {
            let xy = solver.mesh().centre(solver.mesh().cell(nn - 1, i));
            (
                u_exact(xy[0], xy[1]),
                v_exact(xy[0], xy[1]),
                p_exact(xy[0], xy[1]),
            )
        })
        .collect();
    let zeros = vec![0.0; ns];
    let mut field = PatchField::new(solver.mesh());
    solver.initialize(&mut field, |_, _| (0.0, 0.0));
    solver.stamp_acceptors(&mut field, &acc);
    solver.set_acceptor_correction(&zeros);
    let steady_tol = if convecting { 1e-6 } else { 1e-7 };
    let mut steady = f64::INFINITY;
    let mut steps = 0;
    let mut max_div = 0.0_f64;
    for _ in 0..600_000 {
        let before = (field.u.clone(), field.v.clone());
        let r = solver.advance(&mut field, dt).await?;
        assert!(
            r.poisson_converged,
            "pressure solve did not converge: {r:?}"
        );
        solver.stamp_acceptors(&mut field, &acc);
        let flux_scale: f64 = field.flux.iter().map(|f| f.abs()).sum::<f64>().max(1e-300);
        max_div = max_div.max(r.max_divergence / flux_scale);
        steps += 1;
        let change = field
            .u
            .iter()
            .zip(&before.0)
            .chain(field.v.iter().zip(&before.1))
            .map(|(a, b)| (a - b).abs())
            .fold(0.0, f64::max);
        steady = change / dt;
        if steady < steady_tol {
            break;
        }
    }
    assert!(steady < steady_tol, "no steady state: {steady:.3e}");
    let mesh = solver.mesh();
    let (mut sq, mut vol) = (0.0, 0.0);
    for c in 0..mesh.cell_count() {
        if solver.is_acceptor(c) {
            continue;
        }
        let xy = mesh.centre(c);
        let eu = field.u[c] - u_exact(xy[0], xy[1]);
        let ev = field.v[c] - v_exact(xy[0], xy[1]);
        sq += (eu * eu + ev * ev) * mesh.area(c);
        vol += mesh.area(c);
    }
    Ok(((sq / vol).sqrt(), steps, max_div))
}

#[tokio::test]
async fn cylinder_flag_patch_keeps_the_p0_orders() -> CfdResult<()> {
    // Measured with the fillet fixed at 5 mm: Stokes 2.17 / 2.13; upwind 2.00
    // / 1.86 — the patch's cells are so small against the field (cell Péclet
    // ≈ 0.1 at ν = 0.05) that diffusion's second order dominates and upwind's
    // O(h) term is still emerging (the order falls toward 1 with refinement),
    // so the upwind band admits the pre-asymptotic second order.
    // P4 step 2 gate (iii): TVD (van Albada) below upwind at every rung,
    // orders in [1.5, 2.6]. `RTX_CF_SCHEME=none|upwind|tvd` runs one scheme.
    let only = std::env::var("RTX_CF_SCHEME").ok();
    for (convection, gate) in [
        (PatchConvection::None, 1.8..2.6),
        (PatchConvection::Upwind, 0.7..2.4),
        (PatchConvection::TvdVanAlbada, 1.5..2.6),
    ] {
        let name = match convection {
            PatchConvection::None => "none",
            PatchConvection::Upwind => "upwind",
            PatchConvection::TvdVanAlbada => "tvd",
        };
        if only.as_deref().is_some_and(|o| o != name) {
            continue;
        }
        let upwind = [5.134660e-4, 2.244819e-4, 1.336112e-4];
        let mut errs = Vec::new();
        let mut hs = Vec::new();
        for (&ny, &u) in [41usize, 62, 82].iter().zip(&upwind) {
            let (l2, steps, max_div) = march(ny, convection).await?;
            println!(
                "  cylinder-flag {convection:?} ny={ny}: L2 {l2:.6e}, max div {max_div:.2e}, {steps} steps"
            );
            if convection == PatchConvection::TvdVanAlbada {
                println!("    tvd / upwind at ny={ny}: {:.3}", l2 / u);
            }
            errs.push(l2);
            hs.push(0.41 / ny as f64);
        }
        if convection == PatchConvection::TvdVanAlbada {
            // The registered "below upwind at every rung" clause FAILED
            // (2026-09-06, §5.11): ratios 1.061 / 1.069 / 1.035 with orders
            // 1.98 / 1.97. At cell Péclet ≈ 0.1 the Stokes floor (4.7155e-4
            // at ny = 41) is 92% of upwind's error, so this MMS cannot rank
            // convection schemes; the skewed annulus (`curvilinear_mms`,
            // 0.24× upwind at ns = 128) is the discriminating gate. The
            // ratios are recorded here, the order band is the assertion.
            let ratios: Vec<f64> = errs.iter().zip(&upwind).map(|(e, u)| e / u).collect();
            println!("  cylinder-flag TvdVanAlbada / upwind ratios {ratios:?}");
        }
        let o: Vec<f64> = errs
            .windows(2)
            .zip(hs.windows(2))
            .map(|(e, h)| (e[0] / e[1]).ln() / (h[0] / h[1]).ln())
            .collect();
        println!("  cylinder-flag {convection:?} orders {o:?}");
        assert!(
            o.iter().all(|x| gate.contains(x)),
            "{convection:?} orders {o:?} outside {gate:?}"
        );
    }
    Ok(())
}