//! Decaying Taylor–Green vortex: the transient benchmark for PISO.
//!
//! With wavenumber `k = pi` on the unit square,
//!
//! ```text
//!     u =  A(t) sin(pi x) cos(pi y)         A(t) = e^(-2 nu pi^2 t)
//!     v = -A(t) cos(pi x) sin(pi y)
//!     p = -(rho A^2 / 4)(cos 2pi x + cos 2pi y)
//! ```
//!
//! is an exact unsteady Navier–Stokes solution with **zero body force**: the
//! convective term is balanced identically by the pressure gradient and the
//! decay comes from viscosity alone. The normal velocity vanishes on all
//! four walls for all time (`u(0,y) = u(1,y) = 0`, `v(x,0) = v(x,1) = 0`),
//! so the closed staggered box fits it exactly; only the *tangential* wall
//! velocity decays in time, and it reaches the solver by re-setting the
//! wall-velocity hook with the current amplitude before every step.
//!
//! This is the same spatial field as `tests/mms_piso.rs`, which pinned the
//! steady spatial discretisation with a manufactured source. What this adds
//! is exactly the transient machinery that test cannot see: the time
//! derivative, the unsteady pressure–velocity coupling and the projection's
//! splitting error, exercised with no source hook at all. Two checks:
//!
//! 1. The L2 velocity error at `T` falls under simultaneous space–time
//!    refinement (`dt ~ h^2`, matching first-order upwind's spatial error
//!    against explicit Euler's temporal one).
//! 2. The kinetic-energy decay rate matches the closed form
//!    `E(T)/E(0) = e^(-4 nu pi^2 T)` — a scalar with an exact answer, and a
//!    check no steady measurement can make at all.

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())
}

struct Measurement {
    l2_velocity: f64,
    energy_ratio: f64,
    max_div: f64,
    steps: usize,
}

async fn measure(n: usize) -> CfdResult<Measurement> {
    let dx = 1.0 / n as f64;

    // Explicit predictor: dt under the diffusion limit, so dt ~ h^2 and the
    // temporal error refines together with the spatial one.
    let dt = 0.4 * dx * dx / (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);
    // Enough correctors that every step is driven to the divergence
    // tolerance the corrector loop itself measures — 2 is not enough on the
    // finer grids, where 400 Gauss-Seidel sweeps per projection leave a
    // residual the next corrector must mop up.
    let params = PisoParameters {
        corrector_steps: 60,
        time_step: dt,
        tolerance: 1e-9,
        ..PisoParameters::default()
    };
    let mut solver = PisoSolver::new(config, params)?;

    let mut field = FlowField::new(n, n, dx, dx)?;

    // Exact initial condition on every face, boundary faces included (the
    // normal boundary values are zero and stay zero).
    for j in 0..n {
        let y = (j as f64 + 0.5) * dx;
        for i in 0..=n {
            field.u[(j, i)] = u_exact(i as f64 * dx, y, 0.0);
        }
    }
    for j in 0..=n {
        let y = j as f64 * dx;
        for i in 0..n {
            field.v[(j, i)] = v_exact((i as f64 + 0.5) * dx, y, 0.0);
        }
    }
    for j in 0..n {
        for i in 0..n {
            field.p[(j, i)] = p_exact((i as f64 + 0.5) * dx, (j as f64 + 0.5) * dx, 0.0);
        }
    }

    let initial_energy = kinetic_energy(&field, n, dx);

    let empty = BoundaryConditions::new();
    for step in 0..steps {
        // The tangential wall velocity decays with the solution; the
        // predictor differentiates the state at t_n, so the wall belongs to
        // t_n as well.
        let t = step as f64 * dt;
        solver.set_wall_velocity(move |x, y| (u_exact(x, y, t), v_exact(x, y, t)));
        let result = solver.solve_time_step(&mut field, &empty, dt).await?;
        assert!(
            result.solver_result.converged,
            "step {step}: projection left mass residual {:.3e}",
            result.solver_result.final_residual
        );
    }

    let mut squared = 0.0;
    let mut volume = 0.0;
    for j in 0..n {
        let y = (j as f64 + 0.5) * dx;
        for i in 1..n {
            let e = field.u[(j, i)] - u_exact(i as f64 * dx, y, T_END);
            squared += e * e * dx * dx;
            volume += dx * dx;
        }
    }
    for j in 1..n {
        let y = j as f64 * dx;
        for i in 0..n {
            let e = field.v[(j, i)] - v_exact((i as f64 + 0.5) * dx, y, T_END);
            squared += e * e * dx * dx;
            volume += dx * dx;
        }
    }

    let mut max_div: f64 = 0.0;
    for j in 0..n {
        for i in 0..n {
            let div = (field.u[(j, i + 1)] - field.u[(j, i)]) / dx
                + (field.v[(j + 1, i)] - field.v[(j, i)]) / dx;
            max_div = max_div.max(div.abs());
        }
    }

    Ok(Measurement {
        l2_velocity: squared.sqrt() / volume.sqrt(),
        energy_ratio: kinetic_energy(&field, n, dx) / initial_energy,
        max_div,
        steps,
    })
}

/// Discrete kinetic energy over the interior faces.
fn kinetic_energy(field: &FlowField, n: usize, dx: f64) -> f64 {
    let mut energy = 0.0;
    for j in 0..n {
        for i in 1..n {
            energy += 0.5 * RHO * field.u[(j, i)] * field.u[(j, i)] * dx * dx;
        }
    }
    for j in 1..n {
        for i in 0..n {
            energy += 0.5 * RHO * field.v[(j, i)] * field.v[(j, i)] * dx * dx;
        }
    }
    energy
}

#[tokio::test]
async fn taylor_green_decays_at_the_exact_rate() -> CfdResult<()> {
    let resolutions = [16usize, 32, 64];
    let exact_ratio = (-4.0 * NU * PI * PI * T_END).exp();

    let mut measurements = Vec::new();
    for &n in &resolutions {
        measurements.push(measure(n).await?);
    }

    let errors: Vec<f64> = measurements.iter().map(|m| m.l2_velocity).collect();
    let rates: Vec<f64> = errors
        .windows(2)
        .map(|pair| (pair[0] / pair[1]).log2())
        .collect();

    for (i, &n) in resolutions.iter().enumerate() {
        let rate = if i == 0 {
            String::from("    -")
        } else {
            format!("{:5.2}", rates[i - 1])
        };
        println!(
            "  n = {n:3} ({:4} steps)  L2 = {:.6e}  order = {rate}  E(T)/E(0) = {:.5} \
             (exact {exact_ratio:.5})  max div = {:.2e}",
            measurements[i].steps, errors[i], measurements[i].energy_ratio, measurements[i].max_div
        );
    }

    assert!(
        errors.windows(2).all(|pair| pair[1] < pair[0]),
        "the error must fall under refinement; got {errors:?}"
    );

    // Measured: L2 = 2.267e-2, 1.153e-2, 5.841e-3 — orders 0.97 and 0.98,
    // first-order upwind's rate, with dt ~ h^2 keeping the temporal error
    // subordinate.
    for (i, &rate) in rates.iter().enumerate() {
        assert!(
            (0.85..1.5).contains(&rate),
            "refinement {} -> {}: observed order {rate:.3}, expected ~1 from \
             first-order upwind. Errors: {errors:?}",
            resolutions[i],
            resolutions[i + 1]
        );
    }

    // The projection must keep every step divergence-free. Measured 1e-10,
    // 7e-9, 4.5e-7 with the SOR inner solve; the 400-sweep Gauss-Seidel this
    // test originally ran against left 1e-2 here, growing with mesh size.
    for m in &measurements {
        assert!(m.max_div < 1e-6, "max divergence {:.3e}", m.max_div);
    }

    // Energy decay: the deficit against the exact ratio is upwind's excess
    // numerical dissipation and must halve per refinement. Measured deficits
    // 0.0690, 0.0360, 0.0185 (ratios 1.92, 1.95); the finest mesh sits
    // within 2.3% of the closed form.
    let deficits: Vec<f64> = measurements
        .iter()
        .map(|m| exact_ratio - m.energy_ratio)
        .collect();
    for pair in deficits.windows(2) {
        let ratio = pair[0] / pair[1];
        assert!(
            (1.6..2.4).contains(&ratio),
            "energy-deficit refinement ratio {ratio:.2}, expected ~2; \
             deficits {deficits:?}"
        );
    }
    assert!(
        deficits[deficits.len() - 1] < 0.03 * exact_ratio,
        "finest-mesh energy ratio {:.5} is more than 3% from the exact \
         {exact_ratio:.5}",
        measurements[measurements.len() - 1].energy_ratio
    );

    Ok(())
}