//! R8-a: the first coupled 3D FSI — the embedded3 cut-cell fluid (device
//! path) coupled to the 2D Turek–Hron flag (rtx-fea 35×2 Quad8 SVK, total
//! Lagrangian, Newmark γ 0.7) under a SPAN-UNIFORM deformation: the
//! structure's centreline drives the 3D body's polyline each coupled step,
//! the fluid's operator-route wall load, integrated over the span per unit
//! width, loads the structure's wetted nodes. Partitioned: per coupled
//! step the fluid re-runs the same step from a device snapshot for each
//! Aitken subiteration (the overset harness's pattern; `DeviceStep::
//! snapshot / restore`).
//!
//! The body is the embedded flag test's capsule (a semicircular tip, apex
//! on A): its 2D counterpart is the overset's SEMICIRCLE line (P5-2 ny 62:
//! 94.8 mm at 1.914 Hz), not the flat/1.25 mm-corner reference line.
//!
//! Knobs `RTX_E3FSI_*`: `NY` (62), `NZ` (4 = the periodic slab; 0 = the
//! full 0.41 m duct with slip sides), `T_RIGID` (3.0 s of rigid flag),
//! `T_END` (13.0), `RIGID_ONLY` (1 = stop after the rigid phase: target 1),
//! `RTOL` (1e-3), `FLOOR` (1.5e-7 on the centreline vector), `MAX_SUBIT`
//! (12), `STALL_ACCEPT` (5), `GAMMA` (0.7), `SPEED` (3.0 m/s, the band's
//! surface-speed bound), `CSV` (per-step series), `TRACE` (steps whose
//! passes are printed); R8-e: `KIN` (`centre` = the mid-plane nodes drive
//! the fluid's centreline; `surface` = the mean of the bottom and top edge
//! nodes, the nodes the loads land on — the transfer-consistent pair that
//! removes the small-dt instability), `LINE_CSV` (the committed centreline
//! and the per-station mid/bottom/top y displacement per step).
//!
//! `RTX_E3FSI_NY=62 RTX_E3FSI_CSV=<path> RTX_CUDA_ARCH=sm_120 cargo test --release -p rtx-fsi \
//!  --features cuda --test fsi2_embedded3 -- --ignored --nocapture`
#![cfg(feature = "cuda")]

#[path = "fsi2_embedded3/fluid.rs"]
mod fluid;
#[path = "fsi2_harness/mod.rs"]
mod fsi2_harness;
#[path = "fsi2_embedded3/state.rs"]
mod state;

use std::cell::RefCell;
use std::io::Write as _;

use fluid::{CX, CY, Contribution, E3Fluid, HALF, Line, R_CYL};
use fsi2_harness::{FSI2, Interface, clamp_left, flag_mesh, median, mid_amp};
use nalgebra::{DVector, Vector3};
use rtx_fea::analysis::{
    AnalysisConfig, ConvergenceCriteria, DynamicState, NonlinearDynamicAnalysis,
};
use rtx_fea::materials::{LinearElastic, MaterialDatabase};
use rtx_fea::mesh::{MaterialId, Mesh, NodeId};
use rtx_fsi::Subiterated;

pub fn env_f(name: &str, default: f64) -> f64 {
    std::env::var(name)
        .ok()
        .and_then(|v| v.parse().ok())
        .unwrap_or(default)
}

const X0: f64 = 0.25;
/// Centreline stations: the element corners at x 0.25, 0.26, …, 0.59 (the
/// last 10 mm is the capsule's cap).
const STATIONS: usize = 35;

/// R8-h: the flat tip (`RTX_E3FSI_TIP=flat`, corner radius
/// `RTX_E3FSI_TIP_CORNER`, default 0.00125 m): the centreline gains the
/// tip node A (x 0.60) as its last station and the body's tip is flat
/// through it (`DeviceSdf::tip_corner`); unset = the capsule.
pub fn flat_tip() -> Option<f64> {
    match std::env::var("RTX_E3FSI_TIP").as_deref() {
        Err(_) | Ok("capsule") => None,
        Ok("flat") => {
            let rc = env_f("RTX_E3FSI_TIP_CORNER", 0.00125);
            assert!((0.0..=HALF).contains(&rc), "RTX_E3FSI_TIP_CORNER {rc}");
            Some(rc)
        }
        Ok(v) => panic!("RTX_E3FSI_TIP={v}: flat or capsule"),
    }
}

/// The centreline's stations (35, or 36 with the flat tip's A).
fn stations() -> usize {
    STATIONS + usize::from(flat_tip().is_some())
}

/// The flag's structure-side bookkeeping: the centreline nodes, the wetted
/// edges with their reference coordinates, the tip node A.
struct Flag {
    centre: Vec<NodeId>,
    interface: Interface,
    /// (reference x, wetted index or None for the clamp corner), ascending.
    bottom: Vec<(f64, Option<usize>)>,
    top: Vec<(f64, Option<usize>)>,
    /// (reference y, wetted index), ascending (corners included).
    tip: Vec<(f64, usize)>,
    a_node: NodeId,
    /// R8-e: the bottom and top edge nodes at each station's x (the wetted
    /// surface the loads land on).
    surface: Vec<[NodeId; 2]>,
}

impl Flag {
    fn build(mesh: &Mesh) -> Self {
        let find = |x: f64, y: f64| -> NodeId {
            *mesh
                .nodes
                .iter()
                .find(|(_, n)| {
                    let p = n.position();
                    (p.x - x).abs() < 1e-9 && (p.y - y).abs() < 1e-9
                })
                .expect("node")
                .0
        };
        let centre = (0..stations())
            .map(|k| find(X0 + 0.01 * k as f64, 0.2))
            .collect();
        let interface = Interface::build(mesh);
        let edge = |idx: &[usize]| -> Vec<(f64, Option<usize>)> {
            let mut v: Vec<(f64, Option<usize>)> = vec![(X0, None)];
            v.extend(idx.iter().map(|&k| (interface.reference[k].0, Some(k))));
            v.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap());
            v
        };
        let bottom = edge(&interface.bottom);
        let top = edge(&interface.top);
        let mut tip: Vec<(f64, usize)> = interface
            .tip
            .iter()
            .map(|&k| (interface.reference[k].1, k))
            .collect();
        tip.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap());
        Self {
            centre,
            interface,
            bottom,
            top,
            tip,
            a_node: find(0.6, 0.2),
            // Over stations(): with the flat tip, station A (x 0.60) takes the mean of the tip's two corner nodes (merge of R8-e with R8-h).
            surface: (0..stations())
                .map(|k| {
                    let x = X0 + 0.01 * k as f64;
                    [find(x, 0.19), find(x, 0.21)]
                })
                .collect(),
        }
    }
}

/// Linear split of `value` at `s` over sorted stations; the share of a
/// `None` station (the clamp) is dropped.
fn split<T: Copy>(stations: &[(f64, T)], s: f64, value: f64, mut add: impl FnMut(T, f64)) {
    let n = stations.len();
    let s = s.clamp(stations[0].0, stations[n - 1].0);
    let m = stations.partition_point(|st| st.0 < s).clamp(1, n - 1);
    let (a, b) = (stations[m - 1], stations[m]);
    let u = if b.0 > a.0 {
        (s - a.0) / (b.0 - a.0)
    } else {
        0.0
    };
    add(a.1, (1.0 - u) * value);
    add(b.1, u * value);
}

/// The operator route's contributions onto the wetted nodes (per unit
/// span): each is placed by its closest point on the fluid's centreline —
/// the top or bottom edge at the same arc position, or the tip edge by its
/// lateral offset beyond the last station. Returns the nodal loads, the
/// flag's (fx, fy) and the cylinder's (fx, fy) per span, and the flag's fy
/// by part (pressure, shear, diffusive and convective exchange).
/// `RTX_E3FSI_PARTS` (bit mask, default 15 = all) keeps parts off the
/// structure (a diagnostic; the reported loads keep every part).
fn distribute(
    flag: &Flag,
    line: &Line,
    contributions: &[Contribution],
    width: f64,
) -> (Vec<(NodeId, Vector3<f64>)>, [f64; 2], [f64; 2], [f64; 4]) {
    let nw = flag.interface.wetted.len();
    let mut f = vec![[0.0f64; 2]; nw];
    let (mut on_flag, mut on_cyl) = ([0.0f64; 2], [0.0f64; 2]);
    let mut parts_y = [0.0f64; 4];
    let pts = &line.pts;
    let mut cum = vec![0.0; pts.len()];
    for m in 1..pts.len() {
        cum[m] = cum[m - 1]
            + ((pts[m][0] - pts[m - 1][0]).powi(2) + (pts[m][1] - pts[m - 1][1]).powi(2)).sqrt();
    }
    let parts_on = env_f("RTX_E3FSI_PARTS", 15.0) as usize;
    for &(pos, c, part, v) in contributions {
        if c > 1 {
            continue;
        }
        let v = v / width;
        let (x, y) = (pos[0], pos[1]);
        let mut best = (f64::INFINITY, 0usize, 0.0f64);
        for m in 0..pts.len() - 1 {
            let (a, b) = (pts[m], pts[m + 1]);
            let (ex, ey) = (b[0] - a[0], b[1] - a[1]);
            let l2 = ex * ex + ey * ey;
            let u = (((x - a[0]) * ex + (y - a[1]) * ey) / l2).clamp(0.0, 1.0);
            let d = ((x - a[0] - u * ex).powi(2) + (y - a[1] - u * ey).powi(2)).sqrt();
            if d < best.0 {
                best = (d, m, u);
            }
        }
        let d_cyl = ((x - CX).powi(2) + (y - CY).powi(2)).sqrt() - R_CYL;
        if d_cyl < best.0 - HALF {
            on_cyl[c] += v;
            continue;
        }
        on_flag[c] += v;
        if c == 1 {
            parts_y[part] += v;
        }
        if parts_on & (1 << part) == 0 {
            continue;
        }
        let (_, m, u) = best;
        let (a, b) = (pts[m], pts[m + 1]);
        let (ex, ey) = (b[0] - a[0], b[1] - a[1]);
        let len = (ex * ex + ey * ey).sqrt();
        let (tx, ty) = (ex / len, ey / len);
        let (px, py) = (a[0] + u * ex, a[1] + u * ey);
        // Lateral offset: + on the upper side of the centreline.
        let eta = tx * (y - py) - ty * (x - px);
        let along = tx * (x - px) + ty * (y - py);
        let mut add = |k: usize, w: f64| f[k][c] += w;
        if m + 2 == pts.len() && u >= 1.0 && along > 0.0 {
            let yr = 0.2 + eta.clamp(-HALF, HALF);
            split(&flag.tip, yr, v, &mut add);
        } else {
            let xr = X0 + cum[m] + u * len;
            let edge = if eta >= 0.0 { &flag.top } else { &flag.bottom };
            split(edge, xr, v, |k, w| {
                if let Some(k) = k {
                    add(k, w);
                }
            });
        }
    }
    let nodal = flag
        .interface
        .wetted
        .iter()
        .zip(&f)
        .map(|(&id, v)| (id, Vector3::new(v[0], v[1], 0.0)))
        .collect();
    (nodal, on_flag, on_cyl, parts_y)
}

/// The fluid's centreline from the structure's centreline displacement `c`
/// (2 per station) with velocities `(c − c_prev) / dt`.
fn line_of(t: f64, c: &[f64], c_prev: &[f64], dt: f64) -> Line {
    let pts = (0..stations())
        .map(|k| [X0 + 0.01 * k as f64 + c[2 * k], 0.2 + c[2 * k + 1]])
        .collect();
    let vel = (0..stations())
        .map(|k| {
            [
                (c[2 * k] - c_prev[2 * k]) / dt,
                (c[2 * k + 1] - c_prev[2 * k + 1]) / dt,
            ]
        })
        .collect();
    Line { t, pts, vel }
}

#[test]
#[ignore = "R8-a: the coupled FSI2 on embedded3 (GPU; minutes to hours)"]
fn fsi2_on_embedded3() {
    let ny = env_f("RTX_E3FSI_NY", 62.0) as usize;
    let nz = env_f("RTX_E3FSI_NZ", 4.0) as usize;
    let t_rigid = env_f("RTX_E3FSI_T_RIGID", 3.0);
    let t_end = env_f("RTX_E3FSI_T_END", 13.0);
    let rigid_only = env_f("RTX_E3FSI_RIGID_ONLY", 0.0) > 0.5;
    let rtol = env_f("RTX_E3FSI_RTOL", 1e-3);
    let floor = env_f("RTX_E3FSI_FLOOR", 1.5e-7);
    let max_subit = env_f("RTX_E3FSI_MAX_SUBIT", 12.0) as usize;
    let stall_accept = env_f("RTX_E3FSI_STALL_ACCEPT", 5.0);
    let gamma = env_f("RTX_E3FSI_GAMMA", 0.7);
    let speed = env_f("RTX_E3FSI_SPEED", 3.0);
    let trace = env_f("RTX_E3FSI_TRACE", 0.0) as usize;
    let csv_path = std::env::var("RTX_E3FSI_CSV").ok();
    let restore_check = env_f("RTX_E3FSI_RESTORE_CHECK", 0.0) as usize;
    let case = FSI2;

    let mesh = flag_mesh(35, 2);
    let flag_geo = Flag::build(&mesh);
    let zero_c = vec![0.0; 2 * stations()];
    let rest = line_of(0.0, &zero_c, &zero_c, 1.0);
    // `RTX_E3FSI_LOAD=<dir>`: continue from a saved coupled state (extruded
    // onto the full duct when the saved nz differs); `RTX_E3FSI_SAVE=<dir>`
    // saves the state every `RTX_E3FSI_SAVE_EVERY` coupled steps (2000) and at the end.
    let saved = std::env::var("RTX_E3FSI_LOAD")
        .ok()
        .map(|d| state::Saved::load(&d).expect("load the saved state"));
    let save_dir = std::env::var("RTX_E3FSI_SAVE").ok();
    let save_every = env_f("RTX_E3FSI_SAVE_EVERY", 2000.0) as usize;
    let mut fl = E3Fluid::build(ny, nz, speed, rest.clone(), saved.as_ref());
    if let Some(s) = &saved {
        println!(
            "  loaded the coupled state at t {:.4} from {} ({}×{}×{} → nz {})",
            s.t,
            std::env::var("RTX_E3FSI_LOAD").unwrap(),
            s.dims[0],
            s.dims[1],
            s.dims[2],
            fl.grid.nz
        );
    }
    let dt = fl.dt;
    println!(
        "  R8-a FSI2 on embedded3: rigid to {t_rigid} s, coupled to {t_end} s; Aitken rtol {rtol:.1e} floor {floor:.1e} max {max_subit} stall accept {stall_accept}; Newmark γ {gamma}; the body's 2D counterpart = the overset SEMICIRCLE line (ny 62: 94.8 mm, 1.914 Hz)"
    );
    let mut csv = csv_path.as_ref().map(|p| {
        let mut f = std::fs::File::create(p).expect("csv");
        writeln!(
            f,
            "t,phase,ux,uy,drag,lift,drag_flag,lift_flag,subit,dres,residual,cg,fresh,sink_dx,sink_dy,lift_p,lift_shear,lift_xdiff,lift_xconv"
        )
        .unwrap();
        f
    });
    // R8-e `RTX_E3FSI_LINE_CSV`: the committed centreline per coupled step
    // (t, then x, y displacement per station) — the unstable mode's shape.
    let mut line_csv = std::env::var("RTX_E3FSI_LINE_CSV")
        .ok()
        .map(|p| std::fs::File::create(p).expect("line csv"));
    let start = std::time::Instant::now();

    // Phase 1: the rigid flag (target 1: the rest state vs CFD2 136.7 / 10.53).
    let rigid_steps = if saved.is_some() {
        0
    } else {
        (t_rigid / dt).round() as usize
    };
    let mut last = ([0.0; 3], Vec::new());
    for step in 0..rigid_steps {
        let r = fl.step();
        assert!(r.final_residual.is_finite(), "rigid death at step {step}");
        if (step + 1) % 50 == 0 || step + 1 == rigid_steps {
            last = fl.loads();
            let (tot, contrib) = &last;
            let (_, on_flag, on_cyl, _) = distribute(&flag_geo, &rest, contrib, fl.load_width);
            let t = fl.time();
            if let Some(f) = csv.as_mut() {
                writeln!(
                    f,
                    "{t:.6},rigid,0,0,{:.5},{:.5},{:.5},{:.5},0,0,{:.3e},{},{},{:.3e},{:.3e},,,,",
                    tot[0],
                    tot[1],
                    on_flag[0],
                    on_flag[1],
                    r.final_residual,
                    r.poisson_iterations,
                    r.fresh_cells,
                    on_flag[0] + on_cyl[0] - tot[0],
                    on_flag[1] + on_cyl[1] - tot[1]
                )
                .unwrap();
            }
            if (step + 1) % 500 == 0 || step + 1 == rigid_steps {
                println!(
                    "  rigid t {t:.3}: drag/span {:.2} lift/span {:+.2} (flag {:.2} {:+.2}, cylinder {:.2} {:+.2}); residual {:.1e}, CG {}; [{:.0} s]",
                    tot[0],
                    tot[1],
                    on_flag[0],
                    on_flag[1],
                    on_cyl[0],
                    on_cyl[1],
                    r.final_residual,
                    r.poisson_iterations,
                    start.elapsed().as_secs_f64()
                );
            }
        }
    }
    let (tot0, contrib0) = last;
    println!(
        "  RIGID ny {ny} nz {nz} at t {:.3}: drag/span {:.2} (CFD2 136.7), lift/span {:+.2} (CFD2 10.53); {rigid_steps} steps in {:.0} s",
        fl.time(),
        tot0[0],
        tot0[1],
        start.elapsed().as_secs_f64()
    );
    if rigid_only && saved.is_none() {
        return;
    }

    // The structure: FSI2's flag at the coupled step dt.
    let mut db = MaterialDatabase::new();
    db.add_material(
        MaterialId(0),
        LinearElastic::new(case.e_s, case.nu_s).with_density(case.rho_s),
        None,
    );
    let beta = (gamma + 0.5).powi(2) / 4.0;
    let analysis = NonlinearDynamicAnalysis::new(
        mesh.clone(),
        db,
        clamp_left(&mesh),
        dt,
        1,
        AnalysisConfig::default(),
    )
    .with_total_lagrangian()
    .with_convergence_criteria(ConvergenceCriteria {
        max_iterations: 60,
        ..ConvergenceCriteria::default()
    })
    .with_newmark_parameters(gamma, beta);
    let flag = RefCell::new(analysis.stepper().unwrap());
    let centre_dofs: Vec<[usize; 2]> = flag_geo
        .centre
        .iter()
        .map(|&id| {
            let d = flag.borrow().node_dofs(id);
            [d[0], d[1]]
        })
        .collect();
    let a_dofs = flag.borrow().node_dofs(flag_geo.a_node);
    // R8-e `RTX_E3FSI_KIN=surface`: the fluid's centreline is the mean of the
    // bottom and top edge nodes at each station (the nodes the loads land on)
    // instead of the mid-plane nodes (default `centre`).
    let kin_surface = std::env::var("RTX_E3FSI_KIN").is_ok_and(|v| v == "surface");
    let surface_dofs: Vec<[usize; 4]> = flag_geo
        .surface
        .iter()
        .map(|[b, t]| {
            let (db, dt_) = (flag.borrow().node_dofs(*b), flag.borrow().node_dofs(*t));
            [db[0], db[1], dt_[0], dt_[1]]
        })
        .collect();
    let extract = |s: &DynamicState| -> Vec<f64> {
        let mut c = vec![0.0; 2 * stations()];
        if kin_surface {
            for (k, d) in surface_dofs.iter().enumerate() {
                c[2 * k] = 0.5 * (s.displacement[d[0]] + s.displacement[d[2]]);
                c[2 * k + 1] = 0.5 * (s.displacement[d[1]] + s.displacement[d[3]]);
            }
            return c;
        }
        for (k, d) in centre_dofs.iter().enumerate() {
            c[2 * k] = s.displacement[d[0]];
            c[2 * k + 1] = s.displacement[d[1]];
        }
        c
    };

    let (mut flag_state, mut committed_nodal, mut line_n, mut c_fluid_n) = match &saved {
        Some(s) => {
            let nodal: Vec<(NodeId, Vector3<f64>)> = flag_geo
                .interface
                .wetted
                .iter()
                .enumerate()
                .map(|(k, &id)| (id, Vector3::new(s.nodal[2 * k], s.nodal[2 * k + 1], 0.0)))
                .collect();
            let state = DynamicState {
                displacement: DVector::from_vec(s.disp.clone()),
                velocity: DVector::from_vec(s.vel.clone()),
                acceleration: DVector::from_vec(s.acc.clone()),
            };
            (state, nodal, s.line.clone(), s.c_fluid.clone())
        }
        None => {
            let (nodal0, _, _, _) = distribute(&flag_geo, &rest, &contrib0, fl.load_width);
            flag.borrow_mut().set_nodal_forces(&nodal0);
            let state = flag.borrow_mut().rest_state().unwrap();
            // The fluid's own previous line and centreline (its geometry's history).
            let line = Line {
                t: fl.time(),
                ..rest.clone()
            };
            (state, nodal0, line, zero_c.clone())
        }
    };

    let coupled_steps = ((t_end - fl.time()) / dt).round() as usize;
    let fl = RefCell::new(fl);
    let (mut times, mut uy_series, mut ux_series) = (Vec::new(), Vec::new(), Vec::new());
    let (mut drag_s, mut lift_s) = (Vec::new(), Vec::new());
    let (mut total_subit, mut max_seen, mut stalled) = (0usize, 0usize, 0usize);
    let mut death: Option<String> = None;
    let t_fluid = std::cell::Cell::new(0.0f64);
    let t_restore = std::cell::Cell::new(0.0f64);
    let t_loads = std::cell::Cell::new(0.0f64);
    let t_struct = std::cell::Cell::new(0.0f64);
    let phase_start = std::time::Instant::now();
    for step in 0..coupled_steps {
        let t_old = fl.borrow().time();
        let t_new = t_old + dt;
        let predicted = {
            flag.borrow_mut().set_nodal_forces(&committed_nodal);
            let (p, _) = flag.borrow_mut().step(&flag_state).unwrap();
            extract(&p)
        };
        let c_struct_n = extract(&flag_state);
        let snap = fl.borrow_mut().snapshot();
        let dirty = std::cell::Cell::new(false);
        type Pass = (
            DynamicState,
            Vec<(NodeId, Vector3<f64>)>,
            Line,
            Vec<f64>,
            [f64; 3],
            [f64; 2],
            fluid_step::Stats,
        );
        let latest: RefCell<Option<Pass>> = RefCell::new(None);
        let pass = |cand: &[f64]| -> Vec<f64> {
            let line = line_of(t_new, cand, &c_fluid_n, dt);
            let mut f = fl.borrow_mut();
            let tr = std::time::Instant::now();
            f.set_lines(line_n.clone(), line.clone());
            if dirty.get() {
                f.restore(&snap);
            }
            dirty.set(true);
            t_restore.set(t_restore.get() + tr.elapsed().as_secs_f64());
            let tf = std::time::Instant::now();
            let r = f.step();
            t_fluid.set(t_fluid.get() + tf.elapsed().as_secs_f64());
            if !r.final_residual.is_finite() {
                return vec![f64::NAN; cand.len()];
            }
            let tl = std::time::Instant::now();
            let (tot, contrib) = f.loads();
            let (nodal, on_flag, on_cyl, parts_y) =
                distribute(&flag_geo, &line, &contrib, f.load_width);
            t_loads.set(t_loads.get() + tl.elapsed().as_secs_f64());
            let ts = std::time::Instant::now();
            let mut st = flag.borrow_mut();
            st.set_nodal_forces(&nodal);
            let (new_state, _) = st.step(&flag_state).unwrap();
            t_struct.set(t_struct.get() + ts.elapsed().as_secs_f64());
            let out = extract(&new_state);
            if step < trace {
                let res: f64 = out
                    .iter()
                    .zip(cand)
                    .map(|(a, b)| (a - b).powi(2))
                    .sum::<f64>()
                    .sqrt();
                println!(
                    "    step {step} pass: |c_new − c_cand| {res:.3e}, tip cand ({:+.4e}, {:+.4e}), load flag ({:+.3}, {:+.3}) cyl ({:+.3}, {:+.3}) total ({:+.3}, {:+.3}), residual {:.1e}, fresh {}",
                    cand[2 * stations() - 2],
                    cand[2 * stations() - 1],
                    on_flag[0],
                    on_flag[1],
                    on_cyl[0],
                    on_cyl[1],
                    tot[0],
                    tot[1],
                    r.final_residual,
                    r.fresh_cells
                );
            }
            let stats = fluid_step::Stats {
                residual: r.final_residual,
                cg: r.poisson_iterations,
                fresh: r.fresh_cells,
                sink: [
                    on_flag[0] + on_cyl[0] - tot[0],
                    on_flag[1] + on_cyl[1] - tot[1],
                ],
                parts_y,
            };
            *latest.borrow_mut() =
                Some((new_state, nodal, line, cand.to_vec(), tot, on_flag, stats));
            out
        };
        let increment: f64 = predicted
            .iter()
            .zip(&c_struct_n)
            .map(|(a, b)| (a - b).powi(2))
            .sum::<f64>()
            .sqrt();
        let tol = floor.max(rtol * increment);
        let acceptable = (stall_accept * tol).max(0.1 * increment);
        let outcome = Subiterated::aitken(max_subit, tol)
            .unwrap()
            .solve(&predicted, &pass);
        // `RTX_E3FSI_RESTORE_CHECK=N`: on the first N steps, the accepted
        // candidate re-run from the snapshot (the restore's repeatability).
        if step < restore_check && outcome.is_ok() {
            let (cand, tot_a, out_a) = {
                let l = latest.borrow();
                let l = l.as_ref().expect("a pass ran");
                (l.3.clone(), l.4, extract(&l.0))
            };
            let out_b = pass(&cand);
            let tot_b = latest.borrow().as_ref().expect("re-run").4;
            let dout: f64 = out_a
                .iter()
                .zip(&out_b)
                .map(|(a, b)| (a - b).powi(2))
                .sum::<f64>()
                .sqrt();
            println!(
                "  RESTORE CHECK step {step}: load a ({:+.9e}, {:+.9e}) b ({:+.9e}, {:+.9e}), |Δ structure output| {dout:.3e}",
                tot_a[0], tot_a[1], tot_b[0], tot_b[1]
            );
        }
        let (iters, dres) = match outcome {
            Ok(c) => (c.iterations, c.residual),
            Err(
                rtx_fsi::FsiError::CouplingNotConverged {
                    iterations,
                    residual,
                    ..
                }
                | rtx_fsi::FsiError::CouplingDiverged {
                    iterations,
                    residual,
                },
            ) if residual < acceptable => {
                stalled += 1;
                (iterations, residual)
            }
            Err(e) => {
                println!(
                    "  R8-a DEATH at coupled step {step} t = {t_new:.4}: {e:?} (increment {increment:.3e}, tol {tol:.3e}, acceptable {acceptable:.3e})"
                );
                death = Some(format!("step {step} t {t_new:.4}: {e:?}"));
                break;
            }
        };
        total_subit += iters;
        max_seen = max_seen.max(iters);
        let (new_state, nodal, line, cand, tot, on_flag, stats) =
            latest.borrow_mut().take().expect("a pass ran");
        flag_state = new_state;
        committed_nodal = nodal;
        line_n = line;
        c_fluid_n = cand;
        let ux = flag_state.displacement[a_dofs[0]];
        let uy = flag_state.displacement[a_dofs[1]];
        times.push(t_new);
        ux_series.push(ux);
        uy_series.push(uy);
        drag_s.push(tot[0]);
        lift_s.push(tot[1]);
        if let Some(f) = csv.as_mut() {
            writeln!(
                f,
                "{t_new:.6},coupled,{ux:.6e},{uy:.6e},{:.5},{:.5},{:.5},{:.5},{iters},{dres:.3e},{:.3e},{},{},{:.3e},{:.3e},{:.4},{:.4},{:.4},{:.4}",
                tot[0],
                tot[1],
                on_flag[0],
                on_flag[1],
                stats.residual,
                stats.cg,
                stats.fresh,
                stats.sink[0],
                stats.sink[1],
                stats.parts_y[0],
                stats.parts_y[1],
                stats.parts_y[2],
                stats.parts_y[3]
            )
            .unwrap();
        }
        if let Some(f) = line_csv.as_mut() {
            // Then the y displacement of the mid-plane, bottom and top nodes per station.
            let mut c: Vec<String> = c_fluid_n.iter().map(|v| format!("{v:.6e}")).collect();
            for (k, d) in surface_dofs.iter().enumerate() {
                let yc = flag_state.displacement[centre_dofs[k][1]];
                let (yb, yt) = (flag_state.displacement[d[1]], flag_state.displacement[d[3]]);
                c.extend([yc, yb, yt].map(|v| format!("{v:.6e}")));
            }
            writeln!(f, "{t_new:.6},{}", c.join(",")).unwrap();
        }
        let at_end = step + 1 == coupled_steps;
        if let Some(dir) = save_dir
            .as_ref()
            .filter(|_| at_end || (step + 1) % save_every == 0)
        {
            let f = fl.borrow();
            let g = f.grid;
            state::Saved {
                t: t_new,
                dims: [g.nx, g.ny, g.nz],
                h: f.h,
                u: f.field.u.clone(),
                v: f.field.v.clone(),
                w: f.field.w.clone(),
                p: f.field.p.clone(),
                disp: flag_state.displacement.as_slice().to_vec(),
                vel: flag_state.velocity.as_slice().to_vec(),
                acc: flag_state.acceleration.as_slice().to_vec(),
                line: line_n.clone(),
                c_fluid: c_fluid_n.clone(),
                nodal: committed_nodal
                    .iter()
                    .flat_map(|(_, v)| [v.x, v.y])
                    .collect(),
            }
            .save(dir)
            .expect("save the coupled state");
        }
        if (step + 1) % 250 == 0 {
            let w = &uy_series[uy_series.len().saturating_sub(600)..];
            let (mid, amp) = mid_amp(w);
            println!(
                "  t {t_new:.3} ({} steps): uy(A) {uy:+.4e} ux {ux:+.4e} (last ~1 period mid {mid:+.3e} amp {amp:.3e}), drag {:.1} lift {:+.1}; {:.2} subit/step (max {max_seen}, stalled {stalled}); fluid {:.0} s restore {:.0} s loads {:.0} s structure {:.0} s of {:.0} s",
                step + 1,
                tot[0],
                tot[1],
                total_subit as f64 / (step + 1) as f64,
                t_fluid.get(),
                t_restore.get(),
                t_loads.get(),
                t_struct.get(),
                phase_start.elapsed().as_secs_f64()
            );
        }
    }
    // Summary: the last two seconds (or what there is).
    let n = times.len();
    if n > 0 {
        let t_last = times[n - 1];
        let from = times.partition_point(|&t| t < t_last - 2.0);
        let (mid, amp) = mid_amp(&uy_series[from..]);
        let (uxm, uxa) = mid_amp(&ux_series[from..]);
        let f = fsi2_harness::crossing_frequency(&times[from..], &uy_series[from..]);
        let mut d: Vec<f64> = drag_s[from..].to_vec();
        let (dm, _) = mid_amp(&d);
        let (lm, la) = mid_amp(&lift_s[from..]);
        let dmed = median(&mut d);
        println!(
            "  FINAL R8-a ny {ny} nz {nz}: {n} coupled steps to t {t_last:.3}; last 2 s: uy(A) {:.2} ± {:.2} mm, ux(A) {:.2} ± {:.2} mm, f {}, drag mid {dm:.1} (median {dmed:.1}), lift {lm:+.1} ± {la:.1}; {:.2} subit/step (max {max_seen}, stalled {stalled}); death {}; wall {:.0} s",
            1e3 * mid,
            1e3 * amp,
            1e3 * uxm,
            1e3 * uxa,
            f.map_or("n/a".into(), |f| format!("{f:.4} Hz")),
            total_subit as f64 / n as f64,
            death.as_deref().unwrap_or("none"),
            start.elapsed().as_secs_f64()
        );
    }
}

mod fluid_step {
    pub struct Stats {
        pub residual: f64,
        pub cg: usize,
        pub fresh: usize,
        pub sink: [f64; 2],
        pub parts_y: [f64; 4],
    }
}