//! Physics on the moving mesh: the ALE solver must reproduce fixed-grid
//! results when the mesh does not move, and must keep them when it does.
//!
//! Two claims, on the decaying Taylor–Green vortex (`tests/taylor_green.rs`
//! has the closed form and the reasoning; wavenumber `pi` on the unit box,
//! zero body force, normal velocities exactly zero on the fixed boundary):
//!
//! 1. **Degeneracy**: with zero mesh motion on a uniform grid, the
//!    conservative ALE update is algebraically identical to the fixed-grid
//!    PISO scheme — same fluxes, same projection, same inner solve — so the
//!    two solvers must agree step for step to rounding, not to truncation.
//!    This pins every geometric generalisation (non-uniform spacings, swept
//!    volumes, half-face v-weights) to the verified PISO implementation.
//!
//! 2. **Invariance under mesh motion**: the interior mesh lines wiggling
//!    (same arbitrary motion as the DGCL test) must not change what the
//!    scheme converges to. The L2 error against the exact solution still
//!    falls at first order under space–time refinement, and the
//!    kinetic-energy decay still approaches `e^(-4 nu pi^2 T)`. The mesh
//!    motion is fixed in physical space while the grid refines, so finer
//!    meshes resolve the *same* moving-mesh problem.

use rtx_cfd::solvers::incompressible::ale::{
    AleField, AleParameters, AlePisoSolver, SweptFaceRule,
};
use rtx_cfd::solvers::incompressible::{
    BoundaryConditions, FlowField, IncompressibleSolver, PisoParameters, PisoSolver,
};
use rtx_cfd::{CfdConfig, CfdResult};
use std::f64::consts::PI;

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

fn amplitude(t: f64) -> f64 {
    (-2.0 * NU * PI * PI * t).exp()
}

fn u_exact(x: f64, y: f64, t: f64) -> f64 {
    amplitude(t) * (PI * x).sin() * (PI * y).cos()
}

fn v_exact(x: f64, y: f64, t: f64) -> f64 {
    -amplitude(t) * (PI * x).cos() * (PI * y).sin()
}

fn p_exact(x: f64, y: f64, t: f64) -> f64 {
    let a = amplitude(t);
    -RHO * a * a / 4.0 * ((2.0 * PI * x).cos() + (2.0 * PI * y).cos())
}

/// The DGCL test's interior mesh motion, on the unit square: smooth,
/// boundary-fixed, lines out of phase, displacement gradient below 1.
fn moved(reference: f64, t: f64, rate: f64, phase: f64) -> f64 {
    reference + 0.06 * (PI * reference).sin() * (rate * t + phase * reference).sin()
}

fn config() -> CfdConfig {
    CfdConfig::new()
        .with_density(RHO)
        .with_viscosity(RHO * NU)
        .with_reference_velocity(1.0)
        .with_reference_length(1.0)
}

fn ale_solver(tolerance: f64) -> CfdResult<AlePisoSolver> {
    let params = AleParameters {
        corrector_steps: 60,
        tolerance,
        swept_face_rule: SweptFaceRule::Trapezoidal,
        ..AleParameters::default()
    };
    let mut solver = AlePisoSolver::new(config(), params)?;
    solver.set_boundary_velocity(|x, y, t| (u_exact(x, y, t), v_exact(x, y, t)));
    Ok(solver)
}

fn tg_field(n: usize) -> CfdResult<AleField> {
    let mut field = AleField::uniform(n, n, 1.0, 1.0)?;
    let h = 1.0 / n as f64;
    for j in 0..n {
        let y = (j as f64 + 0.5) * h;
        for i in 0..=n {
            field.u[(j, i)] = u_exact(i as f64 * h, y, 0.0);
        }
    }
    for j in 0..=n {
        let y = j as f64 * h;
        for i in 0..n {
            field.v[(j, i)] = v_exact((i as f64 + 0.5) * h, y, 0.0);
        }
    }
    for j in 0..n {
        for i in 0..n {
            field.p[(j, i)] = p_exact((i as f64 + 0.5) * h, (j as f64 + 0.5) * h, 0.0);
        }
    }
    Ok(field)
}

/// Volume-weighted L2 velocity error and kinetic energy on the current
/// (possibly non-uniform) geometry.
fn l2_error_and_energy(field: &AleField, t: f64) -> (f64, f64) {
    let n = field.nx;
    let xc: Vec<f64> = field.x.windows(2).map(|w| 0.5 * (w[0] + w[1])).collect();
    let yc: Vec<f64> = field.y.windows(2).map(|w| 0.5 * (w[0] + w[1])).collect();

    let mut squared = 0.0;
    let mut volume = 0.0;
    let mut energy = 0.0;
    for j in 0..n {
        let h = field.y[j + 1] - field.y[j];
        for i in 1..n {
            let w = xc[i] - xc[i - 1];
            let e = field.u[(j, i)] - u_exact(field.x[i], yc[j], t);
            squared += e * e * w * h;
            volume += w * h;
            energy += 0.5 * RHO * field.u[(j, i)] * field.u[(j, i)] * w * h;
        }
    }
    for j in 1..n {
        let h = yc[j] - yc[j - 1];
        for i in 0..n {
            let w = field.x[i + 1] - field.x[i];
            let e = field.v[(j, i)] - v_exact(xc[i], field.y[j], t);
            squared += e * e * w * h;
            volume += w * h;
            energy += 0.5 * RHO * field.v[(j, i)] * field.v[(j, i)] * w * h;
        }
    }
    (squared.sqrt() / volume.sqrt(), energy)
}

/// March Taylor–Green to `T_END` on a mesh that wiggles when `moving`.
async fn measure(n: usize, moving: bool) -> CfdResult<(f64, f64)> {
    let dt = 0.4 * (1.0 / n as f64).powi(2) / (4.0 * NU);
    let steps = (T_END / dt).ceil() as usize;
    let dt = T_END / steps as f64;

    let mut solver = ale_solver(1e-9)?;
    let mut field = tg_field(n)?;
    let (_, initial_energy) = l2_error_and_energy(&field, 0.0);

    let rx: Vec<f64> = (0..=n).map(|i| i as f64 / n as f64).collect();
    for step in 0..steps {
        let t_new = (step + 1) as f64 * dt;
        let (new_x, new_y): (Vec<f64>, Vec<f64>) = if moving {
            (
                rx.iter().map(|&x| moved(x, t_new, 2.9, 3.0)).collect(),
                rx.iter().map(|&y| moved(y, t_new, 4.3, 2.0)).collect(),
            )
        } else {
            (rx.clone(), rx.clone())
        };
        let result = solver.advance(&mut field, &new_x, &new_y, dt).await?;
        assert!(
            result.solver_result.converged,
            "n = {n} moving = {moving} step {step}: mass residual {:.3e}",
            result.solver_result.final_residual
        );
    }

    let (l2, final_energy) = l2_error_and_energy(&field, T_END);
    Ok((l2, final_energy / initial_energy))
}

#[tokio::test]
async fn zero_motion_on_a_uniform_grid_reproduces_piso_to_rounding() -> CfdResult<()> {
    let n = 16;
    let h = 1.0 / n as f64;
    let dt = 2e-3;
    let steps = 5;

    let mut ale = ale_solver(1e-9)?;
    let mut ale_field = tg_field(n)?;

    let piso_params = PisoParameters {
        corrector_steps: 60,
        time_step: dt,
        tolerance: 1e-9,
        ..PisoParameters::default()
    };
    let mut piso = PisoSolver::new(config(), piso_params)?;
    let mut piso_field = FlowField::new(n, n, h, h)?;
    piso_field.u.copy_from(&ale_field.u);
    piso_field.v.copy_from(&ale_field.v);
    piso_field.p.copy_from(&ale_field.p);

    let lines: Vec<f64> = (0..=n).map(|i| i as f64 / n as f64).collect();
    let empty = BoundaryConditions::new();
    for step in 0..steps {
        let t = step as f64 * dt;
        piso.set_wall_velocity(move |x, y| (u_exact(x, y, t), v_exact(x, y, t)));
        piso.solve_time_step(&mut piso_field, &empty, dt).await?;
        ale.advance(&mut ale_field, &lines, &lines, dt).await?;
    }

    let mut worst: f64 = 0.0;
    for (a, b) in ale_field.u.iter().zip(piso_field.u.iter()) {
        worst = worst.max((a - b).abs());
    }
    for (a, b) in ale_field.v.iter().zip(piso_field.v.iter()) {
        worst = worst.max((a - b).abs());
    }
    println!("  ALE vs PISO after {steps} steps: max |difference| = {worst:.3e}");

    // Same discretisation, different code paths: agreement to rounding.
    // (Not bit-identical — the ALE path forms spacings as differences of
    // node coordinates — but far below any truncation scale.)
    // Measured: 2.2e-16 — one ulp of the velocity scale.
    assert!(
        worst < 1e-12,
        "ALE with zero mesh motion diverged from PISO by {worst:.3e}"
    );
    Ok(())
}

#[tokio::test]
async fn taylor_green_survives_arbitrary_mesh_motion() -> CfdResult<()> {
    let exact_ratio = (-4.0 * NU * PI * PI * T_END).exp();

    let (err_fixed, ratio_fixed) = measure(32, false).await?;
    let (err_coarse, _) = measure(16, true).await?;
    let (err_moving, ratio_moving) = measure(32, true).await?;

    let order = (err_coarse / err_moving).log2();
    println!(
        "  fixed  n = 32: L2 = {err_fixed:.4e}  E(T)/E(0) = {ratio_fixed:.5} (exact {exact_ratio:.5})"
    );
    println!(
        "  moving n = 16: L2 = {err_coarse:.4e}\n  moving n = 32: L2 = {err_moving:.4e}  \
         order = {order:.2}  E(T)/E(0) = {ratio_moving:.5}"
    );

    // Measured: fixed n=32 L2 = 1.1532e-2 (PISO's published Taylor-Green
    // value to four digits); moving 2.4218e-2 -> 1.0729e-2, order 1.17;
    // energy ratios 0.78490 fixed / 0.78986 moving against exact 0.82087 —
    // upwind's dissipation deficit, unchanged by the motion.
    //
    // The moving mesh must not change what the scheme converges to: the
    // error still falls at ~first order (upwind) under refinement...
    assert!(
        (0.85..1.5).contains(&order),
        "moving-mesh refinement 16 -> 32: observed order {order:.3}, expected ~1; \
         errors {err_coarse:.3e} -> {err_moving:.3e}"
    );
    // ...and stays commensurate with the fixed-mesh error at equal
    // resolution — mesh motion may cost accuracy but not the solution.
    assert!(
        err_moving < 2.0 * err_fixed,
        "mesh motion inflated the L2 error {err_fixed:.3e} -> {err_moving:.3e}"
    );
    // Energy decay stays quantitative on the moving mesh: within upwind's
    // dissipation deficit of the closed form at this resolution.
    assert!(
        (ratio_moving - exact_ratio).abs() < 0.05 * exact_ratio,
        "moving-mesh energy ratio {ratio_moving:.5} vs exact {exact_ratio:.5}"
    );
    Ok(())
}