diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/closure.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/closure.rs
index ca39135..0ec824c 100644
--- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/closure.rs
+++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/closure.rs
@@ -160,7 +160,17 @@ impl Mask {
         let f = self.lattice().face(c, q).expect("open face");
         let h = [self.grid.dx, self.grid.dy, self.grid.dz];
         let d1 = cv.distance;
-        let d2 = d1 + h[d] * n[d].abs();
+        // S2-7: the second point's value lives at ITS OWN open-part centroid
+        // (its wall distance), not one lattice step along the axis from
+        // this face's centroid — the two differ by ½h(1 − α)n_t.
+        let d2 = if self.wall_order2_centroid {
+            self.cv_geometry(c, q).distance
+        } else {
+            d1 + h[d] * n[d].abs()
+        };
+        if d2 <= d1 {
+            return linear;
+        }
         (d2 / (d1 * (d2 - d1)), -d1 / (d2 * (d2 - d1)), Some(f))
     }
 }
diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cut.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cut.rs
index e85d64d..414c119 100644
--- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cut.rs
+++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cut.rs
@@ -90,29 +90,7 @@ impl CutGeometry {
         // (0, 0) to its (1, 1) corner in the face's own (a, b) order — the
         // Kuhn split's diagonals: for an x-face (y, z), a y-face (x, z), a
         // z-face (x, y); the same triangles seen from either cell.
-        let tri_area_fraction = |p0: f64, p1: f64, p2: f64| -> f64 {
-            let v = [p0, p1, p2];
-            let pos = v.iter().filter(|&&q| q >= 0.0).count();
-            match pos {
-                0 => 0.0,
-                3 => 1.0,
-                1 => {
-                    let a = v.iter().position(|&q| q >= 0.0).unwrap();
-                    let (b, c) = ((a + 1) % 3, (a + 2) % 3);
-                    (v[a] / (v[a] - v[b])) * (v[a] / (v[a] - v[c]))
-                }
-                _ => {
-                    let a = v.iter().position(|&q| q < 0.0).unwrap();
-                    let (b, c) = ((a + 1) % 3, (a + 2) % 3);
-                    1.0 - (v[a] / (v[a] - v[b])) * (v[a] / (v[a] - v[c]))
-                }
-            }
-        };
-        // Quad corners in (a, b) order: q00, q10, q01, q11; triangles
-        // (q00, q10, q11) and (q00, q11, q01).
-        let quad_fraction = |q00: f64, q10: f64, q01: f64, q11: f64| -> f64 {
-            0.5 * (tri_area_fraction(q00, q10, q11) + tri_area_fraction(q00, q11, q01))
-        };
+        let quad_fraction = quad_fraction;
         let mut a_u = vec![0.0; (nx + 1) * ny * nz];
         let mut a_v = vec![0.0; nx * (ny + 1) * nz];
         let mut a_w = vec![0.0; nx * ny * (nz + 1)];
@@ -221,6 +199,13 @@ impl CutGeometry {
         }
     }
 
+    /// φ at the corner `(k, j, i)` of the corner lattice.
+    #[inline]
+    #[must_use]
+    pub fn corner_phi(&self, k: usize, j: usize, i: usize) -> f64 {
+        self.phi[Self::node(self.grid, k, j, i)]
+    }
+
     /// Total fluid volume.
     #[must_use]
     pub fn fluid_volume(&self) -> f64 {
@@ -244,6 +229,34 @@ impl CutGeometry {
     }
 }
 
+/// The fluid fraction of a triangle from its three corner values of φ
+/// (the linear interpolant; fluid where φ ≥ 0).
+pub(super) fn tri_area_fraction(p0: f64, p1: f64, p2: f64) -> f64 {
+    let v = [p0, p1, p2];
+    let pos = v.iter().filter(|&&q| q >= 0.0).count();
+    match pos {
+        0 => 0.0,
+        3 => 1.0,
+        1 => {
+            let a = v.iter().position(|&q| q >= 0.0).unwrap();
+            let (b, c) = ((a + 1) % 3, (a + 2) % 3);
+            (v[a] / (v[a] - v[b])) * (v[a] / (v[a] - v[c]))
+        }
+        _ => {
+            let a = v.iter().position(|&q| q < 0.0).unwrap();
+            let (b, c) = ((a + 1) % 3, (a + 2) % 3);
+            1.0 - (v[a] / (v[a] - v[b])) * (v[a] / (v[a] - v[c]))
+        }
+    }
+}
+
+/// The fluid fraction of a quad from its corner values in (a, b) order
+/// (q00, q10, q01, q11): the two triangles along the (0, 0)–(1, 1)
+/// diagonal (the Kuhn split's).
+pub(super) fn quad_fraction(q00: f64, q10: f64, q01: f64, q11: f64) -> f64 {
+    0.5 * (tri_area_fraction(q00, q10, q11) + tri_area_fraction(q00, q11, q01))
+}
+
 fn det3(a: [f64; 3], b: [f64; 3], c: [f64; 3]) -> f64 {
     a[0] * (b[1] * c[2] - b[2] * c[1]) - a[1] * (b[0] * c[2] - b[2] * c[0])
         + a[2] * (b[0] * c[1] - b[1] * c[0])
diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs
index 544687c..0d66dd5 100644
--- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs
+++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cutwall.rs
@@ -219,6 +219,8 @@ impl Mask {
             wall_advancing: false,
             exchange_convection_off: false,
             wall_exchange_axis: false,
+            cv_sides_exact: false,
+            wall_order2_centroid: false,
             grad_weights: None,
             diffusion_centroid: false,
             face_shifts: None,
@@ -375,6 +377,12 @@ impl Mask {
                 ap[d] = [minus, plus];
             }
         }
+        // S2-7: the sides' own apertures instead of the whole-face averages.
+        if self.cv_sides_exact {
+            if let Some(exact) = self.cut.as_ref().and_then(|cut| self.exact_cv_sides(cut, c, p)) {
+                ap = exact;
+            }
+        }
         let mut wall = [0.0; 3];
         for d in 0..3 {
             wall[d] = -(ap[d][1] - ap[d][0]) * area[d];
@@ -415,6 +423,86 @@ impl Mask {
         }
     }
 
+    /// S2-7: the control volume's side apertures from the interpolant on the
+    /// sides' OWN corners. An unknown face's control volume is the tile
+    /// between the two adjacent cells' centres: across `c` its sides are
+    /// two HALF faces (the far half of `cell_minus`'s face, the near half
+    /// of `cell_plus`'s), in the own direction the two cells' centre
+    /// planes. φ is linear along every edge, so the mid-edge values are
+    /// exact for the interpolant; each half face / centre plane is a quad
+    /// through the faces' own `quad_fraction`. The averages of whole-face
+    /// apertures the default takes are wrong by O(1) wherever the wall
+    /// crosses a side (the in-plane momentum residual on oblique walls).
+    /// `None` at a domain side (the default stays).
+    fn exact_cv_sides(&self, cut: &CutGeometry, c: usize, p: [i64; 3]) -> Option<[[f64; 2]; 3]> {
+        let lat = self.lattice();
+        let g = self.grid;
+        let mut pm = p;
+        pm[c] -= 1;
+        let cells = [g.kji(lat.cell(pm)?), g.kji(lat.cell(p)?)];
+        // The corner of cell (k, j, i) at unit offsets `o = [di, dj, dk]`.
+        let corner = |cell: (usize, usize, usize), o: [usize; 3]| {
+            cut.corner_phi(cell.0 + o[2], cell.1 + o[1], cell.2 + o[0])
+        };
+        // The value at a cell's corner or, with `half`, at the mid-point of
+        // its edge along `c` (the interpolant's mean of the two corners).
+        let value = |cell: (usize, usize, usize), mut o: [usize; 3], half: bool| -> f64 {
+            if half {
+                o[c] = 0;
+                let a = corner(cell, o);
+                o[c] = 1;
+                0.5 * (a + corner(cell, o))
+            } else {
+                corner(cell, o)
+            }
+        };
+        let mut ap = [[1.0; 2]; 3];
+        for d in 0..3 {
+            if d == c {
+                // The two cells' centre planes across `c`: corners at the
+                // mid-points of the cells' `c` edges.
+                let (d1, d2) = ((c + 1) % 3, (c + 2) % 3);
+                for (side, cell) in cells.iter().enumerate() {
+                    let mid = |o1: usize, o2: usize| {
+                        let mut o = [0; 3];
+                        o[d1] = o1;
+                        o[d2] = o2;
+                        value(*cell, o, true)
+                    };
+                    ap[d][side] = super::cut::quad_fraction(mid(0, 0), mid(1, 0), mid(0, 1), mid(1, 1));
+                }
+            } else {
+                let e = 3 - c - d;
+                for side in 0..2 {
+                    // cell_minus's `d` face at `side`, its half nearer the
+                    // unknown face (c from ½ to 1); cell_plus's, c from 0 to ½.
+                    let half = |cell: (usize, usize, usize), far: bool| -> f64 {
+                        let at = |oc: u8, oe: usize| {
+                            let mut o = [0; 3];
+                            o[d] = side;
+                            o[e] = oe;
+                            match oc {
+                                0 => value(cell, o, false),
+                                1 => value(cell, o, true),
+                                _ => {
+                                    o[c] = 1;
+                                    value(cell, o, false)
+                                }
+                            }
+                        };
+                        if far {
+                            super::cut::quad_fraction(at(1, 0), at(2, 0), at(1, 1), at(2, 1))
+                        } else {
+                            super::cut::quad_fraction(at(0, 0), at(1, 0), at(0, 1), at(1, 1))
+                        }
+                    };
+                    ap[d][side] = 0.5 * (half(cells[0], true) + half(cells[1], false));
+                }
+            }
+        }
+        Some(ap)
+    }
+
     /// The surface velocity component `c` at the foot of the normal from
     /// the face centre `x`. With a cut geometry the signed distance and
     /// the normal come from the geometry's own corner values (the
diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device.rs
index 83e4ddf..4271931 100644
--- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device.rs
+++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device.rs
@@ -173,6 +173,14 @@ impl DeviceStep {
             !solver.params.momentum_volume_tiled,
             "the tiled momentum volume (A1-b) is a host prototype: the device kernels do not carry it"
         );
+        assert!(
+            !solver.params.wall_order2_centroid,
+            "the order-2 second point at the neighbour's centroid (S2-7) is a host prototype: the device kernels do not carry it"
+        );
+        assert!(
+            !solver.params.cv_sides_exact,
+            "the exact control-volume sides (S2-7) are 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;
diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs
index fee560f..789c115 100644
--- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs
+++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs
@@ -147,6 +147,19 @@ pub struct Parameters {
     /// layer — the sharpened S2-7 hypothesis. The device path refuses it.
     /// `RTX_E3_MOMENTUM_VOLUME=tiled`.
     pub momentum_volume_tiled: bool,
+    /// HOST PROTOTYPE (S2-7, 2026-09-19): the momentum control volumes' side
+    /// apertures from the interpolant on the sides' own corners (the half
+    /// faces across the face's direction, the cell-centre planes along it)
+    /// instead of the averages of whole-face apertures, which are wrong by
+    /// O(1) wherever the wall crosses a side. Static bodies; the device
+    /// refuses it. `RTX_E3_CV_SIDES=exact`.
+    pub cv_sides_exact: bool,
+    /// HOST PROTOTYPE (S2-7): the quadratic wall gradient's (order 2) second
+    /// point at the neighbour face's OWN centroid distance instead of one
+    /// lattice step along the dominant axis from this face's centroid (the
+    /// two differ by ½h(1 − α)n_t, O(1) of the step on cut faces). The
+    /// device refuses it. `RTX_E3_WALL_ORDER2=centroid`.
+    pub wall_order2_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
@@ -191,6 +204,8 @@ impl Default for Parameters {
             pressure_centroid: std::env::var("RTX_E3_PRESSURE_CENTROID").is_ok_and(|v| v == "1"),
             momentum_volume_tiled: std::env::var("RTX_E3_MOMENTUM_VOLUME")
                 .is_ok_and(|v| v == "tiled"),
+            cv_sides_exact: std::env::var("RTX_E3_CV_SIDES").is_ok_and(|v| v == "exact"),
+            wall_order2_centroid: std::env::var("RTX_E3_WALL_ORDER2").is_ok_and(|v| v == "centroid"),
             // ON by default since S2-5 (`=0` reproduces the records before it).
             diffusion_centroid: std::env::var("RTX_E3_DIFFUSION_CENTROID")
                 .map_or(true, |v| v != "0"),
@@ -330,6 +345,8 @@ impl Solver {
             m.distance_floor_fine = self.params.distance_floor_fine;
             m.wall_advancing = self.params.wall_advancing;
             m.exchange_convection_off = self.params.exchange_convection_off;
+            m.cv_sides_exact = self.params.cv_sides_exact;
+            m.wall_order2_centroid = self.params.wall_order2_centroid;
             m.diffusion_centroid = self.params.diffusion_centroid;
             if self.params.diffusion_centroid {
                 m.compute_face_shifts();
diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs
index bd7a30a..e47407f 100644
--- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs
+++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/wall.rs
@@ -108,6 +108,12 @@ pub struct Mask {
     pub(super) exchange_convection_off: bool,
     /// The axis-distance implicit wall exchange (S2-5).
     pub(super) wall_exchange_axis: bool,
+    /// S2-7: the momentum control volumes' side apertures from the
+    /// interpolant on the sides' own corners (host prototype).
+    pub(super) cv_sides_exact: bool,
+    /// S2-7: the quadratic wall gradient's second point at the neighbour's
+    /// own centroid distance (host prototype).
+    pub(super) wall_order2_centroid: 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,
@@ -534,6 +540,8 @@ impl Mask {
             wall_advancing: false,
             exchange_convection_off: false,
             wall_exchange_axis: false,
+            cv_sides_exact: false,
+            wall_order2_centroid: false,
             grad_weights: None,
             diffusion_centroid: false,
             face_shifts: None,
diff --git a/crates/specialized/rtx-cfd/tests/embedded3_wall_position_oblique.rs b/crates/specialized/rtx-cfd/tests/embedded3_wall_position_oblique.rs
index 3b766d4..4abe6b8 100644
--- a/crates/specialized/rtx-cfd/tests/embedded3_wall_position_oblique.rs
+++ b/crates/specialized/rtx-cfd/tests/embedded3_wall_position_oblique.rs
@@ -13,7 +13,7 @@
 //! printed).
 use rtx_cfd::solvers::incompressible::ConvectionScheme;
 use rtx_cfd::solvers::incompressible::embedded3::{
-    Body, Boundaries, Field, Fluid, Grid, Parameters, Side, Solver, WallScheme,
+    Body, Boundaries, FaceKind, Field, Fluid, Grid, Parameters, Side, Solver, WallScheme,
 };
 
 const MU: f64 = 0.1;
@@ -56,6 +56,9 @@ struct Reading {
     /// RMS of the pressure about its linear fit, over F·G: full cells, cut cells.
     p_full: f64,
     p_cut: f64,
+    /// Cut cells (fraction < 1) and virtually merged cells in the mask.
+    cut_cells: usize,
+    merged: usize,
 }
 
 /// Least squares of `y = a − b x`: returns (a, b).
@@ -161,6 +164,10 @@ fn reading(n: usize, slope: f64, c0: f64, along_z: bool, s26: Option<bool>) -> R
         solver.advance(&mut field, dt);
     }
     let mask = solver.mask().expect("mask");
+    let cut_cells = (0..g.cells())
+        .filter(|&c| mask.cell_active(c) && mask.vol(c) < 1.0)
+        .count();
+    let merged = mask.merged_cells();
     // The sides pin the flow RATE (the exact profile), so a displaced wall
     // appears as a streamwise pressure slope: F_eff = F − dp/ds, and on the
     // full faces of a central window `u/t_x = F_eff/(2μ) (G_eff²/4 − r²)`
@@ -231,6 +238,8 @@ fn reading(n: usize, slope: f64, c0: f64, along_z: bool, s26: Option<bool>) -> R
         force_p: 1.0 + minus_slope / F,
         p_full: rms(&pf),
         p_cut: rms(&pc),
+        cut_cells,
+        merged,
     }
 }
 
@@ -318,13 +327,31 @@ fn oblique_wall_linear_exactness() {
 #[test]
 #[ignore = "S2-6 instrument: the oblique cut wall's effective position with in-plane flow (minutes on the host)"]
 fn oblique_wall_effective_position() {
-    for (slope, c0) in [
+    let list = |name: &str| -> Option<Vec<f64>> {
+        std::env::var(name)
+            .ok()
+            .map(|v| v.split(',').filter_map(|t| t.trim().parse().ok()).collect())
+    };
+    // S2-7: `RTX_E3_OBLIQUE_SLOPES=0.25,0.5,1` restricts the slopes;
+    // `RTX_E3_OBLIQUE_C0_SHIFTS=0,0.25,0.5,0.75` sweeps the registration
+    // (c0 = 0.53 + shift · h at each n; a shift of 1 is the same
+    // registration); `RTX_E3_OBLIQUE_ZFLOW=0` skips the z-flow rows.
+    let slopes = list("RTX_E3_OBLIQUE_SLOPES");
+    let shifts = list("RTX_E3_OBLIQUE_C0_SHIFTS").unwrap_or_else(|| vec![0.0]);
+    let zflow = std::env::var("RTX_E3_OBLIQUE_ZFLOW").map_or(true, |v| v != "0");
+    for (slope, c0_base) in [
         (0.0, 0.53),
         (0.0, 0.77),
         (0.25, 0.53),
         (0.5, 0.53),
         (1.0, 0.53),
     ] {
+        if slopes
+            .as_ref()
+            .is_some_and(|l| !l.iter().any(|s| (s - slope).abs() < 1e-9))
+        {
+            continue;
+        }
         // `RTX_E3_OBLIQUE_N=64` adds a finer rung to the in-plane mode.
         let extra: Option<usize> = std::env::var("RTX_E3_OBLIQUE_N")
             .ok()
@@ -333,22 +360,401 @@ fn oblique_wall_effective_position() {
         if let Some(n) = extra {
             runs = vec![(n, false)];
         }
+        if !zflow {
+            runs.retain(|r| !r.1);
+        }
         for (n, along_z) in runs {
-            let r = reading(n, slope, c0, along_z, None);
-            println!(
-                "  slope {slope:.2} c0 {c0} n {n} {}: wall offset {:+.4} h (u faces) {:+.4} h (v faces), positive = inside the fluid; F_eff/F {:.5} (profile) {:.5} (pressure slope); pressure about its fit: {:.2e} full cells, {:.2e} near-wall cells (of F·G)",
-                if along_z {
-                    "z-flow (w faces)"
+            for &shift in &shifts {
+                let c0 = c0_base + shift / n as f64;
+                let r = reading(n, slope, c0, along_z, None);
+                println!(
+                    "  slope {slope:.2} c0 {c0:.5} n {n} {}: wall offset {:+.4} h (u faces) {:+.4} h (v faces), positive = inside the fluid; F_eff/F {:.5} (profile) {:.5} (pressure slope); pressure about its fit: {:.2e} full cells, {:.2e} near-wall cells (of F·G); cut cells {} merged {}",
+                    if along_z {
+                        "z-flow (w faces)"
+                    } else {
+                        "in-plane"
+                    },
+                    r.off_u,
+                    r.off_v,
+                    r.force_u,
+                    r.force_p,
+                    r.p_full,
+                    r.p_cut,
+                    r.cut_cells,
+                    r.merged
+                );
+            }
+        }
+    }
+}
+
+/// S2-7 Step 2 probe: the discrete operators applied to the EXACT in-plane
+/// Poiseuille field valued at the faces' open-part centroids — after one
+/// step with one corrector: the predictor's acceleration per near-wall
+/// face in units of the source's `F/ρ` (zero for a consistent momentum
+/// operator on the exact field), by aperture band; the cut cells'
+/// divergence of the exact centroid VALUES, of the exact open-part MEANS
+/// (the exact fluxes: must vanish), of the predicted field `u*` and of the
+/// corrected field, over the cell's largest face flux, by fluid-fraction
+/// class; and the spurious pressure one projection creates at the cut
+/// cells (of F·G). `u*` is recovered from the one correction.
+#[allow(clippy::too_many_lines)]
+fn probe(n: usize, slope: f64, c0: f64) {
+    #[allow(non_snake_case)]
+    let F: f64 = std::env::var("RTX_E3_OBLIQUE_F")
+        .ok()
+        .and_then(|v| v.parse().ok())
+        .unwrap_or(1e-4);
+    let h = 1.0 / n as f64;
+    let (nx, ny, nz) = ((LX * n as f64) as usize, (LY * n as f64) as usize, 2);
+    let norm = (1.0 + slope * slope).sqrt();
+    let (tx, ty) = (1.0 / norm, slope / norm);
+    let gap = W / norm;
+    let r_of = move |x: f64, y: f64| ((y - slope * x - c0) - 0.5 * W) / norm;
+    let speed = move |x: f64, y: f64| {
+        let r = r_of(x, y);
+        if r.abs() < 0.5 * gap {
+            F / (2.0 * MU) * (0.25 * gap * gap - r * r)
+        } else {
+            0.0
+        }
+    };
+    let mut params = parameters(None);
+    params.corrector_steps = 1;
+    let rho = 1.0;
+    let mut solver = Solver::new(
+        Fluid {
+            density: rho,
+            viscosity: MU,
+            reference_velocity: 1.0,
+            reference_length: 1.0,
+        },
+        params,
+    );
+    solver.set_boundary_velocity(move |x, y, _z, _t| {
+        let s = speed(x, y);
+        (s * tx, s * ty, 0.0)
+    });
+    solver.set_momentum_source(move |_, _, _, _| (F * tx, F * ty, 0.0));
+    solver.set_body(Body::from_sdf(move |x, y, _z, _t| 0.5 * gap - r_of(x, y).abs()));
+    let g = Grid::cubic(nx, ny, nz, h);
+    let mut field = Field::new(g);
+    solver.initialize(&mut field);
+    let (su, sv) = {
+        let mask = solver.mask().expect("mask");
+        let t = mask
+            .face_shift_tables()
+            .expect("the centroid shift tables (diffusion_centroid on)");
+        (t[0].clone(), t[1].clone())
+    };
+    // The exact field at the open parts' centroids (a full face: its centre).
+    for k in 0..nz {
+        for j in 0..ny {
+            for i in 0..=nx {
+                let f = g.uface(k, j, i);
+                field.u[f] = tx * speed(i as f64 * h + su[3 * f], (j as f64 + 0.5) * h + su[3 * f + 1]);
+            }
+        }
+        for j in 0..=ny {
+            for i in 0..nx {
+                let f = g.vface(k, j, i);
+                field.v[f] = ty * speed((i as f64 + 0.5) * h + sv[3 * f], j as f64 * h + sv[3 * f + 1]);
+            }
+        }
+    }
+    // The ghost faces (inside the body, in the imposition band) carry the
+    // solver's own reconstruction from the exact fluid field, as in a march.
+    {
+        let (body, mask) = (solver.body().expect("body"), solver.mask().expect("mask"));
+        mask.impose(body, &mut field.u, &mut field.v, &mut field.w, solver.time());
+    }
+    let dt = 0.5 * h * h / (6.0 * MU);
+    solver.advance(&mut field, dt);
+    let mask = solver.mask().expect("mask");
+    let pp = &field.p_prime;
+    // A full cell's fraction is not exactly 1 (the interpolant's volume).
+    let is_cut = |c: usize| mask.vol(c) < 1.0 - 1e-9;
+    let (mut vmin, mut vmax) = (f64::INFINITY, 0.0f64);
+    let window = |x: f64| (x - 0.5 * LX).abs() < 0.6;
+    // The predictor's value of a face: the corrected value plus the one
+    // correction (unknown faces), the face's value otherwise.
+    let star_u = |j: usize, i: usize| -> f64 {
+        let f = g.uface(0, j, i);
+        if i == 0 || i == nx || mask.u_kind(f) != FaceKind::Fluid {
+            return field.u[f];
+        }
+        let (cm, cp) = (g.cell(0, j, i - 1), g.cell(0, j, i));
+        field.u[f] + (dt / rho) * mask.grad_weight(0, f) * (pp[cp] - pp[cm]) / h
+    };
+    let star_v = |j: usize, i: usize| -> f64 {
+        let f = g.vface(0, j, i);
+        if j == 0 || j == ny || mask.v_kind(f) != FaceKind::Fluid {
+            return field.v[f];
+        }
+        let (cm, cp) = (g.cell(0, j - 1, i), g.cell(0, j, i));
+        field.v[f] + (dt / rho) * mask.grad_weight(1, f) * (pp[cp] - pp[cm]) / h
+    };
+    // (a) the predictor's acceleration on near-wall faces, of F/ρ.
+    let bin_of = |a: f64, near: bool| -> Option<usize> {
+        if a < 1.0 {
+            Some(((a * 4.0).floor() as usize).min(3))
+        } else if near {
+            Some(4)
+        } else {
+            None
+        }
+    };
+    let mut acc: [Vec<f64>; 5] = Default::default();
+    for j in 0..ny {
+        for i in 1..nx {
+            let f = g.uface(0, j, i);
+            if !window(i as f64 * h) || mask.u_kind(f) != FaceKind::Fluid {
+                continue;
+            }
+            let (cm, cp) = (g.cell(0, j, i - 1), g.cell(0, j, i));
+            let near = is_cut(cm) || is_cut(cp);
+            let Some(b) = bin_of(mask.a_u(f), near) else { continue };
+            acc[b].push(rho * (star_u(j, i) - field.u_old[f]) / dt / (F * tx));
+        }
+    }
+    if slope > 0.0 {
+        for j in 1..ny {
+            for i in 0..nx {
+                let f = g.vface(0, j, i);
+                if !window((i as f64 + 0.5) * h) || mask.v_kind(f) != FaceKind::Fluid {
+                    continue;
+                }
+                let (cm, cp) = (g.cell(0, j - 1, i), g.cell(0, j, i));
+                let near = is_cut(cm) || is_cut(cp);
+                let Some(b) = bin_of(mask.a_v(f), near) else { continue };
+                acc[b].push(rho * (star_v(j, i) - field.v_old[f]) / dt / (F * ty));
+            }
+        }
+    }
+    // (b) the cut cells' divergence under four valuations, of the cell's
+    // largest face flux; (c) the spurious pressure at cut cells.
+    let sigma_u = |a: f64| a * a * h * h / (12.0 * norm * norm);
+    let sigma_v = |a: f64| slope * slope * a * a * h * h / (12.0 * norm * norm);
+    let class_of = |vol: f64| -> usize {
+        if vol < 0.1 {
+            0
+        } else if vol < 0.5 {
+            1
+        } else if vol < 1.0 {
+            2
+        } else {
+            3
+        }
+    };
+    let mut div: [Vec<[f64; 4]>; 4] = Default::default();
+    let (mut p_cut, mut p_full) = (Vec::new(), Vec::new());
+    let area = h * h;
+    for j in 1..ny - 1 {
+        for i in 1..nx - 1 {
+            let c = g.cell(0, j, i);
+            if !window((i as f64 + 0.5) * h) || !mask.cell_active(c) {
+                continue;
+            }
+            let vol = mask.vol(c);
+            let near = is_cut(c)
+                || [g.cell(0, j, i - 1), g.cell(0, j, i + 1), g.cell(0, j - 1, i), g.cell(0, j + 1, i)]
+                    .iter()
+                    .any(|&q| is_cut(q));
+            if !near {
+                if (r_of((i as f64 + 0.5) * h, (j as f64 + 0.5) * h)).abs() < 0.5 * gap {
+                    p_full.push(pp[c]);
+                    vmin = vmin.min(vol);
+                    vmax = vmax.max(vol);
+                }
+                continue;
+            }
+            let vol = if is_cut(c) { vol } else { 1.0 };
+            let mut d = [0.0; 4];
+            let mut scale: f64 = 0.0;
+            for (comp, fi, fj, sign) in [(0usize, i + 1, j, 1.0), (0, i, j, -1.0), (1, i, j + 1, 1.0), (1, i, j, -1.0)] {
+                let (a, exact, cur, star, sig, t) = if comp == 0 {
+                    let f = g.uface(0, fj, fi);
+                    (mask.a_u(f), field.u_old[f], field.u[f], star_u(fj, fi), sigma_u(mask.a_u(f)), tx)
                 } else {
-                    "in-plane"
-                },
-                r.off_u,
-                r.off_v,
-                r.force_u,
-                r.force_p,
-                r.p_full,
-                r.p_cut
-            );
+                    let f = g.vface(0, fj, fi);
+                    (mask.a_v(f), field.v_old[f], field.v[f], star_v(fj, fi), sigma_v(mask.a_v(f)), ty)
+                };
+                if a <= 0.0 {
+                    continue;
+                }
+                let mean = exact - F / (2.0 * MU) * sig * t;
+                d[0] += sign * a * area * exact;
+                d[1] += sign * a * area * mean;
+                d[2] += sign * a * area * star;
+                d[3] += sign * a * area * cur;
+                scale = scale.max((a * area * exact).abs());
+            }
+            if scale > 0.0 {
+                div[class_of(vol)].push([d[0] / scale, d[1] / scale, d[2] / scale, d[3] / scale]);
+            }
+            if is_cut(c) {
+                p_cut.push(pp[c]);
+            }
+        }
+    }
+    let rms = |v: &[f64]| (v.iter().map(|x| x * x).sum::<f64>() / v.len().max(1) as f64).sqrt();
+    let maxabs = |v: &[f64]| v.iter().fold(0.0f64, |m, x| m.max(x.abs()));
+    let p0 = p_full.iter().sum::<f64>() / p_full.len().max(1) as f64;
+    let p_cut: Vec<f64> = p_cut.iter().map(|p| (p - p0) / (F * gap)).collect();
+    let p_full: Vec<f64> = p_full.iter().map(|p| (p - p0) / (F * gap)).collect();
+    println!(
+        "  probe slope {slope:.2} c0 {c0:.5} n {n}: merged {}; full cells' fraction {vmin:.3e}–{vmax:.3e}",
+        mask.merged_cells()
+    );
+    let names = ["α<¼", "¼–½", "½–¾", "¾–1", "full, next to a cut cell"];
+    for (b, name) in names.iter().enumerate() {
+        println!(
+            "    predictor acceleration of F/ρ on faces {name}: {} faces, rms {:.3e}, max {:.3e}",
+            acc[b].len(),
+            rms(&acc[b]),
+            maxabs(&acc[b])
+        );
+    }
+    let classes = ["vol<0.1 (merged class)", "0.1–0.5", "0.5–1", "full, next to a cut cell"];
+    for (k, name) in classes.iter().enumerate() {
+        let col = |m: usize| div[k].iter().map(|d| d[m]).collect::<Vec<_>>();
+        println!(
+            "    divergence / largest face flux, cells {name}: {} cells; exact centroid values rms {:.3e}, exact means rms {:.3e}, predicted u* rms {:.3e} max {:.3e}, corrected rms {:.3e}",
+            div[k].len(),
+            rms(&col(0)),
+            rms(&col(1)),
+            rms(&col(2)),
+            maxabs(&col(2)),
+            rms(&col(3))
+        );
+    }
+    println!(
+        "    spurious pressure after one projection (of F·G): cut cells rms {:.3e} max {:.3e}; full cells rms {:.3e}",
+        rms(&p_cut),
+        maxabs(&p_cut),
+        rms(&p_full)
+    );
+}
+
+/// The z-flow control of the Step 2 probe: the same channel with the exact
+/// Poiseuille `w(x, y)` at the w faces' open-part centroids (no pressure,
+/// no own-direction variation): the predictor's acceleration on the cut
+/// w faces in units of F/ρ, by aperture band.
+fn probe_z(n: usize, slope: f64, c0: f64) {
+    #[allow(non_snake_case)]
+    let F: f64 = std::env::var("RTX_E3_OBLIQUE_F")
+        .ok()
+        .and_then(|v| v.parse().ok())
+        .unwrap_or(1e-4);
+    let h = 1.0 / n as f64;
+    let (nx, ny, nz) = ((LX * n as f64) as usize, (LY * n as f64) as usize, 2);
+    let norm = (1.0 + slope * slope).sqrt();
+    let gap = W / norm;
+    let r_of = move |x: f64, y: f64| ((y - slope * x - c0) - 0.5 * W) / norm;
+    let speed = move |x: f64, y: f64| {
+        let r = r_of(x, y);
+        if r.abs() < 0.5 * gap {
+            F / (2.0 * MU) * (0.25 * gap * gap - r * r)
+        } else {
+            0.0
+        }
+    };
+    let mut params = parameters(None);
+    params.corrector_steps = 1;
+    let rho = 1.0;
+    let mut solver = Solver::new(
+        Fluid {
+            density: rho,
+            viscosity: MU,
+            reference_velocity: 1.0,
+            reference_length: 1.0,
+        },
+        params,
+    );
+    solver.set_boundary_velocity(move |x, y, _z, _t| (0.0, 0.0, speed(x, y)));
+    solver.set_momentum_source(move |_, _, _, _| (0.0, 0.0, F));
+    solver.set_body(Body::from_sdf(move |x, y, _z, _t| 0.5 * gap - r_of(x, y).abs()));
+    let g = Grid::cubic(nx, ny, nz, h);
+    let mut field = Field::new(g);
+    solver.initialize(&mut field);
+    let sw = solver
+        .mask()
+        .expect("mask")
+        .face_shift_tables()
+        .expect("shift tables")[2]
+        .clone();
+    for k in 0..=nz {
+        for j in 0..ny {
+            for i in 0..nx {
+                let f = g.wface(k, j, i);
+                field.w[f] = speed((i as f64 + 0.5) * h + sw[3 * f], (j as f64 + 0.5) * h + sw[3 * f + 1]);
+            }
+        }
+    }
+    {
+        let (body, mask) = (solver.body().expect("body"), solver.mask().expect("mask"));
+        mask.impose(body, &mut field.u, &mut field.v, &mut field.w, solver.time());
+    }
+    let dt = 0.5 * h * h / (6.0 * MU);
+    solver.advance(&mut field, dt);
+    let mask = solver.mask().expect("mask");
+    let window = |x: f64| (x - 0.5 * LX).abs() < 0.6;
+    let is_cut = |c: usize| mask.vol(c) < 1.0 - 1e-9;
+    let mut acc: [Vec<f64>; 5] = Default::default();
+    for j in 0..ny {
+        for i in 0..nx {
+            let f = g.wface(0, j, i);
+            if !window((i as f64 + 0.5) * h) || mask.w_kind(f) != FaceKind::Fluid {
+                continue;
+            }
+            let a = mask.a_w(f);
+            let c = g.cell(0, j, i);
+            let b = if a < 1.0 {
+                ((a * 4.0).floor() as usize).min(3)
+            } else if is_cut(c) {
+                4
+            } else {
+                continue;
+            };
+            // z-uniform: the projection leaves w alone (dp'/dz = 0).
+            acc[b].push(rho * (field.w[f] - field.w_old[f]) / dt / F);
+        }
+    }
+    let rms = |v: &[f64]| (v.iter().map(|x| x * x).sum::<f64>() / v.len().max(1) as f64).sqrt();
+    let maxabs = |v: &[f64]| v.iter().fold(0.0f64, |m, x| m.max(x.abs()));
+    println!("  probe-z slope {slope:.2} c0 {c0:.5} n {n}:");
+    let names = ["α<¼", "¼–½", "½–¾", "¾–1", "full, in a cut cell"];
+    for (b, name) in names.iter().enumerate() {
+        println!(
+            "    predictor acceleration of F/ρ on w faces {name}: {} faces, rms {:.3e}, max {:.3e}",
+            acc[b].len(),
+            rms(&acc[b]),
+            maxabs(&acc[b])
+        );
+    }
+}
+
+#[test]
+#[ignore = "S2-7 Step 2 probe: the discrete operators on the exact field (seconds per rung on the host)"]
+fn oblique_operator_probe() {
+    let list = |name: &str| -> Option<Vec<f64>> {
+        std::env::var(name)
+            .ok()
+            .map(|v| v.split(',').filter_map(|t| t.trim().parse().ok()).collect())
+    };
+    let slopes = list("RTX_E3_OBLIQUE_SLOPES").unwrap_or_else(|| vec![0.25, 0.5, 1.0]);
+    let shifts = list("RTX_E3_OBLIQUE_C0_SHIFTS").unwrap_or_else(|| vec![0.0]);
+    let ns: Vec<usize> = list("RTX_E3_OBLIQUE_NS")
+        .map(|l| l.iter().map(|&x| x as usize).collect())
+        .unwrap_or_else(|| vec![16, 32]);
+    for &slope in &slopes {
+        for &n in &ns {
+            for &shift in &shifts {
+                probe(n, slope, 0.53 + shift / n as f64);
+                if std::env::var("RTX_E3_OBLIQUE_ZFLOW").map_or(true, |v| v != "0") {
+                    probe_z(n, slope, 0.53 + shift / n as f64);
+                }
+            }
         }
     }
 }