//! S2-6 instrument: the cut wall's effective position on an OBLIQUE wall
//! with the flow IN the plane of the cut. Poiseuille flow (body force `F`
//! along the tangent) in a channel between two embedded parallel planes
//! `y = s x + c0` and `y = s x + c0 + w`; the sides carry the exact
//! solution, z is periodic. The exact field is `U(r) t`, `U = F/(2μ)
//! (G²/4 − r²)`, constant pressure, and the 5-point Laplacian is exact on
//! it, so on full faces of a central window `u/t_x + F r²/(2μ)` is the
//! constant `F G_eff²/(8μ)`: the effective gap from the u faces and from
//! the v faces separately, each an effective wall offset per wall in
//! units of h. Unlike the z-directed flat-wall test this one exercises the
//! own-direction coupling of cut faces, the convective terms' cut-face
//! values and the projection next to the wall (the spurious pressure is
//! printed).
use rtx_cfd::solvers::incompressible::ConvectionScheme;
use rtx_cfd::solvers::incompressible::embedded3::{
    Body, Boundaries, FaceKind, Field, Fluid, Grid, Parameters, Side, Solver, WallScheme,
};

const MU: f64 = 0.1;
const LX: f64 = 2.0;
const LY: f64 = 2.0;
const W: f64 = 0.5;

/// The cut wall's parameters: the environment's (`None`), or the S2-6
/// closures forced on / off for the gate.
fn parameters(s26: Option<bool>) -> Parameters {
    let mut p = 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()
    };
    if let Some(on) = s26 {
        p.diffusion_centroid = true;
        p.wall_distance_oblique = on;
        p.diffusion_transverse = on;
        p.distance_floor_fine = on;
    }
    p
}

struct Reading {
    /// Effective wall offset per wall from the fitted profile of the u / v faces.
    off_u: f64,
    off_v: f64,
    /// The driving force the profile's curvature implies, over F.
    force_u: f64,
    /// 1 − (the fitted streamwise pressure slope)/F: must equal `force_u`.
    force_p: f64,
    /// 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).
fn fit(points: &[(f64, f64)]) -> (f64, f64) {
    let n = points.len() as f64;
    let (sx, sy) = points
        .iter()
        .fold((0.0, 0.0), |s, p| (s.0 + p.0, s.1 + p.1));
    let (mx, my) = (sx / n, sy / n);
    let (sxx, sxy) = points.iter().fold((0.0, 0.0), |s, p| {
        (s.0 + (p.0 - mx) * (p.0 - mx), s.1 + (p.0 - mx) * (p.1 - my))
    });
    let slope = sxy / sxx;
    (my - slope * mx, -slope)
}

fn reading(n: usize, slope: f64, c0: f64, along_z: bool, s26: Option<bool>) -> Reading {
    // The driving force (`RTX_E3_OBLIQUE_F`): the problem is linear in it
    // but for the convective terms, so a small value switches them off.
    #[allow(non_snake_case)]
    let F: f64 = std::env::var("RTX_E3_OBLIQUE_F")
        .ok()
        .and_then(|v| v.parse().ok())
        .unwrap_or(1.0);
    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 solver = Solver::new(
        Fluid {
            density: 1.0,
            viscosity: MU,
            reference_velocity: 1.0,
            reference_length: 1.0,
        },
        parameters(s26),
    );
    // `along_z`: the same channel with the flow along the periodic z (the
    // cross-direction diffusion of w alone: no pressure, no convection).
    solver.set_boundary_velocity(move |x, y, _z, _t| {
        let s = speed(x, y);
        if along_z {
            (0.0, 0.0, s)
        } else {
            (s * tx, s * ty, 0.0)
        }
    });
    // `RTX_E3_OBLIQUE_DRIVE=pressure`: no body force — the sides' exact
    // profile pins the flow rate and the pressure slope alone drives it.
    let by_pressure = std::env::var("RTX_E3_OBLIQUE_DRIVE").is_ok_and(|v| v == "pressure");
    solver.set_momentum_source(move |_, _, _, _| {
        if by_pressure {
            (0.0, 0.0, 0.0)
        } else if along_z {
            (0.0, 0.0, F)
        } else {
            (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);
    for k in 0..=nz {
        for j in 0..ny {
            for i in 0..nx {
                if along_z {
                    field.w[g.wface(k, j, i)] = speed((i as f64 + 0.5) * h, (j as f64 + 0.5) * h);
                }
            }
        }
    }
    for k in 0..nz {
        if along_z {
            break;
        }
        for j in 0..ny {
            for i in 0..=nx {
                field.u[g.uface(k, j, i)] = tx * speed(i as f64 * h, (j as f64 + 0.5) * h);
            }
        }
        for j in 0..=ny {
            for i in 0..nx {
                field.v[g.vface(k, j, i)] = ty * speed((i as f64 + 0.5) * h, j as f64 * h);
            }
        }
    }
    solver.initialize(&mut field);
    let dt = 0.5 * h * h / (6.0 * MU);
    let steps = (2.0 / dt).ceil() as usize;
    for _ in 0..steps {
        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²)`
    // exactly (the 5-point Laplacian is exact on it). Fit both constants.
    let window = |x: f64| (x - 0.5 * LX).abs() < 0.3;
    let (mut pu, mut pv) = (Vec::new(), Vec::new());
    let (mut pf, mut pc) = (Vec::new(), Vec::new());
    for j in 0..ny {
        for i in 0..nx {
            let (xu, yu) = (i as f64 * h, (j as f64 + 0.5) * h);
            let fu = g.uface(0, j, i);
            let (xw, yw) = ((i as f64 + 0.5) * h, (j as f64 + 0.5) * h);
            let fw = g.wface(0, j, i);
            if along_z {
                if window(xw) && r_of(xw, yw).abs() < 0.4 * gap && mask.a_w(fw) >= 1.0 {
                    pu.push((r_of(xw, yw).powi(2), field.w[fw]));
                }
            } else if window(xu) && r_of(xu, yu).abs() < 0.4 * gap && mask.a_u(fu) >= 1.0 {
                pu.push((r_of(xu, yu).powi(2), field.u[fu] / tx));
            }
            let (xv, yv) = ((i as f64 + 0.5) * h, j as f64 * h);
            let fv = g.vface(0, j, i);
            if !along_z
                && slope > 0.0
                && window(xv)
                && r_of(xv, yv).abs() < 0.4 * gap
                && mask.a_v(fv) >= 1.0
            {
                pv.push((r_of(xv, yv).powi(2), field.v[fv] / ty));
            }
            let c = g.cell(0, j, i);
            let (xc, yc) = ((i as f64 + 0.5) * h, (j as f64 + 0.5) * h);
            if window(xc) && mask.cell_active(c) && r_of(xc, yc).abs() < 0.5 * gap + h {
                let s_along = xc * tx + yc * ty;
                let full = r_of(xc, yc).abs() < 0.5 * gap - 1.5 * h;
                if full {
                    pf.push((s_along, field.p[c]));
                } else {
                    pc.push((s_along, field.p[c]));
                }
            }
        }
    }
    let profile = |points: &[(f64, f64)]| {
        if points.is_empty() {
            return (f64::NAN, f64::NAN);
        }
        let (a, b) = fit(points);
        let gap_eff = 2.0 * (a / b).sqrt();
        (0.5 * (gap - gap_eff) / h, 2.0 * MU * b / F)
    };
    let (off_u, force_u) = profile(&pu);
    let (off_v, _) = profile(&pv);
    let (p0, minus_slope) = fit(&pf);
    let rms = |points: &[(f64, f64)]| {
        (points
            .iter()
            .map(|(s, p)| (p - (p0 - minus_slope * s)).powi(2))
            .sum::<f64>()
            / points.len().max(1) as f64)
            .sqrt()
            / (F * gap)
    };
    Reading {
        off_u,
        off_v,
        force_u,
        force_p: 1.0 + minus_slope / F,
        p_full: rms(&pf),
        p_cut: rms(&pc),
        cut_cells,
        merged,
    }
}

/// The linear-exactness mode: in-plane Couette flow `u = K dist t` over ONE
/// embedded oblique wall (no force, constant pressure, the sides carry the
/// exact field). A scheme exact on linear fields returns the wall position
/// to round-off; the fitted zero of the profile on full faces is the offset.
fn couette(n: usize, slope: f64, c0: f64, s26: Option<bool>) -> (f64, f64) {
    const K: f64 = 1.0;
    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 dist = move |x: f64, y: f64| (y - slope * x - c0) / norm;
    let speed = move |x: f64, y: f64| K * dist(x, y).max(0.0);
    let mut solver = Solver::new(
        Fluid {
            density: 1.0,
            viscosity: MU,
            reference_velocity: 1.0,
            reference_length: 1.0,
        },
        parameters(s26),
    );
    solver.set_boundary_velocity(move |x, y, _z, _t| {
        let s = speed(x, y);
        (s * tx, s * ty, 0.0)
    });
    solver.set_body(Body::from_sdf(move |x, y, _z, _t| dist(x, y)));
    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 {
                field.u[g.uface(k, j, i)] = tx * speed(i as f64 * h, (j as f64 + 0.5) * h);
            }
        }
        for j in 0..=ny {
            for i in 0..nx {
                field.v[g.vface(k, j, i)] = ty * speed((i as f64 + 0.5) * h, j as f64 * h);
            }
        }
    }
    solver.initialize(&mut field);
    let dt = 0.5 * h * h / (6.0 * MU);
    let steps = (2.0 / dt).ceil() as usize;
    for _ in 0..steps {
        solver.advance(&mut field, dt);
    }
    let mask = solver.mask().expect("mask");
    // u/t_x = K (dist − δ): fit on full faces of the central window, two to
    // six cells off the wall.
    let mut points = Vec::new();
    let mut worst: f64 = 0.0;
    for j in 0..ny {
        for i in 0..nx {
            let (x, y) = (i as f64 * h, (j as f64 + 0.5) * h);
            let f = g.uface(0, j, i);
            let dd = dist(x, y);
            if (x - 0.5 * LX).abs() < 0.3 && dd > 2.0 * h && dd < 6.0 * h && mask.a_u(f) >= 1.0 {
                points.push((dd, field.u[f] / tx));
                worst = worst.max((field.u[f] / tx - K * dd).abs() / (K * h));
            }
        }
    }
    let (a, minus_b) = fit(&points);
    // y = a − (−b) x with b the slope: the zero sits at dist = −a / b.
    let b = -minus_b;
    (-a / b / h, worst)
}

#[test]
#[ignore = "S2-6 instrument: linear exactness of the cut wall on an oblique wall (in-plane Couette; a minute on the host)"]
fn oblique_wall_linear_exactness() {
    for slope in [0.0, 0.25, 0.5, 1.0] {
        for n in [16usize, 32] {
            let (off, worst) = couette(n, slope, 0.53, None);
            println!(
                "  couette slope {slope:.2} n {n}: wall offset {off:+.4} h (negative = inside the body); worst full-face error {worst:.4} of K·h"
            );
        }
    }
}

#[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() {
    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()
            .and_then(|v| v.parse().ok());
        let mut runs = vec![(16usize, false), (32, false), (16, true), (32, true)];
        if let Some(n) = extra {
            runs = vec![(n, false)];
        }
        if !zflow {
            runs.retain(|r| !r.1);
        }
        for (n, along_z) in runs {
            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 {
                    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);
                }
            }
        }
    }
}

/// The gate (S2-6): with the oblique distance, the transverse centroid
/// correction and the fine floor the cut wall is linear-exact to 0.02 h on
/// an oblique wall (without them it sits 0.05–0.07 h inside the body at
/// every h), and the z-directed Poiseuille offset halves per rung at
/// slope ½ (without them: −0.087 → −0.083 h).
#[test]
fn oblique_wall_position_is_second_order() {
    for slope in [0.5, 1.0] {
        let (fixed, _) = couette(32, slope, 0.53, Some(true));
        let (before, _) = couette(32, slope, 0.53, Some(false));
        println!(
            "  couette slope {slope}: offset {fixed:+.4} h with the S2-6 closures, {before:+.4} h without"
        );
        assert!(fixed.abs() < 0.02, "slope {slope}: {fixed}");
        assert!(
            before.abs() > 2.0 * fixed.abs(),
            "the instrument lost its contrast"
        );
    }
    let coarse = reading(16, 0.5, 0.53, true, Some(true)).off_u;
    let fine = reading(32, 0.5, 0.53, true, Some(true)).off_u;
    println!("  z-flow slope 0.5: offset {coarse:+.4} h at n 16, {fine:+.4} h at n 32");
    assert!(coarse.abs() < 0.04, "n 16 offset {coarse}");
    assert!(
        fine.abs() < 0.65 * coarse.abs(),
        "the offset does not halve: {coarse} → {fine}"
    );
}