//! Turek–Hron CFD1: steady laminar flow (Re = 20) past the rigid cylinder
//! with the rigid flag attached, on the embedded-boundary PISO solver.
//!
//! Geometry and parameters from the FEATFLOW benchmark definition (sourced
//! 2026-08-20, see omni-cortex `docs/turek_hron_geometry_decision.md`):
//! channel `[0, 2.5] x [0, 0.41]`, cylinder centre (0.2, 0.2) radius 0.05,
//! flag `[0.25, 0.6] x [0.19, 0.21]`, `rho = 1000`, `nu = 1e-3`, parabolic
//! inflow with mean `U = 0.2` (max 0.3), no-slip walls, outlet at the right.
//! Reference (level 6): **drag 14.2929, lift 1.11905** on cylinder + flag.
//!
//! The flag is modelled as `[0.20, 0.6] x [0.19, 0.21]`: its left 5 cm lie
//! inside the cylinder, which removes the two 1 mm fluid wedges the literal
//! corners (0.25, 0.19 ± 0.01) would leave between bar and circle — below
//! the grid scale here, and filled in the benchmark's own meshes.
//!
//! This is the first quantitative claim of the embedded solver against an
//! external reference, and the refinement study that was falsifier 4 of the
//! geometry decision. The projection runs the multigrid-preconditioned CG
//! solver; with the SOR projection the in-suite resolution was bounded at
//! h = 10 mm (~0.1 s/step, 609 s for the run below; an hour at 5 mm).
//!
//! Measured, dev profile, each run settled (control-volume drag stagnant to
//! `1e-4` relative over 200 steps after one flow-through time, 12.5 s):
//!
//! | ny | h (mm) | cells | dt (s)  | steps | wall s | surface drag / lift (skipped) | CV drag / lift | CV drag vs ref |
//! |----|--------|-------|---------|-------|--------|-------------------------------|----------------|----------------|
//! | 41 | 10.0   | 10250 | 1.92e-3 | 6500  | 21     | 15.7126 / 0.9355 (3)          | 15.6156 / 1.0785 | +9.25%       |
//! | 62 | 6.61   | 23436 | 1.10e-3 | 11400 | 64     | 15.3450 / 0.7818 (2)          | 15.2829 / 0.9195 | +6.93%       |
//! | 82 | 5.0    | 41000 | 7.35e-4 | 17000 | 173    | 15.3944 / 1.0178 (3)          | 15.0988 / 1.0673 | +5.64%       |
//!
//! Reference drag 14.2929, lift 1.11905. The SOR projection at ny = 41 gave
//! exactly the same four digits (surface 15.7126 / 0.9355, CV 15.6156 /
//! 1.0785) in 609 s — the same discrete system, a different inner solver —
//! and the multigrid run is required to reproduce it.
//!
//! What the three points say: the control-volume drag error falls
//! monotonically with h at an apparent order of 0.71 against the reference
//! (0.70 from the 10 -> 6.6 mm pair, 0.74 from 6.6 -> 5 mm; 0.57 from the
//! reference-free three-grid estimate, whose Richardson extrapolate is
//! 14.04, 1.8% under the reference). Sub-first-order is what a sharp
//! embedded boundary sampled on a Cartesian grid delivers for a blunt body
//! whose cut cells change with every h; the drag has not reached the
//! asymptotic range at 5 mm. The surface-integral route agrees with the
//! control volume to 0.6 / 0.4 / 2.0% but is not monotone (it samples the
//! pressure half a cell off the body and skips the junction samples), and
//! the lift — a 1 N difference of two 100 N-scale pressure integrals over a
//! flag that is 2 / 3 / 4 cells thick — is not monotone either (-3.6%,
//! -17.8%, -4.6% on the control-volume route). The assertions are the
//! measured bands: routes within 3%, CV drag error strictly decreasing with
//! apparent order above 0.5, the finest CV drag within 7% and CV lift
//! within 10% of the reference, and the ny = 41 run reproducing the SOR
//! loads to four digits. `RTX_CFD1_NY=41,62` (comma list) overrides the
//! resolution list for studies; the reference bands apply only when the
//! finest grid is at least ny = 82.

use rtx_cfd::solvers::incompressible::{
    AleBoundaries, EmbeddedBody, EmbeddedParameters, EmbeddedPisoSolver, FlowField,
    PoissonSolverKind, SideBoundary,
};
use rtx_cfd::{CfdConfig, CfdResult};

const L: f64 = 2.5;
const H: f64 = 0.41;
const RHO: f64 = 1000.0;
const NU: f64 = 1e-3;
const U_MEAN: f64 = 0.2;
const REF_DRAG: f64 = 14.2929;
const REF_LIFT: f64 = 1.11905;

fn inflow(y: f64) -> f64 {
    1.5 * U_MEAN * y * (H - y) / (0.5 * H).powi(2)
}

fn body() -> EmbeddedBody {
    EmbeddedBody::union(
        EmbeddedBody::circle(0.2, 0.2, 0.05),
        EmbeddedBody::rectangle(0.20, 0.19, 0.6, 0.21),
    )
}

struct Cfd1 {
    drag_surface: f64,
    lift_surface: f64,
    skipped: usize,
    drag_cv: f64,
    lift_cv: f64,
    steps: usize,
    seconds: f64,
}

/// March CFD1 to a steady state on a grid of `ny` cells across the channel.
/// Steady means the control-volume drag has stopped moving: its relative
/// change over the last 200 steps below `1e-4`, after at least one flow-
/// through time — "the answer stopped moving", not a residual.
async fn run_cfd1(ny: usize) -> CfdResult<Cfd1> {
    let h = H / ny as f64;
    let nx = (L / h).round() as usize;
    let mu = RHO * NU;
    // Explicit predictor: the COMBINED criterion — convective Courant numbers
    // in both directions plus the diffusion number must stay below one —
    // with the local peak velocity taken as 1.5x the inflow peak for the
    // 24% blockage. (Taking 0.4 of the smaller single limit, as the MMS
    // tests do, went NaN here at t ~ 7 s: both limits are active at once.)
    let u_peak = 1.5 * 1.5 * U_MEAN;
    let dt = 0.25 / (2.0 * u_peak / h + 4.0 * NU / (h * h));

    let config = CfdConfig::new()
        .with_density(RHO)
        .with_viscosity(mu)
        .with_reference_velocity(U_MEAN)
        .with_reference_length(0.1);
    let params = EmbeddedParameters {
        corrector_steps: 2,
        tolerance: 1e-7,
        boundaries: AleBoundaries {
            left: SideBoundary::Velocity,
            right: SideBoundary::PressureOutlet,
            bottom: SideBoundary::Velocity,
            top: SideBoundary::Velocity,
        },
        poisson_solver: PoissonSolverKind::Multigrid,
        poisson_precision: rtx_cfd::solvers::incompressible::MgPrecision::F64,
        ..EmbeddedParameters::default()
    };
    let mut solver = EmbeddedPisoSolver::new(config, params)?;
    solver.set_boundary_velocity(|x, y, _| {
        if x <= 0.0 {
            (inflow(y), 0.0)
        } else {
            (0.0, 0.0)
        }
    });
    solver.set_body(body());

    let mut field = FlowField::new(nx, ny, h, h)?;
    // Start from the inflow profile everywhere (the body's faces are
    // overwritten by the mask at initialisation).
    for j in 0..ny {
        let u0 = inflow((j as f64 + 0.5) * h);
        for i in 0..=nx {
            field.u[(j, i)] = u0;
        }
    }
    solver.initialize(&mut field)?;

    // Control volume for the momentum balance: whole cells, in the fluid
    // on its boundary, enclosing cylinder and flag.
    let cv = (
        (0.10 / h).round() as usize,
        (0.75 / h).round() as usize,
        (0.05 / h).round() as usize,
        (0.36 / h).round() as usize,
    );
    let cv_force =
        |field: &FlowField, mask: &rtx_cfd::solvers::incompressible::EmbeddedMask, dt: f64| {
            mask.control_volume_force(
                &field.u,
                &field.v,
                &field.p,
                &field.u_old,
                &field.v_old,
                dt,
                RHO,
                mu,
                None,
                cv,
            )
        };

    let start = std::time::Instant::now();
    let flow_through = L / U_MEAN;
    let min_steps = (flow_through / dt).ceil() as usize;
    let mut history: Vec<f64> = Vec::new();
    let mut steps = 0;
    loop {
        let result = solver.advance(&mut field, dt).await?;
        steps += 1;
        if steps % 50 == 0 {
            let (fx, _) = cv_force(&field, solver.mask().unwrap(), dt);
            history.push(fx);
            // Diagnostics: where is the velocity largest, did the projection
            // converge, how big was the ghost correction.
            let (mut umax, mut at) = (0.0f64, (0usize, 0usize));
            for j in 0..ny {
                for i in 0..=nx {
                    let a = field.u[(j, i)].abs();
                    if a > umax {
                        umax = a;
                        at = (j, i);
                    }
                }
            }
            if steps % 250 == 0 || umax > 3.0 * 1.5 * U_MEAN || !umax.is_finite() {
                println!(
                    "    ny = {ny}: step {steps} t = {:.2} s  drag_cv = {fx:.4}  max|u| = {umax:.4} at (x={:.3}, y={:.3})  \
                     projection: converged {} residual {:.2e} passes {}  ghost corr {:.2e}  [{:.0} s wall]",
                    solver.time(),
                    at.1 as f64 * h,
                    (at.0 as f64 + 0.5) * h,
                    result.solver_result.converged,
                    result.solver_result.final_residual,
                    result.corrector_steps_performed,
                    result.ghost_correction,
                    start.elapsed().as_secs_f64()
                );
            }
            assert!(
                umax.is_finite(),
                "velocity became non-finite at step {steps}"
            );
            if steps >= min_steps && history.len() > 4 {
                let now = history[history.len() - 1];
                let then = history[history.len() - 5];
                if ((now - then) / now).abs() < 1e-4 {
                    break;
                }
            }
        }
        assert!(steps < 400_000, "CFD1 at ny = {ny} did not settle");
    }
    let seconds = start.elapsed().as_secs_f64();

    let mask = solver.mask().unwrap();
    let surface = mask.surface_force(
        solver.body().unwrap(),
        &field.u,
        &field.v,
        &field.p,
        mu,
        solver.time(),
        0.5 * h,
    );
    let (drag_cv, lift_cv) = cv_force(&field, mask, dt);
    Ok(Cfd1 {
        drag_surface: surface.fx,
        lift_surface: surface.fy,
        skipped: surface.skipped,
        drag_cv,
        lift_cv,
        steps,
        seconds,
    })
}

#[tokio::test]
async fn cfd1_drag_and_lift_against_the_featflow_reference() -> CfdResult<()> {
    // h = 10, 6.6, 5 mm: 21 + 64 + 173 s in the dev profile, sequential.
    // `RTX_CFD1_NY` (comma-separated ny list) overrides for studies.
    let resolutions: Vec<usize> = std::env::var("RTX_CFD1_NY").ok().map_or_else(
        || vec![41usize, 62, 82],
        |list| {
            list.split(',')
                .map(|t| {
                    t.trim()
                        .parse()
                        .expect("RTX_CFD1_NY: comma-separated ny list")
                })
                .collect()
        },
    );
    let mut results = Vec::new();
    for &ny in &resolutions {
        let r = run_cfd1(ny).await?;
        println!(
            "  ny = {ny:3} (h = {:.4})   surface: drag {:.4} lift {:.4} (skipped {})   \
             control volume: drag {:.4} lift {:.4}   [{} steps, {:.0} s]   reference drag {REF_DRAG} lift {REF_LIFT}",
            H / ny as f64,
            r.drag_surface,
            r.lift_surface,
            r.skipped,
            r.drag_cv,
            r.lift_cv,
            r.steps,
            r.seconds
        );
        results.push(r);
    }
    let rel = |a: f64, b: f64| ((a - b) / b).abs();

    // The two load routes agree at every resolution (measured 0.6 / 0.4 /
    // 2.0%): a surface integral that samples the wrong side of the body or
    // a control volume that drops a flux term moves one route and not the
    // other.
    for (ny, r) in resolutions.iter().zip(&results) {
        assert!(
            rel(r.drag_surface, r.drag_cv) < 3e-2,
            "ny = {ny}: the two drag routes disagree: surface {:.4} vs control volume {:.4}",
            r.drag_surface,
            r.drag_cv
        );
    }

    // The same discrete system, a different inner solver: the settled
    // loads at ny = 41 must reproduce the SOR-projection values (see the
    // module docs) to four significant digits. A projection solving a
    // different system — wrong coefficient on the outlet column, wrong
    // anchor, a stop that is not the true residual — moves the drag far
    // more.
    const SOR_DRAG_CV: f64 = 15.6156;
    const SOR_DRAG_SURFACE: f64 = 15.7126;
    if let Some(coarse) = resolutions
        .iter()
        .position(|&ny| ny == 41)
        .map(|k| &results[k])
    {
        assert!(
            rel(coarse.drag_cv, SOR_DRAG_CV) < 5e-4
                && rel(coarse.drag_surface, SOR_DRAG_SURFACE) < 5e-4,
            "multigrid projection does not reproduce the SOR-projection loads: control volume {:.4} \
             vs {SOR_DRAG_CV}, surface {:.4} vs {SOR_DRAG_SURFACE}",
            coarse.drag_cv,
            coarse.drag_surface
        );
        // Cost: SOR took 609 s for 6500 steps (0.094 s/step) on the
        // development machine in the dev profile. Wall time is machine- and
        // profile-bound, so it is reported rather than asserted; the
        // multigrid run measured 0.003 s/step.
        println!(
            "  ny = 41 multigrid projection: {:.4} s/step ({} steps, {:.0} s); SOR baseline \
             0.0937 s/step (6500 steps, 609 s)",
            coarse.seconds / coarse.steps as f64,
            coarse.steps,
            coarse.seconds
        );
    }

    // Refinement: the control-volume drag error against the reference
    // falls strictly with h (measured +9.25%, +6.93%, +5.64%), at an
    // apparent order above 0.5 between the coarsest and finest grids
    // (measured 0.71). A discretisation that is not converging — a ghost
    // correction with the wrong sign, a body sampled on the wrong side —
    // keeps the error flat or growing.
    let errors: Vec<f64> = results.iter().map(|r| rel(r.drag_cv, REF_DRAG)).collect();
    for w in errors.windows(2) {
        assert!(
            w[1] < w[0],
            "control-volume drag error is not falling with h: {errors:?} (resolutions {resolutions:?})"
        );
    }
    if resolutions.len() > 1 {
        let (n0, n1) = (
            resolutions[0] as f64,
            resolutions[resolutions.len() - 1] as f64,
        );
        let order = (errors[0] / errors[errors.len() - 1]).ln() / (n1 / n0).ln();
        println!(
            "  control-volume drag error vs reference: {:?}  apparent order {order:.2} \
             (ny {n0} -> {n1})",
            errors
                .iter()
                .map(|e| format!("{e:+.4}"))
                .collect::<Vec<_>>()
        );
        assert!(
            order > 0.5,
            "apparent order of the control-volume drag error is {order:.2} (errors {errors:?})"
        );
    }

    // The finest grid against the reference: the measured bands at h = 5 mm
    // (CV drag +5.64%, CV lift -4.62%), applied only when the study reaches
    // that grid.
    if let Some((ny, fine)) = resolutions.iter().zip(&results).next_back() {
        if *ny >= 82 {
            assert!(
                rel(fine.drag_cv, REF_DRAG) < 0.07,
                "ny = {ny}: control-volume drag {:.4} vs reference {REF_DRAG}",
                fine.drag_cv
            );
            assert!(
                rel(fine.lift_cv, REF_LIFT) < 0.10,
                "ny = {ny}: control-volume lift {:.4} vs reference {REF_LIFT}",
                fine.lift_cv
            );
            // Surface-route lift at the finest grid: measured -9.05% at
            // ny = 82 (the flag is four cells thick; lift is a ~1 N
            // difference of ~100 N-scale integrals) — the measured band.
            assert!(
                rel(fine.lift_surface, REF_LIFT) < 0.12,
                "ny = {ny}: surface-route lift {:.4} vs reference {REF_LIFT}",
                fine.lift_surface
            );
        }
    }
    Ok(())
}