From 74fac32a91e1ae1cf94276876d769c0bc555400a Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Fri, 25 Sep 2026 22:11:31 -0500 Subject: [PATCH] R8-e: the coupled 3D FSI's small-dt instability located in the load/kinematics transfer, not a lagged fluid response MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - tests/fsi2_embedded3_impulse.rs (new, ignored, cuda): the added-mass timing instrument — the R8-a fluid with a prescribed centreline velocity step (a one-step acceleration delta), baseline and impulse passes from one device snapshot, the load difference per step. Pitch shape on the slab ny 62: generalized added mass 39.11 / 39.09 / 39.07 kg/m at dt x 1 / 0.75 / 0.5, the step after carries 1.7 / 1.3 / 0.9 % above the steady tail — the cut-cell fluid answers in the same step (the lagged-added-mass hypothesis refuted). - fsi2_embedded3.rs: `RTX_E3FSI_LINE_CSV` (per-step centreline + mid/bottom/top y per station) showed the growing 2dt mode is a through-thickness mode — the edge nodes move against the mid-plane nodes (~4:1): the fluid reads the mid-plane, the loads land on the edges, so the added mass acts on that mode with the wrong sign (a negative effective added mass; the scalar Newmark model is unstable for it at omega dt below ~1-2 with gamma 0.7, stable at larger dt or gamma 0.9 — R8-a's observed pattern). `RTX_E3FSI_KIN=surface` (default off, byte-identical when unset): the fluid's centreline = the mean of the bottom and top edge nodes (the nodes the loads land on) — dt x 0.75 / 0.5 / 0.25 stable at gamma 0.7. Co-Authored-By: Claude Opus 5.5 (1M context) --- .../rtx-fsi/tests/fsi2_embedded3.rs | 49 +++- .../rtx-fsi/tests/fsi2_embedded3_impulse.rs | 220 ++++++++++++++++++ 2 files changed, 268 insertions(+), 1 deletion(-) create mode 100644 crates/specialized/rtx-fsi/tests/fsi2_embedded3_impulse.rs diff --git a/crates/specialized/rtx-fsi/tests/fsi2_embedded3.rs b/crates/specialized/rtx-fsi/tests/fsi2_embedded3.rs index 99cf106..33b40a0 100644 --- a/crates/specialized/rtx-fsi/tests/fsi2_embedded3.rs +++ b/crates/specialized/rtx-fsi/tests/fsi2_embedded3.rs @@ -18,7 +18,11 @@ //! `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). +//! 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= RTX_CUDA_ARCH=sm_120 cargo test --release -p rtx-fsi \ //! --features cuda --test fsi2_embedded3 -- --ignored --nocapture` @@ -67,6 +71,9 @@ struct Flag { /// (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 { @@ -107,6 +114,12 @@ impl Flag { top, tip, a_node: find(0.6, 0.2), + surface: (0..STATIONS) + .map(|k| { + let x = X0 + 0.01 * k as f64; + [find(x, 0.19), find(x, 0.21)] + }) + .collect(), } } } @@ -287,6 +300,11 @@ fn fsi2_on_embedded3() { .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). @@ -380,8 +398,27 @@ fn fsi2_on_embedded3() { }) .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 { 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]]; @@ -603,6 +640,16 @@ fn fsi2_on_embedded3() { ) .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 = 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() diff --git a/crates/specialized/rtx-fsi/tests/fsi2_embedded3_impulse.rs b/crates/specialized/rtx-fsi/tests/fsi2_embedded3_impulse.rs new file mode 100644 index 0000000..843cb15 --- /dev/null +++ b/crates/specialized/rtx-fsi/tests/fsi2_embedded3_impulse.rs @@ -0,0 +1,220 @@ +//! R8-e: the added-mass timing instrument for the coupled 3D FSI's small-dt +//! instability. The R8-a fluid (the embedded3 cut-cell device path, the +//! slab by default) with the flag's centreline PRESCRIBED: rigid to `T0`, +//! then from one snapshot two passes of `N` steps each — the baseline (the +//! flag at rest) and the impulse (the centreline's velocity steps from 0 to +//! `V` × shape at step `N_PRE`: its wall acceleration is a single-step +//! delta V/dt; `MODE=pulse` returns the velocity to 0 one step later, a +//! +V/dt, −V/dt pair). The difference of the two passes' loads, step by +//! step, is the fluid's response to the impulse: an added-mass response +//! lives in the impulse's own step; a lagged part shows one step later. +//! +//! Knobs `RTX_E3IMP_*`: `NY` (62), `NZ` (4), `T0` (0.2 s), `N_PRE` (4), +//! `N` (24), `V` (0.05 m/s at the tip), `SHAPE` (`pitch` about the root, +//! 1 at A; `heave`; `zig`, `tip`), `MODE` (`step`, `pulse`), `CSV` (per-step series). +//! The fluid's own knobs as the R8-a harness (`RTX_E3FSI_DT_SCALE`, …). +#![cfg(feature = "cuda")] + +#[path = "fsi2_embedded3/fluid.rs"] +mod fluid; +#[path = "fsi2_embedded3/state.rs"] +mod state; + +use std::io::Write as _; + +use fluid::{CX, CY, Contribution, E3Fluid, HALF, Line, R_CYL}; + +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; +const STATIONS: usize = 35; +const CHORD: f64 = 0.35; + +fn line(t: f64, c: &[f64], v: &[f64]) -> Line { + Line { + t, + pts: (0..STATIONS) + .map(|k| [X0 + 0.01 * k as f64, 0.2 + c[k]]) + .collect(), + vel: (0..STATIONS).map(|k| [0.0, v[k]]).collect(), + } +} + +/// The flag's loads from the route's contributions (per unit span): fy on +/// the flag, fy by part, and the generalized force Σ fy φ(x) with the +/// prescribed shape's φ. +fn flag_loads( + contrib: &[Contribution], + pts: &[[f64; 2]], + width: f64, + phi: impl Fn(f64) -> f64, +) -> (f64, f64, [f64; 4]) { + let (mut fy, mut q, mut parts) = (0.0, 0.0, [0.0; 4]); + for &(pos, c, part, v) in contrib { + if c != 1 { + continue; + } + let (x, y) = (pos[0], pos[1]); + let mut best = f64::INFINITY; + 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 u = (((x - a[0]) * ex + (y - a[1]) * ey) / (ex * ex + ey * ey)).clamp(0.0, 1.0); + best = best.min(((x - a[0] - u * ex).powi(2) + (y - a[1] - u * ey).powi(2)).sqrt()); + } + let d_cyl = ((x - CX).powi(2) + (y - CY).powi(2)).sqrt() - R_CYL; + if d_cyl < best - HALF { + continue; + } + let v = v / width; + fy += v; + q += v * phi(x); + parts[part] += v * phi(x); + } + (fy, q, parts) +} + +#[test] +#[ignore = "R8-e: the added-mass timing instrument (GPU, a minute)"] +fn impulse_response() { + let ny = env_f("RTX_E3IMP_NY", 62.0) as usize; + let nz = env_f("RTX_E3IMP_NZ", 4.0) as usize; + let t0 = env_f("RTX_E3IMP_T0", 0.2); + let n_pre = env_f("RTX_E3IMP_N_PRE", 4.0) as usize; + let n = env_f("RTX_E3IMP_N", 24.0) as usize; + let vel = env_f("RTX_E3IMP_V", 0.05); + let shape = std::env::var("RTX_E3IMP_SHAPE").unwrap_or_else(|_| "pitch".into()); + let mode = std::env::var("RTX_E3IMP_MODE").unwrap_or_else(|_| "step".into()); + // The shape per station (`zig`: alternate stations ±, the root clamped; + // `tip`: the last 5 stations), linear between stations. + let phis: Vec = (0..STATIONS) + .map(|k| { + let x = X0 + 0.01 * k as f64; + match shape.as_str() { + "heave" => 1.0, + "zig" => { + if k == 0 { + 0.0 + } else if k % 2 == 0 { + 1.0 + } else { + -1.0 + } + } + "tip" => { + if k + 5 >= STATIONS { + 1.0 + } else { + 0.0 + } + } + _ => ((x - X0) / CHORD).clamp(0.0, 1.0), + } + }) + .collect(); + let phi = |x: f64| -> f64 { + let s = ((x - X0) / 0.01).clamp(0.0, (STATIONS - 1) as f64); + let m = (s.floor() as usize).min(STATIONS - 2); + let u = s - m as f64; + (1.0 - u) * phis[m] + u * phis[m + 1] + }; + let zero = vec![0.0; STATIONS]; + let rest = line(0.0, &zero, &zero); + let mut fl = E3Fluid::build(ny, nz, 3.0, rest.clone(), None); + let dt = fl.dt; + let rigid = (t0 / dt).round() as usize; + for s in 0..rigid { + let r = fl.step(); + assert!(r.final_residual.is_finite(), "rigid death at {s}"); + } + let t_start = fl.time(); + println!( + " R8-e impulse: ny {ny} nz {nz} dt {dt:.4e}, rigid {rigid} steps to t {t_start:.4}; shape {shape}, mode {mode}, V {vel} m/s (acceleration V/dt {:.2} m/s²) at step {n_pre} of {n}", + vel / dt + ); + let snap = fl.snapshot(); + // Per pass: per step (t, tip velocity, fy, Q, Q parts, lift total). + let mut rows: Vec> = Vec::new(); + for pass in 0..2 { + fl.restore(&snap); + let mut c = zero.clone(); + let mut prev = rest.clone(); + prev.t = t_start; + let mut out = Vec::new(); + for s in 0..n { + let t_new = t_start + (s + 1) as f64 * dt; + // The step's centreline velocity (constant over the step). + let v_tip = if pass == 0 || s < n_pre { + 0.0 + } else if mode == "pulse" { + if s == n_pre { vel } else { 0.0 } + } else { + vel + }; + let v: Vec = phis.iter().map(|p| p * v_tip).collect(); + for k in 0..STATIONS { + c[k] += dt * v[k]; + } + let next = line(t_new, &c, &v); + fl.set_lines(prev.clone(), next.clone()); + let r = fl.step(); + assert!( + r.final_residual.is_finite(), + "death at pass {pass} step {s}" + ); + let (tot, contrib) = fl.loads(); + let (fy, q, parts) = flag_loads(&contrib, &next.pts, fl.load_width, &phi); + out.push([ + t_new, v_tip, fy, q, parts[0], parts[1], parts[2], parts[3], tot[1], + ]); + prev = next; + } + rows.push(out); + } + let mut csv = std::env::var("RTX_E3IMP_CSV") + .ok() + .map(|p| std::fs::File::create(p).expect("csv")); + if let Some(f) = csv.as_mut() { + writeln!( + f, + "k,t,v_tip,fy0,q0,fy,q,dfy,dq,dq_p,dq_shear,dq_xdiff,dq_xconv,dlift" + ) + .unwrap(); + } + let q_imp = rows[1][n_pre][3] - rows[0][n_pre][3]; + for s in 0..n { + let (a, b) = (&rows[0][s], &rows[1][s]); + let k = s as i64 - n_pre as i64; + let d = |i: usize| b[i] - a[i]; + if let Some(f) = csv.as_mut() { + writeln!( + f, + "{k},{:.6},{:.4e},{:.6},{:.6},{:.6},{:.6},{:.6e},{:.6e},{:.6e},{:.6e},{:.6e},{:.6e},{:.6e}", + b[0], b[1], a[2], a[3], b[2], b[3], d(2), d(3), d(4), d(5), d(6), d(7), d(8) + ) + .unwrap(); + } + if (-1..=6).contains(&k) { + println!( + " k {k:+}: v_tip {:.3e} ΔQ {:+.5e} ({:+.4} of the impulse step) parts p {:+.4e} shear {:+.4e} xdiff {:+.4e} xconv {:+.4e} Δfy {:+.5e}", + b[1], + d(3), + d(3) / q_imp, + d(4), + d(5), + d(6), + d(7), + d(2) + ); + } + } + println!( + " IMPULSE ΔQ(k=0) {q_imp:+.5e} N/m → generalized added mass −ΔQ·dt/V {:.4} kg/m", + -q_imp * dt / vel + ); +}