From 83922974d78a39ae722a643d7dd16118906dbe29 Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Mon, 21 Sep 2026 19:33:59 -0500 Subject: [PATCH] =?UTF-8?q?R2-b/R3:=20embedded3=5Fadded=5Fmass=20(the=20cu?= =?UTF-8?q?t-cell=20added-mass=20pin,=20the=20overset=20instrument's=20par?= =?UTF-8?q?ameters:=20n=2064/128=20C=5Fm=201.087/1.127=20vs=20the=20confin?= =?UTF-8?q?ed=20potential-flow=201.071=20+=20Stokes=200.064=20=3D=201.135,?= =?UTF-8?q?=20the=20overset's=201.097/1.127);=20flag=20tests:=20the=20caps?= =?UTF-8?q?ule's=20tip=20inset=20(RTX=5FE3=5FFLAG=5FTIP=5FINSET,=20default?= =?UTF-8?q?=20FLAG=5FHALF=20=E2=80=94=20the=20apex=20ON=20the=20benchmark'?= =?UTF-8?q?s=20tip=20A;=20every=20earlier=20record=20had=20it=2010=20mm=20?= =?UTF-8?q?beyond,=20=3D0=20restores),=20root=20fillet=20knob=20RTX=5FE3?= =?UTF-8?q?=5FFLAG=5FFILLET=20(5=20mm:=20no=20effect),=20the=20tip=20CSV?= =?UTF-8?q?=20column=20=3D=20the=20record's=20tip=20in=20recorded=20mode?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Fable 5.1 --- .../rtx-cfd/tests/embedded3_added_mass.rs | 180 ++++++++++++++++++ .../tests/embedded3_flag_reference_2d.rs | 27 ++- .../rtx-cfd/tests/embedded3_flag_wake.rs | 65 ++++++- 3 files changed, 268 insertions(+), 4 deletions(-) create mode 100644 crates/specialized/rtx-cfd/tests/embedded3_added_mass.rs diff --git a/crates/specialized/rtx-cfd/tests/embedded3_added_mass.rs b/crates/specialized/rtx-cfd/tests/embedded3_added_mass.rs new file mode 100644 index 0000000..45189c4 --- /dev/null +++ b/crates/specialized/rtx-cfd/tests/embedded3_added_mass.rs @@ -0,0 +1,180 @@ +//! R2-b: the cut-cell added-mass pin — the overset's instrument +//! (`overset_added_mass.rs`) rebuilt on embedded3. +//! +//! A no-slip circular cylinder of radius `R` translates through a closed +//! box of still fluid, `x_c(t) = A sin ωt`, on the periodic slab (nz cells). +//! The in-phase force is the added-mass reaction `F = −m_a ẍ` with +//! `m_a = ρ π R²` (unbounded potential flow) plus Stokes's viscous +//! correction `4 / √(π β)`, `β = R² ω / ν`, plus a small blockage term for +//! the box. The operator route's F_x is fitted to `c_s sin ωt + c_c cos ωt + +//! c_0` over the last two of `periods` periods: `C_m = c_s / (ρ π R² A ω²)` +//! (with `x_c = A sin ωt`, `−m_a ẍ = m_a A ω² sin ωt`), the viscous part +//! `c_c / (ρ π R² A ω²)` against `4 / √(π β)`. +//! +//! Knobs: `RTX_E3_AM_N` (64), `RTX_E3_AM_NZ` (4), `RTX_E3_AM_PERIODS` (4), +//! `RTX_E3_AM_SCHEME` (tvd | upwind), `RTX_E3_AM_CSV=`. +use rtx_cfd::solvers::incompressible::ConvectionScheme; +use rtx_cfd::solvers::incompressible::embedded3::{ + Body, Boundaries, Field, Fluid, Grid, Parameters, Side, Solver, WallScheme, +}; +use std::f64::consts::PI; +use std::io::Write as _; + +const RHO: f64 = 1.0; +const NU: f64 = 5e-5; +const R: f64 = 0.1; +const CX: f64 = 0.5; +const CY: f64 = 0.5; +const AMP: f64 = 0.005; +const OMEGA: f64 = 2.0 * PI; + +fn env_f(name: &str, default: f64) -> f64 { + std::env::var(name) + .ok() + .and_then(|v| v.parse().ok()) + .unwrap_or(default) +} + +/// Least squares of `y ≈ c_s sin ωt + c_c cos ωt + c_0`. +fn fit(samples: &[(f64, f64)]) -> (f64, f64, f64) { + let mut m = [[0.0f64; 3]; 3]; + let mut rhs = [0.0f64; 3]; + for &(t, y) in samples { + let b = [(OMEGA * t).sin(), (OMEGA * t).cos(), 1.0]; + for i in 0..3 { + rhs[i] += b[i] * y; + for j in 0..3 { + m[i][j] += b[i] * b[j]; + } + } + } + // Gaussian elimination with partial pivoting. + let mut a = [[0.0f64; 4]; 3]; + for i in 0..3 { + a[i][..3].copy_from_slice(&m[i]); + a[i][3] = rhs[i]; + } + for i in 0..3 { + let piv = (i..3).max_by(|&p, &q| a[p][i].abs().total_cmp(&a[q][i].abs())).unwrap(); + a.swap(i, piv); + for k in 0..3 { + if k != i { + let f = a[k][i] / a[i][i]; + for j in 0..4 { + a[k][j] -= f * a[i][j]; + } + } + } + } + (a[0][3] / a[0][0], a[1][3] / a[1][1], a[2][3] / a[2][2]) +} + +#[test] +#[ignore = "R2-b: the cut-cell added-mass pin (minutes per rung on the host)"] +fn cut_cell_added_mass() { + let n = env_f("RTX_E3_AM_N", 64.0) as usize; + let nz = env_f("RTX_E3_AM_NZ", 4.0) as usize; + let periods = env_f("RTX_E3_AM_PERIODS", 4.0) as usize; + let scheme = match std::env::var("RTX_E3_AM_SCHEME").as_deref() { + Ok("upwind") => ConvectionScheme::Upwind, + _ => ConvectionScheme::TvdVanAlbada, + }; + let h = 1.0 / n as f64; + let lz = nz as f64 * h; + let period = 2.0 * PI / OMEGA; + // The overset pin's step rule (its patch spacing = h here). + let dt_raw = 0.4 * (h * h / (4.0 * NU)).min(h).min(0.2 * h / (AMP * OMEGA)); + let steps_per_period = (period / dt_raw).ceil() as usize; + let dt = period / steps_per_period as f64; + let beta = R * R * OMEGA / NU; + let stokes = 4.0 / (PI * beta).sqrt(); + println!( + " cut-cell added mass: n = {n} (R/h {:.1}), nz {nz}, dt = {dt:.3e} ({steps_per_period} per period), A/R = {:.3}, KC = {:.3}, β = {beta:.0}: Stokes C_m ≈ {:.3} (viscous {:.3})", + R / h, + AMP / R, + 2.0 * PI * AMP / R, + 1.0 + stokes, + stokes + ); + let mut solver = Solver::new( + Fluid { + density: RHO, + viscosity: RHO * NU, + reference_velocity: AMP * OMEGA, + reference_length: 2.0 * R, + }, + Parameters { + corrector_steps: 3, + inner_stop_factor: 1e-3, + tolerance: 1e-8, + convection_scheme: scheme, + wall_scheme: WallScheme::CutCell, + boundaries: Boundaries { + z0: Side::Periodic, + z1: Side::Periodic, + ..Boundaries::default() + }, + max_surface_speed: Some(AMP * OMEGA * 1.05), + ..Parameters::default() + }, + ); + solver.set_boundary_velocity(|_, _, _, _| (0.0, 0.0, 0.0)); + let xc = |t: f64| CX + AMP * (OMEGA * t).sin(); + let uc = |t: f64| AMP * OMEGA * (OMEGA * t).cos(); + let body = Body::from_sdf(move |x, y, _z, t| ((x - xc(t)).powi(2) + (y - CY).powi(2)).sqrt() - R) + .with_surface_velocity(move |_x, _y, _z, t| (uc(t), 0.0, 0.0)); + solver.set_moving_body(body); + let g = Grid::cubic(n, n, nz, h); + let mut field = Field::new(g); + solver.initialize(&mut field); + let mut csv = std::env::var("RTX_E3_AM_CSV").ok().map(|p| { + let mut f = std::fs::File::create(p).expect("csv"); + writeln!(f, "t,xc,fx,fy,fx_rec,residual,fresh").unwrap(); + f + }); + let mut samples: Vec<(f64, f64, f64)> = Vec::new(); + let mut worst = 0.0f64; + let start = std::time::Instant::now(); + for step in 0..periods * steps_per_period { + let r = solver.advance(&mut field, dt); + let t = (step + 1) as f64 * dt; + worst = worst.max(r.final_residual); + assert!(r.final_residual.is_finite(), "death at step {step}"); + let mask = solver.mask().expect("mask"); + let body = solver.body().expect("body"); + let f = mask.cut_wall_force(body, &field, RHO * NU, t).expect("cut wall"); + let fr = mask + .cut_wall_force_reconstructed(body, &field, RHO * NU, t, None) + .expect("reconstructed"); + let (fx, fy, fx_rec) = (f[0] / lz, f[1] / lz, fr[0] / lz); + samples.push((t, fx, fx_rec)); + if let Some(c) = csv.as_mut() { + writeln!(c, "{t:.6},{:.6},{fx:.6e},{fy:.6e},{fx_rec:.6e},{:.3e},{}", xc(t), r.final_residual, r.fresh_cells).unwrap(); + } + if (step + 1) % steps_per_period == 0 { + println!( + " period {}: F_x at the end {fx:+.4e} (reconstructed {fx_rec:+.4e}), residual {:.1e}, {:.0} s", + (step + 1) / steps_per_period, + r.final_residual, + start.elapsed().as_secs_f64() + ); + } + } + let t0 = (periods as f64 - 2.0) * period - 1e-12; + let scale = RHO * PI * R * R * AMP * OMEGA * OMEGA; + let op: Vec<(f64, f64)> = samples.iter().filter(|s| s.0 > t0).map(|s| (s.0, s.1)).collect(); + let rec: Vec<(f64, f64)> = samples.iter().filter(|s| s.0 > t0).map(|s| (s.0, s.2)).collect(); + let (cs, cc, c0) = fit(&op); + let (cs_r, cc_r, _) = fit(&rec); + println!( + " ADDED MASS n = {n}: C_m operator {:.4} (reconstructed {:.4}) vs Stokes {:.3}; viscous {:.4} (reconstructed {:.4}) vs {:.3}; offset {:+.3e}; worst residual {worst:.1e}; {:.0} s", + cs / scale, + cs_r / scale, + 1.0 + stokes, + cc / scale, + cc_r / scale, + stokes, + c0, + start.elapsed().as_secs_f64() + ); +} 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 d46193a..c7406b3 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_flag_reference_2d.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_flag_reference_2d.rs @@ -50,9 +50,34 @@ fn deflection(s: f64, t: f64) -> (f64, f64) { ) } +/// R3 (2026-09-21): the capsule's tip inset (`RTX_E3_FLAG_TIP_INSET`, +/// metres, default `FLAG_HALF`): the last centreline point is pulled back +/// along the last segment so the capsule's apex sits ON the benchmark's +/// tip A. Every record before this date had the apex 10 mm beyond A +/// (`=0` restores them). +fn tip_inset() -> f64 { + std::env::var("RTX_E3_FLAG_TIP_INSET") + .ok() + .and_then(|v| v.parse().ok()) + .unwrap_or(FLAG_HALF) +} + /// 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. +/// mode, or (S2-9) the recorded FSI2 kinematics along its stations; the +/// last point carries the tip inset. fn centreline(m: usize, n: usize, t: f64) -> (f64, f64, f64, f64) { + let inset = tip_inset(); + if m < n || inset <= 0.0 { + return centreline_raw(m, n, t); + } + let (ax, ay, _, _) = centreline_raw(n - 1, n, t); + let (bx, by, bvx, bvy) = centreline_raw(n, n, t); + let len = ((bx - ax).powi(2) + (by - ay).powi(2)).sqrt(); + let f = (1.0 - inset / len).max(0.0); + (ax + f * (bx - ax), ay + f * (by - ay), bvx, bvy) +} + +fn centreline_raw(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! { diff --git a/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs b/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs index a2b4c09..a958c55 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs @@ -80,6 +80,45 @@ fn amplitude() -> f64 { env_f("RTX_E3_FLAG_AMP", AMP) } +/// R3: the root fillet radius (`RTX_E3_FLAG_FILLET`, metres; 0 = the sharp +/// concave corner of the plain union; the overset outline carries 5 mm). +fn root_fillet() -> f64 { + env_f("RTX_E3_FLAG_FILLET", 0.0) +} + +/// R3 (2026-09-21): the capsule's tip inset (`RTX_E3_FLAG_TIP_INSET`, +/// metres, default `FLAG_HALF`): the centreline polyline is shortened by +/// this much along its last segment so the capsule's apex sits ON the +/// benchmark's tip A (the overset outline's semicircle apex). Every flag +/// record before this date had the apex 10 mm beyond A (`=0` restores +/// them): at ny 62 on the recorded motion that was drag 253.3 → 240.0 and +/// the lift swing 1,005 → 836 (the overset's 867). +fn tip_inset() -> f64 { + env_f("RTX_E3_FLAG_TIP_INSET", FLAG_HALF) +} + +/// Smooth union with a concave fillet of radius `r` (the plain `min` at r = 0). +fn fillet_union(d1: f64, d2: f64, r: f64) -> f64 { + if r > 0.0 && d1 < r && d2 < r { + r - ((r - d1).powi(2) + (r - d2).powi(2)).sqrt() + } else { + d1.min(d2) + } +} + +/// Pull the polyline's last point back by `inset` along its last segment. +fn inset_last(pts: &mut [(f64, f64, f64, f64)], inset: f64) { + if inset <= 0.0 || pts.len() < 2 { + return; + } + let n = pts.len(); + let (ax, ay, _, _) = pts[n - 2]; + let (bx, by, bvx, bvy) = pts[n - 1]; + let len = ((bx - ax).powi(2) + (by - ay).powi(2)).sqrt(); + let f = (1.0 - inset / len).max(0.0); + pts[n - 1] = (ax + f * (bx - ax), ay + f * (by - ay), bvx, bvy); +} + /// Signed distance to the deflected flag's cross-section (a capsule /// around the centreline polyline of `n` segments) and the centreline's /// velocity at the closest point (transverse only in the analytic mode). @@ -101,6 +140,7 @@ fn flag_2d_recorded(rec: &Recorded, x: f64, y: f64, t: f64) -> (f64, (f64, f64)) let mut c = cell.borrow_mut(); if c.0.to_bits() != t.to_bits() { c.1 = rec.at(t, CY); + inset_last(&mut c.1, tip_inset()); c.0 = t; } let pts = &c.1; @@ -141,6 +181,14 @@ fn flag_2d_analytic(x: f64, y: f64, t: f64) -> (f64, f64) { let (d, v) = deflection(s, t); *p = (FLAG_X0 + s * FLAG_LEN, CY + d, v); } + let inset = tip_inset(); + if inset > 0.0 { + let (ax, ay, _) = c.1[N - 1]; + let (bx, by, bv) = c.1[N]; + let len = ((bx - ax).powi(2) + (by - ay).powi(2)).sqrt(); + let f = (1.0 - inset / len).max(0.0); + c.1[N] = (ax + f * (bx - ax), ay + f * (by - ay), bv); + } c.0 = t; } c.1 @@ -273,7 +321,10 @@ 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)) + let r_fillet = root_fillet(); + let body = Body::from_sdf(move |x, y, z, t| { + fillet_union(cyl(x, y), flag_3d(x, y, z, t, r_edge).0, r_fillet) + }) .with_surface_velocity(move |x, y, z, t| { let (df, (vx, vy)) = flag_3d(x, y, z, t, r_edge); if df <= cyl(x, y) { @@ -295,7 +346,7 @@ fn flag_wake_on_the_device() { } solver.initialize(&mut field); println!( - " flag wake ny {ny} (span {}; inflow {}, z sides {}): {nx}×{ny}×{nz} = {} cells, h {h:.4e}, dt {dt:.3e}, {periods} periods = {t_end:.3} s, {} steps", + " flag wake ny {ny} (span {}; inflow {}, z sides {}; root fillet {:.4} m, tip inset {:.4} m): {nx}×{ny}×{nz} = {} cells, h {h:.4e}, dt {dt:.3e}, {periods} periods = {t_end:.3} s, {} steps", flag_span(), if slab_nz > 0 || inflow_2d { "2d" } else { "3d" }, if slab_nz > 0 { @@ -305,6 +356,8 @@ fn flag_wake_on_the_device() { } else { "wall" }, + root_fillet(), + tip_inset(), g.cells(), (t_end / dt).ceil() as usize ); @@ -368,7 +421,13 @@ fn flag_wake_on_the_device() { let ft = mask .cut_wall_force(body, &field, RHO * NU, t) .expect("wall"); - let tip = deflection(1.0, t).0; + // The tip's transverse deflection: the record's last station in + // recorded mode (until 2026-09-21 this column held the analytic + // first mode even then — R2's fits use the record directly). + let tip = match recorded() { + Some(rec) => rec.at(t, CY).last().map_or(0.0, |p| p.1 - CY), + None => deflection(1.0, t).0, + }; if sample { println!( " t {t:7.4} (tip {tip:+.4}): drag/span {:.1} lift/span {:+.1} N/m (reconstructed {:.1} {:+.1}); total {:.3} {:+.3} N; residual {:.1e} CG {} fresh {}; [{:.0} s]",