embedded3 S2-1 remedy: the reconstructed wall route on the cut geometry's polygons (two probes along the interpolant normal) — sphere MMS 11.3/8.3 % (the best route); wired into the DFG and flag drivers

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

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

@@ -567,6 +567,90 @@ impl Mask {
         Some([p[0] + s[0], p[1] + s[1], p[2] + s[2]])
     }
 
+    /// The reconstructed wall route (S2-1 remedy): on every wall polygon
+    /// (the cell's closure `W_c`, its centroid taken as the cell centre's
+    /// foot on the interpolant surface) the traction from two probes at
+    /// `h` and `2h` along the interpolant normal — the wall pressure by
+    /// linear extrapolation, the wall shear from the quadratic fit of the
+    /// tangential velocity through the probes and the wall velocity (the
+    /// binary wall's validated sampler, on the cut geometry's own
+    /// polygons; the probes are a cell away, past the merged slivers).
+    /// Force on the body: `Σ_c (p_w W_c + μ ∂ₙu_t |W_c| t)` over the
+    /// planes `k0..k1`; `None` without a cut geometry.
+    pub fn cut_wall_force_reconstructed(
+        &self,
+        body: &Body,
+        f: &Field,
+        mu: f64,
+        t: f64,
+        planes: Option<(usize, usize)>,
+    ) -> Option<[f64; 3]> {
+        let cut = self.cut.as_ref()?;
+        let g = self.grid;
+        let (k0, k1) = planes.unwrap_or((0, g.nz));
+        let h = g.dx.min(g.dy).min(g.dz);
+        let (d1, d2) = (h, 2.0 * h);
+        let mut force = [0.0; 3];
+        for (idx, w) in cut.wall.iter().enumerate() {
+            let area = (w[0] * w[0] + w[1] * w[1] + w[2] * w[2]).sqrt();
+            if area == 0.0 || !self.cell_fluid[idx] {
+                continue;
+            }
+            let (k, j, i) = g.kji(idx);
+            if k < k0 || k >= k1 {
+                continue;
+            }
+            let xc = [
+                (i as f64 + 0.5) * g.dx,
+                (j as f64 + 0.5) * g.dy,
+                (k as f64 + 0.5) * g.dz,
+            ];
+            let (s, n) = self.interpolant_distance_and_normal(cut, xc);
+            // The wall point and the outward (into the fluid) normal.
+            let x = [xc[0] - s * n[0], xc[1] - s * n[1], xc[2] - s * n[2]];
+            let at = |d: f64| [x[0] + d * n[0], x[1] + d * n[1], x[2] + d * n[2]];
+            let (x1, x2) = (at(d1), at(d2));
+            let (Some(p1), Some(p2)) = (
+                self.pressure_at(&f.p, x1[0], x1[1], x1[2]),
+                self.pressure_at(&f.p, x2[0], x2[1], x2[2]),
+            ) else {
+                // No fit: the operator's own pressure on this polygon.
+                for c in 0..3 {
+                    force[c] += f.p[idx] * w[c];
+                }
+                continue;
+            };
+            let p_wall = p1 + (p1 - p2) * d1 / (d2 - d1);
+            for c in 0..3 {
+                force[c] += p_wall * w[c];
+            }
+            let (Some(u1), Some(u2)) = (
+                self.velocity_at(body, f, x1[0], x1[1], x1[2], t),
+                self.velocity_at(body, f, x2[0], x2[1], x2[2], t),
+            ) else {
+                continue;
+            };
+            let us = body.surface_velocity(x[0], x[1], x[2], t);
+            let us = [us.0, us.1, us.2];
+            // Tangential components (the normal removed) and the wall
+            // gradient of the quadratic through 0, d1, d2.
+            let tang = |v: [f64; 3]| {
+                let vn = v[0] * n[0] + v[1] * n[1] + v[2] * n[2];
+                [v[0] - vn * n[0], v[1] - vn * n[1], v[2] - vn * n[2]]
+            };
+            let (t1, t2, ts) = (tang(u1), tang(u2), tang(us));
+            let wall_gradient =
+                |f1: f64, f2: f64| (f1 * d2 * d2 - f2 * d1 * d1) / (d1 * d2 * (d2 - d1));
+            for c in 0..3 {
+                let dn = wall_gradient(t1[c] - ts[c], t2[c] - ts[c]);
+                // Traction on the body = −(fluid stress on the fluid side):
+                // the shear the fluid exerts on the wall along +t.
+                force[c] += mu * dn * area;
+            }
+        }
+        Some(force)
+    }
+
     /// The cut-cell load route restricted to the cells (and faces) of the
     /// planes `k0..k1`, divided by the slab's thickness: the load per unit
     /// span on a body's mid-section.

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

@@ -141,7 +141,15 @@ 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);
+            // The reconstructed wall route (two probes on the cut polygons).
+            let fr = mask
+                .cut_wall_force_reconstructed(body, &field, RHO * NU, t, None)
+                .expect("reconstructed");
+            let fs = rtx_cfd::solvers::incompressible::embedded3::SurfaceForce {
+                f: fr,
+                samples: 0,
+                skipped: 0,
+            };
             let (cd_s, cl_s) = (coef * fs.f[0], coef * fs.f[1]);
             let zc = 0.5 * H;
             let p_front = mask
@@ -153,7 +161,7 @@ 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}, 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]",
+                "  t {t:8.4}: c_D {cd:.4} (CV {cd_cv:.4}, reconstructed {cd_s:.4} skipped {})  c_L {cl:.5} (CV {cl_cv:.5}, reconstructed {cl_s:.5})  Δp {dp:.4}  residual {:.1e} CG {}  [{:.0} s]",
                 fs.skipped,
                 r.final_residual,
                 r.poisson_iterations,
@@ -194,7 +202,7 @@ 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; 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",
+        "  FINAL ny {ny}: c_D {cd:.4} (CV {cd_cv:.4}, routes {:.2e} apart; reconstructed {:.4}, {:.2e} from CV)  c_L {cl:.5} (CV {cl_cv:.5}, reconstructed {:.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(),

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

@@ -60,7 +60,7 @@ fn ladder(resolutions: &[usize], scheme: WallScheme) -> Ladder {
             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}  F_sampler rel {a:.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_reconstructed rel {a:.3e}",
             mm.l2_velocity,
             mm.max_div,
             mm.ghost_correction,

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

@@ -258,7 +258,9 @@ pub fn measure(n: usize, scheme: WallScheme, c: (f64, f64, f64)) -> Measurement
     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,
+        WallScheme::CutCell => mask
+            .cut_wall_force_reconstructed(body, &f, MU, t, None)
+            .expect("reconstructed"),
     };
     Measurement {
         l2_velocity: (sq / vol).sqrt(),