S2-9: the flag on the recorded FSI2 kinematics — fsi2_overset_dump_centreline (rtx-fsi) exports the saved instants' wetted-node d/ḋ as a centreline per knot (clamp + 70 columns); tests/embedded3_flag_kinematics (cubic Hermite through the knots with their velocities, t = 0 at RTX_E3_FLAG_KIN_T0, a 0.5 s ramp from the undeflected line, arc-length interpolation along the stations) drives both the 3D/slab flag (RTX_E3_FLAG_KINEMATICS, both in-plane surface-velocity components, the record's period) and the 2D reference; the analytic paths unchanged (same tuples)

Co-Authored-By: Claude Fable 5.1 <[email protected]>

crates/specialized/rtx-cfd/tests/embedded3_flag_kinematics/mod.rs

@@ -0,0 +1,113 @@
+//! S2-9: the recorded FSI2 kinematics for the flag tests (`RTX_E3_FLAG_KINEMATICS=<csv>`
+//! from `fsi2_overset_dump_centreline`): per knot t and per station (x, y, vx, vy);
+//! between knots the cubic Hermite interpolant through the knots with their
+//! velocities (the overset replay's form). The run's t = 0 maps to the recorded
+//! `RTX_E3_FLAG_KIN_T0` (13.0 s, inside the limit cycle) and the motion ramps in
+//! over `RTX_E3_FLAG_KIN_RAMP` (0.5 s) from the undeflected line at `cy`.
+#![allow(dead_code)]
+
+fn env_f(name: &str, default: f64) -> f64 {
+    std::env::var(name).ok().and_then(|v| v.parse().ok()).unwrap_or(default)
+}
+
+pub struct Recorded {
+    pub times: Vec<f64>,
+    /// Per knot: per station (x, y, vx, vy).
+    pub knots: Vec<Vec<(f64, f64, f64, f64)>>,
+    pub t0: f64,
+    pub ramp: f64,
+}
+
+pub fn recorded() -> Option<&'static Recorded> {
+    static REC: std::sync::OnceLock<Option<Recorded>> = std::sync::OnceLock::new();
+    REC.get_or_init(|| {
+        let path = std::env::var("RTX_E3_FLAG_KINEMATICS").ok()?;
+        let text = std::fs::read_to_string(&path).expect("kinematics csv");
+        let mut times = Vec::new();
+        let mut knots = Vec::new();
+        for line in text.lines().skip(1) {
+            let v: Vec<f64> = line.split(',').map(|x| x.trim().parse().expect("number")).collect();
+            times.push(v[0]);
+            knots.push(v[1..].chunks_exact(4).map(|c| (c[0], c[1], c[2], c[3])).collect());
+        }
+        assert!(times.len() >= 2, "kinematics: at least two knots");
+        Some(Recorded {
+            times,
+            knots,
+            t0: env_f("RTX_E3_FLAG_KIN_T0", 13.0),
+            ramp: env_f("RTX_E3_FLAG_KIN_RAMP", 0.5),
+        })
+    })
+    .as_ref()
+}
+
+impl Recorded {
+    /// The largest recorded station speed (the CFL and band bound).
+    pub fn max_speed(&self) -> f64 {
+        self.knots
+            .iter()
+            .flat_map(|k| k.iter().map(|p| (p.2 * p.2 + p.3 * p.3).sqrt()))
+            .fold(0.0f64, f64::max)
+    }
+
+    /// The stations at the run's time `t`: (x, y, vx, vy), ramped from the
+    /// undeflected line at `cy` (each station at its first-knot x).
+    pub fn at(&self, t: f64, cy: f64) -> Vec<(f64, f64, f64, f64)> {
+        let tr = (self.t0 + t).min(*self.times.last().unwrap());
+        let k = match self.times.binary_search_by(|x| x.partial_cmp(&tr).unwrap()) {
+            Ok(i) => i.min(self.times.len() - 2),
+            Err(i) => i.saturating_sub(1).min(self.times.len() - 2),
+        };
+        let h = self.times[k + 1] - self.times[k];
+        let s = ((tr - self.times[k]) / h).clamp(0.0, 1.0);
+        let (s2, s3) = (s * s, s * s * s);
+        let (h00, h10, h01, h11) = (2.0 * s3 - 3.0 * s2 + 1.0, s3 - 2.0 * s2 + s, -2.0 * s3 + 3.0 * s2, s3 - s2);
+        let (d00, d10, d01, d11) = (
+            (6.0 * s2 - 6.0 * s) / h,
+            (3.0 * s2 - 4.0 * s + 1.0) / h,
+            (-6.0 * s2 + 6.0 * s) / h,
+            (3.0 * s2 - 2.0 * s) / h,
+        );
+        let (r, rd) = if t >= self.ramp {
+            (1.0, 0.0)
+        } else {
+            let u = t / self.ramp;
+            (u * u * (3.0 - 2.0 * u), 6.0 * u * (1.0 - u) / self.ramp)
+        };
+        let (a, b) = (&self.knots[k], &self.knots[k + 1]);
+        let first = &self.knots[0];
+        a.iter()
+            .zip(b)
+            .zip(first)
+            .map(|((p, q), f0)| {
+                let x = h00 * p.0 + h10 * h * p.2 + h01 * q.0 + h11 * h * q.2;
+                let y = h00 * p.1 + h10 * h * p.3 + h01 * q.1 + h11 * h * q.3;
+                let vx = d00 * p.0 + d10 * h * p.2 + d01 * q.0 + d11 * h * q.2;
+                let vy = d00 * p.1 + d10 * h * p.3 + d01 * q.1 + d11 * h * q.3;
+                let (xr, yr) = (f0.0, cy);
+                (
+                    (1.0 - r) * xr + r * x,
+                    (1.0 - r) * yr + r * y,
+                    r * vx + rd * (x - xr),
+                    r * vy + rd * (y - yr),
+                )
+            })
+            .collect()
+    }
+
+    /// The point and velocity at arc fraction `s ∈ [0, 1]` along the
+    /// stations' polyline at `t` (linear between stations by arc length).
+    pub fn along(&self, pts: &[(f64, f64, f64, f64)], s: f64) -> (f64, f64, f64, f64) {
+        let mut cum = vec![0.0; pts.len()];
+        for m in 1..pts.len() {
+            let (dx, dy) = (pts[m].0 - pts[m - 1].0, pts[m].1 - pts[m - 1].1);
+            cum[m] = cum[m - 1] + (dx * dx + dy * dy).sqrt();
+        }
+        let target = s.clamp(0.0, 1.0) * cum[pts.len() - 1];
+        let m = cum.partition_point(|&c| c < target).clamp(1, pts.len() - 1);
+        let seg = cum[m] - cum[m - 1];
+        let u = if seg > 0.0 { ((target - cum[m - 1]) / seg).clamp(0.0, 1.0) } else { 0.0 };
+        let (a, b) = (pts[m - 1], pts[m]);
+        (a.0 + u * (b.0 - a.0), a.1 + u * (b.1 - a.1), a.2 + u * (b.2 - a.2), a.3 + u * (b.3 - a.3))
+    }
+}

crates/specialized/rtx-cfd/tests/embedded3_flag_reference_2d.rs

@@ -8,6 +8,8 @@
 //! median-filtered lift swing.
 //!
 //! `cargo test --release -p rtx-cfd --test embedded3_flag_reference_2d -- --ignored --nocapture`
+mod embedded3_flag_kinematics;
+use embedded3_flag_kinematics::{Recorded, recorded};
 use rtx_cfd::solvers::incompressible::{
     AleBoundaries, ConvectionScheme, EmbeddedBody, EmbeddedParameters, EmbeddedPisoSolver,
     FlowField, MgPrecision, MgSmoother, PoissonSolverKind, SideBoundary,
@@ -48,26 +50,42 @@ fn deflection(s: f64, t: f64) -> (f64, f64) {
     )
 }
 
-fn centreline(m: usize, n: usize, t: f64) -> (f64, f64, f64) {
+/// The centreline point `m` of `n` at `t`: (x, y, vx, vy) — the first
+/// mode, or (S2-9) the recorded FSI2 kinematics along its stations.
+fn centreline(m: usize, n: usize, t: f64) -> (f64, f64, f64, f64) {
     let s = m as f64 / n as f64;
+    if let Some(rec) = recorded() {
+        thread_local! {
+            static PTS: std::cell::RefCell<(f64, Vec<(f64, f64, f64, f64)>)> =
+                const { std::cell::RefCell::new((f64::NAN, Vec::new())) };
+        }
+        return PTS.with(|cell| {
+            let mut c = cell.borrow_mut();
+            if c.0.to_bits() != t.to_bits() {
+                c.1 = rec.at(t, CY);
+                c.0 = t;
+            }
+            rec.along(&c.1, s)
+        });
+    }
     let (d, v) = deflection(s, t);
-    (FLAG_X0 + s * FLAG_LEN, CY + d, v)
+    (FLAG_X0 + s * FLAG_LEN, CY + d, 0.0, v)
 }
 
 /// Distance to the capsule flag and the centreline velocity at the foot.
-fn flag_sdf(x: f64, y: f64, t: f64) -> (f64, f64) {
+fn flag_sdf(x: f64, y: f64, t: f64) -> (f64, (f64, f64)) {
     let n = 40;
     let mut best = f64::INFINITY;
-    let mut v_best = 0.0;
+    let mut v_best = (0.0, 0.0);
     for m in 0..n {
-        let (ax, ay, av) = centreline(m, n, t);
-        let (bx, by, bv) = centreline(m + 1, n, t);
+        let (ax, ay, avx, avy) = centreline(m, n, t);
+        let (bx, by, bvx, bvy) = centreline(m + 1, n, t);
         let (ex, ey) = (bx - ax, by - ay);
         let u = (((x - ax) * ex + (y - ay) * ey) / (ex * ex + ey * ey)).clamp(0.0, 1.0);
         let d = ((x - ax - u * ex).powi(2) + (y - ay - u * ey).powi(2)).sqrt();
         if d < best {
             best = d;
-            v_best = av + u * (bv - av);
+            v_best = (avx + u * (bvx - avx), avy + u * (bvy - avy));
         }
     }
     (best - FLAG_HALF, v_best)
@@ -84,8 +102,8 @@ fn samples(t: f64, ds: f64) -> Vec<(f64, f64, f64, f64, f64)> {
     let mut out = Vec::new();
     let n = ((FLAG_LEN / ds).ceil() as usize).max(8);
     for m in 0..n {
-        let (ax, ay, _) = centreline(m, n, t);
-        let (bx, by, _) = centreline(m + 1, n, t);
+        let (ax, ay, _, _) = centreline(m, n, t);
+        let (bx, by, _, _) = centreline(m + 1, n, t);
         let (ex, ey) = (bx - ax, by - ay);
         let len = (ex * ex + ey * ey).sqrt();
         let (tx, ty) = (ex / len, ey / len);
@@ -99,8 +117,8 @@ fn samples(t: f64, ds: f64) -> Vec<(f64, f64, f64, f64, f64)> {
         }
     }
     // The tip: a semicircle around the last centreline point.
-    let (tx0, ty0, _) = centreline(n, n, t);
-    let (px, py, _) = centreline(n - 1, n, t);
+    let (tx0, ty0, _, _) = centreline(n, n, t);
+    let (px, py, _, _) = centreline(n - 1, n, t);
     let ang0 = (ty0 - py).atan2(tx0 - px);
     let n_arc = ((std::f64::consts::PI * FLAG_HALF / ds).ceil() as usize).max(4);
     for k in 0..n_arc {
@@ -144,8 +162,10 @@ async fn flag_wake_2d_reference() -> CfdResult<()> {
     let h = H / ny as f64;
     let nx = (L / h).round() as usize;
     let dt = 8.817e-4;
-    let period = 1.0 / FREQ;
-    let t_end = 2.0 * period;
+    // S2-9: the record's period and the requested number of periods.
+    let env_f = |k: &str, d: f64| std::env::var(k).ok().and_then(|v| v.parse().ok()).unwrap_or(d);
+    let period = if recorded().is_some() { env_f("RTX_E3_FLAG_KIN_PERIOD", 0.5225) } else { 1.0 / FREQ };
+    let t_end = env_f("RTX_E3_FLAG_PERIODS", 2.0) * period;
     let config = CfdConfig::new()
         .with_density(RHO)
         .with_viscosity(RHO * NU)
@@ -174,7 +194,7 @@ async fn flag_wake_2d_reference() -> CfdResult<()> {
         }
     });
     let flag = EmbeddedBody::from_sdf(|x, y, t| flag_sdf(x, y, t).0)
-        .with_surface_velocity(|x, y, t| (0.0, flag_sdf(x, y, t).1));
+        .with_surface_velocity(|x, y, t| flag_sdf(x, y, t).1);
     solver.set_moving_body(EmbeddedBody::union(
         EmbeddedBody::circle(CX, CY, R_CYL),
         flag,

crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs

@@ -19,6 +19,8 @@
 //!  RTX_CUDA_ARCH=sm_120 cargo test --release -p rtx-cfd --features cuda --test embedded3_flag_wake -- --ignored --nocapture`
 #![cfg(feature = "cuda")]
 
+mod embedded3_flag_kinematics;
+use embedded3_flag_kinematics::{Recorded, recorded};
 use rtx_cfd::solvers::incompressible::ConvectionScheme;
 use rtx_cfd::solvers::incompressible::embedded3::step::device::DeviceStep;
 use rtx_cfd::solvers::incompressible::embedded3::{
@@ -80,8 +82,48 @@ fn amplitude() -> f64 {
 
 /// Signed distance to the deflected flag's cross-section (a capsule
 /// around the centreline polyline of `n` segments) and the centreline's
-/// transverse velocity at the closest point.
-fn flag_2d(x: f64, y: f64, t: f64) -> (f64, f64) {
+/// velocity at the closest point (transverse only in the analytic mode).
+fn flag_2d(x: f64, y: f64, t: f64) -> (f64, (f64, f64)) {
+    if let Some(rec) = recorded() {
+        return flag_2d_recorded(rec, x, y, t);
+    }
+    let (d, v) = flag_2d_analytic(x, y, t);
+    (d, (0.0, v))
+}
+
+/// The recorded centreline's capsule and velocity at the closest point.
+fn flag_2d_recorded(rec: &Recorded, x: f64, y: f64, t: f64) -> (f64, (f64, f64)) {
+    thread_local! {
+        static POLY: std::cell::RefCell<(f64, Vec<(f64, f64, f64, f64)>)> =
+            const { std::cell::RefCell::new((f64::NAN, Vec::new())) };
+    }
+    POLY.with(|cell| {
+        let mut c = cell.borrow_mut();
+        if c.0.to_bits() != t.to_bits() {
+            c.1 = rec.at(t, CY);
+            c.0 = t;
+        }
+        let pts = &c.1;
+        let mut best = f64::INFINITY;
+        let mut v_best = (0.0, 0.0);
+        for m in 0..pts.len() - 1 {
+            let (ax, ay, avx, avy) = pts[m];
+            let (bx, by, bvx, bvy) = pts[m + 1];
+            let (ex, ey) = (bx - ax, by - ay);
+            let l2 = ex * ex + ey * ey;
+            let u = (((x - ax) * ex + (y - ay) * ey) / l2).clamp(0.0, 1.0);
+            let (px, py) = (ax + u * ex, ay + u * ey);
+            let d = ((x - px).powi(2) + (y - py).powi(2)).sqrt();
+            if d < best {
+                best = d;
+                v_best = (avx + u * (bvx - avx), avy + u * (bvy - avy));
+            }
+        }
+        (best - FLAG_HALF, v_best)
+    })
+}
+
+fn flag_2d_analytic(x: f64, y: f64, t: f64) -> (f64, f64) {
     const N: usize = 40;
     // The centreline polyline at `t`, once per thread and time (PERF-3
     // P1-2): the solver asks for the surface velocity at ~10⁶ faces per
@@ -130,7 +172,7 @@ fn flag_span() -> f64 {
 
 /// The flag in 3D: the extruded capsule cut to the span with edges
 /// rounded to radius `r` (no cut at the full width).
-fn flag_3d(x: f64, y: f64, z: f64, t: f64, r: f64) -> (f64, f64) {
+fn flag_3d(x: f64, y: f64, z: f64, t: f64, r: f64) -> (f64, (f64, f64)) {
     let (d2, v) = flag_2d(x, y, t);
     let span = flag_span();
     if span >= H {
@@ -159,10 +201,13 @@ fn flag_wake_on_the_device() {
     let slab_nz = env_f("RTX_E3_FLAG_NZ", 0.0) as usize;
     let nz = if slab_nz > 0 { slab_nz } else { ny };
     let r_edge = h;
-    let dt_cfl = 0.3 * h / (U_M.max(2.0 * std::f64::consts::PI * FREQ * AMP));
+    // S2-9: with a recorded kinematics the period and the speed bound are the record's.
+    let rec_period = env_f("RTX_E3_FLAG_KIN_PERIOD", 0.5225);
+    let rec_speed = recorded().map(Recorded::max_speed);
+    let dt_cfl = 0.3 * h / (U_M.max(rec_speed.unwrap_or(2.0 * std::f64::consts::PI * FREQ * AMP)));
     // `RTX_E3_FLAG_DT_SCALE` scales the step (the dt ladder of the loads).
     let dt = dt_cfl.min(0.5 * h * h / (6.0 * NU)) * env_f("RTX_E3_FLAG_DT_SCALE", 1.0);
-    let period = 1.0 / FREQ;
+    let period = if recorded().is_some() { rec_period } else { 1.0 / FREQ };
     let t_end = periods * period;
     let mut solver = Solver::new(
         Fluid {
@@ -196,7 +241,7 @@ fn flag_wake_on_the_device() {
             },
             // The narrow band: the flag's tip speed bounds the surface motion.
             max_surface_speed: Some(
-                (2.0 * std::f64::consts::PI * FREQ * amplitude() * 1.05).max(1e-3),
+                (rec_speed.unwrap_or(2.0 * std::f64::consts::PI * FREQ * amplitude()) * 1.05).max(1e-3),
             ),
             ..Parameters::default()
         },
@@ -218,9 +263,9 @@ fn flag_wake_on_the_device() {
     let cyl = move |x: f64, y: f64| ((x - CX).powi(2) + (y - CY).powi(2)).sqrt() - R_CYL;
     let body = Body::from_sdf(move |x, y, z, t| cyl(x, y).min(flag_3d(x, y, z, t, r_edge).0))
         .with_surface_velocity(move |x, y, z, t| {
-            let (df, v) = flag_3d(x, y, z, t, r_edge);
+            let (df, (vx, vy)) = flag_3d(x, y, z, t, r_edge);
             if df <= cyl(x, y) {
-                (0.0, v, 0.0)
+                (vx, vy, 0.0)
             } else {
                 (0.0, 0.0, 0.0)
             }

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

@@ -649,3 +649,71 @@ fn fsi2_overset_prescribed_motion() {
         mean(&|w| w.3[2])
     );
 }
+
+/// S2-9: the recorded FSI2 kinematics as a CENTRELINE per instant for the
+/// 3D flag (`RTX_FSI2O_REPLAY=dir`, `RTX_FSI2O_NY`, `RTX_FSI2O_CENTRELINE_CSV=out`):
+/// per knot, t and per station (the clamp, then each bottom/top column
+/// root → tip) the deformed midpoint (x, y) and its velocity (vx, vy).
+#[test]
+#[ignore = "S2-9 export: needs RTX_FSI2O_REPLAY (the saved instants)"]
+fn fsi2_overset_dump_centreline() {
+    use fsi2_harness::overset::OversetFluid;
+    use fsi2_harness::replay::Replay;
+    use std::io::Write as _;
+    let Ok(replay_dir) = std::env::var("RTX_FSI2O_REPLAY") else {
+        println!("  RTX_FSI2O_REPLAY unset");
+        return;
+    };
+    let out = std::env::var("RTX_FSI2O_CENTRELINE_CSV").expect("RTX_FSI2O_CENTRELINE_CSV");
+    let ny: usize = std::env::var("RTX_FSI2O_NY").ok().and_then(|v| v.parse().ok()).unwrap_or(62);
+    let case = case_from_env("FSI2O", FSI2);
+    let fluid = OversetFluid::build_case(case, ny, 35, 100, 3).expect("overset fluid");
+    let itf = &fluid.interface;
+    let replay = Replay::load(&replay_dir);
+    let n = 2 * itf.wetted.len();
+    assert_eq!(replay.d[0].len(), n, "the instants' d does not match this interface");
+    // Columns root → tip: bottom is root → tip, top is tip → root.
+    let cols: Vec<(usize, usize)> = itf
+        .bottom
+        .iter()
+        .copied()
+        .zip(itf.top.iter().rev().copied())
+        .collect();
+    for &(b, t) in &cols {
+        assert!((itf.reference[b].0 - itf.reference[t].0).abs() < 1e-9, "bottom / top columns disagree");
+    }
+    let mut f = std::fs::File::create(&out).expect("csv");
+    write!(f, "t").unwrap();
+    for i in 0..=cols.len() {
+        write!(f, ",x{i},y{i},vx{i},vy{i}").unwrap();
+    }
+    writeln!(f).unwrap();
+    let clamp = (fsi2_harness::FLAG_X0, 0.5 * (fsi2_harness::FLAG_Y0 + fsi2_harness::FLAG_Y1));
+    for (k, &tk) in replay.times.iter().enumerate() {
+        let (d, dd) = (&replay.d[k], &replay.dd[k]);
+        write!(f, "{tk:.6}").unwrap();
+        write!(f, ",{:.9},{:.9},0,0", clamp.0, clamp.1).unwrap();
+        for &(b, t) in &cols {
+            let (xb, yb) = (itf.reference[b].0 + d[2 * b], itf.reference[b].1 + d[2 * b + 1]);
+            let (xt, yt) = (itf.reference[t].0 + d[2 * t], itf.reference[t].1 + d[2 * t + 1]);
+            let (vxb, vyb, vxt, vyt) = (dd[2 * b], dd[2 * b + 1], dd[2 * t], dd[2 * t + 1]);
+            write!(
+                f,
+                ",{:.9},{:.9},{:.9},{:.9}",
+                0.5 * (xb + xt),
+                0.5 * (yb + yt),
+                0.5 * (vxb + vxt),
+                0.5 * (vyb + vyt)
+            )
+            .unwrap();
+        }
+        writeln!(f).unwrap();
+    }
+    println!(
+        "  centreline: {} instants over [{:.4}, {:.4}] s, {} stations → {out}",
+        replay.times.len(),
+        replay.times[0],
+        replay.times.last().unwrap(),
+        cols.len() + 1
+    );
+}