test(rtx-fsi): coupling-level rescue rung A (RTX_FSI{2,3}_CRESCUE, default off) + TRACE_FROM autopsy window
CI / Build CPU-Only (Explicit) (push) Canceled after 0s
Documentation / Build User Guide (push) Canceled after 0s
Performance Benchmarks / Run Benchmarks (push) Canceled after 0s
CI / Format Check (push) Canceled after 0s
CI / Clippy Check (push) Canceled after 0s
CI / Build (macos-latest) (push) Canceled after 0s
CI / Test (macos-latest) (push) Canceled after 0s
CI / Build (ubuntu-latest) (push) Canceled after 0s
CI / Test (ubuntu-latest) (push) Canceled after 0s
Documentation / Build API Documentation (push) Canceled after 0s
CI / Python Bindings (maturin) (macos-latest) (push) Canceled after 0s
CI / Python Bindings (maturin) (ubuntu-latest) (push) Canceled after 0s
CI / WASM Build + Size Check (push) Canceled after 0s
CI / Distributed Training Tests (push) Canceled after 0s
CI / CI Success (push) Canceled after 0s

Coupling-rescue campaign (omni-cortex docs/coupling_rescue_campaign.md).
Every recorded FSI2 s=1 coupling death is the END of a multi-step
runaway (74-march burst scan: green marches never exceed 2.4x their own
p95 tip jump; 27/28 deaths burst 5-37 steps first). The traced autopsy
of the (1.27, 2.0) death (replay digit-identical, CSV byte-identical)
shows a growing period-2 instability of the CONVERGED coupled scheme at
the tip's max-velocity crossing: increment x33 and converged load
1e4 -> 1e6 N over 200 steps with the coupling converging on 111 of the
first 112 steps.

- march.rs: RTX_{prefix}_TRACE_FROM (print-only autopsy window: per-pass
  residual/load + one line per step with increment, tol_step, retry_at,
  acceptable, outcome, committed tip jump); RTX_{prefix}_CRESCUE (default
  off): on a fatal stall or a committed tip jump > 3x the running p95
  (trailing 2000 non-rescued steps), reject the step and repeat the
  interval as 2/4/8/16/32 coupled substeps of dt/n (Mayr-Wall-Gee
  reduced-step repetition, five repetitions); per-rescue record printed,
  MarchResult.coupling_rescues/_failures/rescue_records, rate cap 20
  rescues per second of march (dies loudly).
- rescue.rs (new): the substep ladder — each substep a complete coupled
  step at dt/n with its own predictor, fresh coupler, C1 velocity
  chaining, tolerances at the substep's increment; the march's own step
  path is NOT routed through it (digit identity by construction).
- rtx-fea NonlinearDynamicStepper::step_with_dt (step at an explicit dt;
  step() delegates float-for-float); Fsi2Harness::advance_subcycled_with
  (explicit fluid dt; advance_subcycled delegates).

Verified: FSI2 + FSI3 committed defaults and the noise probe
digit-identical knob-off vs same-day / 2026-08-31 baselines; knob ON
on the anchor u=1.00 r=1: green, zero rescues, CSV byte-identical to
TWIN-1's u1.00.csv; fmt + clippy clean on the touched files.

Verdict of rung A on the provocation set: REFUTED 3/3 by mechanism —
every rescued interval was carried but reproduced the rejected step's
motion (dt/2..dt/32 give the same jump), so the runaway lives in the
coupled LOAD at the crossing, not in the time integration; the
SUBCYCLE=2 closure marches through the same crossing. The knob stays,
default off, as the instrument that measured this.

Co-Authored-By: Claude Fable 5.1 <[email protected]>
Claude-Session: https://claude.ai/code/session_01X2GmJXeQ2njUecEKiJZ1G2
This commit is contained in:
Omar Sobh
2026-09-02 21:36:06 -07:00
co-authored by Claude Fable 5.1
parent 0f578087ce
commit e76271ac67
6 changed files with 584 additions and 22 deletions
@@ -101,6 +101,27 @@ pub struct MarchConfig {
/// Print every pass's interface residual for the first `trace_steps`
/// coupled steps (diagnostics; 0 = off).
pub trace_steps: usize,
/// Trace window for a death autopsy (`RTX_{prefix}_TRACE_FROM`,
/// default `usize::MAX` = off): from this coupled step onward,
/// print every pass like `trace_steps` does PLUS one line per step
/// with the predictor increment, the step tolerance, the two
/// acceptance thresholds, the coupling outcome and the committed
/// tip jump. Reporting-only — the knob touches no float on the
/// solver path, so a traced replay is digit-identical to the
/// untraced march (the property every autopsy so far has leaned
/// on). Motivation: every recorded coupling death is the END of a
/// multi-step runaway (committed tip jumps 1070× the march's own
/// p95 for 537 steps before the panic; `scripts/fsi_burst_scan.py`
/// in omni-cortex), and `trace_steps` only sees the release.
pub trace_from: usize,
/// Coupling-level rescue (`RTX_{prefix}_CRESCUE`, default off =
/// bit-identical): a step whose coupling ends above its acceptance
/// window (today's death), or whose committed tip jump exceeds
/// 3× the march's running p95 of |Δuy| (the runaway signature
/// every recorded death carries for 537 steps before its panic),
/// is repeated as 2, 4, 8, 16, 32 coupled substeps of dt/n — see
/// `rescue.rs` and omni-cortex `docs/coupling_rescue_campaign.md`.
pub coupling_rescue: bool,
/// C^1 interface motion (constant acceleration across the step from
/// the previous end velocity) instead of a constant velocity with a
/// jump at the step boundary. See `Fsi2Harness::advance_subcycled`.
@@ -151,6 +172,11 @@ impl MarchConfig {
ffld_dir: std::env::var(key("FFLD")).ok().or(defaults.ffld_dir),
initial_relaxation: num("OMEGA0", defaults.initial_relaxation),
trace_steps: num("TRACE", defaults.trace_steps as f64) as usize,
trace_from: std::env::var(key("TRACE_FROM"))
.ok()
.and_then(|v| v.parse::<usize>().ok())
.unwrap_or(defaults.trace_from),
coupling_rescue: num("CRESCUE", f64::from(u8::from(defaults.coupling_rescue))) != 0.0,
c1_interface: num("C1", f64::from(u8::from(defaults.c1_interface))) != 0.0,
predictor: std::env::var(key("PREDICTOR")).unwrap_or(defaults.predictor),
quiescent_release: num("QUIESCENT", f64::from(u8::from(defaults.quiescent_release)))
@@ -211,6 +237,13 @@ pub struct MarchResult {
/// defeated the plain SVK Newton, the failure that used to kill the
/// march outright.
pub newton_rescues: (usize, usize),
/// Coupling-level rescues (steps repeated as substeps — see
/// `rescue.rs`); zero with the knob off and on every healthy march.
pub coupling_rescues: usize,
/// Kinematic-trigger rescues whose whole ladder failed (the step
/// was kept as the coupling accepted it).
pub coupling_rescue_failures: usize,
pub rescue_records: Vec<super::rescue::RescueRecord>,
pub final_state_finite: bool,
pub elapsed: f64,
}
@@ -293,6 +326,8 @@ pub fn run_march(case: BenchmarkCase, config: &MarchConfig) -> MarchResult {
ref ffld_dir,
initial_relaxation,
trace_steps,
trace_from,
coupling_rescue,
c1_interface,
ref predictor,
quiescent_release,
@@ -424,6 +459,13 @@ pub fn run_march(case: BenchmarkCase, config: &MarchConfig) -> MarchResult {
// trajectory.
let mut interval_drag: Vec<f64> = Vec::new();
let mut interval_lift: Vec<f64> = Vec::new();
// Coupling-level rescue bookkeeping (all inert with the knob off).
let mut coupling_rescues = 0usize;
let mut coupling_rescue_failures = 0usize;
let mut rescue_records: Vec<super::rescue::RescueRecord> = Vec::new();
let mut rescue_steps: Vec<usize> = Vec::new();
let mut jump_window: std::collections::VecDeque<f64> = std::collections::VecDeque::new();
let mut running_p95: Option<f64> = None;
let mut csv = csv_path
.as_ref()
.map(|p| std::fs::File::create(p).expect("csv path"));
@@ -470,6 +512,9 @@ pub fn run_march(case: BenchmarkCase, config: &MarchConfig) -> MarchResult {
let fluid_saved = solver.borrow().snapshot();
let field_saved = field.borrow().clone();
t_save += save_start.elapsed().as_secs_f64();
// The rescue repeats the interval from the committed start state.
let step_start_state = coupling_rescue.then(|| flag_state.clone());
let step_start_nodal = coupling_rescue.then(|| committed_nodal.clone());
type PassResult = (
FlowField,
DynamicState,
@@ -502,7 +547,7 @@ pub fn run_march(case: BenchmarkCase, config: &MarchConfig) -> MarchResult {
let (nodal, conservation, skipped) =
harness.sample_load(&solver_ref, &trial_field, d_candidate);
t_sample.set(t_sample.get() + sample_start.elapsed().as_secs_f64());
if step < trace_steps {
if step < trace_steps || step >= trace_from {
let load: f64 = nodal.iter().map(|(_, f)| f.norm()).sum();
let peak = nodal.iter().map(|(_, f)| f.norm()).fold(0.0, f64::max);
println!(
@@ -517,7 +562,7 @@ pub fn run_march(case: BenchmarkCase, config: &MarchConfig) -> MarchResult {
let (candidate_state, _) = flag_ref.step(&flag_state).unwrap();
t_structure.set(t_structure.get() + structure_start.elapsed().as_secs_f64());
let d_new = extract(&candidate_state);
if step < trace_steps {
if step < trace_steps || step >= trace_from {
let residual: f64 = d_new
.iter()
.zip(d_candidate)
@@ -586,6 +631,62 @@ pub fn run_march(case: BenchmarkCase, config: &MarchConfig) -> MarchResult {
outcome = iqn_ref.solve(&d_predicted, pass);
}
}
let trace_line = (step >= trace_from).then(|| {
let verdict = match &outcome {
Ok(c) => format!("converged in {} (residual {:.3e})", c.iterations, c.residual),
Err(rtx_fsi::FsiError::CouplingNotConverged {
iterations,
residual,
..
}) => format!(
"NOT converged after {iterations} (residual {residual:.3e}) -> {}",
if *residual < acceptable { "accepted stall" } else { "DEATH" }
),
Err(rtx_fsi::FsiError::CouplingDiverged {
iterations,
residual,
}) => format!(
"DIVERGED after {iterations} (residual {residual:.3e}) -> {}",
if *residual < acceptable { "accepted stall" } else { "DEATH" }
),
Err(e) => format!("error {e:?}"),
};
format!(
" TRACE step {step} t = {:.5}: increment {increment:.3e}, tol_step {tol_step:.3e}, \
retry_at {retry_at:.3e}, acceptable {acceptable:.3e}, retries so far {retried_steps}, \
stalls so far {stalled_steps}, Newton rescues {:?}; {verdict}",
t_release + (step + 1) as f64 * dt,
flag.borrow().rescue_counts(),
)
});
let outcome_ok = outcome.is_ok();
// The coupling-level rescue's ladder over this step's interval
// (never called with the knob off; the saved start state exists
// only with it on).
let run_ladder = || -> Result<super::rescue::RescueOutcome, rtx_fsi::FsiError> {
let iv = super::rescue::Interval {
harness: &harness,
solver: &solver,
field: &field,
flag: &flag,
wetted_dofs: &wetted_dofs,
fluid_saved: &fluid_saved,
field_saved: &field_saved,
start_state: step_start_state.as_ref().expect("rescue start state"),
start_nodal: step_start_nodal.as_deref().expect("rescue start load"),
config,
dt,
};
super::rescue::substep_interval(&iv)
};
let mut pending: Option<(&'static str, f64, super::rescue::RescueOutcome)> = None;
let mut tip_rejected = f64::NAN;
if let (Some(line), Err(_)) = (&trace_line, &outcome) {
// A stall's record prints here (accepted or fatal — a death
// panics below, before the commit); a converged step's
// record prints after the commit with its tip jump.
println!("{line}");
}
match outcome {
Ok(converged) => {
total_subiterations += converged.iterations;
@@ -607,26 +708,153 @@ pub fn run_march(case: BenchmarkCase, config: &MarchConfig) -> MarchResult {
total_subiterations += iterations;
max_subiterations = max_subiterations.max(iterations);
}
Err(e) => panic!(
"{} coupling failed at step {step}: {e:?} (Newton rescues so far: {:?})",
case.name,
flag.borrow().rescue_counts()
),
Err(e) => {
assert!(
coupling_rescue,
"{} coupling failed at step {step}: {e:?} (Newton rescues so far: {:?})",
case.name,
flag.borrow().rescue_counts()
);
let residual = match &e {
rtx_fsi::FsiError::CouplingNotConverged { residual, .. }
| rtx_fsi::FsiError::CouplingDiverged { residual, .. } => *residual,
_ => f64::NAN,
};
match run_ladder() {
Ok(o) => pending = Some(("fatal stall", residual, o)),
Err(e2) => panic!(
"{} coupling failed at step {step}: {e:?}; the coupling-level rescue's \
substep ladder {:?} failed too: {e2:?} (coupling rescues so far {}, \
Newton rescues {:?})",
case.name,
super::rescue::LADDER,
coupling_rescues,
flag.borrow().rescue_counts()
),
}
}
}
let prev_uy: Option<f64> = uy_series.last().copied();
if pending.is_none() {
// `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;
// Kinematic trigger: the coupling accepted this step, but the
// committed tip moved beyond the march's own running statistic.
if let (true, Some(p95), Some(prev)) = (coupling_rescue, running_p95, prev_uy) {
let uy_new = flag_state.displacement[a_dofs[1]];
let jump = (uy_new - prev).abs();
if jump > super::rescue::TIP_JUMP_FACTOR * p95 {
let accepted_solver = solver.borrow().snapshot();
let accepted_field = field.borrow().clone();
match run_ladder() {
Ok(o) => {
tip_rejected = uy_new;
pending = Some(("tip jump", jump, o));
}
Err(e2) => {
// Keep the step the coupling accepted; the
// ladder left the fluid mid-failure.
solver.borrow_mut().restore(&accepted_solver);
*field.borrow_mut() = accepted_field;
coupling_rescue_failures += 1;
println!(
" COUPLING RESCUE FAILED at step {step} t = {:.4}: tip jump \
{jump:.3e} > {} x running p95 {p95:.3e}; ladder {:?} ended \
{e2:?}; keeping the accepted step",
t_release + (step + 1) as f64 * dt,
super::rescue::TIP_JUMP_FACTOR,
super::rescue::LADDER
);
}
}
}
}
} else {
// A fatal stall never commits its last pass.
drop(latest.borrow_mut().take());
}
let was_rescued = pending.is_some();
if let Some((trigger, before, o)) = pending {
flag_state = o.state;
committed_nodal = o.nodal;
worst_conservation = worst_conservation.max(o.worst_conservation);
total_skipped += o.skipped;
total_subiterations += o.passes;
stalled_steps += o.stalls;
worst_stall = worst_stall.max(o.worst_residual);
// The main coupler's secant columns predate the rejected step.
if let Some(iqn_ref) = iqn.as_mut() {
iqn_ref.reset_history();
}
coupling_rescues += 1;
rescue_steps.push(step);
let record = super::rescue::RescueRecord {
step,
t: t_release + (step + 1) as f64 * dt,
trigger,
before,
n: o.n,
passes: o.passes,
tip_rejected,
tip_rescued: flag_state.displacement[a_dofs[1]],
};
println!(
" COUPLING RESCUE at step {step} t = {:.4}: {trigger} ({before:.3e}) carried by \
n = {} substeps in {} passes ({} substep stalls, worst {:.3e}); tip {:+.4e} -> \
{:+.4e}",
record.t,
o.n,
o.passes,
o.stalls,
o.worst_residual,
tip_rejected,
record.tip_rescued
);
rescue_records.push(record);
let per_second = (1.0 / dt).round() as usize;
let recent = rescue_steps
.iter()
.filter(|&&s| step - s < per_second)
.count();
assert!(
recent <= super::rescue::RATE_CAP_PER_SECOND,
"{} coupling rescue rate cap: {recent} rescues within one second of march at \
step {step} — a runaway the substeps only delay is not to be hidden",
case.name
);
}
// `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]];
// Running p95 of the committed tip jump over the trailing window
// of non-rescued steps (the kinematic trigger's yardstick).
if coupling_rescue && !was_rescued {
if let Some(prev) = prev_uy {
jump_window.push_back((uy - prev).abs());
if jump_window.len() > super::rescue::TIP_JUMP_WINDOW {
jump_window.pop_front();
}
if jump_window.len() == super::rescue::TIP_JUMP_WINDOW {
let mut sorted: Vec<f64> = jump_window.iter().copied().collect();
sorted.sort_by(|a, b| a.partial_cmp(b).unwrap());
running_p95 = Some(sorted[(0.95 * sorted.len() as f64) as usize]);
}
}
}
if let Some(line) = trace_line.as_ref().filter(|_| outcome_ok) {
let jump = uy - uy_series.last().copied().unwrap_or(uy);
println!("{line}; committed tip jump {jump:+.3e} (uy {uy:+.4e})");
}
times.push(t);
ux_series.push(ux);
uy_series.push(uy);
@@ -712,6 +940,9 @@ pub fn run_march(case: BenchmarkCase, config: &MarchConfig) -> MarchResult {
skipped: total_skipped,
spiked: harness.spiked_total.get(),
newton_rescues: flag.borrow().rescue_counts(),
coupling_rescues,
coupling_rescue_failures,
rescue_records,
final_state_finite: flag_state.displacement.iter().all(|v| v.is_finite()),
elapsed: start.elapsed().as_secs_f64(),
}
@@ -11,6 +11,7 @@
#![allow(dead_code)] // several test crates share this; each uses a subset
pub mod march;
pub mod rescue;
use std::cell::Cell;
use std::sync::{Arc, RwLock};
@@ -653,7 +654,33 @@ impl Fsi2Harness {
subcycle: usize,
v_n: Option<&[f64]>,
) {
let dt = self.dt_fluid * subcycle as f64;
self.advance_subcycled_with(
solver,
field,
d_n,
d_candidate,
subcycle,
v_n,
self.dt_fluid,
);
}
/// [`Self::advance_subcycled`] at an explicit fluid substep `dt_fluid`
/// (the coupled step is `dt_fluid * subcycle`). The march's own path
/// passes `self.dt_fluid` — same arithmetic, digit for digit; the
/// coupling-level rescue passes `dt_fluid / n` for its substeps.
#[allow(clippy::too_many_arguments)]
pub fn advance_subcycled_with(
&self,
solver: &mut EmbeddedPisoSolver,
field: &mut FlowField,
d_n: &[f64],
d_candidate: &[f64],
subcycle: usize,
v_n: Option<&[f64]>,
dt_fluid: f64,
) {
let dt = dt_fluid * subcycle as f64;
let mean_velocity: Vec<f64> = d_candidate
.iter()
.zip(d_n)
@@ -696,7 +723,7 @@ impl Fsi2Harness {
}
};
self.set_geometry(&d_sub, &ddot_sub);
futures::executor::block_on(solver.advance(field, self.dt_fluid)).unwrap();
futures::executor::block_on(solver.advance(field, dt_fluid)).unwrap();
}
}
}
@@ -0,0 +1,277 @@
//! Coupling-level rescue — rung A of omni-cortex
//! `docs/coupling_rescue_campaign.md`.
//!
//! A coupled step the march would otherwise abandon (a fatal stall:
//! the coupling ended above its acceptance window) or reject (a
//! committed tip jump beyond the march's own running statistic — the
//! runaway signature every recorded death shows for 537 steps before
//! its panic) is REPEATED over the same interval as n = 2, 4, 8, 16,
//! 32 coupled substeps of dt/n: Mayr, Wall and Gee's reduced-step
//! repetition (κ = 0.5, at most five repetitions) applied to the
//! partitioned loop, one layer above the structural Newton rescue.
//!
//! Each substep is a complete coupled step at dt/n: its own predictor
//! from the substep's start state, the fluid subcycled at
//! dt_fluid/n from the fluid state the previous substep committed, the
//! flag stepped at dt/n (`NonlinearDynamicStepper::step_with_dt`, whose
//! Newton and rescue ladder already take dt as a parameter), the C¹
//! interface velocity chained through the substep's own start
//! velocity, and a FRESH coupler (the secant history of the other dt
//! is not this map's). Tolerances follow the march's own rule at the
//! substep's own increment. A substep that ends above its acceptance
//! window fails the rung; the next rung restores the start-of-step
//! fluid snapshot and tries twice as many substeps.
//!
//! **Default OFF.** With the knob off nothing here runs — the march's
//! own step logic is not routed through this module, on purpose: the
//! path the digit-identity protocol certifies stays float-for-float
//! untouched, and this duplicate of it only exists past a trigger.
use std::cell::{Cell, RefCell};
use nalgebra::Vector3;
use rtx_cfd::solvers::incompressible::{EmbeddedPisoSolver, EmbeddedSolverState, FlowField};
use rtx_fea::analysis::{DynamicState, NonlinearDynamicStepper};
use rtx_fea::mesh::NodeId;
use rtx_fsi::{FsiError, IqnIls, Subiterated};
use super::Fsi2Harness;
use super::march::MarchConfig;
/// The substep ladder: five repetitions at κ = 0.5.
pub const LADDER: [usize; 5] = [2, 4, 8, 16, 32];
/// Kinematic trigger: a committed tip jump above this multiple of the
/// march's running p95 of |Δuy| rejects the step. Calibrated on all 74
/// recorded s = 1 marches (`scripts/fsi_burst_scan.py`): every march
/// that reached its horizon stays ≤ 2.4×; 27 of 28 deaths exceed 5×.
pub const TIP_JUMP_FACTOR: f64 = 3.0;
/// The running statistic's window (committed, non-rescued steps); the
/// kinematic trigger is inert until the window is full.
pub const TIP_JUMP_WINDOW: usize = 2000;
/// More rescues than this within one second of march is a runaway the
/// substeps only delay — the march dies loudly instead of hiding it.
pub const RATE_CAP_PER_SECOND: usize = 20;
/// One rescue, as recorded in the march result and printed as it fires.
#[derive(Debug, Clone)]
pub struct RescueRecord {
pub step: usize,
pub t: f64,
/// `"fatal stall"` or `"tip jump"`.
pub trigger: &'static str,
/// The rejected step's residual (fatal stall) or |Δuy| (tip jump).
pub before: f64,
/// Substeps that carried the interval; 0 = every rung failed (the
/// step was then kept as the coupling accepted it, or the march
/// died — see the trigger).
pub n: usize,
/// Coupling passes spent across all rungs tried.
pub passes: usize,
/// Tip uy the rejected step had committed (NaN for a fatal stall,
/// which never commits) and the tip the rescue committed.
pub tip_rejected: f64,
pub tip_rescued: f64,
}
/// What a successful ladder rung hands back to the march.
pub struct RescueOutcome {
pub state: DynamicState,
pub nodal: Vec<(NodeId, Vector3<f64>)>,
pub n: usize,
pub passes: usize,
pub stalls: usize,
pub worst_residual: f64,
pub worst_conservation: f64,
pub skipped: usize,
}
/// Everything one interval repetition needs, borrowed from the march.
pub struct Interval<'a, 'b> {
pub harness: &'a Fsi2Harness,
pub solver: &'a RefCell<EmbeddedPisoSolver>,
pub field: &'a RefCell<FlowField>,
pub flag: &'a RefCell<NonlinearDynamicStepper<'b>>,
pub wetted_dofs: &'a [[usize; 2]],
/// The start-of-step fluid snapshot and field the march already
/// saves for its coupling passes.
pub fluid_saved: &'a EmbeddedSolverState,
pub field_saved: &'a FlowField,
pub start_state: &'a DynamicState,
pub start_nodal: &'a [(NodeId, Vector3<f64>)],
pub config: &'a MarchConfig,
/// The full coupled step the interval spans.
pub dt: f64,
}
fn extract(state: &DynamicState, wetted: &[[usize; 2]]) -> Vec<f64> {
let mut d = vec![0.0; 2 * wetted.len()];
for (k, dofs) in wetted.iter().enumerate() {
d[2 * k] = state.displacement[dofs[0]];
d[2 * k + 1] = state.displacement[dofs[1]];
}
d
}
fn extract_velocity(state: &DynamicState, wetted: &[[usize; 2]]) -> Vec<f64> {
let mut v = vec![0.0; 2 * wetted.len()];
for (k, dofs) in wetted.iter().enumerate() {
v[2 * k] = state.velocity[dofs[0]];
v[2 * k + 1] = state.velocity[dofs[1]];
}
v
}
/// Repeat the interval up the ladder; the first rung whose every
/// substep converges or stalls inside its acceptance window wins.
///
/// # Errors
/// The last rung's coupling error when all five fail. The fluid and
/// field are then left as the last failed rung left them — the caller
/// restores or dies.
pub fn substep_interval(iv: &Interval<'_, '_>) -> Result<RescueOutcome, FsiError> {
let mut last_err = None;
let mut passes_so_far = 0usize;
for &n in &LADDER {
match attempt(iv, n) {
Ok(mut outcome) => {
outcome.passes += passes_so_far;
return Ok(outcome);
}
Err((e, passes)) => {
passes_so_far += passes;
last_err = Some(e);
}
}
}
Err(last_err.expect("ladder is non-empty"))
}
/// One rung: `n` coupled substeps of `dt/n` from the saved start.
fn attempt(iv: &Interval<'_, '_>, n: usize) -> Result<RescueOutcome, (FsiError, usize)> {
let cfg = iv.config;
let dt_sub = iv.dt / n as f64;
let dt_fluid_sub = dt_sub / cfg.subcycle as f64;
iv.solver.borrow_mut().restore(iv.fluid_saved);
*iv.field.borrow_mut() = iv.field_saved.clone();
let mut state = iv.start_state.clone();
let mut nodal: Vec<(NodeId, Vector3<f64>)> = iv.start_nodal.to_vec();
let mut passes = 0usize;
let mut stalls = 0usize;
let mut worst_residual = 0.0f64;
let mut worst_conservation = 0.0f64;
let mut skipped = 0usize;
for _ in 0..n {
let d_n = extract(&state, iv.wetted_dofs);
let v_n: Option<Vec<f64>> = cfg
.c1_interface
.then(|| extract_velocity(&state, iv.wetted_dofs));
let d_predicted: Vec<f64> = if cfg.predictor == "kinematic" {
let v = extract_velocity(&state, iv.wetted_dofs);
d_n.iter().zip(&v).map(|(d, v)| d + dt_sub * v).collect()
} else {
let mut flag = iv.flag.borrow_mut();
flag.set_nodal_forces(&nodal);
let (predicted, _) = flag
.step_with_dt(&state, dt_sub)
.expect("structure-alone predictor at the substep dt");
extract(&predicted, iv.wetted_dofs)
};
let sub_saved = iv.solver.borrow().snapshot();
let sub_field = iv.field.borrow().clone();
type PassResult = (
FlowField,
DynamicState,
Vec<(NodeId, Vector3<f64>)>,
f64,
usize,
);
let latest: RefCell<Option<PassResult>> = RefCell::new(None);
let pass_count = Cell::new(0usize);
let state_ref = &state;
let pass = |d_candidate: &[f64]| -> Vec<f64> {
pass_count.set(pass_count.get() + 1);
let mut solver_ref = iv.solver.borrow_mut();
solver_ref.restore(&sub_saved);
let mut trial = sub_field.clone();
iv.harness.advance_subcycled_with(
&mut solver_ref,
&mut trial,
&d_n,
d_candidate,
cfg.subcycle,
v_n.as_deref(),
dt_fluid_sub,
);
let (nodal_c, conservation, skipped_c) =
iv.harness.sample_load(&solver_ref, &trial, d_candidate);
let mut flag = iv.flag.borrow_mut();
flag.set_nodal_forces(&nodal_c);
let (candidate, _) = flag
.step_with_dt(state_ref, dt_sub)
.expect("flag step at the substep dt");
let d_new = extract(&candidate, iv.wetted_dofs);
*latest.borrow_mut() = Some((trial, candidate, nodal_c, conservation, skipped_c));
d_new
};
let increment: f64 = d_predicted
.iter()
.zip(&d_n)
.map(|(a, b)| (a - b) * (a - b))
.sum::<f64>()
.sqrt();
let tol_step = cfg.tol_floor.max(cfg.rtol * increment);
let acceptable = (cfg.stall_accept * tol_step).max(0.1 * increment);
let outcome = if cfg.coupler == "iqn" {
let mut iqn = IqnIls::new(cfg.max_subiterations, tol_step)
.map_err(|e| (e, passes))?
.with_reuse(cfg.reuse)
.with_initial_relaxation(cfg.initial_relaxation)
.map_err(|e| (e, passes))?;
iqn.solve(&d_predicted, pass)
} else {
Subiterated::aitken(cfg.max_subiterations, tol_step)
.map_err(|e| (e, passes))?
.solve(&d_predicted, pass)
};
passes += pass_count.get();
match outcome {
Ok(c) => worst_residual = worst_residual.max(c.residual),
Err(
FsiError::CouplingNotConverged { residual, .. }
| FsiError::CouplingDiverged { residual, .. },
) if residual < acceptable => {
stalls += 1;
worst_residual = worst_residual.max(residual);
}
Err(e) => return Err((e, passes)),
}
let (new_field, new_state, new_nodal, conservation, skipped_c) =
latest.borrow_mut().take().expect("pass ran");
// As in the march: the solver's state after the last pass IS
// the response to the accepted interface; commit it directly.
*iv.field.borrow_mut() = new_field;
state = new_state;
nodal = new_nodal;
worst_conservation = worst_conservation.max(conservation);
skipped += skipped_c;
}
Ok(RescueOutcome {
state,
nodal,
n,
passes,
stalls,
worst_residual,
worst_conservation,
skipped,
})
}