R2-b/R3: embedded3_added_mass (the cut-cell added-mass pin, the overset instrument's parameters: n 64/128 C_m 1.087/1.127 vs the confined potential-flow 1.071 + Stokes 0.064 = 1.135, the overset's 1.097/1.127); flag tests: the capsule's tip inset (RTX_E3_FLAG_TIP_INSET, default FLAG_HALF — the apex ON the benchmark's tip A; every earlier record had it 10 mm beyond, =0 restores), root fillet knob RTX_E3_FLAG_FILLET (5 mm: no effect), the tip CSV column = the record's tip in recorded mode

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

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=<path>`.
+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()
+    );
+}

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! {

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]",