//! Code verification of the SIMPLE solver by manufactured solution.
//!
//! Cavity benchmarks tell you whether the answer looks like the picture in the
//! paper. This tells you the rate at which the discretisation converges to an
//! exact solution, which is a statement about the code rather than about the
//! flow, and which no amount of tuning can fake.
//!
//! # The manufactured solution
//!
//! ```text
//!     u(x, y) =  sin(pi x) cos(pi y)
//!     v(x, y) = -cos(pi x) sin(pi y)
//!     p(x, y) =  sin(pi x) sin(pi y)
//! ```
//!
//! The velocity field is divergence-free by construction —
//! `du/dx = pi cos cos` and `dv/dy = -pi cos cos` — which it must be, or the
//! pressure-correction equation is being asked to solve an inconsistent
//! problem and the measurement means nothing.
//!
//! It also has zero normal velocity on all four sides of the unit square
//! (`u(0, y) = u(1, y) = 0`, `v(x, 0) = v(x, 1) = 0`), so the domain is a
//! closed box: the net mass flux across the boundary is exactly zero and the
//! pure-Neumann pressure-correction system is compatible. The tangential
//! velocities on those walls are *not* zero — `u(x, 0) = sin(pi x)` — so the
//! walls are moving walls, and getting them to the solver is the point of
//! [`SimpleSolver::set_wall_velocity`].
//!
//! The pressure is deliberately *not* the one that goes with this velocity in
//! a force-free flow. With `p = (cos 2pi x + cos 2pi y)/4` the convective and
//! pressure terms cancel identically at unit density, which would leave the
//! convection discretisation untested. This choice keeps all three terms
//! present in the source.
//!
//! # Grid convention
//!
//! On this staggered layout, with `nx` by `ny` cells over the unit square:
//!
//! ```text
//!     u[(j, i)]  at  ( i dx,          (j + 0.5) dy )     i = 0..=nx
//!     v[(j, i)]  at  ( (i + 0.5) dx,  j dy         )     j = 0..=ny
//!     p[(j, i)]  at  ( (i + 0.5) dx,  (j + 0.5) dy )
//! ```
//!
//! Cell `i` is bounded by u-faces `i` and `i + 1`, which is the convention
//! `compute_mass_source` uses to form the divergence and therefore the one
//! that defines the grid.
//!
//! Every u row sits at `y = (j + 0.5) dy`, which is strictly interior; the only
//! u faces on a domain boundary are `i = 0` and `i = nx`. Symmetrically the
//! only boundary v faces are `j = 0` and `j = ny`. This test prescribes the
//! exact solution on exactly those faces and on nothing else — the near-wall
//! lines are unknowns and are solved for.
//!
//! # What this measures
//!
//! First-order upwind convection with second-order diffusion should give an
//! observed order of 1 in the velocity, and the three quantities asserted
//! below are the ones that a plausible-but-wrong discretisation fails:
//!
//! 1. **Observed order of the velocity error.** 0.85 and 0.91 over
//!    16 -> 32 -> 64, approaching 1 from below. It was 0.48. In the Stokes
//!    limit the same measurement gives 2.05 and 2.06, so the shortfall below 1
//!    is first-order upwind convection, not the discretisation.
//! 2. **Divergence, split between the outer ring of cells and the interior.**
//!    A converged SIMPLE solve must satisfy discrete continuity to solver
//!    tolerance in *every* cell. Splitting the measure is what exposes a
//!    continuity equation that is only enforced on the interior: a ring value
//!    orders of magnitude worse than the interior one is a boundary defect
//!    wearing an interior-looking global norm as a disguise.
//! 3. **Maximum pressure error with the mean removed**, since the pressure is
//!    determined only up to a constant. This is the sharpest of the three: a
//!    velocity field can look nearly right while the pressure it was computed
//!    from is wrong by an O(1) amount, and a pressure error that *grows* under
//!    refinement is proof of an inconsistent discretisation no tolerance
//!    tuning can hide.
//!
//! # History
//!
//! Found by this test, in order:
//!
//! - The diffusion conductances omitted the face area, reading `mu / dx` where
//!   finite volume requires `mu * A / delta`, so viscosity was too large by a
//!   factor of `1/h` — 65 times on a 65x65 mesh. Before that fix the error did
//!   not reduce under refinement at all.
//! - The x-momentum equation dropped the body force entirely: `u_source_term`
//!   was computed and never added to the source, while `v_source_term` was.
//!   The u equation was therefore solving a different problem from the one the
//!   manufactured solution defines, which is why `u` was markedly further from
//!   exact than `v` on the same mesh.
//! - The momentum sweeps froze the near-wall lines (`j = 0`, `j = ny - 1` for
//!   u) and treated them as boundaries, imposing the wall value half a cell
//!   inside the domain; continuity was enforced only on the interior cells, so
//!   the outer ring of cells was never made divergence-free. Together these
//!   held the observed order at 0.5 and made the pressure error *grow* under
//!   refinement.
//!
//! That last fix has a consequence worth stating here, because it is what makes
//! the near-wall lines unknowns: once every cell has a continuity equation,
//! every cell needs at least one face the pressure correction is allowed to
//! move. A boundary condition that pins an *interior* face therefore no longer
//! merely wastes work — it over-determines that cell and the solve diverges.
//! `tests/simple_tests.rs` was doing exactly that, prescribing whole u rows and
//! v columns that lie half a cell inside the domain. On a staggered grid the
//! only velocity components living on a boundary are the normal ones; the
//! tangential no-slip or lid condition is not a stored value at all and belongs
//! in [`SimpleSolver::set_wall_velocity`].

use rtx_cfd::solvers::incompressible::{
    BoundaryConditions, ConvectionScheme, FlowField, SimpleParameters, SimpleSolver,
};
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 p_exact(x: f64, y: f64) -> f64 {
    (PI * x).sin() * (PI * y).sin()
}

/// Momentum source `f = rho (u.grad)u - mu lap(u) + grad p`, derived by hand.
///
/// The convective terms collapse neatly:
/// `(u.grad)u = pi sin(pi x) cos(pi x) = (pi/2) sin(2 pi x)` and likewise
/// `(u.grad)v = (pi/2) sin(2 pi y)`, because `sin^2 + cos^2` factors out.
/// The Laplacians are `lap(u) = -2 pi^2 u` and `lap(v) = -2 pi^2 v`.
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)
}

/// Everything one grid resolution has to say about the discretisation.
struct Measurement {
    /// Volume-weighted L2 error of the velocity over the interior faces.
    l2_velocity: f64,
    /// Largest `|div u|` over the outer ring of cells.
    max_div_ring: f64,
    /// Largest `|div u|` over every cell not in the outer ring.
    max_div_interior: f64,
    /// Largest `|p - p_exact|` after removing the mean of the difference; the
    /// pressure is determined only up to an additive constant.
    max_p_error: f64,
}

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

    let config = CfdConfig::new()
        .with_density(RHO)
        .with_viscosity(MU)
        .with_reference_velocity(1.0)
        .with_reference_length(1.0);
    let params = SimpleParameters::default()
        .with_max_iterations(40000)
        .with_tolerance(1e-9)
        .with_convection_scheme(scheme);
    let mut solver = SimpleSolver::new(config, params)?;
    solver.set_momentum_source(source);

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

    // Impose the exact solution on the genuine boundary faces, and only there.
    //
    // These are the faces that lie on the domain boundary: u faces `i = 0` and
    // `i = n` on the left and right sides, v faces `j = 0` and `j = n` on the
    // bottom and top. The momentum sweeps never touch them and an empty
    // boundary-condition set leaves them alone for the whole solve, so this is
    // exactly the Dirichlet data the manufactured solution defines.
    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);
    }

    // The *tangential* velocity of each wall, which is what the near-wall
    // momentum control volumes need for their half-cell diffusion term. There
    // is no room to store it on this grid — no u node lies on the bottom wall —
    // so it reaches the solver as a function of position instead.
    solver.set_wall_velocity(|x, y| (u_exact(x, y), v_exact(x, y)));

    let empty = BoundaryConditions::new();
    for _ in 0..40000 {
        let (mass, momentum) = solver
            .solve_simple_iteration(&mut field, &empty, 0.01)
            .await?;
        if (mass * mass + momentum * momentum).sqrt() < 1e-9 {
            break;
        }
    }

    // L2 error over the interior faces, weighted by cell volume.
    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;
        }
    }

    // Discrete divergence, formed exactly as `compute_mass_source` forms it:
    // cell `(i, j)` is bounded by u faces `i` and `i + 1`, v faces `j` and
    // `j + 1`.
    let mut max_div_ring: f64 = 0.0;
    let mut max_div_interior: 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;
            let on_ring = i == 0 || i == n - 1 || j == 0 || j == n - 1;
            if on_ring {
                max_div_ring = max_div_ring.max(div.abs());
            } else {
                max_div_interior = max_div_interior.max(div.abs());
            }
        }
    }

    // Pressure error with the mean of the difference removed: the pressure
    // correction equation is pure Neumann, so `p` is only defined up to a
    // constant and comparing it raw would measure the anchor, not the solution.
    let mut mean_offset = 0.0;
    for j in 0..n {
        for i in 0..n {
            mean_offset += field.p[(j, i)] - p_exact((i as f64 + 0.5) * dx, (j as f64 + 0.5) * dy);
        }
    }
    mean_offset /= (n * n) as f64;
    let mut max_p_error: f64 = 0.0;
    for j in 0..n {
        for i in 0..n {
            let e = field.p[(j, i)]
                - p_exact((i as f64 + 0.5) * dx, (j as f64 + 0.5) * dy)
                - mean_offset;
            max_p_error = max_p_error.max(e.abs());
        }
    }

    Ok(Measurement {
        l2_velocity: squared.sqrt() / volume.sqrt(),
        max_div_ring,
        max_div_interior,
        max_p_error,
    })
}

/// The velocity error must fall at the rate the convection scheme dictates.
///
/// First-order upwind gives order 1. The lower bound is set just under what is
/// actually achieved; the upper bound is there because a rate well above 1 on a
/// first-order scheme means the error *measure* is wrong, not that the solver
/// is unexpectedly good.
#[tokio::test]
async fn 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, ConvectionScheme::Upwind).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 |p - p_exact| (mean removed) = {:.6e}",
            errors[i], measurements[i].max_p_error
        );
    }

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

    // 0.80, not 1.0. The measured rates are 0.848 (16 -> 32) and 0.905
    // (32 -> 64), approaching 1 from below — which is exactly what first-order
    // upwind convection gives.
    //
    // The wall treatment is *not* the limiter, contrary to what this comment
    // first claimed. Repeating the measurement in the Stokes limit, where
    // convection is negligible and every remaining operator is second order,
    // gives observed order 2.05 and 2.06:
    //
    //     rho = 1.000  (Re = 20.00)   3.5162e-2  1.9537e-2  1.0375e-2   0.85  0.91
    //     rho = 0.001  (Re =  0.02)   2.2131e-3  5.3510e-4  1.2812e-4   2.05  2.06
    //
    // So the near-wall half-cell term is second-order accurate and the whole
    // Stokes discretisation reaches its nominal rate. The shortfall below 1 at
    // Re = 20 is upwind's `O(h)` numerical viscosity and nothing else, which
    // means a second-order convection scheme is now the thing that moves this
    // number — and that work is no longer blocked behind an unverified
    // discretisation.
    //
    // The bound sits just under the worst measured rate rather than at the
    // theoretical 1 because the rate is still climbing at 64. It still has
    // teeth: before the near-wall rows became unknowns the rates were 0.475 and
    // 0.470, and this fails them by a wide margin.
    for (i, &rate) in rates.iter().enumerate() {
        assert!(
            rate > 0.80,
            "refinement {} -> {}: observed order {rate:.3}, below the order 1 a \
             first-order upwind convection scheme must deliver. Errors: {errors:?}",
            resolutions[i],
            resolutions[i + 1]
        );
        assert!(
            rate < 2.3,
            "refinement {} -> {}: observed order {rate:.3}, above what a \
             first-order convection scheme can deliver — suspect the error \
             measure rather than celebrating. Errors: {errors:?}",
            resolutions[i],
            resolutions[i + 1]
        );
    }

    // The pressure must converge too, and this is the assertion with the most
    // teeth of the three. Before the outer ring of cells was given a continuity
    // equation the error read 4.08e-1, 6.24e-1, 7.56e-1 on these three meshes —
    // O(1) on a field spanning [-1, 1], and *growing* under refinement, which is
    // proof of an inconsistent discretisation. No velocity norm revealed it: the
    // velocity error was falling the whole time.
    //
    // It now reads 9.25e-2, 5.23e-2, 2.80e-2 — falling at about the same rate as
    // the velocity (0.82 then 0.90).
    let p_errors: Vec<f64> = measurements.iter().map(|m| m.max_p_error).collect();
    assert!(
        p_errors.windows(2).all(|pair| pair[1] < pair[0]),
        "the pressure error must fall under refinement; got {p_errors:?}"
    );
    assert!(
        p_errors[0] < 0.12,
        "max pressure error {:.4e} on the coarsest mesh, on a field spanning \
         [-1, 1]. Got: {p_errors:?}",
        p_errors[0]
    );

    Ok(())
}

/// Discrete continuity must hold in *every* cell, not just the interior ones.
///
/// A converged SIMPLE solve satisfies the discrete divergence to solver
/// tolerance by construction — the pressure correction exists for no other
/// reason. Any cell where it does not is a cell whose continuity equation was
/// never assembled, and the ring/interior split is what makes that visible: a
/// global norm averages a bad outer ring away against a good interior.
///
/// Measured on this 32x32 solve: ring 2.50e-10, interior 2.53e-7. Before the
/// ring cells had a continuity equation those read 1.04e1 and 1.20e0 — the
/// threshold below fails that by six orders of magnitude.
///
/// The bound is 1e-5 rather than machine zero because what is left is the inner
/// Gauss-Seidel budget, not the discretisation. `pressure_correction_step`
/// takes at most 200 sweeps per outer iteration; raising that to 3000 drops the
/// interior figure to 1.90e-8, and it would keep falling. A stationary
/// iterative solver's truncation is a solver-effort question, so it is not
/// worth spending the wall-clock here to chase — but it is worth saying out
/// loud rather than reporting 2.53e-7 as if it were zero.
#[tokio::test]
async fn continuity_holds_on_the_outer_ring_as_well_as_the_interior() -> CfdResult<()> {
    let m = measure(32, ConvectionScheme::Upwind).await?;
    println!(
        "  n =  32   max |div u| ring = {:.6e}   interior = {:.6e}",
        m.max_div_ring, m.max_div_interior
    );

    assert!(
        m.max_div_interior < 1e-5,
        "interior divergence {:.6e} is not at solver tolerance",
        m.max_div_interior
    );
    assert!(
        m.max_div_ring < 1e-5,
        "outer-ring divergence {:.6e} against interior {:.6e}: the cells on the \
         ring are not having continuity enforced on them",
        m.max_div_ring,
        m.max_div_interior
    );

    Ok(())
}

/// The deferred-correction TVD scheme must lift the observed order toward 2.
///
/// The Stokes-limit measurement already pinned every non-convective operator
/// at second order, so with a second-order convective flux the whole
/// discretisation should approach 2 — and the error at every resolution must
/// be strictly below upwind's, since the schemes differ only in the
/// convective face values.
///
/// Measured (van Albada, 16 -> 32 -> 64): L2 velocity 1.325e-3, 4.406e-4,
/// 1.232e-4 — observed orders 1.59 and 1.84, climbing toward 2, against
/// upwind's 0.85 and 0.91 on the same meshes. The error is 27x to 84x below
/// upwind's at equal resolution. The shortfall from exactly 2 is the limiter
/// clipping at extrema plus the pure-upwind fallback at faces whose
/// far-upwind node lies outside the domain, both of which shrink with h —
/// which is why the rate climbs. The order bound of 1.5 sits under the worst
/// measured rate and still fails first-order upwind by a wide margin.
#[tokio::test]
async fn observed_order_reaches_two_with_tvd_convection() -> CfdResult<()> {
    let resolutions = [16usize, 32, 64];
    let mut tvd = Vec::new();
    let mut upwind = Vec::new();
    for &n in &resolutions {
        tvd.push(measure(n, ConvectionScheme::TvdVanAlbada).await?);
        upwind.push(measure(n, ConvectionScheme::Upwind).await?);
    }
    let errors: Vec<f64> = tvd.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}   TVD L2 = {:.6e}   order = {rate}   upwind L2 = {:.6e}   \
             max |p - p_exact| = {:.6e}",
            errors[i], upwind[i].l2_velocity, tvd[i].max_p_error
        );
    }

    for (i, m) in tvd.iter().enumerate() {
        assert!(
            m.l2_velocity < upwind[i].l2_velocity,
            "TVD error {:.4e} not below upwind {:.4e} at n = {}",
            m.l2_velocity,
            upwind[i].l2_velocity,
            resolutions[i]
        );
    }

    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 > 1.5,
            "refinement {} -> {}: observed order {rate:.3}, below the 1.59 and 1.84 \
             this scheme measures. Errors: {errors:?}",
            resolutions[i],
            resolutions[i + 1]
        );
        assert!(
            rate < 2.5,
            "refinement {} -> {}: observed order {rate:.3} — too good; suspect the \
             error measure. Errors: {errors:?}",
            resolutions[i],
            resolutions[i + 1]
        );
    }

    // TVD must not corrupt continuity or the pressure: divergence stays at
    // solver tolerance and the pressure error still falls under refinement.
    for m in &tvd {
        assert!(m.max_div_ring < 1e-5);
        assert!(m.max_div_interior < 1e-5);
    }
    let p_errors: Vec<f64> = tvd.iter().map(|m| m.max_p_error).collect();
    assert!(
        p_errors.windows(2).all(|pair| pair[1] < pair[0]),
        "the pressure error must fall under refinement; got {p_errors:?}"
    );

    Ok(())
}

/// Van Leer is the same construction with a different limiter; one resolution
/// pins it as implemented (beats upwind) without doubling the suite's runtime.
#[tokio::test]
async fn van_leer_limiter_also_beats_upwind() -> CfdResult<()> {
    let tvd = measure(32, ConvectionScheme::TvdVanLeer).await?;
    let upwind = measure(32, ConvectionScheme::Upwind).await?;
    println!(
        "  n =  32   van Leer L2 = {:.6e}   upwind L2 = {:.6e}",
        tvd.l2_velocity, upwind.l2_velocity
    );
    assert!(tvd.l2_velocity < upwind.l2_velocity);
    assert!(tvd.max_div_interior < 1e-5 && tvd.max_div_ring < 1e-5);
    Ok(())
}