embedded3 S2-1: the traction sampler on the cut field (pressure fit excludes merged cells; sampler route in the sphere MMS and DFG drivers)

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

crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/loads.rs

@@ -83,7 +83,10 @@ impl Mask {
                 }
             }
         }
-        if nodes.iter().all(|&(idx, _)| self.is_fluid_cell(idx)) {
+        // A virtually merged small cell carries its master's pressure
+        // unknown, not its own: it is dropped from the fit.
+        let usable = |idx: usize| self.is_fluid_cell(idx) && self.master(idx).is_none();
+        if nodes.iter().all(|&(idx, _)| usable(idx)) {
             return Some(nodes.iter().map(|&(idx, w)| w * p[idx]).sum());
         }
         // Weighted by the z-direction weight (a z-invariant field then fits
@@ -96,7 +99,7 @@ impl Mask {
             for _ in 0..dirs(ny, fy).len() * dirs(nx, fx).len() {
                 let &(idx, _) = it.next().expect("node");
                 let _ = dk;
-                if self.is_fluid_cell(idx) {
+                if usable(idx) {
                     let (k, j, i) = g.kji(idx);
                     pts.push((
                         (i as f64 + 0.5) * dx,

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

@@ -75,7 +75,12 @@ fn dfg_3d_2z_on_the_device() {
             (0.0, 0.0, 0.0)
         }
     });
-    solver.set_body(Body::cylinder_z(CX, CY, 0.5 * D));
+    // The cylinder extruded across the width, with samples for the
+    // traction route (S2-1).
+    solver.set_body(Body::extruded(
+        rtx_cfd::solvers::incompressible::EmbeddedBody::circle(CX, CY, 0.5 * D),
+        H,
+    ));
     let g = Grid::cubic(nx, ny, nz, h);
     let mut field = Field::new(g);
     // Start from the inflow profile everywhere (a faster approach to steady).
@@ -99,13 +104,18 @@ fn dfg_3d_2z_on_the_device() {
     let coef = 2.0 / (RHO * U_BAR * U_BAR * D * H);
     let csv = std::env::var("RTX_E3_DFG_CSV").ok().map(|p| {
         let mut f = std::fs::File::create(p).expect("csv");
-        writeln!(f, "t,cd_wall,cl_wall,cd_cv,cl_cv,dp,residual,cg").unwrap();
+        writeln!(
+            f,
+            "t,cd_wall,cl_wall,cd_cv,cl_cv,dp,residual,cg,cd_sampler,cl_sampler"
+        )
+        .unwrap();
         f
     });
     let mut csv = csv;
     let sample_every = (steps / 100).max(1);
     let start = std::time::Instant::now();
     let mut last: Option<(f64, f64, f64, f64, f64)> = None;
+    let mut last_sampler = (f64::NAN, f64::NAN);
     let mut settled = false;
     for step in 0..steps {
         let r = device.advance(dt);
@@ -131,6 +141,8 @@ fn dfg_3d_2z_on_the_device() {
             );
             let fcv =
                 mask.control_volume_force_with_walls(&field, dt, RHO, RHO * NU, None, bx, true);
+            let fs = mask.surface_force(body, &field, RHO * NU, t, 0.5 * h);
+            let (cd_s, cl_s) = (coef * fs.f[0], coef * fs.f[1]);
             let zc = 0.5 * H;
             let p_front = mask
                 .pressure_at(&field.p, CX - 0.5 * D, CY, zc)
@@ -141,7 +153,8 @@ fn dfg_3d_2z_on_the_device() {
             let dp = p_front - p_back;
             let (cd, cl, cd_cv, cl_cv) = (coef * fw[0], coef * fw[1], coef * fcv[0], coef * fcv[1]);
             println!(
-                "  t {t:8.4}: c_D {cd:.4} (CV {cd_cv:.4})  c_L {cl:.5} (CV {cl_cv:.5})  Δp {dp:.4}  residual {:.1e} CG {}  [{:.0} s]",
+                "  t {t:8.4}: c_D {cd:.4} (CV {cd_cv:.4}, sampler {cd_s:.4} skipped {})  c_L {cl:.5} (CV {cl_cv:.5}, sampler {cl_s:.5})  Δp {dp:.4}  residual {:.1e} CG {}  [{:.0} s]",
+                fs.skipped,
                 r.final_residual,
                 r.poisson_iterations,
                 start.elapsed().as_secs_f64()
@@ -149,11 +162,12 @@ fn dfg_3d_2z_on_the_device() {
             if let Some(f) = csv.as_mut() {
                 writeln!(
                     f,
-                    "{t:.5},{cd:.6},{cl:.6},{cd_cv:.6},{cl_cv:.6},{dp:.6},{:.3e},{}",
+                    "{t:.5},{cd:.6},{cl:.6},{cd_cv:.6},{cl_cv:.6},{dp:.6},{:.3e},{},{cd_s:.6},{cl_s:.6}",
                     r.final_residual, r.poisson_iterations
                 )
                 .unwrap();
             }
+            last_sampler = (cd_s, cl_s);
             if let Some((pcd, pcl, _, _, pdp)) = last {
                 let rel = ((cd - pcd) / cd)
                     .abs()
@@ -180,8 +194,11 @@ fn dfg_3d_2z_on_the_device() {
     }
     let (cd, cl, cd_cv, cl_cv, dp) = last.expect("samples");
     println!(
-        "  FINAL ny {ny}: c_D {cd:.4} (CV {cd_cv:.4}, routes {:.2e} apart)  c_L {cl:.5} (CV {cl_cv:.5})  Δp {dp:.4}  — reference c_D 6.05–6.25, c_L 0.008–0.010, Δp 0.165–0.175; {:.0} s",
+        "  FINAL ny {ny}: c_D {cd:.4} (CV {cd_cv:.4}, routes {:.2e} apart; sampler {:.4}, {:.2e} from CV)  c_L {cl:.5} (CV {cl_cv:.5}, sampler {:.5})  Δp {dp:.4}  — reference c_D 6.05–6.25, c_L 0.008–0.010, Δp 0.165–0.175; {:.0} s",
         ((cd - cd_cv) / cd).abs(),
+        last_sampler.0,
+        ((last_sampler.0 - cd_cv) / cd_cv).abs(),
+        last_sampler.1,
         start.elapsed().as_secs_f64()
     );
     if let Some(t) = device.timers() {

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

@@ -54,8 +54,13 @@ fn ladder(resolutions: &[usize], scheme: WallScheme) -> Ladder {
             mm.force_cv[1] - fcv[1],
             mm.force_cv[2] - fcv[2],
         ]) / f_scale;
+        let a = norm([
+            mm.force_sampler[0] - fe[0],
+            mm.force_sampler[1] - fe[1],
+            mm.force_sampler[2] - fe[2],
+        ]) / f_scale;
         println!(
-            "  n = {n:3}  L2 u {:.4e} (order {rate})  max div {:.2e}  ghost corr {:.2e}  F_surface {:.4?} rel {s:.3e} (skipped {})  F_cv {:.4?} rel {c:.3e}",
+            "  n = {n:3}  L2 u {:.4e} (order {rate})  max div {:.2e}  ghost corr {:.2e}  F_surface {:.4?} rel {s:.3e} (skipped {})  F_cv {:.4?} rel {c:.3e}  F_sampler rel {a:.3e}",
             mm.l2_velocity,
             mm.max_div,
             mm.ghost_correction,
@@ -63,6 +68,15 @@ fn ladder(resolutions: &[usize], scheme: WallScheme) -> Ladder {
             mm.skipped,
             mm.force_cv
         );
+        if scheme == WallScheme::CutCell {
+            // S2-1 (reported, not gated): the traction sampler on the cut
+            // field against the box route — 1.15× at n 24 on the first read
+            // (`docs/embedded3_campaign.md`, S2-1).
+            println!(
+                "    sampler / box error ratio at n = {n}: {:.3}",
+                a / c.max(1e-300)
+            );
+        }
         se.push(s);
         ce.push(c);
     }

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

@@ -107,9 +107,14 @@ pub struct Measurement {
     pub l2_velocity: f64,
     pub max_div: f64,
     pub ghost_correction: f64,
+    /// The scheme's wall route (the traction sampler on the binary wall,
+    /// the operator route on the cut wall).
     pub force_surface: [f64; 3],
     pub skipped: usize,
     pub force_cv: [f64; 3],
+    /// The traction sampler on the cut wall (S2-1); the wall route again
+    /// on the binary wall.
+    pub force_sampler: [f64; 3],
 }
 
 /// March the manufactured solution with the sphere at `c` to steady state on grid `n`.
@@ -250,6 +255,10 @@ pub fn measure(n: usize, scheme: WallScheme, c: (f64, f64, f64)) -> Measurement
     let (i0, i1) = (n / 8, n - n / 8);
     let src = |x: f64, y: f64, z: f64| source3(x, y, z);
     let force_cv = mask.control_volume_force(&f, dt, RHO, MU, Some(&src), (i0, i1, i0, i1, i0, i1));
+    let force_sampler = match scheme {
+        WallScheme::GhostBinary => surface.f,
+        WallScheme::CutCell => mask.surface_force(body, &f, MU, t, 0.5 * h).f,
+    };
     Measurement {
         l2_velocity: (sq / vol).sqrt(),
         max_div,
@@ -257,5 +266,6 @@ pub fn measure(n: usize, scheme: WallScheme, c: (f64, f64, f64)) -> Measurement
         force_surface: surface.f,
         skipped: surface.skipped,
         force_cv,
+        force_sampler,
     }
 }