//! Rung F2 of the Turek–Hron ladder: a rigid body MOVING through the fixed
//! grid — per-step mask rebuild, fresh cells, and the falsifier-3
//! measurement (fresh-cell pressure noise) of the geometry decision
//! (omni-cortex `docs/turek_hron_geometry_decision.md`).
//!
//! Two claims, in order:
//!
//! 1. **A stationary body run through the moving path is the static path
//!    to the bit.** The moving path rebuilds the mask every step and
//!    re-imposes ghost values from the previous corrected field; for a
//!    body that happens not to move, both are exactly what the static path
//!    holds, so nothing may differ.
//!
//! 2. **A circle translating through the steady manufactured field leaves
//!    the solution at the static-body error level.** The circle's surface
//!    carries the exact field as its velocity (a "phantom" surface), so
//!    the steady manufactured solution stays exact while the mask sweeps
//!    across the grid: velocity faces flip solid → fluid holding the ghost
//!    reconstruction the previous step left, fresh pressure cells are
//!    refilled from neighbours, and any fresh-cell pressure transient
//!    shows up directly against the KNOWN exact pressure. The measured
//!    time-maxima against the static steady-state levels (L2 u 8.489e-3,
//!    L2 p 2.22e-2 at n = 32, upwind) are the falsifier-3 numbers: spikes
//!    well above the static level would send the method to cut cells.

use rtx_cfd::solvers::incompressible::{
    EmbeddedBody, EmbeddedParameters, EmbeddedPisoSolver, FaceKind, FlowField, PoissonSolverKind,
};
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()
}

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

fn boundary_exact(x: f64, y: f64) -> (f64, f64) {
    let u = if x <= 0.0 || x >= 1.0 {
        0.0
    } else {
        u_exact(x, y)
    };
    let v = if y <= 0.0 || y >= 1.0 {
        0.0
    } else {
        v_exact(x, y)
    };
    (u, v)
}

fn solver(n: usize) -> CfdResult<EmbeddedPisoSolver> {
    let config = CfdConfig::new()
        .with_density(RHO)
        .with_viscosity(MU)
        .with_reference_velocity(1.0)
        .with_reference_length(1.0);
    let mut solver = EmbeddedPisoSolver::new(
        config,
        EmbeddedParameters {
            corrector_steps: 2,
            tolerance: 1e-8,
            poisson_solver: PoissonSolverKind::Multigrid,
            ..EmbeddedParameters::default()
        },
    )?;
    solver.set_momentum_source(|x, y, _| source(x, y));
    solver.set_boundary_velocity(|x, y, _| boundary_exact(x, y));
    let _ = n;
    Ok(solver)
}

fn exact_field(n: usize) -> CfdResult<FlowField> {
    let dx = 1.0 / n as f64;
    let mut field = FlowField::new(n, n, dx, dx)?;
    for j in 0..n {
        for i in 0..=n {
            field.u[(j, i)] = u_exact(i as f64 * dx, (j as f64 + 0.5) * dx);
        }
    }
    for j in 0..=n {
        for i in 0..n {
            field.v[(j, i)] = v_exact((i as f64 + 0.5) * dx, j as f64 * dx);
        }
    }
    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);
        }
    }
    for j in 0..n {
        field.u[(j, 0)] = boundary_exact(0.0, (j as f64 + 0.5) * dx).0;
        field.u[(j, n)] = boundary_exact(1.0, (j as f64 + 0.5) * dx).0;
    }
    for i in 0..n {
        field.v[(0, i)] = boundary_exact((i as f64 + 0.5) * dx, 0.0).1;
        field.v[(n, i)] = boundary_exact((i as f64 + 0.5) * dx, 1.0).1;
    }
    Ok(field)
}

fn time_step(n: usize) -> f64 {
    let dx = 1.0 / n as f64;
    let nu = MU / RHO;
    0.4 * (dx * dx / (4.0 * nu)).min(dx)
}

fn phantom_circle(
    cx: impl Fn(f64) -> f64 + Send + Sync + 'static,
    cy: impl Fn(f64) -> f64 + Send + Sync + 'static,
    r: f64,
) -> EmbeddedBody {
    EmbeddedBody::from_sdf(move |x, y, t| ((x - cx(t)).powi(2) + (y - cy(t)).powi(2)).sqrt() - r)
        .with_surface_velocity(|x, y, _| (u_exact(x, y), v_exact(x, y)))
}

/// Claim 1: stationary body, static path vs moving path, bit for bit.
#[tokio::test]
async fn a_stationary_body_through_the_moving_path_is_bit_identical() -> CfdResult<()> {
    let n = 24;
    let dt = time_step(n);
    let mut fixed = solver(n)?;
    fixed.set_body(phantom_circle(|_| 0.5, |_| 0.45, 0.2));
    let mut moving = solver(n)?;
    moving.set_moving_body(phantom_circle(|_| 0.5, |_| 0.45, 0.2));

    let mut a = exact_field(n)?;
    let mut b = exact_field(n)?;
    for _ in 0..100 {
        let ra = fixed.advance(&mut a, dt).await?;
        let rb = moving.advance(&mut b, dt).await?;
        assert_eq!(rb.fresh_cells, 0, "a stationary body produced fresh cells");
        assert_eq!(ra.ghost_correction, rb.ghost_correction);
    }
    let mut max_diff: f64 = 0.0;
    for (x, y) in a.u.iter().zip(b.u.iter()) {
        max_diff = max_diff.max((x - y).abs());
    }
    for (x, y) in a.v.iter().zip(b.v.iter()) {
        max_diff = max_diff.max((x - y).abs());
    }
    for (x, y) in a.p.iter().zip(b.p.iter()) {
        max_diff = max_diff.max((x - y).abs());
    }
    assert!(
        max_diff == 0.0,
        "moving path with a stationary body differs from the static path by {max_diff:.3e}"
    );
    Ok(())
}

/// The subiteration seam: snapshot the solver + clone the field mid-run of
/// a MOVING body, advance further (a discarded coupling candidate), then
/// restore and advance the same steps again — the re-run must be
/// bit-identical to a run that never diverted. This is what lets an FSI
/// coupling re-run one fluid step under updated interface geometry.
#[tokio::test]
async fn snapshot_restore_rerun_is_bit_identical() -> CfdResult<()> {
    let n = 24;
    let dt = time_step(n);
    let mover = || phantom_circle(|t| 0.42 + 0.30 * t, |t| 0.48 + 0.15 * t, 0.2);

    // Reference: an uninterrupted run of 30 steps.
    let mut reference = solver(n)?;
    reference.set_moving_body(mover());
    let mut ref_field = exact_field(n)?;
    for _ in 0..30 {
        reference.advance(&mut ref_field, dt).await?;
    }

    // Diverted run: 18 steps, snapshot, 12 steps of a discarded candidate,
    // restore, the real 12 steps.
    let mut solver_d = solver(n)?;
    solver_d.set_moving_body(mover());
    let mut field = exact_field(n)?;
    for _ in 0..18 {
        solver_d.advance(&mut field, dt).await?;
    }
    let saved_state = solver_d.snapshot();
    let saved_field = field.clone();
    for _ in 0..12 {
        solver_d.advance(&mut field, dt).await?; // discarded candidate
    }
    solver_d.restore(&saved_state);
    field = saved_field;
    let mut rerun_fresh = 0usize;
    for _ in 0..12 {
        rerun_fresh += solver_d.advance(&mut field, dt).await?.fresh_cells;
    }
    // Cells must actually flip in the re-run window, or the restore of the
    // mask was never exercised against a mask that changes.
    assert!(
        rerun_fresh > 0,
        "no cells flipped after the restore — the test is vacuous"
    );

    assert_eq!(
        solver_d.time().to_bits(),
        reference.time().to_bits(),
        "restored time diverges"
    );
    let mut max_diff: f64 = 0.0;
    for (x, y) in field.u.iter().zip(ref_field.u.iter()) {
        max_diff = max_diff.max((x - y).abs());
    }
    for (x, y) in field.v.iter().zip(ref_field.v.iter()) {
        max_diff = max_diff.max((x - y).abs());
    }
    for (x, y) in field.p.iter().zip(ref_field.p.iter()) {
        max_diff = max_diff.max((x - y).abs());
    }
    assert!(
        max_diff == 0.0,
        "restored re-run differs from the uninterrupted run by {max_diff:.3e}"
    );
    Ok(())
}

/// Mask hysteresis is inert for a body that does not move: the reference
/// classification and the exact classification agree at every cell (a
/// fluid cell has `phi > 0 > -band`, a solid cell `phi <= 0 < band`), so
/// a stationary body with any band is the static path to the bit.
#[tokio::test]
async fn mask_hysteresis_is_inert_for_a_stationary_body() -> CfdResult<()> {
    let n = 24;
    let dt = time_step(n);
    let mut fixed = solver(n)?;
    fixed.set_body(phantom_circle(|_| 0.5, |_| 0.45, 0.2));
    let mut sticky = solver(n)?;
    sticky.set_moving_body(phantom_circle(|_| 0.5, |_| 0.45, 0.2));
    sticky.set_mask_hysteresis(0.5);

    let mut a = exact_field(n)?;
    let mut b = exact_field(n)?;
    for _ in 0..100 {
        fixed.advance(&mut a, dt).await?;
        let rb = sticky.advance(&mut b, dt).await?;
        assert_eq!(rb.fresh_cells, 0, "a stationary body produced fresh cells");
    }
    let mut max_diff: f64 = 0.0;
    for (x, y) in a.u.iter().zip(b.u.iter()) {
        max_diff = max_diff.max((x - y).abs());
    }
    for (x, y) in a.v.iter().zip(b.v.iter()) {
        max_diff = max_diff.max((x - y).abs());
    }
    for (x, y) in a.p.iter().zip(b.p.iter()) {
        max_diff = max_diff.max((x - y).abs());
    }
    assert!(
        max_diff == 0.0,
        "hysteresis on a stationary body differs from the static path by {max_diff:.3e}"
    );
    Ok(())
}

/// Mask hysteresis delays each flip by the band: a translating circle at
/// constant velocity flips its first cell later by ~band / (speed * dt)
/// steps, flips no MORE cells than the plain rebuild over the same
/// traverse, and the delayed timeline is deterministic across re-runs.
#[tokio::test]
async fn mask_hysteresis_delays_flips_by_the_band() -> CfdResult<()> {
    let n = 24;
    let dt = time_step(n);
    let steps = 200;
    fn mover() -> EmbeddedBody {
        phantom_circle(|t| 0.42 + 0.30 * t, |t| 0.48 + 0.15 * t, 0.2)
    }

    let run = |band: f64| async move {
        let mut solver = solver(n)?;
        solver.set_moving_body(mover());
        solver.set_mask_hysteresis(band);
        let mut field = exact_field(n)?;
        let mut timeline = Vec::with_capacity(steps);
        for _ in 0..steps {
            timeline.push(solver.advance(&mut field, dt).await?.fresh_cells);
        }
        CfdResult::Ok(timeline)
    };
    let plain = run(0.0).await?;
    let sticky = run(0.5).await?;
    let rerun = run(0.5).await?;
    assert_eq!(
        sticky, rerun,
        "the sticky flip timeline is not deterministic"
    );

    let first = |t: &[usize]| t.iter().position(|&f| f > 0);
    let first_plain = first(&plain).expect("the plain traverse flips no cells — vacuous");
    let first_sticky = first(&sticky).expect("the sticky traverse flips no cells");
    let total_plain: usize = plain.iter().sum();
    let total_sticky: usize = sticky.iter().sum();
    println!(
        "  first fresh cell: plain step {first_plain}, band 0.5 step {first_sticky}; \
         totals over {steps} steps: plain {total_plain}, band 0.5 {total_sticky}"
    );
    assert!(
        first_sticky > first_plain,
        "the band did not delay the first flip (plain {first_plain}, sticky {first_sticky})"
    );
    assert!(
        total_sticky <= total_plain,
        "hysteresis flipped MORE cells ({total_sticky}) than the plain rebuild ({total_plain})"
    );
    assert!(
        total_sticky * 2 > total_plain,
        "hysteresis suppressed most flips ({total_sticky} of {total_plain}) — the band is \
         acting as a freeze, not a delay"
    );
    Ok(())
}

/// Claim 2: the translating phantom circle. Static steady-state baselines
/// at n = 32 (upwind, from `tests/embedded_mms.rs`): L2 u 8.489e-3,
/// L2 p 2.22e-2.
#[tokio::test]
async fn translating_circle_holds_the_manufactured_field() -> CfdResult<()> {
    let n = 32;
    let dt = time_step(n);
    let dx = 1.0 / n as f64;
    let steps = 300;

    let mut solver = solver(n)?;
    solver.set_moving_body(phantom_circle(
        |t| 0.42 + 0.30 * t,
        |t| 0.48 + 0.15 * t,
        0.2,
    ));
    let mut field = exact_field(n)?;
    solver.initialize(&mut field)?;

    let mut total_fresh = 0usize;
    let mut max_l2_u: f64 = 0.0;
    let mut max_l2_p: f64 = 0.0;
    let mut max_div: f64 = 0.0;
    let mut max_ghost_corr: f64 = 0.0;

    let mut max_residual: f64 = 0.0;
    for _step in 0..steps {
        let result = solver.advance(&mut field, dt).await?;
        total_fresh += result.fresh_cells;
        max_ghost_corr = max_ghost_corr.max(result.ghost_correction.abs());
        max_residual = max_residual.max(result.solver_result.final_residual);

        let mask = solver.mask().expect("mask");
        // L2 velocity error over the current fluid faces.
        let mut squared = 0.0;
        let mut volume = 0.0;
        for j in 0..n {
            for i in 1..n {
                if mask.u_kind(j, i) == FaceKind::Fluid {
                    let e = field.u[(j, i)] - u_exact(i as f64 * dx, (j as f64 + 0.5) * dx);
                    squared += e * e * dx * dx;
                    volume += dx * dx;
                }
            }
        }
        for j in 1..n {
            for i in 0..n {
                if mask.v_kind(j, i) == FaceKind::Fluid {
                    let e = field.v[(j, i)] - v_exact((i as f64 + 0.5) * dx, j as f64 * dx);
                    squared += e * e * dx * dx;
                    volume += dx * dx;
                }
            }
        }
        max_l2_u = max_l2_u.max((squared / volume).sqrt());

        // Mean-shifted L2 pressure error over the current fluid cells, and
        // the divergence.
        let mut diff_sum = 0.0;
        let mut cells = 0usize;
        for j in 0..n {
            for i in 0..n {
                if mask.is_fluid_cell(j, i) {
                    diff_sum +=
                        field.p[(j, i)] - p_exact((i as f64 + 0.5) * dx, (j as f64 + 0.5) * dx);
                    cells += 1;
                }
            }
        }
        let shift = diff_sum / cells as f64;
        let mut p_sq = 0.0;
        for j in 0..n {
            for i in 0..n {
                if mask.is_fluid_cell(j, i) {
                    let e = field.p[(j, i)]
                        - shift
                        - p_exact((i as f64 + 0.5) * dx, (j as f64 + 0.5) * dx);
                    p_sq += e * e;
                    // Bulk divergence: only cells whose four faces are all
                    // fluid unknowns. The end-of-step ghost re-imposition
                    // legitimately changes the PRESCRIBED fluxes of
                    // body-adjacent cells after the projection (the next
                    // projection honours them — the same one-step lag the
                    // static path has); the projection's own residual below
                    // is the continuity claim for those.
                    if mask.u_kind(j, i) == FaceKind::Fluid
                        && mask.u_kind(j, i + 1) == FaceKind::Fluid
                        && mask.v_kind(j, i) == FaceKind::Fluid
                        && mask.v_kind(j + 1, i) == FaceKind::Fluid
                    {
                        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());
                    }
                }
            }
        }
        max_l2_p = max_l2_p.max((p_sq / cells as f64).sqrt());
    }

    println!(
        "  {steps} steps, circle centre moved ({:.3}, {:.3}); fresh cells {total_fresh}; \
         max L2 u {max_l2_u:.4e} (static steady 8.489e-3, ratio {:.2}); \
         max L2 p {max_l2_p:.4e} (static steady 2.22e-2, ratio {:.2}); \
         max bulk |div u| {max_div:.2e}; max projection residual {max_residual:.2e}; \
         max ghost correction {max_ghost_corr:.2e}",
        0.30 * steps as f64 * dt,
        0.15 * steps as f64 * dt,
        max_l2_u / 8.489e-3,
        max_l2_p / 2.22e-2,
    );

    assert!(
        total_fresh > 20,
        "the circle should sweep cells fresh; got {total_fresh} — the test is vacuous"
    );
    assert!(
        max_div < 1e-5,
        "a bulk fluid cell is not divergence-free under motion: {max_div:.3e}"
    );
    assert!(
        max_residual < 1e-6,
        "the projection failed to converge during the sweep: residual {max_residual:.3e}"
    );
    // Falsifier 3: fresh-cell pressure transients must stay at the level of
    // the static discretisation error, not orders above it.
    assert!(
        max_l2_u < 2.0 * 8.489e-3,
        "velocity error under motion {max_l2_u:.3e} vs static steady 8.489e-3"
    );
    assert!(
        max_l2_p < 3.0 * 2.22e-2,
        "pressure error under motion {max_l2_p:.3e} vs static steady 2.22e-2 — fresh-cell \
         spikes; the geometry decision's falsifier 3 fires and cut cells are next"
    );
    Ok(())
}

/// The O(h) price of mask hysteresis, measured: the translating phantom
/// circle with a 0.25h band. The effective wall lags the true surface by
/// up to the band, so the error against the manufactured field must rise
/// above the no-hysteresis moving level — boundedly, at the
/// discretisation's own order, not as a blowup. The printed ratios are
/// the measurement; the asserts exclude a runaway.
#[tokio::test]
async fn translating_circle_with_hysteresis_pays_a_bounded_lag() -> CfdResult<()> {
    let n = 32;
    let dt = time_step(n);
    let dx = 1.0 / n as f64;
    let steps = 300;

    let mut solver = solver(n)?;
    solver.set_moving_body(phantom_circle(
        |t| 0.42 + 0.30 * t,
        |t| 0.48 + 0.15 * t,
        0.2,
    ));
    solver.set_mask_hysteresis(0.25);
    let mut field = exact_field(n)?;
    solver.initialize(&mut field)?;

    let mut total_fresh = 0usize;
    let mut max_l2_u: f64 = 0.0;
    let mut max_l2_p: f64 = 0.0;
    for _step in 0..steps {
        let result = solver.advance(&mut field, dt).await?;
        total_fresh += result.fresh_cells;

        let mask = solver.mask().expect("mask");
        let mut squared = 0.0;
        let mut volume = 0.0;
        for j in 0..n {
            for i in 1..n {
                if mask.u_kind(j, i) == FaceKind::Fluid {
                    let e = field.u[(j, i)] - u_exact(i as f64 * dx, (j as f64 + 0.5) * dx);
                    squared += e * e * dx * dx;
                    volume += dx * dx;
                }
            }
        }
        for j in 1..n {
            for i in 0..n {
                if mask.v_kind(j, i) == FaceKind::Fluid {
                    let e = field.v[(j, i)] - v_exact((i as f64 + 0.5) * dx, j as f64 * dx);
                    squared += e * e * dx * dx;
                    volume += dx * dx;
                }
            }
        }
        max_l2_u = max_l2_u.max((squared / volume).sqrt());

        let mut diff_sum = 0.0;
        let mut cells = 0usize;
        for j in 0..n {
            for i in 0..n {
                if mask.is_fluid_cell(j, i) {
                    diff_sum +=
                        field.p[(j, i)] - p_exact((i as f64 + 0.5) * dx, (j as f64 + 0.5) * dx);
                    cells += 1;
                }
            }
        }
        let shift = diff_sum / cells as f64;
        let mut p_sq = 0.0;
        for j in 0..n {
            for i in 0..n {
                if mask.is_fluid_cell(j, i) {
                    let e = field.p[(j, i)]
                        - shift
                        - p_exact((i as f64 + 0.5) * dx, (j as f64 + 0.5) * dx);
                    p_sq += e * e;
                }
            }
        }
        max_l2_p = max_l2_p.max((p_sq / cells as f64).sqrt());
    }

    println!(
        "  band 0.25h over {steps} steps: fresh cells {total_fresh}; \
         max L2 u {max_l2_u:.4e} (static steady 8.489e-3, ratio {:.2}); \
         max L2 p {max_l2_p:.4e} (static steady 2.22e-2, ratio {:.2})",
        max_l2_u / 8.489e-3,
        max_l2_p / 2.22e-2,
    );

    assert!(
        total_fresh > 20,
        "the circle should sweep cells fresh; got {total_fresh} — the test is vacuous"
    );
    // Measured 2026-08-26: ratios 1.17 (u) and 2.12 (p) — within a percent
    // of the no-hysteresis moving levels (1.16 / 2.11). The band is held to
    // the same bounds as the plain moving test.
    assert!(
        max_l2_u < 2.0 * 8.489e-3,
        "velocity error with a 0.25h band {max_l2_u:.3e} vs static steady 8.489e-3 — \
         the lag is not O(h)-bounded"
    );
    assert!(
        max_l2_p < 3.0 * 2.22e-2,
        "pressure error with a 0.25h band {max_l2_p:.3e} vs static steady 2.22e-2 — \
         the lag is not O(h)-bounded"
    );
    Ok(())
}