embedded3 S2-5: the cut wall sat ½(1−α)h inside the body — cross diffusion over the open-part centroid spacing (RTX_E3_DIFFUSION_CENTROID; host + e3_cut.cu, shift tables, point-implicit excess); flat-wall effective-position instrument; DFG 2D-1 ladder tests (device + host); knobs tried and refuted along the way (oblique distance, axis exchange, centroid pressure gradient)

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

crates/specialized/rtx-cfd/src/kernels/cuda/e3_cut.cu

@@ -26,6 +26,7 @@ struct E3Cut {
     const unsigned int *owner;          /* the master of a merged cell (itself otherwise) */
     const unsigned int *fold_ptr, *fold_idx; /* CSR: the slaves of every cell */
     double *cell_flux;                  /* scratch per cell */
+    const double *s_u, *s_v, *s_w;      /* open-part centroid shifts per face, 3 interleaved (S2-5) */
 };
 
 /* Face index of component c at lattice (i, j, k); −1 outside (z wraps when periodic). */
@@ -95,14 +96,20 @@ __device__ double cut_face_update(const E3Params& g, const E3Ptrs& f, const E3Cu
     for (int d = 0; d < 3; ++d) wall[d] = -(app[d] - apm[d]) * area[d];
     double h_min = fmin(fmin(g.dx, g.dy), g.dz);
     const double* dist = c == 0 ? m.d_u : (c == 1 ? m.d_v : m.d_w);
-    double distance = fmax(dist[fidx] + 0.5 * h[c] * (1.0 - alpha), CUT_DISTANCE_FLOOR * h_min);
+    /* bit 4 of wall_order: the oblique wall distance (the in-plane part of the wall normal) */
+    double n_t = 1.0;
+    if (g.wall_order & 16) {
+        double a_w0 = sqrt(wall[0] * wall[0] + wall[1] * wall[1] + wall[2] * wall[2]);
+        if (a_w0 > 0.0) { double n_c = wall[c] / a_w0; n_t = sqrt(fmax(1.0 - n_c * n_c, 0.0)); }
+    }
+    double distance = fmax(dist[fidx] + 0.5 * h[c] * (1.0 - alpha) * n_t, CUT_DISTANCE_FLOOR * h_min);
     const double* ubt = c == 0 ? m.ub_u : (c == 1 ? m.ub_v : m.ub_w);
     double ub = ubt[fidx];
     /* face position */
     double x[3];
     for (int d = 0; d < 3; ++d) x[d] = (p[d] + (d == c ? 0.0 : 0.5)) * h[d];
 
-    double mass_out = 0.0, conv = 0.0, diff = 0.0;
+    double mass_out = 0.0, conv = 0.0, diff = 0.0, wall_implicit = 0.0, wall_rhs = 0.0;
     for (int d = 0; d < 3; ++d) {
         double a_d = area[d];
         int q[3];
@@ -171,9 +178,41 @@ __device__ double cut_face_update(const E3Params& g, const E3Ptrs& f, const E3Cu
         conv += m_plus * (u_plus + delta_plus) - m_minus * (u_minus + delta_minus);
         /* diffusion */
         double g_minus = apm[d], g_plus = app[d];
-        if (f_up1 >= 0) diff += mu * g_plus * a_d * (up1 - u0) / h[d];
+        /* bit 5 of wall_order: the exchange with a SOLID neighbour over the axis
+           distance to the wall, implicit (S2-5) */
+        double delta_x = h[d];
+        if (g.wall_order & 32) {
+            double a_w0 = sqrt(wall[0] * wall[0] + wall[1] * wall[1] + wall[2] * wall[2]);
+            if (a_w0 > 0.0) {
+                double n_d = fabs(wall[d]) / a_w0;
+                if (n_d >= 1e-12) delta_x = fmin(distance / n_d, h[d]);
+            }
+        }
+        int solid_up = (g.wall_order & 32) && f_up1 >= 0 && cut_ap(m, c)[f_up1] == 0.0;
+        int solid_dn = (g.wall_order & 32) && f_dn1 >= 0 && cut_ap(m, c)[f_dn1] == 0.0;
+        /* bit 6 of wall_order: cross diffusion over the open-part centroid spacing,
+           the part beyond 1/h point-implicit (S2-5) */
+        int centroid = (g.wall_order & 64) && d != c;
+        const double* sh = c == 0 ? m.s_u : (c == 1 ? m.s_v : m.s_w);
+        if (solid_up) { double kx = mu * g_plus * a_d / delta_x; wall_implicit += kx; wall_rhs += kx * up1; }
+        else if (f_up1 >= 0 && centroid) {
+            double dl = h[d] + (sh[3 * f_up1 + d] - sh[3 * fidx + d]);
+            dl = fmin(fmax(dl, 0.25 * h[d]), 2.0 * h[d]);
+            double kx = mu * g_plus * a_d * (1.0 / dl - 1.0 / h[d]);
+            if (kx > 0.0) { wall_implicit += kx; wall_rhs += kx * up1; diff += mu * g_plus * a_d * (up1 - u0) / h[d]; }
+            else diff += (mu * g_plus * a_d / h[d] + kx) * (up1 - u0);
+        }
+        else if (f_up1 >= 0) diff += mu * g_plus * a_d * (up1 - u0) / h[d];
         else if (sides[d][1] == SIDE_VELOCITY) diff += mu * g_plus * a_d * (beyond_p - u0) / (0.5 * h[d]);
-        if (f_dn1 >= 0) diff -= mu * g_minus * a_d * (u0 - dn1) / h[d];
+        if (solid_dn) { double kx = mu * g_minus * a_d / delta_x; wall_implicit += kx; wall_rhs += kx * dn1; }
+        else if (f_dn1 >= 0 && centroid) {
+            double dl = h[d] - (sh[3 * f_dn1 + d] - sh[3 * fidx + d]);
+            dl = fmin(fmax(dl, 0.25 * h[d]), 2.0 * h[d]);
+            double kx = mu * g_minus * a_d * (1.0 / dl - 1.0 / h[d]);
+            if (kx > 0.0) { wall_implicit += kx; wall_rhs += kx * dn1; diff -= mu * g_minus * a_d * (u0 - dn1) / h[d]; }
+            else diff -= (mu * g_minus * a_d / h[d] + kx) * (u0 - dn1);
+        }
+        else if (f_dn1 >= 0) diff -= mu * g_minus * a_d * (u0 - dn1) / h[d];
         else if (sides[d][0] == SIDE_VELOCITY) diff -= mu * g_minus * a_d * (u0 - beyond_m) / (0.5 * h[d]);
     }
     /* The mass fluxes above are volume fluxes: the momentum flux carries rho. */
@@ -191,7 +230,7 @@ __device__ double cut_face_update(const E3Params& g, const E3Ptrs& f, const E3Cu
        open face away from the body along the wall normal's dominant axis
        (order 2; the neighbour's old value explicit). */
     double c1 = 1.0 / distance, shear_explicit = 0.0;
-    if (g.wall_order >= 2 && a_w > 0.0) {
+    if ((g.wall_order & 15) >= 2 && a_w > 0.0) {
         double nw[3] = { wall[0] / a_w, wall[1] / a_w, wall[2] / a_w };
         int d = 0;
         for (int kk = 1; kk < 3; ++kk) if (fabs(nw[kk]) > fabs(nw[d])) d = kk;
@@ -208,7 +247,7 @@ __device__ double cut_face_update(const E3Params& g, const E3Ptrs& f, const E3Cu
     double shear = mu * a_w * c1;
     double v_eff = fmax(alpha, CUT_INERTIA_FLOOR) * h[c] * area[c];
     double inertia = rho * v_eff / g.dt;
-    return (inertia * u0 - conv + diff + pressure + source + shear * ub - shear_explicit) / (inertia + shear);
+    return (inertia * u0 - conv + diff + pressure + source + shear * ub - shear_explicit + wall_rhs) / (inertia + shear + wall_implicit);
 }
 
 /* The predictor on the open interior faces of component c. */

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

@@ -210,6 +210,11 @@ impl Mask {
             scheme: crate::solvers::incompressible::ConvectionScheme::Upwind,
             density: 1.0,
             wall_order: 1,
+            wall_distance_oblique: false,
+            wall_exchange_axis: false,
+            grad_weights: None,
+            diffusion_centroid: false,
+            face_shifts: None,
         };
         mask.compute_merging(None);
         Ok(mask)
@@ -307,6 +312,94 @@ impl Mask {
         self.compute_merging(Some(old));
     }
 
+    /// The shift of a face's open-part centroid from the face centre:
+    /// `½h(1 − α)` along the wall normal's in-plane part, away from the
+    /// body (zero for a full face or without the centroid diffusion). Read
+    /// from the tables of [`Self::compute_face_shifts`].
+    pub(super) fn face_shift(&self, c: usize, p: [i64; 3]) -> [f64; 3] {
+        let (Some(t), Some(f)) = (self.face_shifts.as_ref(), self.lattice().face(c, p)) else {
+            return [0.0; 3];
+        };
+        [t[c][3 * f], t[c][3 * f + 1], t[c][3 * f + 2]]
+    }
+
+    /// The per-face shift tables (three components interleaved).
+    #[must_use]
+    pub fn face_shift_tables(&self) -> Option<&[Vec<f64>; 3]> {
+        self.face_shifts.as_ref()
+    }
+
+    /// Build the open-part centroid shifts of every cut face (S2-5).
+    pub fn compute_face_shifts(&mut self) {
+        let g = self.grid;
+        let h = [g.dx, g.dy, g.dz];
+        let lat = self.lattice();
+        let sizes = [g.n_ufaces(), g.n_vfaces(), g.n_wfaces()];
+        let mut tables = [
+            vec![0.0; 3 * sizes[0]],
+            vec![0.0; 3 * sizes[1]],
+            vec![0.0; 3 * sizes[2]],
+        ];
+        for c in 0..3 {
+            let (ni, nj, nk) = (
+                g.nx + usize::from(c == 0),
+                g.ny + usize::from(c == 1),
+                g.nz + usize::from(c == 2),
+            );
+            for k in 0..nk {
+                for j in 0..nj {
+                    for i in 0..ni {
+                        let p = [i as i64, j as i64, k as i64];
+                        let Some(f) = lat.face(c, p) else { continue };
+                        let Some(alpha) = self.aperture(c, p) else {
+                            continue;
+                        };
+                        if alpha <= 0.0 || alpha >= 1.0 {
+                            continue;
+                        }
+                        // Interior faces only (a control volume needs both cells).
+                        let on_side = p[c] == 0 || p[c] as usize == [g.nx, g.ny, g.nz][c];
+                        if on_side && !(c == 2 && self.periodic_z) {
+                            continue;
+                        }
+                        let cv = self.cv_geometry(c, p);
+                        let mut n = cv.wall;
+                        n[c] = 0.0;
+                        let a = (n[0] * n[0] + n[1] * n[1] + n[2] * n[2]).sqrt();
+                        if a == 0.0 {
+                            continue;
+                        }
+                        for d in 0..3 {
+                            // `wall` points into the body: the open part lies the other way.
+                            tables[c][3 * f + d] = -0.5 * h[d] * (1.0 - alpha) * n[d] / a;
+                        }
+                    }
+                }
+            }
+        }
+        self.face_shifts = Some(tables);
+    }
+
+    /// The distance over which a fluid face exchanges momentum with a solid
+    /// neighbour face along `d`: the full spacing, or (S2-5) the axis
+    /// distance from the open part's centroid to the wall, `min(h, d_f/|n_d|)`.
+    pub(super) fn exchange_delta(&self, cv: &CvGeometry, d: usize) -> f64 {
+        let h = [self.grid.dx, self.grid.dy, self.grid.dz][d];
+        if !self.wall_exchange_axis {
+            return h;
+        }
+        let a_w =
+            (cv.wall[0] * cv.wall[0] + cv.wall[1] * cv.wall[1] + cv.wall[2] * cv.wall[2]).sqrt();
+        if a_w == 0.0 {
+            return h;
+        }
+        let n_d = cv.wall[d].abs() / a_w;
+        if n_d < 1e-12 {
+            return h;
+        }
+        (cv.distance / n_d).min(h)
+    }
+
     /// The wall-gradient coefficients of the unknown face of component
     /// `c` at `p` with control volume `cv`: `u'(0) = c_1 (u_f − U_b) + c_2
     /// (u_n − U_b)` with `u_n` the face returned (one lattice step away
@@ -420,7 +513,21 @@ impl Mask {
                 _ => cut.d_w[f],
             }
         });
-        let distance = (phi_face + 0.5 * h[c] * (1.0 - alpha)).max(DISTANCE_FLOOR * h_min);
+        // The open part's centroid sits ½h(1 − α) from the face centre IN THE
+        // FACE PLANE: its wall distance gains that times the wall normal's
+        // in-plane part (1 for a wall parallel to the face normal).
+        let n_t = if self.wall_distance_oblique {
+            let a_w = (wall[0] * wall[0] + wall[1] * wall[1] + wall[2] * wall[2]).sqrt();
+            if a_w > 0.0 {
+                let n_c = wall[c] / a_w;
+                (1.0 - n_c * n_c).max(0.0).sqrt()
+            } else {
+                1.0
+            }
+        } else {
+            1.0
+        };
+        let distance = (phi_face + 0.5 * h[c] * (1.0 - alpha) * n_t).max(DISTANCE_FLOOR * h_min);
         CvGeometry {
             alpha,
             ap,
@@ -670,6 +777,10 @@ impl Mask {
                 }
             }
         }
+        let gw = self.gradient_weight_force(&f.p, Some((k0, k1)));
+        for c in 0..3 {
+            pressure[c] += gw[c];
+        }
         let lat = self.lattice();
         let values: [&[f64]; 3] = [&f.u, &f.v, &f.w];
         let w_range = if self.periodic_z { 0..nz } else { 1..nz };

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

@@ -112,7 +112,8 @@ impl Mask {
                                     };
                                     let u_face = upwind(m_plus, u0, un) + delta;
                                     let on_fluid = -rho * m_plus * (u_face - u0)
-                                        + mu * cv.ap[d][1] * a_d * (un - u0) / h[d];
+                                        + mu * cv.ap[d][1] * a_d * (un - u0)
+                                            / self.exchange_delta(&cv, d);
                                     force[c] -= on_fluid;
                                 }
                             }
@@ -129,7 +130,8 @@ impl Mask {
                                     };
                                     let u_face = upwind(m_minus, ud, u0) + delta;
                                     let on_fluid = rho * m_minus * (u_face - u0)
-                                        + mu * cv.ap[d][0] * a_d * (ud - u0) / h[d];
+                                        + mu * cv.ap[d][0] * a_d * (ud - u0)
+                                            / self.exchange_delta(&cv, d);
                                     force[c] -= on_fluid;
                                 }
                             }
@@ -210,4 +212,102 @@ impl Mask {
         let x = self.cut_wall_exchange_force(body, &old, mu, rho, t, None)?;
         Some([p[0] + s[0] + x[0], p[1] + s[1] + x[1], p[2] + s[2] + x[2]])
     }
+
+    /// S2-5 host prototype: per-face pressure-gradient weights `ω = h/δ`,
+    /// `δ` the axis distance between the two cells' fluid centroids. A cut
+    /// cell of fluid fraction `v` with unit wall normal `n` (into the body)
+    /// has its centroid shifted by `−½(1 − v) h n` from the cell centre
+    /// (exact for an axis-aligned cut); `δ` is clamped to `[¼h, 2h]`.
+    pub fn compute_gradient_weights(&mut self) {
+        let Some(cut) = self.cut.as_ref() else {
+            return;
+        };
+        let g = self.grid;
+        let h = [g.dx, g.dy, g.dz];
+        let cells = g.cells();
+        let mut shift = vec![[0.0f64; 3]; cells];
+        for idx in 0..cells {
+            if !self.cell_fluid[idx] {
+                continue;
+            }
+            let w = cut.wall[idx];
+            let a = (w[0] * w[0] + w[1] * w[1] + w[2] * w[2]).sqrt();
+            if a == 0.0 {
+                continue;
+            }
+            let v = self.vol(idx);
+            for d in 0..3 {
+                shift[idx][d] = -0.5 * (1.0 - v) * h[d] * w[d] / a;
+            }
+        }
+        let weight = |d: usize, minus: usize, plus: usize| -> f64 {
+            if !(self.cell_fluid[minus] && self.cell_fluid[plus]) {
+                return 1.0;
+            }
+            let delta = h[d] + shift[plus][d] - shift[minus][d];
+            h[d] / delta.clamp(0.25 * h[d], 2.0 * h[d])
+        };
+        let mut wu = vec![1.0; g.n_ufaces()];
+        let mut wv = vec![1.0; g.n_vfaces()];
+        let mut ww = vec![1.0; g.n_wfaces()];
+        for k in 0..g.nz {
+            for j in 0..g.ny {
+                for i in 0..g.nx {
+                    let c = g.cell(k, j, i);
+                    if i > 0 {
+                        wu[g.uface(k, j, i)] = weight(0, g.cell(k, j, i - 1), c);
+                    }
+                    if j > 0 {
+                        wv[g.vface(k, j, i)] = weight(1, g.cell(k, j - 1, i), c);
+                    }
+                    if k > 0 {
+                        ww[g.wface(k, j, i)] = weight(2, g.cell(k - 1, j, i), c);
+                    }
+                }
+            }
+        }
+        self.grad_weights = Some((wu, wv, ww));
+    }
+
+    /// The pressure force the gradient weights add to the closure sum
+    /// `Σ p_c W_c` (a force on the body): `Σ_f α_f A (ω_f − 1)(p_+ − p_−)`.
+    #[must_use]
+    pub fn gradient_weight_force(&self, p: &[f64], planes: Option<(usize, usize)>) -> [f64; 3] {
+        let mut force = [0.0; 3];
+        if self.grad_weights.is_none() {
+            return force;
+        }
+        let g = self.grid;
+        let area = [g.dy * g.dz, g.dx * g.dz, g.dx * g.dy];
+        let (k0, k1) = planes.unwrap_or((0, g.nz));
+        for k in k0..k1 {
+            for j in 0..g.ny {
+                for i in 0..g.nx {
+                    let c = g.cell(k, j, i);
+                    if i > 0 {
+                        let f = g.uface(k, j, i);
+                        force[0] += self.a_u(f)
+                            * area[0]
+                            * (self.grad_weight(0, f) - 1.0)
+                            * (p[c] - p[g.cell(k, j, i - 1)]);
+                    }
+                    if j > 0 {
+                        let f = g.vface(k, j, i);
+                        force[1] += self.a_v(f)
+                            * area[1]
+                            * (self.grad_weight(1, f) - 1.0)
+                            * (p[c] - p[g.cell(k, j - 1, i)]);
+                    }
+                    if k > 0 {
+                        let f = g.wface(k, j, i);
+                        force[2] += self.a_w(f)
+                            * area[2]
+                            * (self.grad_weight(2, f) - 1.0)
+                            * (p[c] - p[g.cell(k - 1, j, i)]);
+                    }
+                }
+            }
+        }
+        force
+    }
 }

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

@@ -117,6 +117,9 @@ impl Solver {
         let mut mass_out = 0.0;
         let mut conv = 0.0;
         let mut diff = 0.0;
+        let shift0 = mask.face_shift(c, p);
+        // The implicit exchange with solid neighbour faces (S2-5).
+        let (mut wall_implicit, mut wall_rhs) = (0.0, 0.0);
         for d in 0..3 {
             let ed = e(d);
             let a_d = area[d];
@@ -187,7 +190,32 @@ impl Solver {
             conv += m_plus * (u_plus + delta_plus) - m_minus * (u_minus + delta_minus);
             // Diffusion through the plus / minus faces.
             let (g_minus, g_plus) = (cv.ap[d][0], cv.ap[d][1]);
+            // The cross-direction spacing between open-part centroids (S2-5).
+            let centroid = d != c && mask.diffusion_centroid;
+            let spacing = |q: [i64; 3], sign: f64| -> f64 {
+                let delta = h[d] + sign * (mask.face_shift(c, q)[d] - shift0[d]);
+                delta.clamp(0.25 * h[d], 2.0 * h[d])
+            };
+            let solid = |q: [i64; 3]| mask.wall_exchange_axis && ap(c, q) == Some(0.0);
             diff += match up1 {
+                Some(un) if solid(add(p, ed, 1)) => {
+                    let k = mu * g_plus * a_d / mask.exchange_delta(&cv, d);
+                    wall_implicit += k;
+                    wall_rhs += k * un;
+                    0.0
+                }
+                Some(un) if centroid => {
+                    // The part of the centroid coupling beyond `1/h` is taken
+                    // point-implicitly (the explicit limit is `h`'s).
+                    let k = mu * g_plus * a_d * (1.0 / spacing(add(p, ed, 1), 1.0) - 1.0 / h[d]);
+                    if k > 0.0 {
+                        wall_implicit += k;
+                        wall_rhs += k * un;
+                        mu * g_plus * a_d * (un - u0) / h[d]
+                    } else {
+                        (mu * g_plus * a_d / h[d] + k) * (un - u0)
+                    }
+                }
                 Some(un) => mu * g_plus * a_d * (un - u0) / h[d],
                 None => {
                     if sides[d][1] == Side::Velocity {
@@ -198,6 +226,22 @@ impl Solver {
                 }
             };
             diff -= match dn1 {
+                Some(ud) if solid(add(p, ed, -1)) => {
+                    let k = mu * g_minus * a_d / mask.exchange_delta(&cv, d);
+                    wall_implicit += k;
+                    wall_rhs += k * ud;
+                    0.0
+                }
+                Some(ud) if centroid => {
+                    let k = mu * g_minus * a_d * (1.0 / spacing(add(p, ed, -1), -1.0) - 1.0 / h[d]);
+                    if k > 0.0 {
+                        wall_implicit += k;
+                        wall_rhs += k * ud;
+                        mu * g_minus * a_d * (u0 - ud) / h[d]
+                    } else {
+                        (mu * g_minus * a_d / h[d] + k) * (u0 - ud)
+                    }
+                }
                 Some(ud) => mu * g_minus * a_d * (u0 - ud) / h[d],
                 None => {
                     if sides[d][0] == Side::Velocity {
@@ -218,7 +262,8 @@ impl Solver {
         let conv = rho * (conv - mass_out * u0);
         let p_plus = lat.cell(cell_plus).map_or(0.0, |ci| field.p[ci]);
         let p_minus = lat.cell(cell_minus).map_or(0.0, |ci| field.p[ci]);
-        let pressure = -(p_plus - p_minus) * cv.alpha * area[c];
+        let idx_f = lat.face(c, p).expect("the face");
+        let pressure = -(p_plus - p_minus) * cv.alpha * area[c] * mask.grad_weight(c, idx_f);
         let v_u = cv.alpha * h[c] * area[c];
         let source = self.momentum_source.as_ref().map_or(0.0, |f| {
             let s = f(x[0], x[1], x[2], t_old);
@@ -238,8 +283,9 @@ impl Solver {
         };
         let v_eff = fraction.max(INERTIA_FLOOR) * h[c] * area[c];
         let inertia = rho * v_eff / dt;
-        let u_star = (inertia * u0 - conv + diff + pressure + source + shear * ub - shear_explicit)
-            / (inertia + shear);
+        let u_star = (inertia * u0 - conv + diff + pressure + source + shear * ub - shear_explicit
+            + wall_rhs)
+            / (inertia + shear + wall_implicit);
         let v_alpha = fraction * h[c] * area[c];
         (u_star, rho * (v_eff - v_alpha) * (u_star - u0) / dt)
     }

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

@@ -163,6 +163,10 @@ impl DeviceStep {
     /// Allocates the device fields for `grid`; the momentum source is
     /// tabulated at `t = 0` (steady sources only in Stage 1).
     pub fn new(solver: Solver, grid: Grid) -> Self {
+        assert!(
+            !solver.params.pressure_centroid,
+            "the centroid pressure gradient (S2-5) is a host prototype: the device kernels do not carry it"
+        );
         let rt = runtime();
         let nu = (grid.nx + 1) * grid.ny * grid.nz;
         let nv = grid.nx * (grid.ny + 1) * grid.nz;
@@ -284,7 +288,11 @@ impl DeviceStep {
             bz0: side_code(b.z0),
             bz1: side_code(b.z1),
             scheme: scheme_code(self.solver.params.convection_scheme),
-            wall_order: i32::from(self.solver.params.wall_order),
+            // Low 4 bits: the shear closure's order; bit 4: the oblique wall distance.
+            wall_order: i32::from(self.solver.params.wall_order)
+                + 16 * i32::from(self.solver.params.wall_distance_oblique)
+                + 32 * i32::from(self.solver.params.wall_exchange_axis)
+                + 64 * i32::from(self.solver.params.diffusion_centroid),
             dx: g.dx,
             dy: g.dy,
             dz: g.dz,

crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/cut.rs

@@ -51,11 +51,11 @@ fn cut_kernels() -> &'static CutKernels {
     })
 }
 
-/// `struct E3Cut` in e3_cut.cu: 18 device pointers.
+/// `struct E3Cut` in e3_cut.cu: 21 device pointers.
 #[repr(C)]
 #[derive(Clone, Copy)]
 struct E3CutPtrs {
-    ptrs: [u64; 18],
+    ptrs: [u64; 21],
 }
 unsafe impl DeviceRepr for E3CutPtrs {}
 unsafe impl ValidAsZeroBits for E3CutPtrs {}
@@ -74,6 +74,8 @@ pub(super) struct DeviceCut {
     fold_ptr: CudaSlice<u32>,
     fold_idx: CudaSlice<u32>,
     cell_flux: CudaSlice<f64>,
+    /// The open-part centroid shifts per face (S2-5; one dummy entry when off).
+    shift: [CudaSlice<f64>; 3],
     pub(super) merged: usize,
 }
 
@@ -228,6 +230,10 @@ impl DeviceCut {
             fold_ptr: up_u(&fold_ptr),
             fold_idx: up_u(&fold_idx),
             cell_flux: rt.stream.alloc_zeros::<f64>(nc).expect("alloc"),
+            shift: match mask.face_shift_tables() {
+                Some(t) => [up_f(&t[0]), up_f(&t[1]), up_f(&t[2])],
+                None => [up_f(&[0.0]), up_f(&[0.0]), up_f(&[0.0])],
+            },
             merged,
         })
     }
@@ -258,6 +264,9 @@ impl DeviceCut {
                 pu(&self.fold_ptr),
                 pu(&self.fold_idx),
                 pf(&self.cell_flux),
+                pf(&self.shift[0]),
+                pf(&self.shift[1]),
+                pf(&self.shift[2]),
             ],
         }
     }

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

@@ -96,6 +96,33 @@ pub struct Parameters {
     /// next open face along the wall normal's dominant axis (S2-4).
     /// `Parameters::default()` reads `RTX_E3_WALL_ORDER` (default 1).
     pub wall_order: u8,
+    /// The cut face's wall distance with the wall's obliquity: `φ + ½h(1 −
+    /// α)|n_t|` (`n_t` the wall normal's part in the face plane) instead
+    /// of `φ + ½h(1 − α)`, which over-reads the distance on oblique walls
+    /// by `½h(1 − α)(1 − |n_t|)` at every h (S2-5). `Parameters::default()`
+    /// reads `RTX_E3_WALL_DISTANCE=oblique` (default: the recorded form).
+    pub wall_distance_oblique: bool,
+    /// The diffusive exchange of a fluid face with a SOLID neighbour face
+    /// over the axis distance to the wall, `δ = min(h, d_f/|n_d|)`, and
+    /// implicit — instead of the full `h`, which places the no-slip value
+    /// deeper in the body than the wall (an effective radius ≈ 0.28 h short
+    /// on the DFG 2D-1 ladder; S2-5). `RTX_E3_WALL_EXCHANGE=axis`.
+    pub wall_exchange_axis: bool,
+    /// HOST PROTOTYPE (S2-5): the pressure gradient across cut cells over
+    /// the distance between the cells' FLUID CENTROIDS (estimated from the
+    /// volume fraction and the wall normal) instead of `h` — a symmetric
+    /// per-face weight `ω = h/δ` in the Poisson coefficient, the velocity
+    /// correction and the predictor's pressure force. The device path
+    /// refuses it. `RTX_E3_PRESSURE_CENTROID=1`.
+    pub pressure_centroid: bool,
+    /// S2-5: the cross-direction diffusion between two faces over the
+    /// distance between their OPEN-PART CENTROIDS (a cut face's velocity
+    /// is its open part's mean, ½h(1 − α) off the face centre along the
+    /// wall's in-plane normal) instead of `h`. With `h` the coupling of a
+    /// cut face to its neighbour is weak by `(1 + α)/2` and the no-slip
+    /// surface sits `½(1 − α) h` inside the body (the flat-wall instrument).
+    /// `RTX_E3_DIFFUSION_CENTROID=1`.
+    pub diffusion_centroid: bool,
 }
 
 impl Default for Parameters {
@@ -116,6 +143,11 @@ impl Default for Parameters {
                 .ok()
                 .and_then(|v| v.parse().ok())
                 .unwrap_or(1),
+            wall_distance_oblique: std::env::var("RTX_E3_WALL_DISTANCE")
+                .is_ok_and(|v| v == "oblique"),
+            wall_exchange_axis: std::env::var("RTX_E3_WALL_EXCHANGE").is_ok_and(|v| v == "axis"),
+            pressure_centroid: std::env::var("RTX_E3_PRESSURE_CENTROID").is_ok_and(|v| v == "1"),
+            diffusion_centroid: std::env::var("RTX_E3_DIFFUSION_CENTROID").is_ok_and(|v| v == "1"),
         }
     }
 }
@@ -241,6 +273,15 @@ impl Solver {
             m.scheme = self.params.convection_scheme;
             m.density = self.fluid.density;
             m.wall_order = self.params.wall_order;
+            m.wall_distance_oblique = self.params.wall_distance_oblique;
+            m.wall_exchange_axis = self.params.wall_exchange_axis;
+            m.diffusion_centroid = self.params.diffusion_centroid;
+            if self.params.diffusion_centroid {
+                m.compute_face_shifts();
+            }
+            if self.params.pressure_centroid {
+                m.compute_gradient_weights();
+            }
             m
         })
         .expect("embedded mask")
@@ -360,6 +401,24 @@ impl Solver {
     // The projection's apertures: step-averaged on a moving cut wall
     // (1 without a cut geometry).
     #[inline]
+    /// The pressure-gradient weight `ω = h/δ` of a u / v / w face (1
+    /// without the centroid prototype).
+    pub(super) fn gu(&self, k: usize, j: usize, i: usize) -> f64 {
+        self.mask
+            .as_ref()
+            .map_or(1.0, |m| m.grad_weight(0, m.grid().uface(k, j, i)))
+    }
+    pub(super) fn gv(&self, k: usize, j: usize, i: usize) -> f64 {
+        self.mask
+            .as_ref()
+            .map_or(1.0, |m| m.grad_weight(1, m.grid().vface(k, j, i)))
+    }
+    pub(super) fn gw(&self, k: usize, j: usize, i: usize) -> f64 {
+        self.mask
+            .as_ref()
+            .map_or(1.0, |m| m.grad_weight(2, m.grid().wface(k, j, i)))
+    }
+
     pub(super) fn au(&self, k: usize, j: usize, i: usize) -> f64 {
         self.mask
             .as_ref()

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

@@ -39,28 +39,28 @@ impl Solver {
                             extra += ae_outlet;
                         }
                     } else if self.u_is_unknown(k, j, i + 1) {
-                        problem.ae[idx] = ae_interior * self.au(k, j, i + 1);
+                        problem.ae[idx] = ae_interior * self.au(k, j, i + 1) * self.gu(k, j, i + 1);
                     }
                     if i == 0 {
                         if b.x0 == outlet {
                             extra += ae_outlet;
                         }
                     } else if self.u_is_unknown(k, j, i) {
-                        problem.aw[idx] = ae_interior * self.au(k, j, i);
+                        problem.aw[idx] = ae_interior * self.au(k, j, i) * self.gu(k, j, i);
                     }
                     if j + 1 == ny {
                         if b.y1 == outlet {
                             extra += an_outlet;
                         }
                     } else if self.v_is_unknown(k, j + 1, i) {
-                        problem.an[idx] = an_interior * self.av(k, j + 1, i);
+                        problem.an[idx] = an_interior * self.av(k, j + 1, i) * self.gv(k, j + 1, i);
                     }
                     if j == 0 {
                         if b.y0 == outlet {
                             extra += an_outlet;
                         }
                     } else if self.v_is_unknown(k, j, i) {
-                        problem.as_[idx] = an_interior * self.av(k, j, i);
+                        problem.as_[idx] = an_interior * self.av(k, j, i) * self.gv(k, j, i);
                     }
                     if k + 1 == nz && !periodic {
                         if b.z1 == outlet {
@@ -434,7 +434,8 @@ impl Solver {
             for j in 0..ny {
                 for i in 1..nx {
                     if self.u_is_unknown(k, j, i) {
-                        let dp_dx = (pp[g.cell(k, j, i)] - pp[g.cell(k, j, i - 1)]) / dx;
+                        let dp_dx =
+                            self.gu(k, j, i) * (pp[g.cell(k, j, i)] - pp[g.cell(k, j, i - 1)]) / dx;
                         let f = g.uface(k, j, i);
                         field.u[f] = field.u_star[f] - (dt / rho) * dp_dx;
                     }
@@ -453,7 +454,8 @@ impl Solver {
             for i in 0..nx {
                 for j in 1..ny {
                     if self.v_is_unknown(k, j, i) {
-                        let dp_dy = (pp[g.cell(k, j, i)] - pp[g.cell(k, j - 1, i)]) / dy;
+                        let dp_dy =
+                            self.gv(k, j, i) * (pp[g.cell(k, j, i)] - pp[g.cell(k, j - 1, i)]) / dy;
                         let f = g.vface(k, j, i);
                         field.v[f] = field.v_star[f] - (dt / rho) * dp_dy;
                     }

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

@@ -99,6 +99,17 @@ pub struct Mask {
     pub(super) density: f64,
     /// The cut wall's shear closure order (S2-4).
     pub(super) wall_order: u8,
+    /// The oblique wall distance of the cut faces (S2-5).
+    pub(super) wall_distance_oblique: bool,
+    /// The axis-distance implicit wall exchange (S2-5).
+    pub(super) wall_exchange_axis: bool,
+    /// The centroid prototype's pressure-gradient weights per u / v / w face.
+    /// The centroid-distance cross diffusion (S2-5).
+    pub(super) diffusion_centroid: bool,
+    /// The open-part centroid shifts per u / v / w face, three components
+    /// interleaved (built with `diffusion_centroid`).
+    pub(super) face_shifts: Option<[Vec<f64>; 3]>,
+    pub(super) grad_weights: Option<(Vec<f64>, Vec<f64>, Vec<f64>)>,
 }
 
 /// The z lattice position of a query: the lower plane index, the upper
@@ -512,6 +523,11 @@ impl Mask {
             scheme: crate::solvers::incompressible::ConvectionScheme::Upwind,
             density: 1.0,
             wall_order: 1,
+            wall_distance_oblique: false,
+            wall_exchange_axis: false,
+            grad_weights: None,
+            diffusion_centroid: false,
+            face_shifts: None,
         })
     }
 
@@ -616,6 +632,15 @@ impl Mask {
     /// Fluid volume fraction of a cell (1 on the binary wall).
     #[inline]
     #[must_use]
+    pub fn grad_weight(&self, c: usize, idx: usize) -> f64 {
+        self.grad_weights.as_ref().map_or(1.0, |w| match c {
+            0 => w.0[idx],
+            1 => w.1[idx],
+            _ => w.2[idx],
+        })
+    }
+    #[inline]
+    #[must_use]
     pub fn vol(&self, idx: usize) -> f64 {
         self.cut.as_ref().map_or(1.0, |c| c.vol[idx])
     }

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

@@ -80,7 +80,11 @@ fn dfg_2d_1_on_the_device() {
     // 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),
+        rtx_cfd::solvers::incompressible::EmbeddedBody::circle(
+            CX,
+            CY,
+            0.5 * D + env_f("RTX_E3_DFG_DR", 0.0) * h,
+        ),
         lz,
     ));
     let g = Grid::cubic(nx, ny, nz, h);

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

@@ -0,0 +1,140 @@
+//! S2-5 host ladder: DFG 2D-1 (Re 20) on a periodic-z slab of 4 cells with
+//! every knob of the cut wall read from the environment — the instrument
+//! for host-only prototypes (`RTX_E3_PRESSURE_CENTROID=1`). Reference
+//! c_D 5.57954, c_L 0.010619, Δp 0.117520. `RTX_E3_DFG_NY` (default 62).
+use rtx_cfd::solvers::incompressible::ConvectionScheme;
+use rtx_cfd::solvers::incompressible::embedded3::{
+    Body, Boundaries, Field, Fluid, Grid, Parameters, Side, Solver, WallScheme,
+};
+
+const H: f64 = 0.41;
+const L: f64 = 2.2;
+const D: f64 = 0.1;
+const CX: f64 = 0.2;
+const CY: f64 = 0.2;
+const U_M: f64 = 0.3;
+const U_BAR: f64 = 2.0 / 3.0 * U_M;
+const RHO: f64 = 1.0;
+const NU: f64 = 1e-3;
+
+fn inflow(y: f64, _z: f64) -> f64 {
+    4.0 * U_M * y * (H - y) / (H * H)
+}
+
+#[test]
+#[ignore = "host DFG at ny 31 with the routes split (about half an hour)"]
+fn dfg_2d1_on_the_host() {
+    let ny: usize = std::env::var("RTX_E3_DFG_NY")
+        .ok()
+        .and_then(|v| v.parse().ok())
+        .unwrap_or(62);
+    let h = H / ny as f64;
+    let nx = (L / h).round() as usize;
+    let nz = 4;
+    let lz = nz as f64 * h;
+    let dt = (0.3 * h / U_M).min(0.5 * h * h / (6.0 * NU));
+    let mut solver = Solver::new(
+        Fluid {
+            density: RHO,
+            viscosity: RHO * NU,
+            reference_velocity: U_BAR,
+            reference_length: D,
+        },
+        Parameters {
+            corrector_steps: 2,
+            tolerance: 1e-8,
+            convection_scheme: ConvectionScheme::TvdVanAlbada,
+            wall_scheme: WallScheme::CutCell,
+            boundaries: Boundaries {
+                x1: Side::PressureOutlet,
+                z0: Side::Periodic,
+                z1: Side::Periodic,
+                ..Boundaries::default()
+            },
+            ..Parameters::default()
+        },
+    );
+    solver.set_boundary_velocity(|x, y, z, _t| {
+        if x <= 0.0 {
+            (inflow(y, z), 0.0, 0.0)
+        } else {
+            (0.0, 0.0, 0.0)
+        }
+    });
+    solver.set_body(Body::extruded(
+        rtx_cfd::solvers::incompressible::EmbeddedBody::circle(CX, CY, 0.5 * D),
+        lz,
+    ));
+    let g = Grid::cubic(nx, ny, nz, h);
+    let mut field = Field::new(g);
+    for k in 0..nz {
+        for j in 0..ny {
+            let u0 = inflow((j as f64 + 0.5) * h, (k as f64 + 0.5) * h);
+            for i in 0..=nx {
+                field.u[g.uface(k, j, i)] = u0;
+            }
+        }
+    }
+    solver.initialize(&mut field);
+    let coef = 2.0 / (RHO * U_BAR * U_BAR * D * lz);
+    let steps = (12.0 / dt).ceil() as usize;
+    let start = std::time::Instant::now();
+    let mut last = (0.0, 0.0);
+    for step in 0..steps {
+        let r = solver.advance(&mut field, dt);
+        if (step + 1) % (steps / 20).max(1) == 0 || step + 1 == steps {
+            let t = solver.time();
+            let mask = solver.mask().unwrap();
+            let body = solver.body().unwrap();
+            let mu = RHO * NU;
+            let (po, so) = mask.cut_wall_force_parts(body, &field, mu, t).unwrap();
+            let ex = mask
+                .cut_wall_exchange_force(body, &field, mu, RHO, t, None)
+                .unwrap();
+            let (pr, sr) = mask
+                .cut_wall_force_reconstructed_parts(body, &field, mu, t, None)
+                .unwrap();
+            let margin = 1.5 * D;
+            let ci = |x: f64| ((x / h).round() as usize).clamp(2, nx - 2);
+            let cj = |y: f64| ((y / h).round() as usize).clamp(2, ny - 2);
+            let bx = (
+                ci(CX - margin),
+                ci(CX + margin),
+                cj(CY - 0.15),
+                cj(CY + 0.15),
+                0,
+                nz,
+            );
+            let fcv = mask.control_volume_force_with_walls(&field, dt, RHO, mu, None, bx, false);
+            let cd = |f: [f64; 3]| coef * f[0];
+            println!(
+                "  t {t:7.3}: c_D operator {:.4} (p {:.4} + s {:.4} + exchange {:.4}) | reconstructed {:.4} (p {:.4} + s {:.4}) | box {:.4} c_L {:.5} Δp {:.5}; residual {:.1e} [{:.0} s]",
+                cd(po) + cd(so) + cd(ex),
+                cd(po),
+                cd(so),
+                cd(ex),
+                cd(pr) + cd(sr),
+                cd(pr),
+                cd(sr),
+                cd(fcv),
+                coef * (po[1] + so[1] + ex[1]),
+                mask.pressure_at(&field.p, CX - 0.5 * D, CY, 0.5 * lz)
+                    .unwrap_or(f64::NAN)
+                    - mask
+                        .pressure_at(&field.p, CX + 0.5 * D, CY, 0.5 * lz)
+                        .unwrap_or(f64::NAN),
+                r.final_residual,
+                start.elapsed().as_secs_f64()
+            );
+            let now = (cd(po) + cd(so) + cd(ex), cd(fcv));
+            if (now.0 - last.0).abs() < 1e-4 * now.0.abs()
+                && (now.1 - last.1).abs() < 1e-4 * now.1.abs()
+                && t > 2.0
+            {
+                println!("  settled");
+                break;
+            }
+            last = now;
+        }
+    }
+}

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

@@ -0,0 +1,102 @@
+//! S2-5 instrument: where does the cut wall sit? Poiseuille flow along z
+//! (periodic, body force `f`) between an EMBEDDED flat wall at `y = y_w`
+//! (cut at a chosen fraction θ of a cell) and the domain's top wall. The
+//! flow rate per unit width is `f (H − y_w)³ / (12 μ)`, so the measured
+//! rate gives the effective wall position `y_eff`; the offset
+//! `(y_eff − y_w)/h` must vanish at second order and is the number the
+//! DFG ladder reads as an effective radius ≈ 0.2 h short.
+use rtx_cfd::solvers::incompressible::ConvectionScheme;
+use rtx_cfd::solvers::incompressible::embedded3::{
+    Body, Boundaries, Field, Fluid, Grid, Parameters, Side, Solver, WallScheme,
+};
+
+const MU: f64 = 0.1;
+const F: f64 = 1.0;
+const HY: f64 = 1.0;
+
+fn offset(ny: usize, theta: f64) -> (f64, Vec<(f64, f64, f64)>) {
+    let h = HY / ny as f64;
+    let (nx, nz) = (3 * ny, 2);
+    let y_w = (ny as f64 / 4.0).floor() * h + theta * h;
+    let exact = move |y: f64| {
+        if y > y_w {
+            F / (2.0 * MU) * (y - y_w) * (HY - y)
+        } else {
+            0.0
+        }
+    };
+    let mut solver = Solver::new(
+        Fluid {
+            density: 1.0,
+            viscosity: MU,
+            reference_velocity: 1.0,
+            reference_length: 1.0,
+        },
+        Parameters {
+            corrector_steps: 2,
+            tolerance: 1e-10,
+            convection_scheme: ConvectionScheme::Upwind,
+            wall_scheme: WallScheme::CutCell,
+            boundaries: Boundaries {
+                z0: Side::Periodic,
+                z1: Side::Periodic,
+                ..Boundaries::default()
+            },
+            ..Parameters::default()
+        },
+    );
+    solver.set_boundary_velocity(move |_x, y, _z, _t| (0.0, 0.0, exact(y)));
+    solver.set_momentum_source(|_, _, _, _| (0.0, 0.0, F));
+    solver.set_body(Body::from_sdf(move |_x, y, _z, _t| y - y_w));
+    let g = Grid::cubic(nx, ny, nz, h);
+    let mut field = Field::new(g);
+    for k in 0..=nz {
+        for j in 0..ny {
+            for i in 0..nx {
+                if k < nz || true {
+                    let idx = g.wface(k.min(nz), j, i);
+                    field.w[idx] = exact((j as f64 + 0.5) * h);
+                }
+            }
+        }
+    }
+    solver.initialize(&mut field);
+    let dt = 0.5 * h * h / (6.0 * MU);
+    let steps = (3.0 / dt).ceil() as usize;
+    for _ in 0..steps {
+        solver.advance(&mut field, dt);
+    }
+    let mask = solver.mask().expect("mask");
+    let i = nx / 2;
+    let mut q = 0.0;
+    let mut profile = Vec::new();
+    for j in 0..ny {
+        let f = g.wface(0, j, i);
+        let a = mask.a_w(f);
+        q += a * field.w[f] * h;
+        let y = (j as f64 + 0.5) * h;
+        if a > 0.0 && profile.len() < 4 {
+            profile.push((a, field.w[f], exact(y)));
+        }
+    }
+    let y_eff = HY - (12.0 * MU * q / F).cbrt();
+    ((y_eff - y_w) / h, profile)
+}
+
+#[test]
+#[ignore = "S2-5 instrument: the cut wall's effective position on a flat wall (a minute on the host)"]
+fn flat_wall_effective_position() {
+    for ny in [16usize, 32] {
+        for theta in [0.05, 0.25, 0.5, 0.75, 0.95] {
+            let (off, profile) = offset(ny, theta);
+            let p: Vec<String> = profile
+                .iter()
+                .map(|(a, w, e)| format!("α {a:.2} w {w:.5} (exact at the face centre {e:.5})"))
+                .collect();
+            println!(
+                "  ny {ny} θ {theta:.2}: effective wall offset {off:+.4} h (positive = the wall sits inside the fluid); first open faces: {}",
+                p.join("; ")
+            );
+        }
+    }
+}