From f1714fb926693010a4516516a03dc4f7e5078a90 Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Sun, 20 Sep 2026 21:27:26 -0500 Subject: [PATCH] =?UTF-8?q?S2-9:=20the=20flag=20on=20the=20recorded=20FSI2?= =?UTF-8?q?=20kinematics=20=E2=80=94=20fsi2=5Foverset=5Fdump=5Fcentreline?= =?UTF-8?q?=20(rtx-fsi)=20exports=20the=20saved=20instants'=20wetted-node?= =?UTF-8?q?=20d/=E1=B8=8B=20as=20a=20centreline=20per=20knot=20(clamp=20+?= =?UTF-8?q?=2070=20columns);=20tests/embedded3=5Fflag=5Fkinematics=20(cubi?= =?UTF-8?q?c=20Hermite=20through=20the=20knots=20with=20their=20velocities?= =?UTF-8?q?,=20t=20=3D=200=20at=20RTX=5FE3=5FFLAG=5FKIN=5FT0,=20a=200.5=20?= =?UTF-8?q?s=20ramp=20from=20the=20undeflected=20line,=20arc-length=20inte?= =?UTF-8?q?rpolation=20along=20the=20stations)=20drives=20both=20the=203D/?= =?UTF-8?q?slab=20flag=20(RTX=5FE3=5FFLAG=5FKINEMATICS,=20both=20in-plane?= =?UTF-8?q?=20surface-velocity=20components,=20the=20record's=20period)=20?= =?UTF-8?q?and=20the=202D=20reference;=20the=20analytic=20paths=20unchange?= =?UTF-8?q?d=20(same=20tuples)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Fable 5.1 --- .../tests/embedded3_flag_kinematics/mod.rs | 113 ++++++++++++++++++ .../tests/embedded3_flag_reference_2d.rs | 48 +++++--- .../rtx-cfd/tests/embedded3_flag_wake.rs | 61 ++++++++-- .../rtx-fsi/tests/turek_hron_fsi2_overset.rs | 68 +++++++++++ 4 files changed, 268 insertions(+), 22 deletions(-) create mode 100644 crates/specialized/rtx-cfd/tests/embedded3_flag_kinematics/mod.rs diff --git a/crates/specialized/rtx-cfd/tests/embedded3_flag_kinematics/mod.rs b/crates/specialized/rtx-cfd/tests/embedded3_flag_kinematics/mod.rs new file mode 100644 index 0000000..6e399d3 --- /dev/null +++ b/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=` +//! 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, + /// Per knot: per station (x, y, vx, vy). + pub knots: Vec>, + pub t0: f64, + pub ramp: f64, +} + +pub fn recorded() -> Option<&'static Recorded> { + static REC: std::sync::OnceLock> = 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 = 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)) + } +} diff --git a/crates/specialized/rtx-cfd/tests/embedded3_flag_reference_2d.rs b/crates/specialized/rtx-cfd/tests/embedded3_flag_reference_2d.rs index c7bcee5..17ae920 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_flag_reference_2d.rs +++ b/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, diff --git a/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs b/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs index d66d94e..51dcec8 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs +++ b/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) } diff --git a/crates/specialized/rtx-fsi/tests/turek_hron_fsi2_overset.rs b/crates/specialized/rtx-fsi/tests/turek_hron_fsi2_overset.rs index e03d442..1f3112b 100644 --- a/crates/specialized/rtx-fsi/tests/turek_hron_fsi2_overset.rs +++ b/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 + ); +}