rtx-fsi: C2 — FSI2's coupled march, and the wrong attractor measured

The unsteady coupling on the piston pattern with the real solvers: per
TIME STEP, Aitken subiterations of (8x-subcycled embedded TVD/multigrid
fluid <-> flag nonlinear-Newmark step), the fluid re-runnable through
snapshot/restore + a field clone, the flag re-runnable because the
stepper commits nothing, the moving polygon carrying the flag's actual
interface velocity. Conservation 8.3e-12 over 9,263 coupled steps; the
rigid phase reproduces this solver's own CFD2 values at both grids
before anything couples.

What the 30 s studies measured (ny = 62 AND 82, three protocols): the
coupled system self-excites into a wake-forced cycle at 3.729/3.728 Hz
with uy(A) +-17.3 mm at both grids - grid-converged, protocol-
independent, kinematically genuine (ux locks 2x uy) - and NOT the
benchmark's mode-2 resonant cycle (1.93 Hz / 81.6 mm; the flag's vacuum
mode 2 is 1.9245 Hz). Diagnosis measured stepwise: the staggered phase
lag starves the resonant channel (a subcycle=2 probe redirected growth
to 1.9 Hz as predicted), and tighter dt_c is blocked because wall-
velocity noise = displacement-tolerance / dt_c while the needed
tolerance sits below the discrete interface noise floor (~1.3e-4 at
full inflow, mask flips through beta dt^2/m). The route to the
benchmark cycle is lowering that floor, not iterating against it.

Robustness, both measured: rare wild tractions (19 samples in 2.4M)
clamped at 20x the sample median - clamped, not dropped: a hard drop
makes the pass discontinuous and the subiteration bounces at step
scale; and Aitken "divergence" verdicts within 5x tolerance accepted as
noise bounces over well-predicted steps (counted, bounded).

The committed default (t_end = 7, ~9.5 min) pins the deterministic
release response (uy 3.773 +- 3.792 mm, band +-35%); study horizons
>= 25 s pin the measured 3.73 Hz attractor so any material change is
loud. Full study record in the module docs and omni-cortex
solver_status.md.

Co-Authored-By: Claude Fable 5 <[email protected]>
Claude-Session: https://claude.ai/code/session_01Lnyrw33Lu6rUhW42E9KHwq

crates/specialized/rtx-fsi/tests/turek_hron_fsi2.rs

@@ -0,0 +1,918 @@
+//! Turek–Hron FSI2: the self-excited flapping flag — rung C2 of the ladder
+//! (omni-cortex `docs/turek_hron_geometry_decision.md`).
+//!
+//! Re = 100 channel flow (`U = 1`) past the rigid cylinder with the elastic
+//! flag at density ratio `rho_s / rho_f = 10`: the flow destabilises the
+//! flag into a large-amplitude limit cycle. Reference (FEATFLOW level 4,
+//! dt = 0.001): `ux(A) = −14.85 ± 12.70 mm [3.86 Hz]`,
+//! `uy(A) = 1.30 ± 81.6 mm [1.93 Hz]`, drag `215.06 ± 77.65`,
+//! lift `0.61 ± 237.8`.
+//!
+//! # The march
+//!
+//! Unlike FSI1's steady fixed point, FSI2 marches in time: per TIME STEP
+//! the fluid step and the flag's nonlinear-Newmark step are subiterated
+//! with Aitken until the end-of-step interface displacement converges —
+//! the coupled piston benchmark's structure, with the real 2-D solvers.
+//! The fluid step is re-runnable inside a subiteration through
+//! [`EmbeddedPisoSolver::snapshot`]/`restore` plus a [`FlowField`] clone;
+//! the flag step is re-runnable because [`NonlinearDynamicStepper::step`]
+//! commits nothing. The moving polygon carries the flag's actual interface
+//! velocity — finite-differenced end-of-step positions over `dt`,
+//! interpolated along the nearest edge (`polygon_interface_velocity`) —
+//! replacing the zero-velocity closure FSI1's steady case allowed. The
+//! fluid keeps TVD convection (shedding physics; limiter chatter is
+//! harmless in time marching) and the multigrid projection.
+//!
+//! # Phases (the validation ladder inside FSI2)
+//!
+//! 1. **Rigid flag** to `t_release`: the ramped inflow over the fixed
+//!    geometry must land near the CFD2 steady state this solver already
+//!    measured (surface drag 121.4 / 123.3 at ny = 62 / 82 vs its CFD2
+//!    runs' ~121 / 122.6) — the harness's fluid configuration is checked
+//!    against a known number before anything couples.
+//! 2. **Release**: the flag starts at rest under the sampled fluid load
+//!    (consistent initial acceleration), and the coupled march runs to
+//!    `t_end`.
+//!
+//! # What the 2026-08-21 study measured (t = 30 s marches, release)
+//!
+//! **The coupled system self-excites at every configuration tried, and at
+//! the loosely-coupled default (8 fluid substeps per coupled step, ~1
+//! subiteration) it lands in a wake-forced cycle at 3.729 / 3.728 Hz with
+//! uy(A) ± 17.3 mm at BOTH ny = 62 and ny = 82 — grid-converged, and
+//! protocol-independent (release at t = 6 and coupled-from-t = 0 reach
+//! the same state). This is NOT the benchmark's cycle** (1.93 Hz,
+//! ± 81.6 mm). The identification is clean: the flag's vacuum mode 2 is
+//! 1.9245 Hz (modal analysis, 35x2 Quad8) — the reference cycle IS mode-2
+//! resonance — while 3.73 Hz matches no structural mode (mode 3 is
+//! 5.26 Hz); the measured state is the heavy flag's off-resonance forced
+//! response at the wake's own shedding frequency, and its ux mean
+//! (−0.8 mm) matches the foreshortening scaling (amp/81.6)^2 x (−14.85).
+//! Loads at ny = 82: drag 141.6 ± 53.9 (ref 215.06 ± 77.65), lift
+//! 49 ± 508 (ref 0.61 ± 237.8) — consistent with the small-amplitude
+//! state.
+//!
+//! Why mode 2 does not win here, measured stepwise: (a) at the default
+//! coupling the motion-load staggered phase lag (~omega dt_c) starves the
+//! resonant channel — a subcycle = 2 probe (lag / 4) redirected early
+//! growth into 1.9 / 2.85 Hz exactly as that predicts; (b) but the probe
+//! then destabilised: at a fixed interface-DISPLACEMENT tolerance the
+//! no-slip closure's wall-velocity noise is tol / dt_c (~0.3 m/s at
+//! 2e-4 / 6.5e-4 — 30% of the mean inflow), and the fluid pumped up and
+//! blew through the coupling. The displacement tolerance a smaller dt_c
+//! needs (~1e-5) sits BELOW the discrete interface noise floor (~1e-4 at
+//! full inflow, from mask flips through Newmark's beta dt^2/m). **The
+//! route to the benchmark cycle is lowering the interface noise floor**
+//! (smoother load sampling / mask transitions, or a vector quasi-Newton
+//! interface solver in place of scalar Aitken), not more iterations
+//! against it.
+//!
+//! Robustness findings, both measured: rare wild tractions from
+//! near-degenerate reconstructions (19 samples in 2.4 million) killed a
+//! t = 25.8 s march through the flag's Newton until the spike CLAMP (20x
+//! the sample median, direction kept — clamping, not dropping: a hard
+//! drop makes the coupling pass discontinuous and the subiteration
+//! bounces at the step scale) and a 60-iteration Newton budget; with
+//! both, the same march runs to t = 30 clean.
+//!
+//! Machinery invariants asserted every run: load-transfer conservation
+//! (partition of unity) at 1e-10 (measured 8e-12 over 9,263 steps),
+//! coupled convergence bookkeeping, finite fields. The committed default
+//! (t_end = 7) pins the deterministic release response; study horizons
+//! (t_end >= 25) pin the measured attractor so any material change is
+//! loud. Full trajectories: the session scratchpad study logs.
+//!
+//! Environment knobs: `RTX_FSI2_NY` (fluid resolution, default 62),
+//! `RTX_FSI2_T_RELEASE` (default 6 s), `RTX_FSI2_T_END` (default 7 s —
+//! the committed onset segment; studies run 30), `RTX_FSI2_SUBCYCLE`
+//! (fluid substeps per coupled step, default 8), `RTX_FSI2_TOL` /
+//! `RTX_FSI2_RTOL` (interface tolerance floor and its
+//! relative-to-increment part), `RTX_FSI2_MAXSUB` (Aitken budget,
+//! default 12), `RTX_FSI2_FLAG_NX` (flag mesh, default 35),
+//! `RTX_FSI2_CSV` (trajectory dump path).
+
+use std::cell::RefCell;
+use std::io::Write as _;
+use std::sync::{Arc, RwLock};
+
+use nalgebra::Vector3;
+use rtx_cfd::CfdConfig;
+use rtx_cfd::solvers::incompressible::{
+    AleBoundaries, ConvectionScheme, EmbeddedBody, EmbeddedParameters, EmbeddedPisoSolver,
+    FlowField, PoissonSolverKind, SideBoundary, polygon_interface_velocity,
+    polygon_signed_distance,
+};
+use rtx_fea::analysis::{
+    AnalysisConfig, ConvergenceCriteria, DynamicState, NonlinearDynamicAnalysis,
+};
+use rtx_fea::assembly::dof_mapping::DofComponent;
+use rtx_fea::boundary::dirichlet::{DirichletBC, DirichletType};
+use rtx_fea::boundary::{BoundaryCondition, BoundaryConditionSet, SpatialFunction};
+use rtx_fea::materials::{LinearElastic, MaterialDatabase};
+use rtx_fea::mesh::{Element, ElementType, MaterialId, Mesh, Node, NodeId};
+use rtx_fsi::{FluidFace, Subiterated, WettedSurface};
+
+const L: f64 = 2.5;
+const H: f64 = 0.41;
+const RHO_F: f64 = 1000.0;
+const NU_F: f64 = 1e-3;
+const U_MEAN: f64 = 1.0;
+const RHO_S: f64 = 10_000.0;
+const E_S: f64 = 1.4e6;
+const NU_S: f64 = 0.4;
+
+const FLAG_X0: f64 = 0.25;
+const FLAG_X1: f64 = 0.6;
+const FLAG_Y0: f64 = 0.19;
+const FLAG_Y1: f64 = 0.21;
+
+// FEATFLOW level-4, dt 0.001 reference values.
+const REF_UY_MEAN: f64 = 1.30e-3;
+const REF_UY_AMP: f64 = 81.6e-3;
+const REF_UY_FREQ: f64 = 1.93;
+const REF_UX_MEAN: f64 = -14.85e-3;
+const REF_UX_AMP: f64 = 12.70e-3;
+const REF_DRAG_MEAN: f64 = 215.06;
+const REF_LIFT_AMP: f64 = 237.8;
+
+fn circle_sdf(x: f64, y: f64) -> f64 {
+    ((x - 0.2).powi(2) + (y - 0.2).powi(2)).sqrt() - 0.05
+}
+
+/// The ramped parabolic inflow of the benchmark definition.
+fn inflow(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)
+}
+
+fn env_or(name: &str, default: f64) -> f64 {
+    std::env::var(name)
+        .map(|v| v.parse().expect(name))
+        .unwrap_or(default)
+}
+
+/// The flag's Quad8 mesh (as in FSI1 and the CSM tests).
+fn flag_mesh(nx: usize, ny: usize) -> Mesh {
+    let mut mesh = Mesh::new(2).unwrap();
+    let (lx, ly) = (2 * nx + 1, 2 * ny + 1);
+    let mut grid = vec![vec![None; ly]; lx];
+    for (i, column) in grid.iter_mut().enumerate() {
+        for (j, slot) in column.iter_mut().enumerate() {
+            if i % 2 == 1 && j % 2 == 1 {
+                continue;
+            }
+            let x = FLAG_X0 + (FLAG_X1 - FLAG_X0) * i as f64 / (2 * nx) as f64;
+            let y = FLAG_Y0 + (FLAG_Y1 - FLAG_Y0) * j as f64 / (2 * ny) as f64;
+            *slot = Some(mesh.add_node(Node::new_2d(x, y)));
+        }
+    }
+    for i in 0..nx {
+        for j in 0..ny {
+            let (a, b) = (2 * i, 2 * j);
+            let nodes = vec![
+                grid[a][b].unwrap(),
+                grid[a + 2][b].unwrap(),
+                grid[a + 2][b + 2].unwrap(),
+                grid[a][b + 2].unwrap(),
+                grid[a + 1][b].unwrap(),
+                grid[a + 2][b + 1].unwrap(),
+                grid[a + 1][b + 2].unwrap(),
+                grid[a][b + 1].unwrap(),
+            ];
+            mesh.add_element(Element::new(ElementType::Quad8, nodes, MaterialId(0)).unwrap())
+                .unwrap();
+        }
+    }
+    mesh
+}
+
+/// The wetted-interface bookkeeping (FSI1's, plus vertex velocities).
+struct Interface {
+    wetted: Vec<NodeId>,
+    reference: Vec<(f64, f64)>,
+    /// Ordered boundary walk: indices into `wetted` (`usize::MAX` marks
+    /// the fixed anchor vertices inside the cylinder / at the clamp).
+    walk: Vec<(usize, (f64, f64))>,
+}
+
+impl Interface {
+    fn build(mesh: &Mesh) -> Self {
+        let eps = 1e-9;
+        let on_bottom = |p: Vector3<f64>| (p.y - FLAG_Y0).abs() < eps;
+        let on_top = |p: Vector3<f64>| (p.y - FLAG_Y1).abs() < eps;
+        let on_tip = |p: Vector3<f64>| (p.x - FLAG_X1).abs() < eps;
+        let clamped = |p: Vector3<f64>| (p.x - FLAG_X0).abs() < eps;
+
+        let mut wetted: Vec<(NodeId, (f64, f64))> = mesh
+            .nodes
+            .iter()
+            .filter(|(_, node)| {
+                let p = node.position();
+                (on_bottom(p) || on_top(p) || on_tip(p)) && !clamped(p)
+            })
+            .map(|(&id, node)| (id, (node.position().x, node.position().y)))
+            .collect();
+        wetted.sort_by_key(|(id, _)| *id);
+        let index_of = |id: NodeId| wetted.iter().position(|(w, _)| *w == id).unwrap();
+
+        let mut bottom: Vec<(NodeId, f64)> = mesh
+            .nodes
+            .iter()
+            .filter(|(_, n)| on_bottom(n.position()) && !clamped(n.position()))
+            .map(|(&id, n)| (id, n.position().x))
+            .collect();
+        bottom.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap());
+        let mut tip: Vec<(NodeId, f64)> = mesh
+            .nodes
+            .iter()
+            .filter(|(_, n)| {
+                let p = n.position();
+                on_tip(p) && !on_bottom(p) && !on_top(p)
+            })
+            .map(|(&id, n)| (id, n.position().y))
+            .collect();
+        tip.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap());
+        let mut top: Vec<(NodeId, f64)> = mesh
+            .nodes
+            .iter()
+            .filter(|(_, n)| on_top(n.position()) && !clamped(n.position()))
+            .map(|(&id, n)| (id, n.position().x))
+            .collect();
+        top.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap());
+
+        let mut walk: Vec<(usize, (f64, f64))> = Vec::new();
+        walk.push((usize::MAX, (0.22, FLAG_Y0)));
+        walk.push((usize::MAX, (FLAG_X0, FLAG_Y0)));
+        for (id, _) in &bottom {
+            walk.push((index_of(*id), (0.0, 0.0)));
+        }
+        for (id, _) in &tip {
+            walk.push((index_of(*id), (0.0, 0.0)));
+        }
+        for (id, _) in &top {
+            walk.push((index_of(*id), (0.0, 0.0)));
+        }
+        walk.push((usize::MAX, (FLAG_X0, FLAG_Y1)));
+        walk.push((usize::MAX, (0.22, FLAG_Y1)));
+
+        let reference = wetted.iter().map(|(_, p)| *p).collect();
+        Self {
+            wetted: wetted.into_iter().map(|(id, _)| id).collect(),
+            reference,
+            walk,
+        }
+    }
+
+    /// Deformed polygon vertices for the interface vector `d`.
+    fn polygon(&self, d: &[f64]) -> Vec<(f64, f64)> {
+        self.walk
+            .iter()
+            .map(|&(k, anchor)| {
+                if k == usize::MAX {
+                    anchor
+                } else {
+                    let (x0, y0) = self.reference[k];
+                    (x0 + d[2 * k], y0 + d[2 * k + 1])
+                }
+            })
+            .collect()
+    }
+
+    /// Per-vertex velocities for the interface velocity vector `ddot`
+    /// (anchors do not move).
+    fn walk_velocities(&self, ddot: &[f64]) -> Vec<(f64, f64)> {
+        self.walk
+            .iter()
+            .map(|&(k, _)| {
+                if k == usize::MAX {
+                    (0.0, 0.0)
+                } else {
+                    (ddot[2 * k], ddot[2 * k + 1])
+                }
+            })
+            .collect()
+    }
+
+    /// Deformed wetted node positions for the transfer.
+    fn deformed_nodes(&self, d: &[f64]) -> Vec<Vector3<f64>> {
+        self.reference
+            .iter()
+            .enumerate()
+            .map(|(k, &(x0, y0))| Vector3::new(x0 + d[2 * k], y0 + d[2 * k + 1], 0.0))
+            .collect()
+    }
+}
+
+fn clamp_left(mesh: &Mesh) -> BoundaryConditionSet {
+    let clamped: Vec<NodeId> = mesh
+        .nodes
+        .iter()
+        .filter(|(_, node)| (node.position().x - FLAG_X0).abs() < 1e-9)
+        .map(|(&id, _)| id)
+        .collect();
+    let mut set = BoundaryConditionSet::new();
+    for component in [DofComponent::DisplacementX, DofComponent::DisplacementY] {
+        set.add_condition(BoundaryCondition::Dirichlet(DirichletBC {
+            nodes: clamped.clone(),
+            components: vec![component],
+            condition_type: DirichletType::Spatial(SpatialFunction(Box::new(|_| 0.0))),
+            time_range: None,
+            ramping_factor: 1.0,
+            gradual_enforcement: false,
+        }));
+    }
+    set
+}
+
+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 crossings of the mean.
+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 {
+            crossings.push(times[k - 1] + (a / (a - b)) * (times[k] - times[k - 1]));
+        }
+    }
+    if crossings.len() < 3 {
+        return None;
+    }
+    Some((crossings.len() - 1) as f64 / (crossings.last().unwrap() - crossings.first().unwrap()))
+}
+
+#[test]
+#[allow(clippy::too_many_lines)]
+fn fsi2_flapping_flag() {
+    let ny = env_or("RTX_FSI2_NY", 62.0) as usize;
+    let t_release = env_or("RTX_FSI2_T_RELEASE", 6.0);
+    let t_end = env_or("RTX_FSI2_T_END", 7.0);
+    // Per-step interface tolerance: max(absolute floor, RTOL x that step's
+    // interface increment). The floor is measured, not wished, and it
+    // RIDES WITH THE LOADS: the interface map carries a noise floor from
+    // discrete mask flips under vanishing geometry changes (each flip's
+    // load jump maps through Newmark's beta dt^2 / m into displacement) —
+    // measured ~5.8e-7 per pass at 2% inflow and ~1.3e-4 at full inflow
+    // on ny = 62. Sub-cell interface accuracy is not the fluid's to
+    // promise. A step that stalls at the floor — including an Aitken
+    // "divergence" verdict within 5x of the tolerance, which is noise
+    // bouncing over a well-predicted (tiny) first residual, not added
+    // mass (the ratio-10 flag's per-node added-mass gain is far below
+    // one) — is ACCEPTED at its last candidate and counted
+    // (`stalled_steps`), never silently retried to a budget: the count
+    // and the worst stalled residual are reported and bounded at the
+    // end. Genuine runaway (residual far beyond the noise scale) still
+    // panics. At the limit cycle the floor is ~1-2% of the per-step
+    // interface increment, so the committed trajectory carries
+    // noise-level interface error each step — recorded, and bounded by
+    // the two-grid amplitude rule before belief.
+    let tol_floor = env_or("RTX_FSI2_TOL", 2e-4);
+    let rtol = env_or("RTX_FSI2_RTOL", 1e-2);
+    let max_subiterations_budget = env_or("RTX_FSI2_MAXSUB", 12.0) as usize;
+    let flag_nx = env_or("RTX_FSI2_FLAG_NX", 35.0) as usize;
+    let csv_path = std::env::var("RTX_FSI2_CSV").ok();
+
+    let h = H / ny as f64;
+    let nx = (L / h).round() as usize;
+    let mu = RHO_F * NU_F;
+    let u_peak = 1.5 * 1.5 * U_MEAN;
+    // The fluid's explicit limit; the coupling (and the flag's Newmark)
+    // run at `subcycle` fluid steps per coupled step — the structure and
+    // the transfer need nowhere near the fluid's dt (CSM3 measured 0.23%
+    // frequency error at dt = 5e-3; dt_c here is ~2.6e-3 at ny = 62), and
+    // the per-pass cost is dominated by the structure solve. Within a
+    // pass the interface geometry is interpolated linearly across the
+    // substeps, so the mask still moves less than a cell per fluid step.
+    let dt_fluid = 0.25 / (2.0 * u_peak / h + 4.0 * NU_F / (h * h));
+    let subcycle = env_or("RTX_FSI2_SUBCYCLE", 8.0) as usize;
+    let dt = dt_fluid * subcycle as f64;
+
+    let mesh = flag_mesh(flag_nx, 2);
+    let interface = Interface::build(&mesh);
+    let a_node = mesh
+        .nodes
+        .iter()
+        .find(|(_, n)| (n.position().x - 0.6).abs() < 1e-9 && (n.position().y - 0.2).abs() < 1e-9)
+        .map(|(&id, _)| id)
+        .expect("point A");
+
+    // The deformable geometry AND its velocity, behind one lock: the
+    // fluid's per-step mask rebuild reads the polygon; the no-slip closure
+    // reads both.
+    let zero_d = vec![0.0; 2 * interface.wetted.len()];
+    let shared = Arc::new(RwLock::new((
+        interface.polygon(&zero_d),
+        interface.walk_velocities(&zero_d),
+    )));
+    let sdf_shared = shared.clone();
+    let vel_shared = shared.clone();
+
+    let config = CfdConfig::new()
+        .with_density(RHO_F)
+        .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,
+        // TVD, deliberately: FSI2 marches in time and needs the shedding
+        // physics upwind's numerical viscosity killed on these grids
+        // (CFD3's finding). The limiter chatter that defeats steady fixed
+        // points (FSI1's finding) is harmless here — each step's fixed
+        // point is the interface displacement of THAT step, not a steady
+        // load.
+        convection_scheme: ConvectionScheme::TvdVanAlbada,
+    };
+    let mut solver = EmbeddedPisoSolver::new(config, params).unwrap();
+    solver.set_boundary_velocity(|x, y, t| {
+        if x <= 0.0 {
+            (inflow(y, t), 0.0)
+        } else {
+            (0.0, 0.0)
+        }
+    });
+    solver.set_moving_body(
+        EmbeddedBody::from_sdf(move |x, y, _| {
+            let geometry = sdf_shared.read().unwrap();
+            circle_sdf(x, y).min(polygon_signed_distance(&geometry.0, x, y))
+        })
+        .with_surface_velocity(move |x, y, _| {
+            let geometry = vel_shared.read().unwrap();
+            if circle_sdf(x, y) <= polygon_signed_distance(&geometry.0, x, y) {
+                (0.0, 0.0)
+            } else {
+                polygon_interface_velocity(&geometry.0, &geometry.1, x, y)
+            }
+        }),
+    );
+
+    // Start at rest; the ramp brings the inflow up from zero.
+    let mut field = FlowField::new(nx, ny, h, h).unwrap();
+    solver.initialize(&mut field).unwrap();
+
+    // Phase 1: rigid flag to t_release.
+    let start = std::time::Instant::now();
+    let rigid_steps = (t_release / dt_fluid).round() as usize;
+    for _ in 0..rigid_steps {
+        futures::executor::block_on(solver.advance(&mut field, dt_fluid)).unwrap();
+    }
+    // Surface drag and lift on cylinder + flag at the current geometry.
+    let measure_force = |solver: &EmbeddedPisoSolver, field: &FlowField| -> (f64, f64) {
+        let mask = solver.mask().unwrap();
+        let body = solver.body().unwrap();
+        let vertices = shared.read().unwrap().0.clone();
+        let mut drag = 0.0;
+        let mut lift = 0.0;
+        let poly_probe = EmbeddedBody::polygon(vertices.clone());
+        for s in poly_probe.surface_samples(0.5 * h) {
+            if circle_sdf(s.x, s.y) < 1e-9 {
+                continue;
+            }
+            if let Some((tx, ty)) = mask.traction_at(
+                body, &field.u, &field.v, &field.p, mu, 0.0, s.x, s.y, s.nx, s.ny,
+            ) {
+                drag += tx * s.ds;
+                lift += ty * s.ds;
+            }
+        }
+        let circle_probe = EmbeddedBody::circle(0.2, 0.2, 0.05);
+        for s in circle_probe.surface_samples(0.5 * h) {
+            if polygon_signed_distance(&vertices, s.x, s.y) < 1e-9 {
+                continue;
+            }
+            if let Some((tx, ty)) = mask.traction_at(
+                body, &field.u, &field.v, &field.p, mu, 0.0, s.x, s.y, s.nx, s.ny,
+            ) {
+                drag += tx * s.ds;
+                lift += ty * s.ds;
+            }
+        }
+        (drag, lift)
+    };
+    // The fluid harness check: surface drag on cylinder + flag near the
+    // CFD2 value this solver measured on this geometry (ny = 62: ~121;
+    // the reference is 136.700 with the boundary layer barely a cell).
+    let (rigid_drag, rigid_lift) = measure_force(&solver, &field);
+    println!(
+        "  rigid phase: {rigid_steps} steps to t = {t_release:.1} s in {:.0} s wall; \
+         surface drag {rigid_drag:.1} (CFD2 ref 136.7, this grid measured ~121), \
+         lift {rigid_lift:.1}",
+        start.elapsed().as_secs_f64()
+    );
+
+    // The flag: nonlinear Newmark stepper at the coupled dt.
+    let mut db = MaterialDatabase::new();
+    db.add_material(
+        MaterialId(0),
+        LinearElastic::new(E_S, NU_S).with_density(RHO_S),
+        None,
+    );
+    // A deep Newton budget: a mid-swing subiteration can hand the flag a
+    // large sudden load change (the coupled lift swings hundreds of N
+    // within a period); typical steps converge in 1-2 iterations, and a
+    // t = 25.8 s failure at the default budget of 25 is what set this.
+    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()
+    });
+    let flag = RefCell::new(analysis.stepper().unwrap());
+    let wetted_dofs: Vec<[usize; 2]> = interface
+        .wetted
+        .iter()
+        .map(|&id| {
+            let dofs = flag.borrow().node_dofs(id);
+            [dofs[0], dofs[1]]
+        })
+        .collect();
+    let a_dofs = flag.borrow().node_dofs(a_node);
+    let extract = |state: &DynamicState| -> Vec<f64> {
+        let mut d = vec![0.0; 2 * wetted_dofs.len()];
+        for (k, dofs) in wetted_dofs.iter().enumerate() {
+            d[2 * k] = state.displacement[dofs[0]];
+            d[2 * k + 1] = state.displacement[dofs[1]];
+        }
+        d
+    };
+
+    // Tractions on the flag's wetted surface for a given geometry, from
+    // the solver's current field/mask; returns the transferred nodal
+    // forces, the conservation defect, and the samples dropped (probe
+    // failures plus spike rejections).
+    let spiked_total = std::cell::Cell::new(0usize);
+    let sample_load = |solver: &EmbeddedPisoSolver,
+                       field: &FlowField,
+                       d: &[f64]|
+     -> (Vec<(NodeId, Vector3<f64>)>, f64, usize) {
+        let vertices = interface.polygon(d);
+        let poly_probe = EmbeddedBody::polygon(vertices);
+        let mask = solver.mask().unwrap();
+        let body = solver.body().unwrap();
+        let mut faces = Vec::new();
+        let mut tractions: Vec<Vector3<f64>> = Vec::new();
+        let mut skipped = 0usize;
+        for s in poly_probe.surface_samples(0.5 * h) {
+            if circle_sdf(s.x, s.y) < 1e-9 {
+                continue; // buried in the cylinder
+            }
+            match mask.traction_at(
+                body, &field.u, &field.v, &field.p, mu, 0.0, s.x, s.y, s.nx, s.ny,
+            ) {
+                Some((tx, ty)) => {
+                    faces.push(FluidFace {
+                        centroid: Vector3::new(s.x, s.y, 0.0),
+                        normal: Vector3::new(s.nx, s.ny, 0.0),
+                        area: s.ds,
+                    });
+                    tractions.push(Vector3::new(tx, ty, 0.0));
+                }
+                None => skipped += 1,
+            }
+        }
+        // Spike guard: a near-degenerate reconstruction can return a
+        // finite but wild traction (the linear-fit condition sits just
+        // above its truncation threshold at concave junctions). CLAMP
+        // samples to 20x the median magnitude, keeping their direction —
+        // the physical load varies smoothly along the surface — and COUNT
+        // them: a non-zero count is a measurement of the pathology, not a
+        // silent repair. Clamping, not dropping: a hard drop threshold
+        // makes the coupling pass discontinuous in the candidate geometry
+        // (a boundary sample flips in/out of the kept set between
+        // subiterations, and the load jumps by the spike magnitude —
+        // measured as a residual bouncing at the scale of the step
+        // increment); the clamp is continuous.
+        let mut magnitudes: Vec<f64> = tractions.iter().map(nalgebra::Vector3::norm).collect();
+        magnitudes.sort_by(|a, b| a.partial_cmp(b).unwrap());
+        let median = magnitudes.get(magnitudes.len() / 2).copied().unwrap_or(0.0);
+        if median > 0.0 {
+            let cap = 20.0 * median;
+            for traction in &mut tractions {
+                let norm = traction.norm();
+                if norm > cap {
+                    *traction *= cap / norm;
+                    spiked_total.set(spiked_total.get() + 1);
+                }
+            }
+        }
+        let nodes_now = interface.deformed_nodes(d);
+        let surface = WettedSurface::build(&faces, &nodes_now).expect("transfer build");
+        let nodal = surface.transfer_load(&faces, &tractions).unwrap();
+        let total_sampled: Vector3<f64> =
+            faces.iter().zip(&tractions).map(|(f, t)| t * f.area).sum();
+        let total_nodal: Vector3<f64> = nodal.iter().sum();
+        let conservation = (total_nodal - total_sampled).norm() / total_sampled.norm().max(1e-30);
+        (
+            interface
+                .wetted
+                .iter()
+                .zip(nodal)
+                .map(|(&id, f)| (id, f))
+                .collect(),
+            conservation,
+            skipped,
+        )
+    };
+
+    // Phase 2: release. The flag starts at rest under the current fluid
+    // load (consistent initial acceleration — the step response about the
+    // steady deflection is the seed perturbation for the instability).
+    let (nodal0, conservation0, _) = sample_load(&solver, &field, &zero_d);
+    flag.borrow_mut().set_nodal_forces(&nodal0);
+    let mut flag_state = flag.borrow_mut().rest_state().unwrap();
+    let mut committed_nodal = nodal0;
+    let mut worst_conservation = conservation0;
+
+    let solver = RefCell::new(solver);
+    let field = RefCell::new(field);
+
+    let coupled_steps = ((t_end - t_release) / dt).round() as usize;
+    let mut times = Vec::with_capacity(coupled_steps);
+    let mut ux_series = Vec::with_capacity(coupled_steps);
+    let mut uy_series = Vec::with_capacity(coupled_steps);
+    let mut total_subiterations = 0usize;
+    let mut max_subiterations = 0usize;
+    let mut total_skipped = 0usize;
+    let mut stalled_steps = 0usize;
+    let mut worst_stall = 0.0f64;
+    let mut force_times: Vec<f64> = Vec::new();
+    let mut drag_series: Vec<f64> = Vec::new();
+    let mut lift_series: Vec<f64> = Vec::new();
+    let mut csv = csv_path.map(|p| std::fs::File::create(p).expect("csv path"));
+
+    let phase_start = std::time::Instant::now();
+    for step in 0..coupled_steps {
+        let d_n = extract(&flag_state);
+
+        // Predictor: the structure alone under the committed load.
+        flag.borrow_mut().set_nodal_forces(&committed_nodal);
+        let (predicted, _) = flag.borrow_mut().step(&flag_state).unwrap();
+        let d_predicted = extract(&predicted);
+
+        let fluid_saved = solver.borrow().snapshot();
+        let field_saved = field.borrow().clone();
+        type PassResult = (
+            FlowField,
+            DynamicState,
+            Vec<(NodeId, Vector3<f64>)>,
+            f64,
+            usize,
+        );
+        let latest: RefCell<Option<PassResult>> = RefCell::new(None);
+
+        let pass = |d_candidate: &[f64]| -> Vec<f64> {
+            // Interface velocity of THIS candidate, constant over the step.
+            let ddot: Vec<f64> = d_candidate
+                .iter()
+                .zip(&d_n)
+                .map(|(new, old)| (new - old) / dt)
+                .collect();
+            // Subcycled fluid steps from the SAME start-of-step state,
+            // geometry interpolated to each substep's end time.
+            let mut solver_ref = solver.borrow_mut();
+            solver_ref.restore(&fluid_saved);
+            let mut trial_field = field_saved.clone();
+            for m in 1..=subcycle {
+                let fraction = m as f64 / subcycle as f64;
+                let d_sub: Vec<f64> = d_n
+                    .iter()
+                    .zip(d_candidate)
+                    .map(|(old, new)| old + fraction * (new - old))
+                    .collect();
+                {
+                    let mut geometry = shared.write().unwrap();
+                    geometry.0 = interface.polygon(&d_sub);
+                    geometry.1 = interface.walk_velocities(&ddot);
+                }
+                futures::executor::block_on(solver_ref.advance(&mut trial_field, dt_fluid))
+                    .unwrap();
+            }
+            // Load on the candidate geometry, flag answers from the
+            // committed state.
+            let (nodal, conservation, skipped) =
+                sample_load(&solver_ref, &trial_field, d_candidate);
+            let mut flag_ref = flag.borrow_mut();
+            flag_ref.set_nodal_forces(&nodal);
+            let (candidate_state, _) = flag_ref.step(&flag_state).unwrap();
+            let d_new = extract(&candidate_state);
+            *latest.borrow_mut() =
+                Some((trial_field, candidate_state, nodal, conservation, skipped));
+            d_new
+        };
+
+        let increment: f64 = d_predicted
+            .iter()
+            .zip(&d_n)
+            .map(|(a, b)| (a - b) * (a - b))
+            .sum::<f64>()
+            .sqrt();
+        let tol_step = tol_floor.max(rtol * increment);
+        let mut scheme = Subiterated::aitken(max_subiterations_budget, tol_step).unwrap();
+        match scheme.solve(&d_predicted, pass) {
+            Ok(converged) => {
+                total_subiterations += converged.iterations;
+                max_subiterations = max_subiterations.max(converged.iterations);
+            }
+            Err(
+                rtx_fsi::FsiError::CouplingNotConverged {
+                    iterations,
+                    residual,
+                    ..
+                }
+                | rtx_fsi::FsiError::CouplingDiverged {
+                    iterations,
+                    residual,
+                },
+            ) if residual < 5.0 * tol_step => {
+                // The noise floor, not divergence: accept the last
+                // candidate, count it, and bound it at the end.
+                stalled_steps += 1;
+                worst_stall = worst_stall.max(residual);
+                total_subiterations += iterations;
+                max_subiterations = max_subiterations.max(iterations);
+            }
+            Err(e) => panic!("coupling failed at step {step}: {e:?}"),
+        }
+        // `latest` holds the response to the accepted interface (the last
+        // pass) — commit it directly; the fluid, mask and flag are
+        // consistent with that interface without an extra pass.
+        let (new_field, new_flag_state, nodal, conservation, skipped) =
+            latest.borrow_mut().take().expect("pass ran");
+        *field.borrow_mut() = new_field;
+        flag_state = new_flag_state;
+        committed_nodal = nodal;
+        worst_conservation = worst_conservation.max(conservation);
+        total_skipped += skipped;
+
+        let t = t_release + (step + 1) as f64 * dt;
+        let ux = flag_state.displacement[a_dofs[0]];
+        let uy = flag_state.displacement[a_dofs[1]];
+        times.push(t);
+        ux_series.push(ux);
+        uy_series.push(uy);
+        if (step + 1) % 10 == 0 {
+            let (drag, lift) = measure_force(&solver.borrow(), &field.borrow());
+            force_times.push(t);
+            drag_series.push(drag);
+            lift_series.push(lift);
+            if let Some(file) = csv.as_mut() {
+                writeln!(file, "{t:.6},{ux:.6e},{uy:.6e},{drag:.6e},{lift:.6e}").unwrap();
+            }
+        } else if let Some(file) = csv.as_mut() {
+            writeln!(file, "{t:.6},{ux:.6e},{uy:.6e},,").unwrap();
+        }
+        if (step + 1) % 1000 == 0 {
+            let window = &uy_series[uy_series.len().saturating_sub(1000)..];
+            let (w_mid, w_amp) = mid_amp(window);
+            println!(
+                "  t = {t:.3} s ({step} steps): uy(A) = {uy:.3e} (window mid {w_mid:.3e} \
+                 amp {w_amp:.3e}), {:.1} subit/step, {:.0} s wall",
+                total_subiterations as f64 / (step + 1) as f64,
+                phase_start.elapsed().as_secs_f64()
+            );
+        }
+    }
+
+    let elapsed = start.elapsed().as_secs_f64();
+    let mean_subiterations = total_subiterations as f64 / coupled_steps.max(1) as f64;
+
+    // Measure over the last three seconds (or the last half, if shorter).
+    let window_start = times
+        .iter()
+        .position(|&t| t >= t_end - 3.0)
+        .unwrap_or(times.len() / 2);
+    let uy_window = &uy_series[window_start..];
+    let ux_window = &ux_series[window_start..];
+    let t_window = &times[window_start..];
+    let (uy_mid, uy_amp) = mid_amp(uy_window);
+    let (ux_mid, ux_amp) = mid_amp(ux_window);
+    let frequency = crossing_frequency(t_window, uy_window);
+
+    // Onset: amplitude of the first quarter of the coupled march vs the
+    // last quarter.
+    let quarter = uy_series.len() / 4;
+    let (_, amp_early) = mid_amp(&uy_series[..quarter.max(1)]);
+    let (_, amp_late) = mid_amp(&uy_series[uy_series.len() - quarter.max(1)..]);
+
+    // Loads over the same window.
+    let force_start = force_times
+        .iter()
+        .position(|&t| t >= t_end - 3.0)
+        .unwrap_or(force_times.len() / 2);
+    let (drag_mid, drag_amp) = mid_amp(&drag_series[force_start..]);
+    let (lift_mid, lift_amp) = mid_amp(&lift_series[force_start..]);
+    println!(
+        "  loads over the window: drag {drag_mid:.2} ± {drag_amp:.2} (ref {REF_DRAG_MEAN} ± \
+         77.65), lift {lift_mid:.2} ± {lift_amp:.2} (ref 0.61 ± {REF_LIFT_AMP})"
+    );
+
+    println!(
+        "  FSI2 (fluid ny = {ny}, flag {flag_nx}x2 Quad8, dt = {dt:.2e}): coupled {coupled_steps} \
+         steps in {:.0} s wall total; {mean_subiterations:.1} subit/step (max \
+         {max_subiterations}); {stalled_steps} stalled steps (worst residual \
+         {worst_stall:.2e}); worst conservation {worst_conservation:.2e}; skipped samples \
+         {total_skipped} (of which {} spike-rejected)\n  measured over [{:.1}, {t_end:.1}] s: uy(A) = {:.4} ± {:.4} mm \
+         (ref {:.2} ± {:.1}), ux(A) = {:.4} ± {:.4} mm (ref {:.2} ± {:.2}), f = {} Hz \
+         (ref {REF_UY_FREQ}); onset amp {:.3e} -> {:.3e} m",
+        elapsed,
+        spiked_total.get(),
+        t_window.first().unwrap_or(&t_release),
+        uy_mid * 1e3,
+        uy_amp * 1e3,
+        REF_UY_MEAN * 1e3,
+        REF_UY_AMP * 1e3,
+        ux_mid * 1e3,
+        ux_amp * 1e3,
+        REF_UX_MEAN * 1e3,
+        REF_UX_AMP * 1e3,
+        frequency.map_or("n/a".to_string(), |f| format!("{f:.3}")),
+        amp_early,
+        amp_late,
+    );
+    // Machinery invariants — asserted at every resolution.
+    assert!(
+        worst_conservation < 1e-10,
+        "load transfer lost force: {worst_conservation:.3e}"
+    );
+    assert!(
+        flag_state.displacement.iter().all(|v| v.is_finite()),
+        "flag state went non-finite"
+    );
+    assert!(
+        mean_subiterations < 10.0,
+        "coupling is grinding: {mean_subiterations:.1} subiterations/step"
+    );
+    // Stalls at the noise floor are tolerated but must stay the exception;
+    // a coupling stalling on most steps is not converging, it is drifting.
+    assert!(
+        stalled_steps * 5 < coupled_steps.max(1),
+        "coupling stalled on {stalled_steps} of {coupled_steps} steps \
+         (worst residual {worst_stall:.2e})"
+    );
+
+    let _ = (amp_early, amp_late);
+
+    // Physics bands, by horizon. The march is deterministic, so short
+    // horizons carry tight regression bands; long horizons pin the
+    // MEASURED loosely-coupled attractor — not benchmark agreement (see
+    // the module docs: the reference's mode-2 resonant cycle at 1.93 Hz /
+    // 81.6 mm is not reached by this coupling; the measured state is the
+    // wake-forced 3.73 Hz / ±17.3 mm cycle at BOTH grids). If a change
+    // moves these numbers, that is a finding either way and must be loud.
+    if ny == 62 && flag_nx == 35 && (t_end - 7.0).abs() < 1e-9 && (t_release - 6.0).abs() < 1e-9 {
+        // The committed default: the release response over [6, 7] s,
+        // measured 2026-08-21 as uy mid 3.773 mm, amp 3.792 mm. The band
+        // is ±35% for cross-platform floating-point drift in a growing
+        // transient, not an accuracy claim.
+        assert!(
+            (2.4e-3..5.2e-3).contains(&uy_mid),
+            "uy release-response mid {uy_mid:.4e} outside the measured band \
+             [2.4e-3, 5.2e-3]"
+        );
+        assert!(
+            (2.4e-3..5.2e-3).contains(&uy_amp),
+            "uy release-response amp {uy_amp:.4e} outside the measured band \
+             [2.4e-3, 5.2e-3]"
+        );
+    } else if t_end >= 25.0 {
+        // Study horizons: the measured attractor of the loosely-coupled
+        // (subcycle 8) march — f = 3.729 / 3.728 Hz and uy amp 17.3 mm at
+        // ny = 62 / 82 (2026-08-21).
+        if let Some(f) = frequency {
+            assert!(
+                (f - 3.73).abs() / 3.73 < 0.10,
+                "uy frequency {f:.3} left the measured 3.73 Hz attractor \
+                 (benchmark reference {REF_UY_FREQ}) — a material change"
+            );
+        }
+        assert!(
+            (12e-3..24e-3).contains(&uy_amp),
+            "uy amplitude {uy_amp:.4e} left the measured ±17.3 mm attractor \
+             band [12e-3, 24e-3]"
+        );
+    }
+}