//! Code verification of the PISO solver by manufactured solution.
//!
//! The manufactured field, its momentum source and the grid convention are
//! exactly those of `tests/mms_navier_stokes.rs` — see that file for the
//! derivation. PISO is a transient stepper, so instead of iterating an outer
//! loop it is marched in time under the steady forcing until the field stops
//! changing; the steady state it lands on satisfies the same spatial
//! discretisation (first-order upwind convection, second-order diffusion,
//! half-cell wall treatment), so the observed order should match SIMPLE's:
//! approaching 1, limited by upwind's `O(h)` numerical viscosity.
//!
//! Until this file existed PISO had no verification of any kind — not a unit
//! test, not a benchmark. The first run of this measurement, against the old
//! implementation, is what confirmed the inverted pressure-correction sign
//! and the frozen near-wall lines recorded in `piso.rs`'s module docs.

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 MU: f64 = 0.05;

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

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

fn source(x: f64, y: f64) -> (f64, f64) {
    let fx = RHO * 0.5 * PI * (2.0 * PI * x).sin()
        + 2.0 * PI * PI * MU * u_exact(x, y)
        + PI * (PI * x).cos() * (PI * y).sin();
    let fy = RHO * 0.5 * PI * (2.0 * PI * y).sin()
        + 2.0 * PI * PI * MU * v_exact(x, y)
        + PI * (PI * x).sin() * (PI * y).cos();
    (fx, fy)
}

struct Measurement {
    l2_velocity: f64,
    max_div: f64,
}

/// March the manufactured problem on an `n` by `n` grid to steady state.
async fn measure(n: usize) -> CfdResult<Measurement> {
    let dx = 1.0 / n as f64;
    let dy = dx;

    // Explicit predictor: dt must respect the diffusion limit `dx^2 / (4 nu)`
    // (the binding one here, with nu = 0.05 and |u| <= 1).
    let nu = MU / RHO;
    let dt = 0.4 * (dx * dx / (4.0 * nu)).min(dx);

    let config = CfdConfig::new()
        .with_density(RHO)
        .with_viscosity(MU)
        .with_reference_velocity(1.0)
        .with_reference_length(1.0);
    let params = PisoParameters {
        corrector_steps: 2,
        time_step: dt,
        tolerance: 1e-8,
        ..PisoParameters::default()
    };
    let mut solver = PisoSolver::new(config, params)?;
    solver.set_momentum_source(source);
    solver.set_wall_velocity(|x, y| (u_exact(x, y), v_exact(x, y)));

    let mut field = FlowField::new(n, n, dx, dy)?;
    for j in 0..n {
        let y = (j as f64 + 0.5) * dy;
        field.u[(j, 0)] = u_exact(0.0, y);
        field.u[(j, n)] = u_exact(1.0, y);
    }
    for i in 0..n {
        let x = (i as f64 + 0.5) * dx;
        field.v[(0, i)] = v_exact(x, 0.0);
        field.v[(n, i)] = v_exact(x, 1.0);
    }

    // March to steady state: stop when the field stops moving, measured as
    // `max |u^{n+1} - u^n| / dt`, the discrete time derivative.
    let empty = BoundaryConditions::new();
    let mut steady_residual = f64::INFINITY;
    for _step in 0..200_000 {
        let u_before = field.u.clone();
        let v_before = field.v.clone();

        solver.solve_time_step(&mut field, &empty, dt).await?;

        let mut max_change: f64 = 0.0;
        for (a, b) in field.u.iter().zip(u_before.iter()) {
            max_change = max_change.max((a - b).abs());
        }
        for (a, b) in field.v.iter().zip(v_before.iter()) {
            max_change = max_change.max((a - b).abs());
        }
        steady_residual = max_change / dt;
        // 1e-6, not tighter: each step's projection is converged to the
        // solver's mass tolerance, not to machine zero, and the leftover
        // per-step noise floors |du/dt| just below 1e-6. The L2 errors being
        // measured are 1e-2 to 1e-3, so a 1e-6 stationarity floor
        // contributes nothing to them.
        if steady_residual < 1e-6 {
            break;
        }
    }
    assert!(
        steady_residual < 1e-6,
        "PISO did not reach a steady state: |du/dt| = {steady_residual:.3e}"
    );

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

    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)]) / dy;
            max_div = max_div.max(div.abs());
        }
    }

    Ok(Measurement {
        l2_velocity: squared.sqrt() / volume.sqrt(),
        max_div,
    })
}

/// The steady state PISO marches to must converge to the exact solution at
/// the rate the spatial discretisation dictates — order approaching 1 for
/// first-order upwind — and must be divergence-free in every cell.
///
/// Measured (16 -> 32 -> 64): L2 velocity 3.516214e-2, 1.953750e-2,
/// 1.037512e-2, orders 0.85 and 0.91, max |div u| ~ 1e-9 everywhere. The
/// errors agree with SIMPLE's on the same meshes (3.516212e-2, 1.953751e-2,
/// 1.037523e-2) to six or seven significant figures: two different
/// algorithms — implicit under-relaxed outer iteration against explicit time
/// marching with projection — land on the same discrete steady solution,
/// which is exactly what sharing a spatial discretisation must produce and
/// is very hard for two independently wrong solvers to fake.
#[tokio::test]
async fn piso_observed_order_matches_the_convection_scheme() -> CfdResult<()> {
    let resolutions = [16usize, 32, 64];
    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}   L2 velocity error = {:.6e}   observed order = {rate}   \
             max |div u| = {:.6e}",
            errors[i], measurements[i].max_div
        );
    }

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

    for (i, &rate) in rates.iter().enumerate() {
        assert!(
            rate > 0.75,
            "refinement {} -> {}: observed order {rate:.3}, below the order 1 \
             first-order upwind must deliver. Errors: {errors:?}",
            resolutions[i],
            resolutions[i + 1]
        );
        assert!(
            rate < 2.3,
            "refinement {} -> {}: observed order {rate:.3}, above what this \
             scheme can deliver — suspect the error measure. Errors: {errors:?}",
            resolutions[i],
            resolutions[i + 1]
        );
    }

    // Every cell, outer ring included, must satisfy continuity: the
    // projection exists for no other reason.
    for (m, &n) in measurements.iter().zip(&resolutions) {
        assert!(
            m.max_div < 1e-5,
            "max |div u| = {:.3e} at n = {n}: the projection is not removing \
             the divergence",
            m.max_div
        );
    }

    Ok(())
}