embedded3 S2-7: oblique instrument gains registration sweep / slope subset / z-flow skip / merged-cell count (RTX_E3_OBLIQUE_{C0_SHIFTS,SLOPES,ZFLOW,NS}) and the operator probe (oblique_operator_probe: predictor acceleration on the exact centroid-valued field by aperture band, cut-cell divergence under four valuations, spurious pressure; z-flow control); two host prototypes, default off, device refuses them: RTX_E3_CV_SIDES=exact (the momentum CV's side apertures from the interpolant on the half faces / cell-centre planes) and RTX_E3_WALL_ORDER2=centroid (the quadratic wall gradient's second point at the neighbour's centroid distance); quad_fraction / tri_area_fraction lifted to module fns (default path digit-identical: oblique record, host suite 21/21, device cut tests)

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

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);
+                }
+            }
         }
     }
 }