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

//! P0 step 1: the patch mesh's metrics (`docs/overset_metal_campaign.md`
//! §5.3). A Cartesian patch reproduces `dx, dy` exactly; every cell of a
//! skewed periodic annulus is closed (`Σ sign·S_f = 0`) to rounding, seam
//! included; the annulus area converges at second order; a folded patch
//! is refused.

use rtx_cfd::mesh::patch_gen::{annulus_skewed, cartesian, channel_sheared};
use rtx_cfd::mesh::{PatchMesh, PatchSide};
use std::f64::consts::PI;

fn closure_defect(mesh: &PatchMesh) -> f64 {
    let mut worst = 0.0_f64;
    for c in 0..mesh.cell_count() {
        let mut sum = [0.0, 0.0];
        for (f, sign) in mesh.cell_faces(c) {
            let s = mesh.faces()[f].s;
            sum[0] += sign * s[0];
            sum[1] += sign * s[1];
        }
        worst = worst.max(sum[0].abs()).max(sum[1].abs());
    }
    worst
}

#[test]
fn cartesian_patch_metrics_are_exact() {
    let (nx, ny, lx, ly) = (8, 5, 2.0, 1.0);
    let mesh = cartesian(nx, ny, lx, ly, false).unwrap();
    let (dx, dy) = (lx / nx as f64, ly / ny as f64);
    for c in 0..mesh.cell_count() {
        assert!((mesh.area(c) - dx * dy).abs() < 1e-15);
        let (k, i) = mesh.cell_ki(c);
        assert!((mesh.centre(c)[0] - (i as f64 + 0.5) * dx).abs() < 1e-14);
        assert!((mesh.centre(c)[1] - (k as f64 + 0.5) * dy).abs() < 1e-14);
    }
    for (f, face) in mesh.faces().iter().enumerate() {
        let expect_s = if mesh.is_sface(f) {
            [dy, 0.0]
        } else {
            [0.0, dx]
        };
        assert!((face.s[0] - expect_s[0]).abs() < 1e-15 && (face.s[1] - expect_s[1]).abs() < 1e-15);
        if face.owner.is_some() && face.neigh.is_some() {
            assert!((face.w - 0.5).abs() < 1e-14);
            let expect = if mesh.is_sface(f) {
                [dx, 0.0]
            } else {
                [0.0, dy]
            };
            assert!((face.d[0] - expect[0]).abs() < 1e-14 && (face.d[1] - expect[1]).abs() < 1e-14);
        }
    }
    assert_eq!(closure_defect(&mesh), 0.0);
    assert_eq!(mesh.faces().len(), ny * (nx + 1) + (ny + 1) * nx);
    let sides: Vec<_> = (0..mesh.faces().len())
        .filter_map(|f| mesh.side(f))
        .collect();
    assert_eq!(sides.iter().filter(|&&s| s == PatchSide::Inner).count(), nx);
    assert_eq!(sides.iter().filter(|&&s| s == PatchSide::Outer).count(), nx);
    assert_eq!(
        sides.iter().filter(|&&s| s == PatchSide::SStart).count(),
        ny
    );
    assert_eq!(sides.iter().filter(|&&s| s == PatchSide::SEnd).count(), ny);
    mesh.validate(60.0).unwrap();
}

#[test]
fn skewed_annulus_cells_are_closed_seam_included() {
    let mesh = annulus_skewed([0.3, -0.2], 0.5, 1.5, 24, 6, 0.3, 3.0).unwrap();
    assert!(mesh.periodic().is_some());
    assert_eq!(mesh.faces().len(), 6 * 24 + 7 * 24);
    let defect = closure_defect(&mesh);
    assert!(defect < 1e-15, "closure defect {defect:.3e}");
    // The seam face is shared: cells ns-1 and 0 of every row see it with
    // opposite signs and no other face joins them.
    for k in 0..mesh.nn() {
        let west_of_0 = mesh.cell_faces(mesh.cell(k, 0))[0];
        let east_of_last = mesh.cell_faces(mesh.cell(k, mesh.ns() - 1))[1];
        assert_eq!(west_of_0.0, east_of_last.0);
        assert_eq!(west_of_0.1, -east_of_last.1);
        let f = &mesh.faces()[west_of_0.0];
        assert_eq!(f.owner, Some(mesh.cell(k, mesh.ns() - 1)));
        assert_eq!(f.neigh, Some(mesh.cell(k, 0)));
        assert!(f.w > 0.0 && f.w < 1.0);
    }
    // No boundary faces along s; nn faces per column across.
    assert!(
        !(0..mesh.faces().len())
            .any(|f| { matches!(mesh.side(f), Some(PatchSide::SStart | PatchSide::SEnd)) })
    );
    // Seam node cells wrap.
    let cells = mesh.node_cells(3, 0);
    assert_eq!(cells.len(), 4);
    assert!(cells.contains(&mesh.cell(2, mesh.ns() - 1)) && cells.contains(&mesh.cell(3, 0)));
    // Stretching: the across-patch spacing grows by the requested factor.
    let eta = rtx_cfd::mesh::patch_gen::stretched_fractions(6, 3.0);
    let ratio = (eta[6] - eta[5]) / (eta[1] - eta[0]);
    assert!((ratio - 3.0).abs() < 1e-12, "stretch ratio {ratio}");
    assert!((eta[6] - 1.0).abs() < 1e-15 && eta[0] == 0.0);
    mesh.validate(60.0).unwrap();
}

#[test]
fn annulus_area_converges_at_second_order() {
    let exact = PI * (1.5_f64.powi(2) - 0.5_f64.powi(2));
    let errs: Vec<f64> = [16usize, 32, 64]
        .iter()
        .map(|&ns| {
            let mesh = annulus_skewed([0.0, 0.0], 0.5, 1.5, ns, ns / 4, 0.0, 1.0).unwrap();
            let area: f64 = (0..mesh.cell_count()).map(|c| mesh.area(c)).sum();
            (area - exact).abs()
        })
        .collect();
    let orders: Vec<f64> = errs.windows(2).map(|p| (p[0] / p[1]).log2()).collect();
    println!("annulus area errors {errs:?}, orders {orders:?}");
    assert!(
        orders.iter().all(|&o| o > 1.9 && o < 2.1),
        "orders {orders:?}"
    );
}

#[test]
fn sheared_channel_is_periodic_with_a_shift_and_congruent() {
    let mesh = channel_sheared(2.0, 1.0, 8, 4, 0.4, true).unwrap();
    assert_eq!(mesh.periodic(), Some([2.0, 0.0]));
    let a0 = mesh.area(0);
    assert!((0..mesh.cell_count()).all(|c| (mesh.area(c) - a0).abs() < 1e-14));
    // n-faces stay horizontal; s-faces are tilted.
    for (f, face) in mesh.faces().iter().enumerate() {
        if mesh.is_sface(f) {
            assert!(face.s[1] != 0.0);
        } else {
            assert_eq!(face.s[0], 0.0);
        }
    }
    // The seam face's d spans one cell width (with the shift applied).
    let seam = &mesh.faces()[mesh.sface(1, 0)];
    assert!(
        (seam.d[0] - 0.25).abs() < 1e-14 && seam.d[1].abs() < 1e-14,
        "seam d {:?}",
        seam.d
    );
    assert!(closure_defect(&mesh) < 1e-15);
    mesh.validate(60.0).unwrap();
}

#[test]
fn a_folded_patch_is_refused() {
    let mut inner = Vec::new();
    let mut outer = Vec::new();
    for i in 0..=8 {
        let th = -2.0 * PI * (i % 8) as f64 / 8.0; // clockwise, as the generator does
        inner.push([th.cos(), th.sin()]);
        outer.push([0.5 * th.cos(), 0.5 * th.sin()]); // inside the inner ring: inverted cells
    }
    let err = rtx_cfd::mesh::patch_gen::transfinite(&inner, &outer, 3, 1.0, Some([0.0, 0.0]));
    assert!(err.is_err());
    // A mismatched periodic column is refused too.
    let bad = PatchMesh::from_nodes(
        4,
        2,
        (0..15).map(|i| i as f64).collect(),
        vec![0.0; 15],
        Some([0.0, 0.0]),
    );
    assert!(bad.is_err());
}