rtx-fsi: C3 — FSI3 opened: the added-mass regime, its impulse artefact, and the coupler hygiene it demanded

The march is extracted from the FSI2 test into a shared
tests/fsi2_harness/march.rs (MarchConfig from env / MarchResult /
window statistics), the harness parameterised by BenchmarkCase (FSI2,
FSI3: inflow, solid density, modulus, rigid-flag drag reference), and
turek_hron_fsi3.rs written on top (Re 200, density ratio 1, E 5.6e6).
The committed FSI2 default is bit-identical through all of it —
every new knob defaults to FSI2's behaviour; re-verified twice to every
printed digit (uy 3.7732 / 3.7920 mm).

FSI3's first contact, traced pass by pass (RTX_FSI3_TRACE), exposed
the added-mass instability in its purest form and, one measured
mechanism at a time, what a partitioned coupling at unit density ratio
needs:

- C^1 interface motion (MarchConfig::c1_interface): a constant
  per-step interface velocity is a velocity JUMP at each step start,
  and the incompressible fluid answers with an impulsive added-mass
  load ~ rho L dv / dt_fluid (8x the physical reaction under
  subcycling): 1,600 N at release, 48,000 N and a 59 mm response one
  step later, the flag's Newton dead the pass after. Constant
  acceleration from the previous end velocity to 2 dd/dt - v_n removes
  the impulse (loads 1,700-2,400 N).
- IQN first-pass relaxation as a knob (initial_relaxation, FSI3 0.05):
  |1 - omega (1 + g)| must contract; 0.5 diverges past gain 3.
- The divergence verdict waits for the secant (IqnIls): the exploratory
  second pass on a high-gain map legitimately overshoots 10x before the
  first secant column exists. Pinned by a gain-40 model test.
- Kinematic predictor (predictor: "kinematic", velocity only): the
  structure-alone predictor ignores an added mass comparable to the
  flag's and overshoots 2-5x, drawing 5-6x loads every first pass; and
  NOT with the acceleration — Newmark average acceleration carries an
  inconsistent initial acceleration as a sign-alternating mode
  (d + dt v + dt^2 a / 2 predicted 22 mm at release; converged 0.14).
- Quiescent release (quiescent_release): the structure-alone
  "consistent" initial acceleration M^-1 F ignores the added mass.
- Coupler hygiene (IqnIls): a stalled or diverged step's secant columns
  are no longer retained (a bistable mask flip's columns extrapolated a
  30 mm interface jump on the next step); two-window stagnation
  detection reports a plateau early instead of bouncing to the budget
  (a single-window test misjudged a slowly converging step and is
  recorded as such); trust region tightened to 10x the residual. A
  noise-column filter at the tolerance was measured to HURT (stalled a
  converging step at 5.5e-4) and is disabled (threshold 0).
- The floor measured, not borrowed (fsi2_interface_noise.rs gains
  RTX_NOISE_CASE=fsi3): flip jumps 3.6e-4 (12x FSI2's), the subcycle-8
  release map stalling near 1e-6, the subcycle-2 map converging to
  9e-10; and through the release transient (flag at ~0.3 m/s) the
  subcycle-8 floor rides up to ~1e-3 with the motion — which moved the
  FSI3 default to subcycle 2 at a 3e-5 floor.

FSI3 committed default (ny 62, release t = 4, t_end 4.2, subcycle 2):
581 steps, 3.8 subit/step (max 8), 0 stalls, 0 retries, conservation
4.3e-13; uy 10.1 +- 27.5 mm and ux -3.1 +- 3.2 mm over the first 0.2 s
(reference cycle 1.45 +- 34.90, -2.86 +- 2.70). Machinery invariants
pinned; physics bands await the study record. Rigid-flag drag 426.9 vs
CFD3's 439.45 (-2.9%).

The FSI2 mode-2 study pin (IQN / subcycle 2, t_end >= 16) was measured
with the pre-hygiene coupler; re-verify on the next s = 2 study run.

48 lib tests green, clippy clean.

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

crates/specialized/rtx-fsi/tests/fsi2_harness/mod.rs

@@ -8,7 +8,9 @@
 //! Everything here is code motion from the tenth-session FSI2 test —
 //! the physics and defaults are unchanged unless a test says otherwise.
 
-#![allow(dead_code)] // two test crates share this; each uses a subset
+#![allow(dead_code)] // several test crates share this; each uses a subset
+
+pub mod march;
 
 use std::cell::Cell;
 use std::sync::{Arc, RwLock};
@@ -40,18 +42,61 @@ pub const FLAG_X1: f64 = 0.6;
 pub const FLAG_Y0: f64 = 0.19;
 pub const FLAG_Y1: f64 = 0.21;
 
+/// The parameters that distinguish the self-excited Turek–Hron cases on
+/// the shared geometry: mean inflow (Re = 100 U), solid density and
+/// Young's modulus. Everything else — channel, cylinder, flag, fluid —
+/// is common.
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct BenchmarkCase {
+    pub name: &'static str,
+    pub u_mean: f64,
+    pub rho_s: f64,
+    pub e_s: f64,
+    pub nu_s: f64,
+    /// The rigid-flag CFD drag on this geometry at this Re (the fluid
+    /// harness check before anything couples): CFD2 / CFD3 means.
+    pub rigid_drag_reference: f64,
+}
+
+/// FSI2: Re 100, density ratio 10 — the heavy flag's resonant flapping.
+pub const FSI2: BenchmarkCase = BenchmarkCase {
+    name: "FSI2",
+    u_mean: 1.0,
+    rho_s: 10_000.0,
+    e_s: 1.4e6,
+    nu_s: 0.4,
+    rigid_drag_reference: 136.7,
+};
+
+/// FSI3: Re 200, density ratio 1 (mu_s = 2e6 → E = 5.6e6) — the
+/// added-mass regime.
+pub const FSI3: BenchmarkCase = BenchmarkCase {
+    name: "FSI3",
+    u_mean: 2.0,
+    rho_s: 1_000.0,
+    e_s: 5.6e6,
+    nu_s: 0.4,
+    rigid_drag_reference: 439.45,
+};
+
 pub 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.
-pub fn inflow(y: f64, t: f64) -> f64 {
+/// The ramped parabolic inflow of the benchmark definition, for a mean
+/// inflow `u_mean`.
+pub fn inflow_for(u_mean: f64, y: f64, t: f64) -> f64 {
     let ramp = if t < 2.0 {
         0.5 * (1.0 - (std::f64::consts::PI * t / 2.0).cos())
     } else {
         1.0
     };
-    ramp * 1.5 * U_MEAN * y * (H - y) / (0.5 * H).powi(2)
+    ramp * 1.5 * u_mean * y * (H - y) / (0.5 * H).powi(2)
+}
+
+/// FSI2's inflow.
+pub fn inflow(y: f64, t: f64) -> f64 {
+    inflow_for(U_MEAN, y, t)
 }
 
 pub fn env_or(name: &str, default: f64) -> f64 {
@@ -260,6 +305,7 @@ pub fn crossing_frequency(times: &[f64], series: &[f64]) -> Option<f64> {
 /// lock, and the load sampling with its spike clamp and (optional)
 /// surface smoothing.
 pub struct Fsi2Harness {
+    pub case: BenchmarkCase,
     pub mesh: Mesh,
     pub interface: Interface,
     pub a_node: NodeId,
@@ -284,11 +330,23 @@ impl Fsi2Harness {
         ny: usize,
         flag_nx: usize,
         smooth_in_h: f64,
+    ) -> (Self, EmbeddedPisoSolver, FlowField) {
+        Self::build_case(FSI2, ny, flag_nx, smooth_in_h)
+    }
+
+    /// Build the harness for a benchmark case (FSI2 or FSI3 parameters
+    /// on the shared geometry).
+    pub fn build_case(
+        case: BenchmarkCase,
+        ny: usize,
+        flag_nx: usize,
+        smooth_in_h: f64,
     ) -> (Self, EmbeddedPisoSolver, FlowField) {
         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;
+        let u_mean = case.u_mean;
+        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.
         let dt_fluid = 0.25 / (2.0 * u_peak / h + 4.0 * NU_F / (h * h));
@@ -315,7 +373,7 @@ impl Fsi2Harness {
         let config = CfdConfig::new()
             .with_density(RHO_F)
             .with_viscosity(mu)
-            .with_reference_velocity(U_MEAN)
+            .with_reference_velocity(u_mean)
             .with_reference_length(0.1);
         let params = EmbeddedParameters {
             corrector_steps: 2,
@@ -336,9 +394,9 @@ impl Fsi2Harness {
             convection_scheme: ConvectionScheme::TvdVanAlbada,
         };
         let mut solver = EmbeddedPisoSolver::new(config, params).unwrap();
-        solver.set_boundary_velocity(|x, y, t| {
+        solver.set_boundary_velocity(move |x, y, t| {
             if x <= 0.0 {
-                (inflow(y, t), 0.0)
+                (inflow_for(u_mean, y, t), 0.0)
             } else {
                 (0.0, 0.0)
             }
@@ -363,6 +421,7 @@ impl Fsi2Harness {
         solver.initialize(&mut field).unwrap();
 
         let harness = Self {
+            case,
             mesh,
             interface,
             a_node,
@@ -518,21 +577,51 @@ impl Fsi2Harness {
         d_n: &[f64],
         d_candidate: &[f64],
         subcycle: usize,
+        v_n: Option<&[f64]>,
     ) {
         let dt = self.dt_fluid * subcycle as f64;
-        let ddot: Vec<f64> = d_candidate
+        let mean_velocity: Vec<f64> = d_candidate
             .iter()
             .zip(d_n)
             .map(|(new, old)| (new - old) / dt)
             .collect();
         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();
-            self.set_geometry(&d_sub, &ddot);
+            let (d_sub, ddot_sub): (Vec<f64>, Vec<f64>) = match v_n {
+                // Constant velocity over the step: the geometry moves
+                // linearly and the wall velocity JUMPS at the step
+                // boundary — harmless for a heavy flag, but the
+                // incompressible fluid answers a velocity jump with an
+                // impulsive added-mass load ~ rho L dv / dt_fluid, which
+                // at unit density ratio destroyed the flag in one step.
+                None => (
+                    d_n.iter()
+                        .zip(d_candidate)
+                        .map(|(old, new)| old + fraction * (new - old))
+                        .collect(),
+                    mean_velocity.clone(),
+                ),
+                // C^1 interface motion: constant acceleration across the
+                // step from the previous end-of-step velocity to the
+                // trapezoidal end velocity 2 dd/dt - v_n (Newmark
+                // average acceleration's own kinematics), so the wall
+                // velocity is continuous at the step boundary and the
+                // impulse is gone. The end-of-substep velocity goes with
+                // the end-of-substep geometry.
+                Some(v_start) => {
+                    let mut d_sub = Vec::with_capacity(d_n.len());
+                    let mut ddot_sub = Vec::with_capacity(d_n.len());
+                    for k in 0..d_n.len() {
+                        let v_end = 2.0 * mean_velocity[k] - v_start[k];
+                        let accel = (v_end - v_start[k]) / dt;
+                        let tau = fraction * dt;
+                        d_sub.push(d_n[k] + v_start[k] * tau + 0.5 * accel * tau * tau);
+                        ddot_sub.push(v_start[k] + accel * tau);
+                    }
+                    (d_sub, ddot_sub)
+                }
+            };
+            self.set_geometry(&d_sub, &ddot_sub);
             futures::executor::block_on(solver.advance(field, self.dt_fluid)).unwrap();
         }
     }