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

//! P1 gates (`docs/overset_metal_campaign.md` §2.2 P1, `next_session.md`
//! 2026-09-04): the curvilinear patch moves and deforms under an exact
//! discrete geometric conservation law.
//!
//! 1. The 2-D identity `Σ_f sign_f S̄_f·δc_f = V^{n+1} − V^n` holds to
//!    rounding on general quadrilaterals under the trapezoidal face rule
//!    and fails visibly under end-of-step areas.
//! 2. Uniform flow on an arbitrarily wiggling AND bending annulus (whole-
//!    patch translation, rotation, quadratic bending, short-wave interior
//!    wiggle) is preserved to ≤ 1e-12 over 400 steps at the `ale_dgcl.rs`
//!    standard (pressure stop 1e-13); the end-of-step negative control
//!    must FAIL visibly.
//! 3. A stationary mesh pushed through the moving path is bit-identical to
//!    the static path (both diffusion variants).
//! 4. Snapshot/restore on a moving mesh re-runs bit-identically.
//! 5. Taylor–Green order is unchanged under a general (non-tensor) interior
//!    mesh motion (`ale_taylor_green.rs` pattern).

use rtx_cfd::mesh::PatchMesh;
use rtx_cfd::mesh::patch_gen::{annulus_skewed, cartesian};
use rtx_cfd::solvers::incompressible::{
    CurvilinearParameters, CurvilinearPisoSolver, NormalDiffusion, PatchConvection, PatchField,
    StepGeometry, SweptFaceRule,
};
use rtx_cfd::{CfdConfig, CfdResult};
use std::f64::consts::PI;

const U_UNIFORM: f64 = 0.7;
const V_UNIFORM: f64 = -0.4;

/// The annulus every moving-mesh test starts from (the P0 skewed,
/// stretched, periodic O-grid).
fn base_annulus(ns: usize, nn: usize) -> CfdResult<PatchMesh> {
    annulus_skewed([0.0, 0.0], 0.5, 1.5, ns, nn, 0.3, 3.0)
}

/// Arbitrary smooth motion of the whole annulus: quadratic bending, a
/// short-wavelength interior wiggle, a rigid rotation and a translation —
/// nothing is fixed, no symmetry survives, and the displacement gradient
/// stays below one so no cell folds. A function of the reference position
/// only, so the periodic seam (column `ns` a bitwise copy of column 0)
/// stays exact.
fn deform(x: f64, y: f64, t: f64) -> [f64; 2] {
    let x1 = x + 0.08 * (1.7 * t).sin() * y * y;
    let y1 = y + 0.06 * (2.3 * t + 1.0).sin() * x * x;
    let x2 = x1 + 0.03 * (4.0 * y1 + 2.9 * t).sin();
    let y2 = y1 + 0.03 * (4.0 * x1 + 4.3 * t + 1.0).sin();
    let th = 0.3 * (1.1 * t).sin();
    let (c, s) = (th.cos(), th.sin());
    [
        c * x2 - s * y2 + 0.2 * (0.9 * t).sin(),
        s * x2 + c * y2 + 0.15 * ((1.3 * t).cos() - 1.0),
    ]
}

/// `base` with every node moved by `motion(x_ref, y_ref, t)`.
fn moved(
    base: &PatchMesh,
    t: f64,
    motion: &dyn Fn(f64, f64, f64) -> [f64; 2],
) -> CfdResult<PatchMesh> {
    let (ns, nn) = (base.ns(), base.nn());
    let count = (ns + 1) * (nn + 1);
    let mut x = Vec::with_capacity(count);
    let mut y = Vec::with_capacity(count);
    for n in 0..count {
        let r = base.node_xy(n);
        let m = motion(r[0], r[1], t);
        x.push(m[0]);
        y.push(m[1]);
    }
    PatchMesh::from_nodes(ns, nn, x, y, base.periodic())
}

fn max_abs_diff(a: &[f64], b: &[f64]) -> f64 {
    a.iter()
        .zip(b)
        .map(|(x, y)| (x - y).abs())
        .fold(0.0, f64::max)
}

fn fields_identical(a: &PatchField, b: &PatchField) -> bool {
    let bits = |v: &[f64], w: &[f64]| v.iter().zip(w).all(|(x, y)| x.to_bits() == y.to_bits());
    bits(&a.u, &b.u) && bits(&a.v, &b.v) && bits(&a.p, &b.p) && bits(&a.flux, &b.flux)
}

// ---------------------------------------------------------------- gate 1

#[test]
fn swept_volumes_close_the_gcl_identity_on_general_quads() -> CfdResult<()> {
    let base = base_annulus(32, 8)?;
    let mut worst_trap = 0.0_f64;
    let mut worst_end = 0.0_f64;
    for step in 0..40 {
        let (t0, t1) = (0.05 * step as f64, 0.05 * (step + 1) as f64);
        let old = moved(&base, t0, &deform)?;
        let new = moved(&base, t1, &deform)?;
        let trap = StepGeometry::new(&old, &new, SweptFaceRule::Trapezoidal);
        let end = StepGeometry::new(&old, &new, SweptFaceRule::EndOfStep);
        worst_trap = worst_trap.max(trap.gcl_defect(&old, &new));
        worst_end = worst_end.max(end.gcl_defect(&old, &new));
    }
    println!(
        "  GCL defect (relative to cell volume): trapezoidal {worst_trap:.3e}, end-of-step {worst_end:.3e}"
    );
    assert!(
        worst_trap < 1e-13,
        "trapezoidal swept volumes do not close the volume increment: {worst_trap:.3e}"
    );
    assert!(
        worst_end > 1e-4,
        "the end-of-step control only reaches {worst_end:.3e}: motion too tame"
    );
    Ok(())
}

// ------------------------------------------------------------- gates 2, 3

async fn uniform_flow_deviation(
    rule: SweptFaceRule,
    steps: usize,
    dt: f64,
) -> CfdResult<(f64, f64, usize)> {
    let config = CfdConfig::new()
        .with_density(1.0)
        .with_viscosity(0.05)
        .with_reference_velocity(1.0)
        .with_reference_length(1.0);
    // Pressure stop at the rounding floor (the ale_dgcl.rs lesson): with an
    // engineering tolerance a partial p' accumulates and re-perturbs u.
    let params = CurvilinearParameters {
        corrector_steps: 2,
        tolerance: 1e-13,
        swept_face_rule: rule,
        ..CurvilinearParameters::default()
    };
    let base = base_annulus(32, 8)?;
    let mut solver = CurvilinearPisoSolver::new(config, params, moved(&base, 0.0, &deform)?)?;
    solver.set_boundary_velocity(|_x, _y, _t| (U_UNIFORM, V_UNIFORM));
    let mut field = PatchField::new(solver.mesh());
    solver.initialize(&mut field, |_, _| (U_UNIFORM, V_UNIFORM));

    let mut worst = 0.0_f64;
    let mut iterations = 0;
    for step in 0..steps {
        let t_new = (step + 1) as f64 * dt;
        solver.set_mesh(moved(&base, t_new, &deform)?)?;
        let r = solver.advance(&mut field, dt).await?;
        iterations += r.poisson_iterations;
        for &u in &field.u {
            worst = worst.max((u - U_UNIFORM).abs());
        }
        for &v in &field.v {
            worst = worst.max((v - V_UNIFORM).abs());
        }
    }
    let max_p = field.p.iter().fold(0.0_f64, |m, &p| m.max(p.abs()));
    Ok((worst, max_p, iterations))
}

#[tokio::test]
async fn uniform_flow_stays_exactly_uniform_on_a_wiggling_and_bending_annulus() -> CfdResult<()> {
    let (worst, max_p, iters) =
        uniform_flow_deviation(SweptFaceRule::Trapezoidal, 400, 1e-3).await?;
    println!(
        "  DGCL on the annulus: max |u - U| over 400 steps = {worst:.3e}, final max |p| = {max_p:.3e}, \
         pressure iterations {iters}"
    );
    let scale = U_UNIFORM.abs().max(V_UNIFORM.abs());
    assert!(
        worst <= 1e-12 * scale,
        "DGCL violated: uniform flow deviated by {worst:.3e} on the moving annulus"
    );
    assert!(max_p < 1e-11, "pressure moved off constant: {max_p:.3e}");
    Ok(())
}

#[tokio::test]
async fn end_of_step_face_areas_violate_the_gcl_visibly() -> CfdResult<()> {
    let (worst, max_p, _) = uniform_flow_deviation(SweptFaceRule::EndOfStep, 400, 1e-3).await?;
    println!(
        "  GCL-violating control: max |u - U| over 400 steps = {worst:.3e}, max |p| = {max_p:.3e}"
    );
    let scale = U_UNIFORM.abs().max(V_UNIFORM.abs());
    assert!(
        worst > 1e-6 * scale,
        "the GCL-violating rule deviated only {worst:.3e}: the DGCL test has lost its teeth"
    );
    Ok(())
}

// ------------------------------------------------------------- gates 4, 5

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

fn u_mms(x: f64, y: f64) -> f64 {
    (PI * x).sin() * (PI * y).cos()
}
fn v_mms(x: f64, y: f64) -> f64 {
    -(PI * x).cos() * (PI * y).sin()
}
fn source_mms(x: f64, y: f64) -> (f64, f64) {
    let conv = RHO * 0.5 * PI;
    (
        conv * (2.0 * PI * x).sin()
            + 2.0 * PI * PI * MU * u_mms(x, y)
            + PI * (PI * x).cos() * (PI * y).sin(),
        conv * (2.0 * PI * y).sin()
            + 2.0 * PI * PI * MU * v_mms(x, y)
            + PI * (PI * x).sin() * (PI * y).cos(),
    )
}

fn mms_solver(
    mesh: PatchMesh,
    diffusion: NormalDiffusion,
) -> CfdResult<(CurvilinearPisoSolver, PatchField)> {
    let config = CfdConfig::new().with_density(RHO).with_viscosity(MU);
    let params = CurvilinearParameters {
        normal_diffusion: diffusion,
        ..CurvilinearParameters::default()
    };
    let mut solver = CurvilinearPisoSolver::new(config, params, mesh)?;
    solver.set_boundary_velocity(|x, y, t| (u_mms(x, y) * (1.0 + 0.1 * t), v_mms(x, y)));
    solver.set_momentum_source(|x, y, _| source_mms(x, y));
    let mut field = PatchField::new(solver.mesh());
    solver.initialize(&mut field, |x, y| (u_mms(x, y), v_mms(x, y)));
    Ok((solver, field))
}

#[tokio::test]
async fn stationary_mesh_through_the_moving_path_is_bit_identical() -> CfdResult<()> {
    for diffusion in [NormalDiffusion::Explicit, NormalDiffusion::LineImplicit] {
        let mesh = base_annulus(24, 6)?;
        let (mut plain, mut f_plain) = mms_solver(mesh.clone(), diffusion)?;
        let (mut moving, mut f_moving) = mms_solver(mesh.clone(), diffusion)?;
        let dt = 1e-3;
        for _ in 0..20 {
            plain.advance(&mut f_plain, dt).await?;
            moving.set_mesh(mesh.clone())?;
            moving.advance(&mut f_moving, dt).await?;
        }
        assert_eq!(plain.time().to_bits(), moving.time().to_bits());
        assert!(
            fields_identical(&f_plain, &f_moving),
            "{diffusion:?}: the moving path with a stationary mesh differs from the static path by \
             u {:.3e} v {:.3e} p {:.3e} flux {:.3e}",
            max_abs_diff(&f_plain.u, &f_moving.u),
            max_abs_diff(&f_plain.v, &f_moving.v),
            max_abs_diff(&f_plain.p, &f_moving.p),
            max_abs_diff(&f_plain.flux, &f_moving.flux),
        );
        println!(
            "  {diffusion:?}: stationary mesh through the moving path — bit-identical over 20 steps"
        );
    }
    Ok(())
}

#[tokio::test]
async fn snapshot_restore_rerun_on_a_moving_mesh_is_bit_identical() -> CfdResult<()> {
    let base = base_annulus(24, 6)?;
    let (mut solver, mut field) =
        mms_solver(moved(&base, 0.0, &deform)?, NormalDiffusion::LineImplicit)?;
    let dt = 1e-3;
    let mut step = 0usize;
    for _ in 0..5 {
        step += 1;
        solver.set_mesh(moved(&base, step as f64 * dt, &deform)?)?;
        solver.advance(&mut field, dt).await?;
    }
    // Name the next mesh BEFORE the snapshot so the pending mesh is part of
    // what restore must bring back.
    solver.set_mesh(moved(&base, (step + 1) as f64 * dt, &deform)?)?;
    let saved = solver.snapshot();
    let field_saved = field.clone();
    let step_saved = step;
    for _ in 0..10 {
        step += 1;
        if step > step_saved + 1 {
            solver.set_mesh(moved(&base, step as f64 * dt, &deform)?)?;
        }
        solver.advance(&mut field, dt).await?;
    }
    let reference = (field.clone(), solver.time());

    solver.restore(&saved);
    field = field_saved;
    step = step_saved;
    for _ in 0..10 {
        step += 1;
        if step > step_saved + 1 {
            solver.set_mesh(moved(&base, step as f64 * dt, &deform)?)?;
        }
        solver.advance(&mut field, dt).await?;
    }
    assert_eq!(solver.time().to_bits(), reference.1.to_bits());
    assert!(
        fields_identical(&field, &reference.0),
        "re-run on the moving mesh differs: u {:.3e} v {:.3e} p {:.3e} flux {:.3e}",
        max_abs_diff(&field.u, &reference.0.u),
        max_abs_diff(&field.v, &reference.0.v),
        max_abs_diff(&field.p, &reference.0.p),
        max_abs_diff(&field.flux, &reference.0.flux),
    );
    println!("  snapshot/restore on the moving mesh — bit-identical re-run over 10 steps");
    Ok(())
}

#[test]
fn set_mesh_refuses_a_topology_change() -> CfdResult<()> {
    let config = CfdConfig::new().with_density(RHO).with_viscosity(MU);
    let mut solver = CurvilinearPisoSolver::new(
        config,
        CurvilinearParameters::default(),
        base_annulus(24, 6)?,
    )?;
    assert!(solver.set_mesh(base_annulus(24, 8)?).is_err());
    assert!(solver.set_mesh(cartesian(24, 6, 1.0, 1.0, false)?).is_err());
    assert!(solver.set_mesh(base_annulus(24, 6)?).is_ok());
    Ok(())
}

// ------------------------------------------------- linear-field falsifier

/// A linear divergence-free field `u = (U + G y, V + G x)` is an exact steady
/// Stokes solution with flat pressure on any moving mesh. `L_f` is exact on
/// linear fields and the pressure is zero, so the only discrete error is the
/// upwind mesh-flux term, `O(h ∂w ∂u)` per unit time — it must FALL with
/// refinement. An `O(displacement · G)` error that does not fall means the
/// field is being carried with the mesh.
#[tokio::test]
async fn linear_field_on_the_moving_annulus_converges() -> CfdResult<()> {
    const G: f64 = 0.3;
    let lin = |x: f64, y: f64| (U_UNIFORM + G * y, V_UNIFORM + G * x);
    let dt = 1e-3;
    let steps = 400;
    let mut errs = Vec::new();
    for ns in [32usize, 64, 128] {
        let config = CfdConfig::new().with_density(1.0).with_viscosity(0.05);
        let params = CurvilinearParameters {
            tolerance: 1e-10,
            convection: PatchConvection::None,
            ..CurvilinearParameters::default()
        };
        let base = base_annulus(ns, ns / 4)?;
        let mut solver = CurvilinearPisoSolver::new(config, params, moved(&base, 0.0, &deform)?)?;
        solver.set_boundary_velocity(move |x, y, _t| lin(x, y));
        let mut field = PatchField::new(solver.mesh());
        solver.initialize(&mut field, lin);
        let mut worst = 0.0_f64;
        for step in 0..steps {
            let t_new = (step + 1) as f64 * dt;
            solver.set_mesh(moved(&base, t_new, &deform)?)?;
            solver.advance(&mut field, dt).await?;
            let mesh = solver.mesh();
            for c in 0..mesh.cell_count() {
                let xy = mesh.centre(c);
                let (ue, ve) = lin(xy[0], xy[1]);
                worst = worst.max((field.u[c] - ue).abs().max((field.v[c] - ve).abs()));
            }
        }
        let max_p = field.p.iter().fold(0.0_f64, |m, &p| m.max(p.abs()));
        println!(
            "  linear field, ns = {ns}: max |u - u_lin| over {steps} steps = {worst:.3e}, max |p| = {max_p:.3e}"
        );
        errs.push(worst);
    }
    let o: Vec<f64> = errs.windows(2).map(|w| (w[0] / w[1]).log2()).collect();
    println!("  linear-field orders {o:?}");
    assert!(
        o.iter().all(|&x| x > 0.7),
        "the moving-mesh error on a linear field does not fall with refinement: {errs:?}"
    );
    Ok(())
}

/// Diagnostic runner: march `exact(x, y, t)` (its own Dirichlet data) in
/// the Stokes limit on `base` moved by `motion` (or fixed), returning the
/// L2 velocity error at `t_end` on the final geometry and max |p|.
async fn stokes_march(
    base: &PatchMesh,
    motion: Option<&dyn Fn(f64, f64, f64) -> [f64; 2]>,
    exact: fn(f64, f64, f64) -> (f64, f64),
    nu: f64,
    dt: f64,
    steps: usize,
) -> CfdResult<(f64, f64)> {
    let config = CfdConfig::new().with_density(1.0).with_viscosity(nu);
    let params = CurvilinearParameters {
        tolerance: 1e-8,
        convection: PatchConvection::None,
        ..CurvilinearParameters::default()
    };
    let start = match motion {
        Some(m) => moved(base, 0.0, m)?,
        None => base.clone(),
    };
    let mut solver = CurvilinearPisoSolver::new(config, params, start)?;
    solver.set_boundary_velocity(exact);
    let mut field = PatchField::new(solver.mesh());
    solver.initialize(&mut field, |x, y| exact(x, y, 0.0));
    for step in 0..steps {
        if let Some(m) = motion {
            solver.set_mesh(moved(base, (step + 1) as f64 * dt, m)?)?;
        }
        solver.advance(&mut field, dt).await?;
    }
    let t_end = steps as f64 * dt;
    let mesh = solver.mesh();
    let (mut sq, mut vol) = (0.0, 0.0);
    for c in 0..mesh.cell_count() {
        let xy = mesh.centre(c);
        let (ue, ve) = exact(xy[0], xy[1], t_end);
        let (eu, ev) = (field.u[c] - ue, field.v[c] - ve);
        sq += (eu * eu + ev * ev) * mesh.area(c);
        vol += mesh.area(c);
    }
    let max_p = field.p.iter().fold(0.0_f64, |m, &p| m.max(p.abs()));
    Ok(((sq / vol).sqrt(), max_p))
}

/// The linear field on the moving UNIT SQUARE (boundary nodes sliding along
/// fixed walls, non-periodic sides, four corners): the error must fall with
/// refinement (measured 2.99e-4 / 1.39e-4 / 6.80e-5, first order).
#[tokio::test]
async fn linear_field_on_the_moving_square_converges() -> CfdResult<()> {
    fn lin(x: f64, y: f64, _t: f64) -> (f64, f64) {
        (U_UNIFORM + 0.3 * y, V_UNIFORM + 0.3 * x)
    }
    let mut errs = Vec::new();
    for n in [16usize, 32, 64] {
        let base = cartesian(n, n, 1.0, 1.0, false)?;
        let (fixed, _) = stokes_march(&base, None, lin, 0.02, 1e-3, 250).await?;
        let (moving, max_p) = stokes_march(&base, Some(&tg_motion), lin, 0.02, 1e-3, 250).await?;
        println!(
            "  linear on the square, n = {n}: fixed L2 {fixed:.3e}, moving L2 {moving:.3e}, max |p| {max_p:.3e}"
        );
        errs.push(moving);
    }
    let o: Vec<f64> = errs.windows(2).map(|w| (w[0] / w[1]).log2()).collect();
    println!("  orders {o:?}");
    assert!(
        o.iter().all(|&x| x > 0.7),
        "linear field on the moving square: orders {o:?}"
    );
    Ok(())
}

/// Taylor–Green (Stokes limit) on the moving ANNULUS — the boundary moves
/// with the mesh and the patch is periodic: the second order of the P0
/// annulus gate must survive the motion.
#[tokio::test]
async fn taylor_green_stokes_on_the_moving_annulus_is_second_order() -> CfdResult<()> {
    fn tg(x: f64, y: f64, t: f64) -> (f64, f64) {
        (tg_u(x, y, t), tg_v(x, y, t))
    }
    let mut errs = Vec::new();
    for ns in [32usize, 64, 128] {
        let base = base_annulus(ns, ns / 4)?;
        let mut h = f64::INFINITY;
        for c in 0..base.cell_count() {
            for (f, _) in base.cell_faces(c) {
                let d = base.faces()[f].d;
                h = h.min((d[0] * d[0] + d[1] * d[1]).sqrt());
            }
        }
        let dt = 0.4 * h * h / (4.0 * NU);
        let steps = (T_END / dt).ceil() as usize;
        let dt = T_END / steps as f64;
        let (fixed, _) = stokes_march(&base, None, tg, NU, dt, steps).await?;
        let (moving, max_p) = stokes_march(&base, Some(&deform), tg, NU, dt, steps).await?;
        println!(
            "  TG Stokes on the annulus, ns = {ns}: fixed L2 {fixed:.3e}, moving L2 {moving:.3e}, max |p| {max_p:.3e}, {steps} steps"
        );
        errs.push(moving);
    }
    let o: Vec<f64> = errs.windows(2).map(|w| (w[0] / w[1]).log2()).collect();
    println!("  orders {o:?}");
    assert!(
        o.iter().all(|&x| x >= 1.8),
        "TG Stokes on the moving annulus: orders {o:?}"
    );
    Ok(())
}

// ---------------------------------------------------------------- gate 6

const NU: f64 = 0.02;
const T_END: f64 = 0.25;

fn tg_amp(t: f64) -> f64 {
    (-2.0 * NU * PI * PI * t).exp()
}
fn tg_u(x: f64, y: f64, t: f64) -> f64 {
    tg_amp(t) * (PI * x).sin() * (PI * y).cos()
}
fn tg_v(x: f64, y: f64, t: f64) -> f64 {
    -tg_amp(t) * (PI * x).cos() * (PI * y).sin()
}
fn tg_p(x: f64, y: f64, t: f64) -> f64 {
    let a = tg_amp(t);
    -RHO * a * a / 4.0 * ((2.0 * PI * x).cos() + (2.0 * PI * y).cos())
}

/// Interior motion of the unit square: the `ale_dgcl.rs` tensor wiggle plus
/// a non-tensor term that skews the cells; boundary-fixed (every term
/// carries `sin(πx)` or `sin(πy)`); displacement gradient below one.
fn tg_motion(x: f64, y: f64, t: f64) -> [f64; 2] {
    let (sx, sy) = ((PI * x).sin(), (PI * y).sin());
    let t = t * env_f64("RTX_CURV_TG_RATE", 1.0);
    let a = env_f64("RTX_CURV_TG_AMP", 1.0);
    let b = env_f64("RTX_CURV_TG_SKEW", 1.0);
    [
        x + a * 0.06 * sx * (2.9 * t + 3.0 * x).sin()
            + b * 0.03 * sx * sy * (4.3 * t + 2.0 * y).sin(),
        y + a * 0.06 * sy * (4.3 * t + 2.0 * y).sin()
            + b * 0.03 * sx * sy * (2.9 * t + 3.0 * x + 1.0).sin(),
    ]
}

/// Print-only diagnostic knobs (`RTX_CURV_*`), never gates.
fn env_f64(k: &str, d: f64) -> f64 {
    std::env::var(k)
        .ok()
        .and_then(|v| v.parse().ok())
        .unwrap_or(d)
}

/// March Taylor–Green on the patch to `T_END`, returning the L2 velocity
/// error on the final geometry and the kinetic-energy ratio.
async fn tg_measure(n: usize, moving: bool, convection: PatchConvection) -> CfdResult<(f64, f64)> {
    let h = 1.0 / n as f64;
    let dt = env_f64("RTX_CURV_DTFRAC", 0.4) * h * h / (4.0 * NU);
    let steps = (T_END / dt).ceil() as usize;
    let dt = T_END / steps as f64;

    let config = CfdConfig::new()
        .with_density(RHO)
        .with_viscosity(RHO * NU)
        .with_reference_velocity(1.0)
        .with_reference_length(1.0);
    let params = CurvilinearParameters {
        tolerance: 1e-6,
        convection,
        normal_diffusion: if std::env::var("RTX_CURV_LINE").is_ok() {
            NormalDiffusion::LineImplicit
        } else {
            NormalDiffusion::Explicit
        },
        ..CurvilinearParameters::default()
    };
    let base = cartesian(n, n, 1.0, 1.0, false)?;
    // Diagnostic: freeze the "moving" run on the deformed mesh at this
    // time and never move it (separates static skew from skewing motion).
    let freeze = std::env::var("RTX_CURV_TG_FREEZE")
        .ok()
        .and_then(|v| v.parse::<f64>().ok());
    // The moving run STARTS on the deformed mesh at t = 0: the motion is not
    // the identity there (its phase terms), and starting on the Cartesian
    // mesh made the first step sweep 2–4 cells at once — a mesh CFL far
    // above one for the explicit upwinded mesh flux, which imprinted an
    // O(displacement) error that no refinement removed (measured: L2
    // 1.14e-2 / 1.17e-2 / 1.51e-2 at n = 16/32/64 in the Stokes limit
    // against a second-order fixed ladder).
    let start = match (moving, freeze) {
        (true, Some(tf)) => moved(&base, tf, &tg_motion)?,
        (true, None) => moved(&base, 0.0, &tg_motion)?,
        (false, _) => base.clone(),
    };
    let moving = moving && freeze.is_none();
    let mut solver = CurvilinearPisoSolver::new(config, params, start)?;
    solver.set_boundary_velocity(|x, y, t| (tg_u(x, y, t), tg_v(x, y, t)));
    let mut field = PatchField::new(solver.mesh());
    solver.initialize(&mut field, |x, y| (tg_u(x, y, 0.0), tg_v(x, y, 0.0)));
    for c in 0..base.cell_count() {
        let xy = solver.mesh().centre(c);
        field.p[c] = tg_p(xy[0], xy[1], 0.0);
    }
    let energy = |field: &PatchField, mesh: &PatchMesh| -> f64 {
        (0..mesh.cell_count())
            .map(|c| 0.5 * RHO * (field.u[c] * field.u[c] + field.v[c] * field.v[c]) * mesh.area(c))
            .sum()
    };
    let e0 = energy(&field, solver.mesh());

    for step in 0..steps {
        if moving {
            let t_new = (step + 1) as f64 * dt;
            solver.set_mesh(moved(&base, t_new, &tg_motion)?)?;
        }
        let r = solver.advance(&mut field, dt).await?;
        assert!(
            r.poisson_converged,
            "n = {n} moving = {moving} step {step}: {r:?}"
        );
    }

    let mesh = solver.mesh();
    let (mut sq, mut vol) = (0.0, 0.0);
    let (mut worst, mut worst_xy) = (0.0_f64, [0.0; 2]);
    for c in 0..mesh.cell_count() {
        let xy = mesh.centre(c);
        let eu = field.u[c] - tg_u(xy[0], xy[1], T_END);
        let ev = field.v[c] - tg_v(xy[0], xy[1], T_END);
        sq += (eu * eu + ev * ev) * mesh.area(c);
        vol += mesh.area(c);
        let e = (eu * eu + ev * ev).sqrt();
        if e > worst {
            worst = e;
            worst_xy = xy;
        }
    }
    if std::env::var("RTX_CURV_TRACE").is_ok() {
        // Against the exact field at the REFERENCE (undeformed) centres: if
        // this is much smaller, the field is being carried with the mesh.
        let (mut sq_ref, mut sq_lag) = (0.0, 0.0);
        for c in 0..mesh.cell_count() {
            let r = base.centre(c);
            let eu = field.u[c] - tg_u(r[0], r[1], T_END);
            let ev = field.v[c] - tg_v(r[0], r[1], T_END);
            sq_ref += (eu * eu + ev * ev) * mesh.area(c);
            let xy = mesh.centre(c);
            let du = tg_u(xy[0], xy[1], T_END) - tg_u(r[0], r[1], T_END);
            let dv = tg_v(xy[0], xy[1], T_END) - tg_v(r[0], r[1], T_END);
            sq_lag += (du * du + dv * dv) * mesh.area(c);
        }
        println!(
            "    n = {n} moving = {moving}: max pointwise error {worst:.3e} at ({:.3}, {:.3}); \
             L2 vs exact at reference centres {:.3e}; full-Lagrangian L2 would be {:.3e}",
            worst_xy[0],
            worst_xy[1],
            (sq_ref / vol).sqrt(),
            (sq_lag / vol).sqrt()
        );
    }
    Ok(((sq / vol).sqrt(), energy(&field, mesh) / e0))
}

/// Diagnostic (print-only): the worst mesh quality the Taylor–Green motion
/// produces at resolution `n` over `[0, T_END]` — largest interior
/// non-orthogonality angle, smallest/largest cell area relative to `h²`,
/// and the largest face aspect ratio.
#[test]
fn taylor_green_motion_mesh_quality() -> CfdResult<()> {
    let n: usize = env_f64("RTX_CURV_N", 64.0) as usize;
    let base = cartesian(n, n, 1.0, 1.0, false)?;
    let h2 = 1.0 / (n * n) as f64;
    let (mut worst_angle, mut min_area, mut max_area, mut max_aspect) =
        (0.0_f64, f64::INFINITY, 0.0_f64, 0.0_f64);
    for k in 0..=50 {
        let t = T_END * k as f64 / 50.0;
        let m = moved(&base, t, &tg_motion)?;
        for f in m.faces() {
            if f.owner.is_some() && f.neigh.is_some() {
                let len = (f.s[0] * f.s[0] + f.s[1] * f.s[1]).sqrt();
                let dl = (f.d[0] * f.d[0] + f.d[1] * f.d[1]).sqrt();
                let cos = ((f.s[0] * f.d[0] + f.s[1] * f.d[1]) / (len * dl)).clamp(-1.0, 1.0);
                worst_angle = worst_angle.max(cos.acos().to_degrees());
            }
        }
        for c in 0..m.cell_count() {
            min_area = min_area.min(m.area(c) / h2);
            max_area = max_area.max(m.area(c) / h2);
            let faces = m.cell_faces(c);
            let len = |f: usize| {
                let s = m.faces()[f].s;
                (s[0] * s[0] + s[1] * s[1]).sqrt()
            };
            let (a, b) = (len(faces[0].0), len(faces[2].0));
            max_aspect = max_aspect.max((a / b).max(b / a));
        }
    }
    println!(
        "  TG motion at n = {n}: max non-orthogonality {worst_angle:.1} deg, cell area / h^2 in [{min_area:.3}, {max_area:.3}], max aspect {max_aspect:.2}"
    );
    let mut worst_gcl = 0.0_f64;
    for k in 0..50 {
        let (t0, t1) = (T_END * k as f64 / 50.0, T_END * (k + 1) as f64 / 50.0);
        let (old, new) = (moved(&base, t0, &tg_motion)?, moved(&base, t1, &tg_motion)?);
        worst_gcl = worst_gcl
            .max(StepGeometry::new(&old, &new, SweptFaceRule::Trapezoidal).gcl_defect(&old, &new));
    }
    println!("  TG motion at n = {n}: GCL defect {worst_gcl:.3e}");
    Ok(())
}

/// Gate: the observed order under mesh motion equals the fixed-mesh order,
/// in both convection modes — upwind (≈ 1, the physical scheme) and the
/// Stokes limit (≥ 1.8, the P0 second-order gate) — and the moving-mesh
/// error stays commensurate with the fixed one at equal resolution
/// (measured 1.05× upwind, 2.4× Stokes: the centred mesh flux costs a
/// constant, not an order).
#[tokio::test]
async fn taylor_green_order_is_unchanged_under_mesh_motion() -> CfdResult<()> {
    let exact_ratio = (-4.0 * NU * PI * PI * T_END).exp();
    let ns = [16usize, 32, 64, 128];
    let only = std::env::var("RTX_CURV_N")
        .ok()
        .and_then(|v| v.parse::<usize>().ok());
    for convection in [PatchConvection::Upwind, PatchConvection::None] {
        let (band, error_factor): (std::ops::Range<f64>, f64) = match convection {
            PatchConvection::Upwind => (0.7..1.6, 2.0),
            PatchConvection::None => (1.8..2.4, 3.0),
        };
        let mut fixed = Vec::new();
        let mut moving = Vec::new();
        for &n in &ns {
            if only.is_some_and(|o| o != n) || (only.is_none() && n == 128) {
                continue;
            }
            let (ef, rf) = tg_measure(n, false, convection).await?;
            let (em, rm) = tg_measure(n, true, convection).await?;
            println!(
                "  {convection:?} n = {n}: fixed L2 {ef:.4e} E(T)/E(0) {rf:.5} | moving L2 {em:.4e} \
                 E(T)/E(0) {rm:.5} (exact {exact_ratio:.5}, ratio moving/fixed {:.3})",
                em / ef
            );
            assert!(
                em < error_factor * ef,
                "{convection:?} n = {n}: mesh motion inflated the L2 error {ef:.3e} -> {em:.3e}"
            );
            // The energy decay is the scheme's business, not the motion's:
            // upwind misses the closed form by 8.6% at n = 16 on the FIXED
            // mesh (the staggered test's 5% gate was set at n = 32), so the
            // claim is that motion leaves the decay where the fixed mesh
            // puts it (measured: 0.3% upwind, 0.03% Stokes).
            assert!(
                (rm - rf).abs() < 0.01 * rf,
                "{convection:?} n = {n}: mesh motion changed the energy decay {rf:.5} -> {rm:.5}"
            );
            fixed.push(ef);
            moving.push(em);
        }
        let order =
            |e: &[f64]| -> Vec<f64> { e.windows(2).map(|w| (w[0] / w[1]).log2()).collect() };
        let (of, om) = (order(&fixed), order(&moving));
        println!("  {convection:?} orders: fixed {of:?}, moving {om:?}");
        for (a, b) in of.iter().zip(&om) {
            assert!(
                band.contains(a),
                "{convection:?} fixed order {a:.3} outside {band:?}"
            );
            assert!(
                band.contains(b),
                "{convection:?} moving order {b:.3} outside {band:?}"
            );
            assert!(
                (a - b).abs() < 0.3,
                "{convection:?}: mesh motion changed the observed order: fixed {a:.3}, moving {b:.3}"
            );
        }
    }
    Ok(())
}