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

//! A-P2, step S2 (`docs/overset_metal_campaign.md` §5.9): the overlap map
//! alone. Classification of the background against the P2 MMS geometry
//! (unit square, phantom circle at (0.6, 0.45) r = 0.2, annulus to
//! r = 0.32) with the donor invariants; linear fields reproduced to
//! rounding in both directions; quadratic fields at second order under
//! refinement; the inverse bilinear on every dual quad of the skewed annulus.

use rtx_cfd::CfdResult;
use rtx_cfd::mesh::PatchMesh;
use rtx_cfd::mesh::patch_gen::annulus_skewed;
use rtx_cfd::solvers::incompressible::overset::overlap::{DEFAULT_OVERLAP_ROWS, inverse_bilinear};
use rtx_cfd::solvers::incompressible::{CellClass, FlowField, OverlapMap};

const CX: f64 = 0.6;
const CY: f64 = 0.45;
const R0: f64 = 0.2;
/// Chosen so the OUTER patch spacing is exactly h at `nn = n/4` with the
/// 3× geometric stretch: span = h · total(nn)/3 = 4.92 h at n = 32.
const R1: f64 = 0.354;

/// The P2 patch at background resolution `n`: `ns = 9n/4` (arc ≈ h at the
/// outer ring), `nn = n/4`, skew 0.3, stretch 3.
fn patch(n: usize) -> CfdResult<PatchMesh> {
    annulus_skewed([CX, CY], R0, R1, 9 * n / 4, n / 4, 0.3, 3.0)
}

fn overlap(n: usize) -> CfdResult<(OverlapMap, PatchMesh)> {
    let p = patch(n)?;
    let h = 1.0 / n as f64;
    let map = OverlapMap::build(&p, n, n, h, h, DEFAULT_OVERLAP_ROWS)?;
    Ok((map, p))
}

fn lattice_field(n: usize, f: &dyn Fn(f64, f64) -> f64) -> FlowField {
    let h = 1.0 / n as f64;
    let mut field = FlowField::new(n, n, h, h).unwrap();
    for j in 0..n {
        for i in 0..=n {
            field.u[(j, i)] = f(i as f64 * h, (j as f64 + 0.5) * h);
        }
    }
    for j in 0..=n {
        for i in 0..n {
            field.v[(j, i)] = f((i as f64 + 0.5) * h, j as f64 * h);
        }
    }
    for j in 0..n {
        for i in 0..n {
            field.p[(j, i)] = f((i as f64 + 0.5) * h, (j as f64 + 0.5) * h);
        }
    }
    field
}

fn patch_field(p: &PatchMesh, f: &dyn Fn(f64, f64) -> f64) -> Vec<f64> {
    (0..p.cell_count())
        .map(|c| {
            let xy = p.centre(c);
            f(xy[0], xy[1])
        })
        .collect()
}

/// Largest interpolation error of `f` at the fringe cells, fringe faces and
/// acceptors.
fn interpolation_errors(n: usize, f: &dyn Fn(f64, f64) -> f64) -> CfdResult<(f64, f64)> {
    let (map, p) = overlap(n)?;
    let h = 1.0 / n as f64;
    let vals = patch_field(&p, f);
    let cells = map.fringe_cell_values(&vals);
    let mut worst_fringe = 0.0_f64;
    for (e, v) in map.fringe_cells.iter().zip(&cells) {
        let exact = f((e.i as f64 + 0.5) * h, (e.j as f64 + 0.5) * h);
        worst_fringe = worst_fringe.max((v - exact).abs());
    }
    let mut field = FlowField::new(n, n, h, h).unwrap();
    map.stamp_fringe_faces(&mut field, &vals, &vals);
    for e in &map.fringe_u {
        let exact = f(e.i as f64 * h, (e.j as f64 + 0.5) * h);
        worst_fringe = worst_fringe.max((field.u[(e.j, e.i)] - exact).abs());
    }
    for e in &map.fringe_v {
        let exact = f((e.i as f64 + 0.5) * h, e.j as f64 * h);
        worst_fringe = worst_fringe.max((field.v[(e.j, e.i)] - exact).abs());
    }
    let bg = lattice_field(n, f);
    let acc = map.acceptor_values(&bg, &bg.p);
    let mut worst_acc = 0.0_f64;
    for (a, (u, v, pp)) in map.acceptors.iter().zip(&acc) {
        let xy = p.centre(a.cell);
        let exact = f(xy[0], xy[1]);
        worst_acc = worst_acc.max(
            (u - exact)
                .abs()
                .max((v - exact).abs())
                .max((pp - exact).abs()),
        );
    }
    Ok((worst_fringe, worst_acc))
}

#[test]
fn classification_on_the_mms_geometry_holds_the_donor_invariants() -> CfdResult<()> {
    for n in [32usize, 64, 128] {
        let (map, p) = overlap(n)?;
        let h = 1.0 / n as f64;
        let mut counts = [0usize; 3];
        for j in 0..n {
            for i in 0..n {
                counts[match map.class(j, i) {
                    CellClass::Active => 0,
                    CellClass::Fringe => 1,
                    CellClass::Hole => 2,
                }] += 1;
            }
        }
        println!(
            "  n = {n} (patch {}x{}): active {}, fringe {}, hole {}; prescribed u {} v {}; acceptors {}; donor rows {:?}",
            p.ns(),
            p.nn(),
            counts[0],
            counts[1],
            counts[2],
            map.fringe_u.len(),
            map.fringe_v.len(),
            map.acceptors.len(),
            map.donor_rows
        );
        // The hole covers the body and the inner part of the patch: more
        // than the circle's area, less than the patch's outer disc.
        let circle = std::f64::consts::PI * R0 * R0 / (h * h);
        let disc = std::f64::consts::PI * R1 * R1 / (h * h);
        assert!(
            (counts[2] as f64) > circle && (counts[2] as f64) < disc,
            "hole count {} outside ({circle:.0}, {disc:.0})",
            counts[2]
        );
        assert!(counts[1] > 0 && map.acceptors.len() == p.ns());
        // Every fringe cell has all its non-active-side faces prescribed.
        for e in &map.fringe_cells {
            let (j, i) = (e.j, e.i);
            for (jj, ii, is_u) in [
                (j, i, true),
                (j, i + 1, true),
                (j, i, false),
                (j + 1, i, false),
            ] {
                let other_active = if is_u {
                    let left = ii > 0 && map.class(jj, ii - 1) == CellClass::Active;
                    let right = ii < n && map.class(jj, ii) == CellClass::Active;
                    left || right
                } else {
                    let below = jj > 0 && map.class(jj - 1, ii) == CellClass::Active;
                    let above = jj < n && map.class(jj, ii) == CellClass::Active;
                    below || above
                };
                if !other_active {
                    let listed = if is_u {
                        map.fringe_u.iter().any(|f| f.j == jj && f.i == ii)
                    } else {
                        map.fringe_v.iter().any(|f| f.j == jj && f.i == ii)
                    };
                    assert!(
                        listed,
                        "fringe cell ({j}, {i}) face ({jj}, {ii}, u = {is_u}) not prescribed"
                    );
                }
            }
        }
    }
    Ok(())
}

#[test]
fn linear_fields_are_reproduced_to_rounding_in_both_directions() -> CfdResult<()> {
    let lin = |x: f64, y: f64| 0.3 + 1.7 * x - 0.9 * y;
    for n in [32usize, 64] {
        let (fringe, acc) = interpolation_errors(n, &lin)?;
        println!("  n = {n}: linear field — fringe {fringe:.3e}, acceptors {acc:.3e}");
        assert!(
            fringe < 1e-13,
            "fringe interpolation not linear-exact: {fringe:.3e}"
        );
        assert!(
            acc < 1e-13,
            "acceptor interpolation not linear-exact: {acc:.3e}"
        );
    }
    Ok(())
}

#[test]
fn quadratic_fields_are_interpolated_at_second_order() -> CfdResult<()> {
    let quad = |x: f64, y: f64| x * x + x * y - y * y;
    let mut fringe = Vec::new();
    let mut acc = Vec::new();
    for n in [32usize, 64, 128] {
        let (f, a) = interpolation_errors(n, &quad)?;
        println!("  n = {n}: quadratic field — fringe {f:.3e}, acceptors {a:.3e}");
        fringe.push(f);
        acc.push(a);
    }
    let order = |e: &[f64]| -> Vec<f64> { e.windows(2).map(|w| (w[0] / w[1]).log2()).collect() };
    let (of, oa) = (order(&fringe), order(&acc));
    println!("  orders: fringe {of:?}, acceptors {oa:?}");
    // Measured: fringe 1.40 / 1.91, acceptors 1.96 / 1.99. The fringe's
    // first rung is pre-asymptotic — at nn = 8 the dual quads of the 3×
    // stretched annulus are markedly non-parallelogram (per-row ratio
    // 1.17 against 1.08 at nn = 16), which sets the bilinear map's error
    // constant; the asymptotic pair is second order.
    assert!(
        of.iter().all(|&o| o > 1.0) && of.last().is_some_and(|&o| o > 1.8),
        "fringe orders {of:?}"
    );
    assert!(oa.iter().all(|&o| o > 1.8), "acceptor orders {oa:?}");
    Ok(())
}

#[test]
fn inverse_bilinear_converges_on_every_dual_quad_of_the_skewed_annulus() -> CfdResult<()> {
    let p = patch(64)?;
    let (ns, nn) = (p.ns(), p.nn());
    let mut worst = 0.0_f64;
    for k in 0..nn - 1 {
        for i in 0..ns {
            let i1 = (i + 1) % ns;
            let q = [
                p.centre(p.cell(k, i)),
                p.centre(p.cell(k, i1)),
                p.centre(p.cell(k + 1, i1)),
                p.centre(p.cell(k + 1, i)),
            ];
            for (s, t) in [(0.5, 0.5), (0.1, 0.9), (0.85, 0.2), (0.01, 0.01)] {
                let n = [(1.0 - s) * (1.0 - t), s * (1.0 - t), s * t, (1.0 - s) * t];
                let x = (0..4).map(|a| n[a] * q[a][0]).sum::<f64>();
                let y = (0..4).map(|a| n[a] * q[a][1]).sum::<f64>();
                let w = inverse_bilinear(&q, x, y).unwrap_or_else(|| {
                    panic!("no convergence in dual quad ({k}, {i}) at ({s}, {t})")
                });
                for a in 0..4 {
                    worst = worst.max((w[a] - n[a]).abs());
                }
            }
        }
    }
    println!("  inverse bilinear: worst weight error over every dual quad {worst:.3e}");
    assert!(worst < 1e-12);
    Ok(())
}

/// The classification must follow a translating patch: cells flip as the
/// hole boundary sweeps over their centres.
#[test]
fn classification_follows_a_translating_patch() -> CfdResult<()> {
    let n = 32;
    let h = 1.0 / n as f64;
    let at = |cx: f64| -> CfdResult<OverlapMap> {
        let p = annulus_skewed([cx, CY], R0, R1, 9 * n / 4, n / 4, 0.3, 3.0)?;
        OverlapMap::build(&p, n, n, h, h, DEFAULT_OVERLAP_ROWS)
    };
    let base = at(0.6)?;
    for cx in [0.6 - 0.25 * h, 0.6 - 0.8 * h, 0.6 - 2.0 * h] {
        let moved = at(cx)?;
        let changed = (0..n)
            .flat_map(|j| (0..n).map(move |i| (j, i)))
            .filter(|&(j, i)| base.class(j, i) != moved.class(j, i))
            .count();
        println!(
            "  patch translated by {:.2} h: {changed} background cells reclassified",
            (0.6 - cx) / h
        );
        assert!(
            changed > 0,
            "no cell reclassified after a {:.2} h translation",
            (0.6 - cx) / h
        );
    }
    Ok(())
}