//! Turek–Hron FSI3: the added-mass flapping flag — rung C3 of the ladder
//! (omni-cortex `docs/turek_hron_geometry_decision.md`).
//!
//! Re = 200 channel flow (`U = 2`) past the rigid cylinder with the
//! elastic flag at density ratio `rho_s / rho_f = 1` and `E = 5.6e6`:
//! the flag is as light as the fluid it displaces, so the fluid's
//! added mass is comparable to the structural mass — the regime in
//! which a staggered partitioned coupling is repulsive and the
//! quasi-Newton interface solver is genuinely needed per step (the
//! coupled piston benchmark's territory). Reference (FEATFLOW level 4,
//! dt = 0.0005): `ux(A) = −2.86 ± 2.70 mm [10.92 Hz]`,
//! `uy(A) = 1.45 ± 34.90 mm [5.46 Hz]`, drag `460.2 ± 27.47`,
//! lift `2.37 ± 153.75`.
//!
//! The march is FSI2's (`fsi2_harness::march`) with the FSI3 case
//! parameters; the rigid-flag phase is checked against this solver's
//! own CFD3 drag on the same geometry. The vacuum modal frequencies do
//! NOT identify this cycle the way they identified FSI2's: with the
//! fluid as heavy as the flag, the wetted frequency is set by the added
//! mass (the reference 5.46 Hz sits far below the vacuum mode 2 of
//! ~12 Hz that E/rho alone would give).
//!
//! # What is pinned
//!
//! Machinery invariants every run (conservation, finite state, coupling
//! bookkeeping). Physics bands arrive with the measurements — this
//! file's first commits carry the harness through the release response
//! only; the study record lives in omni-cortex `solver_status.md`.
//!
//! Environment knobs: `RTX_FSI3_*` — the same set as `RTX_FSI2_*` (see
//! `turek_hron_fsi2.rs`) plus `RTX_FSI3_OMEGA0`, `RTX_FSI3_TRACE`,
//! `RTX_FSI3_C1`, `RTX_FSI3_PREDICTOR`, `RTX_FSI3_QUIESCENT`; defaults:
//! ny 62, t_release 4, t_end 4.2, subcycle 2, tol 6e-5, rtol 1e-2,
//! maxsub 30, coupler `iqn`, reuse 2, omega0 0.05, C^1 interface,
//! kinematic predictor, quiescent release — every one of them a
//! measured necessity at unit density ratio (see the config comments).

mod fsi2_harness;

use fsi2_harness::FSI3;
use fsi2_harness::march::{MarchConfig, run_march};

// FEATFLOW level-4, dt 0.0005 reference values.
const REF_UY_MEAN: f64 = 1.45e-3;
const REF_UY_AMP: f64 = 34.90e-3;
const REF_UY_FREQ: f64 = 5.46;
const REF_UX_MEAN: f64 = -2.86e-3;
const REF_UX_AMP: f64 = 2.70e-3;
const REF_DRAG_MEAN: f64 = 460.2;
const REF_DRAG_AMP: f64 = 27.47;
const REF_LIFT_MEAN: f64 = 2.37;
const REF_LIFT_AMP: f64 = 153.75;

#[test]
fn fsi3_added_mass_flag() {
    // IQN by default: at density ratio 1 the per-step coupling gain is
    // high enough that a scalar relaxation is the wrong tool. The
    // tolerance floor is FSI2's loose-coupling value until FSI3's own is
    // measured; the increment-relative part does the work in the
    // meantime.
    let config = MarchConfig::from_env(
        "FSI3",
        MarchConfig {
            ny: 62,
            flag_nx: 35,
            t_release: 4.0,
            t_end: 4.2,
            // Subcycle 2, not FSI2's default 8: at unit density ratio the
            // interface noise floor rides with the MOTION as well as the
            // loads (mask flips per step ∝ velocity, each ~250 N — a
            // step-scale displacement on a light flag). Traced at
            // subcycle 8 through the release transient (flag at ~0.3
            // m/s): a step bouncing between 4e-4 and 1e-3 for twelve
            // passes, landing on the 13th by luck once and not at all
            // the next run. Subcycle 2 moves the interface 4x less per
            // step; the probe's subcycle-2 release map converges to
            // 9e-10.
            subcycle: 2,
            // The MEASURED FSI3 floor at this coupling (noise probe,
            // RTX_NOISE_CASE=fsi3, t = 4, ny = 62): flip-scan jumps
            // 3.6e-4 (12x FSI2's — 4x the dynamic pressure into a 10x
            // more compliant flag), the s = 8 release map stalling near
            // 1e-6 in 12 passes, and a traced step-1 stall at ~5e-5 from
            // ~250 N load flips; at subcycle 8 through the release
            // transient the floor rode up to ~1e-3 with the motion, which
            // is what moved the default to subcycle 2, budgeted at 3e-5
            // (wall-velocity noise 0.09 m/s, 4.5% of U). It must NOT be
            // FSI2's 2e-4 borrowed blindly with the constant-velocity
            // closure: a step accepted at one pass then left an interface
            // VELOCITY jump the incompressible fluid answered with an
            // impulsive added-mass load (1,600 N at release, 48,000 N
            // and 59 mm one step later, the flag's Newton dead by the
            // pass after) — the C^1 interface below removed that
            // channel.
            // 6e-5, not the 3e-5 the rest-state probe suggested: the
            // floor RIDES WITH THE MOTION, and through the cycle (flag
            // at ~1 m/s) a subcycle-2 step stalled at 1.9e-4 against a
            // 5x-floor window of 1.5e-4 at a turning point, where the
            // increment-relative part of the window collapses while the
            // flip noise does not (t = 6.18 s of the first study). The
            // window at 5 x 6e-5 = 3e-4 sits above the cycle's measured
            // stall level; wall-velocity noise 0.17 m/s against cycle
            // velocities of ~1 m/s.
            tol_floor: 6e-5,
            rtol: 1e-2,
            max_subiterations: 30,
            coupler: "iqn".into(),
            reuse: 2,
            smooth_in_h: 0.0,
            csv_path: None,
            // The blind first relaxation must not explode the structure
            // on a map whose per-pass gain is in the hundreds; 0.05 keeps
            // the exploratory pass bounded (the secant takes over from
            // pass 3, and the divergence verdict waits for it).
            initial_relaxation: 0.05,
            trace_steps: 0,
            // Continuous interface velocity across step boundaries: the
            // constant-velocity closure's jump at each step start drew an
            // impulsive added-mass load the unit-ratio flag could not
            // survive.
            c1_interface: true,
            // Extrapolate the converged interface motion: the
            // structure-alone predictor overshoots 2-5x at unit density
            // ratio and each step's first pass became a violent
            // excursion (traced at the t = 4 release).
            predictor: "kinematic".into(),
            // The structure-alone consistent initial acceleration ignores
            // the added mass; at unit density ratio it is wildly wrong
            // and Newmark carries it as a sign-alternating mode.
            quiescent_release: true,
        },
    );
    let MarchConfig {
        ny, flag_nx, t_end, ..
    } = config;
    let result = run_march(FSI3, &config);
    let w = result.window(2.0);

    println!(
        "  loads over the window: drag {:.2} ± {:.2} (ref {REF_DRAG_MEAN} ± {REF_DRAG_AMP}), \
         lift {:.2} ± {:.2} (ref {REF_LIFT_MEAN} ± {REF_LIFT_AMP})",
        w.drag_mid, w.drag_amp, w.lift_mid, w.lift_amp
    );
    println!(
        "  FSI3 (fluid ny = {ny}, flag {flag_nx}x2 Quad8, dt = {:.2e}, {}): coupled {} \
         steps in {:.0} s wall total; {:.1} subit/step (max {}); {} stalled steps, {} \
         history-reset retries (worst residual {:.2e}); worst conservation {:.2e}; \
         skipped samples {} (of which {} spike-clamped)\n  measured over [{:.1}, {t_end:.1}] s: \
         uy(A) = {:.4} ± {:.4} mm (ref {:.2} ± {:.2}), ux(A) = {:.4} ± {:.4} mm (ref {:.2} ± \
         {:.2}), f = {} Hz (ref {REF_UY_FREQ}); onset amp {:.3e} -> {:.3e} m",
        result.dt,
        config.coupler,
        result.coupled_steps,
        result.elapsed,
        result.mean_subiterations,
        result.max_subiterations,
        result.stalled_steps,
        result.retried_steps,
        result.worst_stall,
        result.worst_conservation,
        result.skipped,
        result.spiked,
        w.t_start,
        w.uy_mid * 1e3,
        w.uy_amp * 1e3,
        REF_UY_MEAN * 1e3,
        REF_UY_AMP * 1e3,
        w.ux_mid * 1e3,
        w.ux_amp * 1e3,
        REF_UX_MEAN * 1e3,
        REF_UX_AMP * 1e3,
        w.frequency.map_or("n/a".to_string(), |f| format!("{f:.3}")),
        w.amp_early,
        w.amp_late,
    );

    // Machinery invariants — asserted at every resolution.
    assert!(
        result.worst_conservation < 1e-10,
        "load transfer lost force: {:.3e}",
        result.worst_conservation
    );
    assert!(result.final_state_finite, "flag state went non-finite");
    assert!(
        result.mean_subiterations < 15.0,
        "coupling is grinding: {:.1} subiterations/step",
        result.mean_subiterations
    );
    assert!(
        result.stalled_steps * 5 < result.coupled_steps.max(1),
        "coupling stalled on {} of {} steps (worst residual {:.2e})",
        result.stalled_steps,
        result.coupled_steps,
        result.worst_stall
    );
}