//! Turek–Hron CFD2 (steady, Re = 100) and CFD3 (periodic vortex shedding,
//! Re = 200) past the rigid cylinder + flag, on the embedded-boundary PISO
//! solver with the multigrid projection.
//!
//! Geometry, parameters and reference values from the FEATFLOW benchmark
//! tables (sourced 2026-08-20, omni-cortex
//! `docs/turek_hron_geometry_decision.md`); the body model and conventions
//! are `tests/turek_hron_cfd.rs`'s (flag extended into the cylinder, loads
//! by the surface-stress and control-volume routes):
//!
//! - CFD2: `U = 1`, Re = 100, steady. Reference (level 6): **drag 136.700,
//!   lift 10.5343**.
//! - CFD3: `U = 2`, Re = 200, periodic. Reference (level 4, dt 0.005):
//!   **drag 439.45 ± 5.6183, lift −11.893 ± 437.81, frequency 4.3956 Hz**,
//!   with the benchmark's inflow ramp `(1 − cos(pi t / 2)) / 2` for
//!   `t < 2 s`.
//!
//! Both cases run the inflow ramp (it is part of CFD3's definition and a
//! gentler start for CFD2's explicit march). CFD2 is settled the way CFD1
//! is: the control-volume drag stagnant to 1e-4 relative over 200 steps
//! after one flow-through time. CFD3 marches to `t = 9 s` and measures over
//! `t in [6, 9]` (~13 shedding periods): mean and amplitude as
//! `(max + min)/2 ± (max − min)/2` of the control-volume series, the
//! frequency from linearly-interpolated upward zero crossings of the lift
//! about its mean, and a periodicity check that the two halves of the
//! window agree on the lift amplitude.
//!
//! Measured (TVD van Albada + multigrid, dev profile; surface route
//! primary, control volume printed as the diagnostic — its central-
//! difference evaluation truncation grows with the convective flux and the
//! two routes differ ~15–25% here where they agreed to 0.6% at Re 20):
//!
//! | case | ny | surface drag | surface lift | f (Hz) | wall |
//! |------|----|--------------|--------------|--------|------|
//! | CFD2 | 41 | 119.9 ± 0.000 (−12.3%) | −3.4 | steady | 57 s |
//! | CFD2 | 62 | 121.4 ± 0.000 (−11.2%) | +30.2 | steady | 176 s |
//! | CFD3 | 41 | 409.0 ± 8.2 (−6.9%) | −184 ± 438.0 (amp +0.05%) | 4.2746 (−2.8%) | 184 s |
//! | CFD3 | 62 | 413.0 ± 11.9 (−6.0%) | +160 ± 555.6 (amp +27%) | 4.3400 (−1.3%) | 618 s |
//! | CFD2 | 82 | 122.6 ± 0.000 (−10.3%) | +8.4 | steady | 496 s |
//! | CFD3 | 82 | 394.2 ± 9.5 (−10.3%) | −2.6 ± 557.2 (amp +27%) | 4.3939 (**−0.04%**) | 1506 s |
//!
//! References: CFD2 drag 136.700, lift 10.5343; CFD3 drag 439.45 ± 5.62,
//! lift −11.893 ± 437.81, f 4.3956. What holds and what does not: the
//! shedding frequency converges cleanly (−2.8% → −1.3% → **−0.04%** at
//! h = 5 mm) and the lift MEAN collapses onto the reference (−184 → +160 →
//! −2.6 vs −11.9); the CFD2 control-volume drag converges (152.4 → 143.3 →
//! 139.4, +2.0% at 5 mm) while its surface drag sits ~−10% (the Re 100–200
//! boundary layer is ~5–10 mm — barely a cell); the CFD3 lift amplitude
//! reads +27% at both 6.6 and 5 mm, unconverged (the flag is 3 / 4 cells
//! thick, and the reference itself needed their level 4). Pre-asymptotic
//! numbers are recorded, not asserted. Suite defaults: CFD2 at ny = 62,
//! CFD3 at ny = 41 (its cost); `RTX_CFD2_NY` / `RTX_CFD3_NY` override.

use rtx_cfd::solvers::incompressible::{
    AleBoundaries, ConvectionScheme, 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 CFD2_U: f64 = 1.0;
const CFD2_REF_DRAG: f64 = 136.700;
const CFD2_REF_LIFT: f64 = 10.5343;

const CFD3_U: f64 = 2.0;
const CFD3_REF_DRAG_MEAN: f64 = 439.45;
const CFD3_REF_DRAG_AMP: f64 = 5.6183;
const CFD3_REF_LIFT_MEAN: f64 = -11.893;
const CFD3_REF_LIFT_AMP: f64 = 437.81;
const CFD3_REF_FREQUENCY: f64 = 4.3956;

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

/// The ramped parabolic inflow of the benchmark definition.
fn inflow(u_mean: f64, y: f64, t: f64) -> f64 {
    let ramp = if t < 2.0 {
        0.5 * (1.0 - (std::f64::consts::PI * t / 2.0).cos())
    } else {
        1.0
    };
    ramp * 1.5 * u_mean * y * (H - y) / (0.5 * H).powi(2)
}

struct Runner {
    solver: EmbeddedPisoSolver,
    field: FlowField,
    nx: usize,
    ny: usize,
    h: f64,
    dt: f64,
    mu: f64,
    cv: (usize, usize, usize, usize),
}

impl Runner {
    fn new(u_mean: f64, ny: usize) -> CfdResult<Self> {
        let h = H / ny as f64;
        let nx = (L / h).round() as usize;
        let mu = RHO * NU;
        // Combined explicit criterion with the blockage's local peak
        // (see turek_hron_cfd.rs for the failure that taught it).
        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,
            // Upwind's numerical viscosity (|u| h / 2 ~ 10x the physical nu
            // on these grids) suppressed CFD3's vortex shedding entirely:
            // the ny = 41 upwind run produced ONE lift zero-crossing in
            // three seconds. The limited scheme restores the physics.
            convection_scheme: ConvectionScheme::TvdVanAlbada,
        };
        let mut solver = EmbeddedPisoSolver::new(config, params)?;
        solver.set_boundary_velocity(move |x, y, t| {
            if x <= 0.0 {
                (inflow(u_mean, y, t), 0.0)
            } else {
                (0.0, 0.0)
            }
        });
        solver.set_body(body());

        // Start at rest: the ramp brings the inflow up from zero.
        let mut field = FlowField::new(nx, ny, h, h)?;
        solver.initialize(&mut field)?;

        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,
        );
        Ok(Self {
            solver,
            field,
            nx,
            ny,
            h,
            dt,
            mu,
            cv,
        })
    }

    fn cv_force(&self) -> (f64, f64) {
        self.solver.mask().unwrap().control_volume_force(
            &self.field.u,
            &self.field.v,
            &self.field.p,
            &self.field.u_old,
            &self.field.v_old,
            self.dt,
            RHO,
            self.mu,
            None,
            self.cv,
        )
    }

    fn surface_force(&self) -> rtx_cfd::solvers::incompressible::SurfaceForce {
        self.solver.mask().unwrap().surface_force(
            self.solver.body().unwrap(),
            &self.field.u,
            &self.field.v,
            &self.field.p,
            self.mu,
            self.solver.time(),
            0.5 * self.h,
        )
    }

    async fn step(&mut self) -> CfdResult<()> {
        self.solver.advance(&mut self.field, self.dt).await?;
        let umax = self
            .field
            .u
            .iter()
            .fold(0.0f64, |acc, &value| acc.max(value.abs()));
        assert!(
            umax.is_finite(),
            "velocity became non-finite at t = {:.3}",
            self.solver.time()
        );
        Ok(())
    }
}

/// One sampled series of both load routes.
struct Series {
    times: Vec<f64>,
    surface_drag: Vec<f64>,
    surface_lift: Vec<f64>,
    skipped_max: usize,
    cv_drag: Vec<f64>,
    cv_lift: Vec<f64>,
    steps: usize,
    seconds: f64,
}

/// March to `t_end`, sampling both load routes every 25 steps once
/// `t >= t_start`. With the TVD convection even the nominally steady CFD2
/// oscillates a little on coarse grids (the upwind run was steady only
/// because its numerical viscosity was ten times the physical one), so
/// every case is measured the same way: time statistics over a window,
/// never a single snapshot.
async fn run_sampled(u_mean: f64, ny: usize, t_start: f64, t_end: f64) -> CfdResult<Series> {
    let mut runner = Runner::new(u_mean, ny)?;
    let start = std::time::Instant::now();
    let mut series = Series {
        times: Vec::new(),
        surface_drag: Vec::new(),
        surface_lift: Vec::new(),
        skipped_max: 0,
        cv_drag: Vec::new(),
        cv_lift: Vec::new(),
        steps: 0,
        seconds: 0.0,
    };
    while runner.solver.time() < t_end {
        runner.step().await?;
        series.steps += 1;
        if series.steps % 25 == 0 && runner.solver.time() >= t_start {
            let surface = runner.surface_force();
            let (cx, cy) = runner.cv_force();
            series.times.push(runner.solver.time());
            series.surface_drag.push(surface.fx);
            series.surface_lift.push(surface.fy);
            series.skipped_max = series.skipped_max.max(surface.skipped);
            series.cv_drag.push(cx);
            series.cv_lift.push(cy);
        }
    }
    series.seconds = start.elapsed().as_secs_f64();
    assert!(series.times.len() > 50, "too few samples in the window");
    Ok(series)
}

/// Mid-range mean and half-range amplitude of a series.
fn mid_amp(series: &[f64]) -> (f64, f64) {
    let max = series.iter().copied().fold(f64::MIN, f64::max);
    let min = series.iter().copied().fold(f64::MAX, f64::min);
    (0.5 * (max + min), 0.5 * (max - min))
}

/// Frequency from linearly-interpolated upward zero crossings about the
/// mean; `None` with fewer than four crossings.
fn crossing_frequency(times: &[f64], series: &[f64]) -> Option<f64> {
    let (mean, _) = mid_amp(series);
    let mut crossings: Vec<f64> = Vec::new();
    for k in 1..series.len() {
        let (a, b) = (series[k - 1] - mean, series[k] - mean);
        if a < 0.0 && b >= 0.0 {
            let frac = a / (a - b);
            crossings.push(times[k - 1] + frac * (times[k] - times[k - 1]));
        }
    }
    (crossings.len() >= 4).then(|| {
        (crossings.len() - 1) as f64 / (crossings.last().unwrap() - crossings.first().unwrap())
    })
}

fn ny_list(var: &str, default: &[usize]) -> Vec<usize> {
    std::env::var(var)
        .ok()
        .map(|s| {
            s.split(',')
                .map(|t| t.trim().parse().expect("integer ny"))
                .collect()
        })
        .unwrap_or_else(|| default.to_vec())
}

#[tokio::test]
async fn cfd2_steady_drag_and_lift() -> CfdResult<()> {
    let resolutions = ny_list("RTX_CFD2_NY", &[62]);
    let rel = |a: f64, b: f64| ((a - b) / b).abs();
    for &ny in &resolutions {
        let r = run_sampled(CFD2_U, ny, 8.0, 10.0).await?;
        let (drag_s, drag_s_amp) = mid_amp(&r.surface_drag);
        let (lift_s, lift_s_amp) = mid_amp(&r.surface_lift);
        let (drag_c, _) = mid_amp(&r.cv_drag);
        let (lift_c, _) = mid_amp(&r.cv_lift);
        println!(
            "  CFD2 ny = {ny:3} (h = {:.4})   surface: drag {drag_s:.3} ± {drag_s_amp:.3} lift {lift_s:.3} ± {lift_s_amp:.3} (skipped ≤ {})   \
             control volume means: drag {drag_c:.3} lift {lift_c:.3}   [{} steps, {:.0} s]   reference drag {CFD2_REF_DRAG} lift {CFD2_REF_LIFT}",
            H / ny as f64,
            r.skipped_max,
            r.steps,
            r.seconds
        );
        assert!(
            rel(drag_s, CFD2_REF_DRAG) < 0.15,
            "ny = {ny}: surface drag mean {drag_s:.3} vs reference {CFD2_REF_DRAG}"
        );
    }
    Ok(())
}

#[tokio::test]
async fn cfd3_shedding_frequency_and_loads() -> CfdResult<()> {
    let resolutions = ny_list("RTX_CFD3_NY", &[41]);
    let rel = |a: f64, b: f64| ((a - b) / b).abs();
    for &ny in &resolutions {
        let r = run_sampled(CFD3_U, ny, 6.0, 9.0).await?;
        let (drag_mean, drag_amp) = mid_amp(&r.surface_drag);
        let (lift_mean, lift_amp) = mid_amp(&r.surface_lift);
        let frequency = crossing_frequency(&r.times, &r.surface_lift);
        let half = r.surface_lift.len() / 2;
        let (_, amp_first) = mid_amp(&r.surface_lift[..half]);
        let (_, amp_second) = mid_amp(&r.surface_lift[half..]);
        let (cv_drag_mean, cv_drag_amp) = mid_amp(&r.cv_drag);
        let (cv_lift_mean, cv_lift_amp) = mid_amp(&r.cv_lift);
        println!(
            "  CFD3 ny = {ny:3} (h = {:.4})   surface: drag {drag_mean:.2} ± {drag_amp:.2}, lift {lift_mean:.2} ± {lift_amp:.2}, f = {frequency:?} Hz (skipped ≤ {})   \
             CV: drag {cv_drag_mean:.2} ± {cv_drag_amp:.2}, lift {cv_lift_mean:.2} ± {cv_lift_amp:.2}   half-window lift amps {amp_first:.2}/{amp_second:.2}   \
             [{} steps, {:.0} s]   reference drag {CFD3_REF_DRAG_MEAN} ± {CFD3_REF_DRAG_AMP}, lift {CFD3_REF_LIFT_MEAN} ± {CFD3_REF_LIFT_AMP}, f {CFD3_REF_FREQUENCY}",
            H / ny as f64,
            r.skipped_max,
            r.steps,
            r.seconds
        );
        let frequency = frequency.expect("the wake must shed: fewer than four lift zero-crossings");
        assert!(
            (amp_first - amp_second).abs() < 0.15 * amp_second.max(1e-9),
            "ny = {ny}: lift amplitude still drifting: halves {amp_first:.2} / {amp_second:.2}"
        );
        assert!(
            rel(frequency, CFD3_REF_FREQUENCY) < 0.10,
            "ny = {ny}: shedding frequency {frequency:.4} vs reference {CFD3_REF_FREQUENCY}"
        );
        assert!(
            rel(drag_mean, CFD3_REF_DRAG_MEAN) < 0.15,
            "ny = {ny}: mean surface drag {drag_mean:.2} vs reference {CFD3_REF_DRAG_MEAN}"
        );
        assert!(
            rel(lift_amp, CFD3_REF_LIFT_AMP) < 0.35,
            "ny = {ny}: surface lift amplitude {lift_amp:.2} vs reference {CFD3_REF_LIFT_AMP}"
        );
    }
    Ok(())
}