//! R7-2a: the composite MAC projection on a coarse grid with one nested
//! ratio-2 patch (`embedded3::composite::MacProjection`, a host prototype
//! nothing else calls).
//!
//! Gates (thresholds asserted, registered before the runs):
//!
//! - P1 (`projection_random`): a random field projected — `D G p = b − A p`
//!   to ≤ 1e-12 (relative to max |A p|) for a random `p`; after the
//!   projection max |D u| ≤ 1e-11 × max |D u*| in every class of cell (coarse
//!   interface, fine interface, the rest); the coarse interface fluxes equal
//!   the sums of their fine fluxes to ≤ 1e-13; the projection idempotent
//!   (max |P P u − P u| ≤ 1e-10 × max |P u|);
//! - P2 (`projection_mms_ladder`, ignored): `u* = u_div + ∇q` on the unit
//!   cube, patch = the middle half, coarse n = 16/32/64: the projected
//!   velocity against `u_div` at the face centres, split into fine interface
//!   faces, coarse interface faces, fine interior faces, coarse faces and the
//!   domain boundary faces; L2 and L∞ orders ≥ 1.8 on the last pair for the
//!   four non-boundary classes (Quadratic interface; Octree / Direct
//!   reported);
//! - P3 (`projection_cut_sphere`, ignored): a Neumann sphere cut into the
//!   patch (apertures on the fine faces): the div gate of P1, a translating
//!   body (uniform `u*` + its wall flux) left unchanged, the orders of P2
//!   away from the body (interface classes ≥ 1.8 in L2), and the composite
//!   against the uniformly fine projection with the same cut;
//! - P4 (`projection_cost`, ignored): wall time and unknowns against the
//!   uniformly fine MAC projection solved by embedded3's production PCG.
//!
//! CSVs go to `$R7_OUT` when set.

use rtx_cfd::solvers::incompressible::MultigridParameters;
use rtx_cfd::solvers::incompressible::embedded3::Grid;
use rtx_cfd::solvers::incompressible::embedded3::composite::{
    CompositeSpec, FaceClass, Interface, MacField, MacProjection, Sdf, uniform_problem,
};
use rtx_cfd::solvers::incompressible::embedded3::poisson::{Problem, solve_pcg};
use std::f64::consts::PI;
use std::io::Write;
use std::sync::Arc;

type Field3 = fn(usize, [f64; 3]) -> f64;
type Scalar = fn(f64, f64, f64) -> f64;

fn zero(_x: f64, _y: f64, _z: f64) -> f64 {
    0.0
}

fn out_file(name: &str) -> Option<std::fs::File> {
    let dir = std::env::var("R7_OUT").ok()?;
    std::fs::create_dir_all(&dir).ok()?;
    std::fs::File::create(format!("{dir}/{name}")).ok()
}

fn maxabs(v: &[f64]) -> f64 {
    v.iter().fold(0.0f64, |m, x| m.max(x.abs()))
}

struct Lcg(u64);

impl Lcg {
    fn next(&mut self) -> f64 {
        self.0 = self
            .0
            .wrapping_mul(6_364_136_223_846_793_005)
            .wrapping_add(1_442_695_040_888_963_407);
        ((self.0 >> 11) as f64 / (1u64 << 53) as f64) * 2.0 - 1.0
    }
}

fn projection(
    g: Grid,
    lo: [usize; 3],
    hi: [usize; 3],
    iface: Interface,
    body: Option<Sdf>,
    dirichlet: Scalar,
) -> MacProjection {
    let spec = CompositeSpec {
        coarse: g,
        lo,
        hi,
        interface: iface,
        body,
        source: &zero,
        dirichlet: &dirichlet,
    };
    MacProjection::new(&spec, 2)
}

// ---------------------------------------------------------------- P1

struct P1 {
    identity: f64,
    div_before: f64,
    /// max |D u| after: coarse interface cells, fine interface cells, rest.
    div_after: [f64; 3],
    conservation: f64,
    idempotence: f64,
    iterations: usize,
    rel: f64,
}

fn random_field(mp: &MacProjection, rng: &mut Lcg) -> MacField {
    let mut u = mp.zeros();
    for d in 0..3 {
        u.coarse[d].iter_mut().for_each(|v| *v = rng.next());
        u.fine[d].iter_mut().for_each(|v| *v = rng.next());
    }
    u
}

fn class_max(mp: &MacProjection, div: &[f64]) -> [f64; 3] {
    let mut m = [0.0f64; 3];
    for (r, v) in div.iter().enumerate() {
        let c = if !mp.comp.at_interface[r] {
            2
        } else if r < mp.comp.n_coarse {
            0
        } else {
            1
        };
        m[c] = m[c].max(v.abs());
    }
    m
}

/// max |coarse flux − Σ fine fluxes| over the coarse interface faces,
/// relative to the largest such flux.
fn conservation(mp: &MacProjection, u: &MacField) -> f64 {
    let hc2 = mp.comp.coarse.dx * mp.comp.coarse.dx;
    let mut sums: [std::collections::HashMap<usize, f64>; 3] = Default::default();
    for f in &mp.comp.iface {
        *sums[f.axis].entry(f.coarse_face).or_insert(0.0) +=
            mp.fine_area[f.axis][f.face] * u.fine[f.axis][f.face];
    }
    let (mut d, mut m) = (0.0f64, 0.0f64);
    for a in 0..3 {
        for (&cf, &s) in &sums[a] {
            assert_eq!(mp.coarse_class[a][cf], FaceClass::CoarseInterface);
            d = d.max((hc2 * u.coarse[a][cf] - s).abs());
            m = m.max(s.abs());
        }
    }
    d / m.max(f64::MIN_POSITIVE)
}

fn run_p1(g: Grid, lo: [usize; 3], hi: [usize; 3], body: Option<Sdf>, seed: u64) -> P1 {
    let mut mp = projection(g, lo, hi, Interface::Quadratic, body, zero);
    let n = mp.comp.unknowns();
    let mut rng = Lcg(seed);
    // D G p = b − A p for a random p.
    let p: Vec<f64> = (0..n).map(|_| rng.next()).collect();
    let gp = mp.gradient(&p);
    let mut dgp = vec![0.0; n];
    mp.divergence(&gp, None, &mut dgp);
    let mut ap = vec![0.0; n];
    mp.comp.a.apply(&p, &mut ap);
    let b = mp.boundary_rhs();
    let mut id = 0.0f64;
    for r in 0..n {
        id = id.max((dgp[r] - (b[r] - ap[r])).abs());
    }
    let identity = id / maxabs(&ap);
    // Project a random field, twice.
    let mut u = random_field(&mp, &mut rng);
    let mut p1 = vec![0.0; n];
    let st = mp.project(&mut u, &mut p1, None, 1e-13, 400);
    let mut div = vec![0.0; n];
    mp.divergence(&u, None, &mut div);
    let div_after = class_max(&mp, &div);
    let cons = conservation(&mp, &u);
    let mut u2 = u.clone();
    let mut p2 = vec![0.0; n];
    let _ = mp.project(&mut u2, &mut p2, None, 1e-13, 400);
    let (mut dd, mut mm) = (0.0f64, 0.0f64);
    for d in 0..3 {
        for (a, b) in u.coarse[d].iter().zip(&u2.coarse[d]) {
            dd = dd.max((a - b).abs());
            mm = mm.max(a.abs());
        }
        for (i, (a, b)) in u.fine[d].iter().zip(&u2.fine[d]).enumerate() {
            if mp.fine_class[d][i] != FaceClass::Closed {
                dd = dd.max((a - b).abs());
                mm = mm.max(a.abs());
            }
        }
    }
    P1 {
        identity,
        div_before: st.div_before,
        div_after,
        conservation: cons,
        idempotence: dd / mm,
        iterations: st.iterations,
        rel: st.rel_residual,
    }
}

fn check_p1(label: &str, r: &P1, csv: &mut Option<std::fs::File>) {
    let rel = r.div_after.map(|v| v / r.div_before);
    eprintln!(
        "P1 {label}: DG=b-A {:.2e}; div before {:.3e}, after/before [c-iface {:.2e}, f-iface {:.2e}, rest {:.2e}]; conservation {:.2e}; idempotence {:.2e}; {} it rel {:.2e}",
        r.identity,
        r.div_before,
        rel[0],
        rel[1],
        rel[2],
        r.conservation,
        r.idempotence,
        r.iterations,
        r.rel
    );
    if let Some(f) = csv.as_mut() {
        writeln!(
            f,
            "{label},{:.3e},{:.3e},{:.3e},{:.3e},{:.3e},{:.3e},{:.3e},{},{:.3e}",
            r.identity,
            r.div_before,
            rel[0],
            rel[1],
            rel[2],
            r.conservation,
            r.idempotence,
            r.iterations,
            r.rel
        )
        .ok();
    }
    assert!(
        r.identity <= 1e-12,
        "{label}: D G != b - A ({:.2e})",
        r.identity
    );
    for (c, v) in rel.iter().enumerate() {
        assert!(*v <= 1e-11, "{label}: class {c} div {v:.2e}");
    }
    assert!(
        r.conservation <= 1e-13,
        "{label}: interface not conservative"
    );
    assert!(
        r.idempotence <= 1e-10,
        "{label}: not idempotent ({:.2e})",
        r.idempotence
    );
}

#[test]
fn projection_smoke() {
    let g = Grid::cubic(12, 12, 12, 1.0 / 12.0);
    let r = run_p1(g, [3; 3], [9; 3], None, 7);
    check_p1("smoke n12", &r, &mut None);
}

#[test]
fn projection_random() {
    let mut csv = out_file("p1_random.csv");
    if let Some(f) = csv.as_mut() {
        writeln!(
            f,
            "case,dg_identity,div_before,div_c_iface_rel,div_f_iface_rel,div_rest_rel,conservation,idempotence,iters,rel_res"
        )
        .ok();
    }
    let cases: [(&str, Grid, [usize; 3], [usize; 3], Option<Sdf>); 3] = [
        (
            "cube24_middle",
            Grid::cubic(24, 24, 24, 1.0 / 24.0),
            [6; 3],
            [18; 3],
            None,
        ),
        (
            "box20x16x12_offcentre",
            Grid::cubic(20, 16, 12, 1.0 / 20.0),
            [4, 3, 3],
            [13, 11, 8],
            None,
        ),
        (
            "cube24_sphere",
            Grid::cubic(24, 24, 24, 1.0 / 24.0),
            [6; 3],
            [18; 3],
            Some(sphere()),
        ),
    ];
    for (label, g, lo, hi, body) in cases {
        let r = run_p1(g, lo, hi, body, 12345);
        check_p1(label, &r, &mut csv);
    }
}

// ---------------------------------------------------------------- P2 / P3

/// `u_div`: each component independent of its own coordinate.
fn udiv(d: usize, x: [f64; 3]) -> f64 {
    match d {
        0 => (1.3 * PI * x[1] + 0.2).cos() * (0.9 * PI * x[2] + 0.6).cos(),
        1 => (1.1 * PI * x[2] + 0.4).cos() * (1.2 * PI * x[0] + 0.1).cos(),
        _ => (0.8 * PI * x[0] + 0.3).cos() * (1.4 * PI * x[1] + 0.1).cos(),
    }
}

/// `q = ½ sin(πx) sin(2πy) sin(πz)`: q = 0 and ∂²q/∂n² = 0 on the cube.
fn q_grad(d: usize, x: [f64; 3]) -> f64 {
    let (a, b, c) = (PI * x[0], 2.0 * PI * x[1], PI * x[2]);
    0.5 * match d {
        0 => PI * a.cos() * b.sin() * c.sin(),
        1 => 2.0 * PI * a.sin() * b.cos() * c.sin(),
        _ => PI * a.sin() * b.sin() * c.cos(),
    }
}

fn ustar_mms(d: usize, x: [f64; 3]) -> f64 {
    udiv(d, x) + q_grad(d, x)
}

/// A localized gradient part inside the patch (the case local refinement
/// is for): `q + exp(−r²/s²)` about the cube's centre, s = 0.06.
const BUMP_S: f64 = 0.06;

fn q_bump(x: f64, y: f64, z: f64) -> f64 {
    let r2 = (x - 0.5).powi(2) + (y - 0.5).powi(2) + (z - 0.5).powi(2);
    (-r2 / (BUMP_S * BUMP_S)).exp()
}

fn ustar_bump(d: usize, x: [f64; 3]) -> f64 {
    let g = -2.0 * (x[d] - 0.5) / (BUMP_S * BUMP_S) * q_bump(x[0], x[1], x[2]);
    ustar_mms(d, x) + g
}

const SPH_C: [f64; 3] = [0.52, 0.49, 0.51];
const SPH_R: f64 = 0.12;
const SPH_A: f64 = 2.0 * PI;

fn sphere() -> Sdf {
    Arc::new(|x, y, z| {
        ((x - SPH_C[0]).powi(2) + (y - SPH_C[1]).powi(2) + (z - SPH_C[2]).powi(2)).sqrt() - SPH_R
    })
}

/// Potential flow past the sphere (U = 1 along x) plus a rigid swirl about
/// the z axis through its centre: divergence-free, u·n = 0 on the sphere.
fn udiv_sph(d: usize, x: [f64; 3]) -> f64 {
    let r = [x[0] - SPH_C[0], x[1] - SPH_C[1], x[2] - SPH_C[2]];
    let r2 = r[0] * r[0] + r[1] * r[1] + r[2] * r[2];
    let rr = r2.sqrt();
    let r5 = r2 * r2 * rr;
    let k = 1.5 * SPH_R.powi(3);
    let pot = match d {
        0 => 1.0 + 0.5 * SPH_R.powi(3) / (r2 * rr) - k * r[0] * r[0] / r5,
        _ => -k * r[0] * r[d] / r5,
    };
    let swirl = match d {
        0 => -0.7 * r[1],
        1 => 0.7 * r[0],
        _ => 0.0,
    };
    pot + swirl
}

/// `q = cos(a (r − R))`: ∂q/∂n = 0 on the sphere; Dirichlet = q outside.
fn q_sph(x: f64, y: f64, z: f64) -> f64 {
    let r = ((x - SPH_C[0]).powi(2) + (y - SPH_C[1]).powi(2) + (z - SPH_C[2]).powi(2)).sqrt();
    (SPH_A * (r - SPH_R)).cos()
}

fn q_sph_grad(d: usize, x: [f64; 3]) -> f64 {
    let r = [x[0] - SPH_C[0], x[1] - SPH_C[1], x[2] - SPH_C[2]];
    let rr = (r[0] * r[0] + r[1] * r[1] + r[2] * r[2]).sqrt();
    -SPH_A * (SPH_A * (rr - SPH_R)).sin() * r[d] / rr
}

fn ustar_sph(d: usize, x: [f64; 3]) -> f64 {
    udiv_sph(d, x) + q_sph_grad(d, x)
}

#[derive(Default, Clone, Copy)]
struct Err {
    s2: f64,
    vol: f64,
    max: f64,
    count: usize,
}

impl Err {
    fn add(&mut self, e: f64, v: f64) {
        self.s2 += e * e * v;
        self.vol += v;
        self.max = self.max.max(e.abs());
        self.count += 1;
    }
    fn l2(&self) -> f64 {
        (self.s2 / self.vol.max(f64::MIN_POSITIVE)).sqrt()
    }
}

const CLASSES: [&str; 6] = ["f_iface", "c_iface", "fine", "coarse", "boundary", "cut"];

/// Face errors against `exact` by class (see `CLASSES`); a fine face is
/// `cut` when it is partly closed or touches a cut cell.
fn face_errors(mp: &MacProjection, u: &MacField, exact: Field3) -> [Err; 6] {
    let mut e = [Err::default(); 6];
    let (hc, hf) = (mp.comp.coarse.dx, mp.comp.fine.dx);
    for d in 0..3 {
        for (i, v) in u.coarse[d].iter().enumerate() {
            let class = match mp.coarse_class[d][i] {
                FaceClass::CoarseInterface => 1,
                FaceClass::Coarse => 3,
                FaceClass::Boundary => 4,
                _ => continue,
            };
            let x = mp.face_centre(false, d, i);
            e[class].add(v - exact(d, x), hc * hc * hc);
        }
        for (i, v) in u.fine[d].iter().enumerate() {
            let class = match mp.fine_class[d][i] {
                FaceClass::FineInterface => 0,
                FaceClass::Fine => {
                    let cut_cell = mp
                        .fine_face_cells(d, i)
                        .iter()
                        .flatten()
                        .any(|&c| mp.comp.cut[c]);
                    if cut_cell || mp.fine_area[d][i] < hf * hf * (1.0 - 1e-12) {
                        5
                    } else {
                        2
                    }
                }
                _ => continue,
            };
            let x = mp.face_centre(true, d, i);
            e[class].add(v - exact(d, x), mp.fine_area[d][i] * hf);
        }
    }
    e
}

fn order(a: f64, b: f64) -> f64 {
    (a / b).log2()
}

struct Rung {
    n: usize,
    err: [Err; 6],
    div_rel: f64,
    iterations: usize,
    unknowns: usize,
}

fn mms_rung(
    n: usize,
    iface: Interface,
    body: Option<Sdf>,
    ustar: Field3,
    exact: Field3,
    dirichlet: Scalar,
) -> (Rung, MacProjection, MacField) {
    let g = Grid::cubic(n, n, n, 1.0 / n as f64);
    let mut mp = projection(g, [n / 4; 3], [3 * n / 4; 3], iface, body, dirichlet);
    let mut u = mp.sample(&ustar);
    let mut p = vec![0.0; mp.comp.unknowns()];
    let st = mp.project(&mut u, &mut p, None, 1e-13, 400);
    assert!(st.converged, "n {n}: {st:?}");
    let err = face_errors(&mp, &u, exact);
    let r = Rung {
        n,
        err,
        div_rel: st.div_after / st.div_before,
        iterations: st.iterations,
        unknowns: mp.comp.unknowns(),
    };
    (r, mp, u)
}

fn ladder_report(
    name: &str,
    rungs: &[Rung],
    csv: &mut Option<std::fs::File>,
) -> Vec<[(f64, f64); 6]> {
    let mut orders = Vec::new();
    for (w, r) in rungs.iter().enumerate() {
        for (c, label) in CLASSES.iter().enumerate() {
            let e = r.err[c];
            if e.count == 0 {
                continue;
            }
            let (o2, oi) = if w > 0 {
                let p = rungs[w - 1].err[c];
                (order(p.l2(), e.l2()), order(p.max, e.max))
            } else {
                (f64::NAN, f64::NAN)
            };
            eprintln!(
                "{name} n {:3} {label:9} faces {:8} L2 {:.3e} Linf {:.3e}  orders {o2:5.2} {oi:5.2}",
                r.n,
                e.count,
                e.l2(),
                e.max
            );
            if let Some(f) = csv.as_mut() {
                writeln!(
                    f,
                    "{name},{},{label},{},{:.6e},{:.6e},{o2:.3},{oi:.3},{:.3e},{},{}",
                    r.n,
                    e.count,
                    e.l2(),
                    e.max,
                    r.div_rel,
                    r.iterations,
                    r.unknowns
                )
                .ok();
            }
        }
        if w > 0 {
            let mut o = [(f64::NAN, f64::NAN); 6];
            for c in 0..6 {
                let (p, e) = (rungs[w - 1].err[c], r.err[c]);
                if e.count > 0 && p.count > 0 {
                    o[c] = (order(p.l2(), e.l2()), order(p.max, e.max));
                }
            }
            orders.push(o);
        }
    }
    orders
}

const CSV_HEAD: &str =
    "case,n,class,faces,l2,linf,order_l2,order_linf,div_after_rel,iters,unknowns";

#[test]
#[ignore = "P2 ladder (minutes)"]
fn projection_mms_ladder() {
    let mut csv = out_file("p2_mms.csv");
    if let Some(f) = csv.as_mut() {
        writeln!(f, "{CSV_HEAD}").ok();
    }
    let mut verdict = Vec::new();
    for (name, iface) in [
        ("quadratic", Interface::Quadratic),
        ("octree", Interface::Octree),
        ("direct", Interface::Direct),
    ] {
        let rungs: Vec<Rung> = [16, 32, 64]
            .iter()
            .map(|&n| mms_rung(n, iface, None, ustar_mms, udiv, zero).0)
            .collect();
        for r in &rungs {
            assert!(
                r.div_rel <= 1e-10,
                "{name} n {}: div {:.2e}",
                r.n,
                r.div_rel
            );
        }
        let orders = ladder_report(name, &rungs, &mut csv);
        if iface == Interface::Quadratic {
            let last = orders.last().expect("pair");
            for c in 0..4 {
                let (o2, oi) = last[c];
                let ok = o2 >= 1.8 && oi >= 1.8;
                eprintln!(
                    "P2 GATE {} L2 {o2:.2} Linf {oi:.2} -> {}",
                    CLASSES[c],
                    if ok { "HELD" } else { "FAILED" }
                );
                verdict.push((CLASSES[c], ok));
            }
        }
    }
    for (c, ok) in verdict {
        assert!(ok, "P2 order gate failed on {c}");
    }
}

// ------------------------------------------------ uniform MAC comparator

/// A uniform-grid MAC projection in the same integrated form (the
/// comparator of P3/P4): open faces are those with a positive coefficient,
/// the Dirichlet boundary as on the composite.
struct UniformMac {
    g: Grid,
    prob: Problem,
    /// Per axis: face area (0 closed) and the Dirichlet value on boundary
    /// faces.
    area: [Vec<f64>; 3],
    pb: [Vec<f64>; 3],
}

fn ufi(g: Grid, d: usize, c: [usize; 3]) -> usize {
    match d {
        0 => g.uface(c[2], c[1], c[0]),
        1 => g.vface(c[2], c[1], c[0]),
        _ => g.wface(c[2], c[1], c[0]),
    }
}

impl UniformMac {
    fn new(n: usize, body: Option<&Sdf>, dirichlet: Scalar) -> Self {
        let g = Grid::cubic(n, n, n, 1.0 / n as f64);
        let (prob, _) = uniform_problem(g, body, &zero, &dirichlet);
        let h = g.dx;
        let dims = [n; 3];
        let mut area: [Vec<f64>; 3] = Default::default();
        let mut pb: [Vec<f64>; 3] = Default::default();
        for d in 0..3 {
            let nf =
                (n + usize::from(d == 0)) * (n + usize::from(d == 1)) * (n + usize::from(d == 2));
            area[d] = vec![0.0; nf];
            pb[d] = vec![0.0; nf];
        }
        for k in 0..n {
            for j in 0..n {
                for i in 0..n {
                    let idx = g.cell(k, j, i);
                    if !prob.active[idx] {
                        continue;
                    }
                    // The plus face of each axis, and the minus face on the
                    // domain boundary.
                    for d in 0..3 {
                        let mut c = [i, j, k];
                        let coef = match d {
                            0 => prob.ae[idx],
                            1 => prob.an[idx],
                            _ => prob.at[idx],
                        };
                        let centre = |c: [usize; 3]| {
                            let mut x = [
                                (c[0] as f64 + 0.5) * h,
                                (c[1] as f64 + 0.5) * h,
                                (c[2] as f64 + 0.5) * h,
                            ];
                            x[d] = c[d] as f64 * h;
                            x
                        };
                        if c[d] == 0 {
                            let f = ufi(g, d, c);
                            area[d][f] = h * h;
                            let x = centre(c);
                            pb[d][f] = dirichlet(x[0], x[1], x[2]);
                        }
                        c[d] += 1;
                        let f = ufi(g, d, c);
                        if c[d] == dims[d] {
                            area[d][f] = h * h;
                            let x = centre(c);
                            pb[d][f] = dirichlet(x[0], x[1], x[2]);
                        } else if coef > 0.0 {
                            area[d][f] = coef * h;
                        }
                    }
                }
            }
        }
        Self { g, prob, area, pb }
    }

    fn sample(&self, f: Field3) -> [Vec<f64>; 3] {
        let (g, h) = (self.g, self.g.dx);
        std::array::from_fn(|d| {
            (0..self.area[d].len())
                .map(|idx| {
                    let nx = g.nx + usize::from(d == 0);
                    let ny = g.ny + usize::from(d == 1);
                    let c = [idx % nx, (idx / nx) % ny, idx / (nx * ny)];
                    let mut x = [
                        (c[0] as f64 + 0.5) * h,
                        (c[1] as f64 + 0.5) * h,
                        (c[2] as f64 + 0.5) * h,
                    ];
                    x[d] = c[d] as f64 * h;
                    f(d, x)
                })
                .collect()
        })
    }

    fn divergence(&self, u: &[Vec<f64>; 3], out: &mut [f64]) {
        let g = self.g;
        for k in 0..g.nz {
            for j in 0..g.ny {
                for i in 0..g.nx {
                    let idx = g.cell(k, j, i);
                    if !self.prob.active[idx] {
                        out[idx] = 0.0;
                        continue;
                    }
                    let mut s = 0.0;
                    for d in 0..3 {
                        let mut c = [i, j, k];
                        let m = ufi(g, d, c);
                        c[d] += 1;
                        let p = ufi(g, d, c);
                        s += self.area[d][p] * u[d][p] - self.area[d][m] * u[d][m];
                    }
                    out[idx] = s;
                }
            }
        }
    }

    fn subtract_gradient(&self, p: &[f64], u: &mut [Vec<f64>; 3]) {
        let (g, h) = (self.g, self.g.dx);
        let n = [g.nx, g.ny, g.nz];
        for d in 0..3 {
            let nx = g.nx + usize::from(d == 0);
            let ny = g.ny + usize::from(d == 1);
            for idx in 0..u[d].len() {
                if self.area[d][idx] == 0.0 {
                    continue;
                }
                let c = [idx % nx, (idx / nx) % ny, idx / (nx * ny)];
                let cell = |c: [usize; 3]| g.cell(c[2], c[1], c[0]);
                let mut m = c;
                if c[d] == 0 {
                    u[d][idx] -= (p[cell(c)] - self.pb[d][idx]) / (0.5 * h);
                } else if c[d] == n[d] {
                    m[d] -= 1;
                    u[d][idx] -= (self.pb[d][idx] - p[cell(m)]) / (0.5 * h);
                } else {
                    m[d] -= 1;
                    u[d][idx] -= (p[cell(c)] - p[cell(m)]) / h;
                }
            }
        }
    }

    /// Returns (setup s, solve s, mac s, iterations, max |div| after).
    fn project(&self, u: &mut [Vec<f64>; 3]) -> (f64, f64, f64, usize, f64) {
        let t0 = std::time::Instant::now();
        let nc = self.g.cells();
        let mut div = vec![0.0; nc];
        self.divergence(u, &mut div);
        let mut prob = self.prob.clone();
        for (r, d) in prob.rhs.iter_mut().zip(&div) {
            *r -= d;
        }
        let mut mac_s = t0.elapsed().as_secs_f64();
        let l1: f64 = prob.rhs.iter().map(|v| v.abs()).sum();
        let params = MultigridParameters {
            max_iterations: 2000,
            ..MultigridParameters::default()
        };
        let mut p = vec![0.0; nc];
        let sol = solve_pcg(&prob, &mut p, &params, 1e-12 * l1, None);
        assert!(sol.converged, "uniform projection did not converge");
        let t1 = std::time::Instant::now();
        self.subtract_gradient(&p, u);
        self.divergence(u, &mut div);
        mac_s += t1.elapsed().as_secs_f64();
        (
            sol.setup_ns as f64 * 1e-9,
            sol.iterate_ns as f64 * 1e-9,
            mac_s,
            sol.iterations,
            maxabs(&div),
        )
    }
}

#[test]
#[ignore = "P3 cut sphere ladder (minutes)"]
fn projection_cut_sphere() {
    let mut csv = out_file("p3_sphere.csv");
    if let Some(f) = csv.as_mut() {
        writeln!(f, "{CSV_HEAD}").ok();
    }
    let mut extra = out_file("p3_sphere_checks.csv");
    if let Some(f) = extra.as_mut() {
        writeln!(
            f,
            "n,div_after_rel,cut_cells,translate_div_before,translate_change,vs_fine_max_diff,vs_fine_exact_err_fine,vs_fine_exact_err_composite"
        )
        .ok();
    }
    let mut rungs = Vec::new();
    for n in [16usize, 32, 64] {
        let (r, mut mp, u) = mms_rung(
            n,
            Interface::Quadratic,
            Some(sphere()),
            ustar_sph,
            udiv_sph,
            q_sph,
        );
        let cut_cells = mp.comp.cut.iter().filter(|&&c| c).count();
        // A translating body: uniform u* plus its wall flux is already
        // divergence-free and must come out unchanged.
        let vel = [0.3, -0.2, 0.5];
        let wall = mp.translating_wall_flux(vel);
        let mut ut = mp.sample(&move |d, _x| vel[d]);
        let u0 = ut.clone();
        let mut pt = vec![0.0; mp.comp.unknowns()];
        let mut div0 = vec![0.0; mp.comp.unknowns()];
        mp.divergence(&ut, Some(&wall), &mut div0);
        let translate_div = maxabs(&div0);
        let st = mp.project(&mut ut, &mut pt, Some(&wall), 1e-12, 400);
        // By linearity P(u0) − u0 − P(0) = G A⁻¹ (D u0 + wall): the change
        // beyond what the Dirichlet data (q on the boundary) drives alone.
        let mut uz = mp.zeros();
        let mut pz = vec![0.0; mp.comp.unknowns()];
        let _ = mp.project(&mut uz, &mut pz, None, 1e-12, 400);
        let mut change = 0.0f64;
        for d in 0..3 {
            for (i, ((a, b), z)) in ut.coarse[d]
                .iter()
                .zip(&u0.coarse[d])
                .zip(&uz.coarse[d])
                .enumerate()
            {
                if mp.coarse_class[d][i] != FaceClass::Covered {
                    change = change.max((a - b - z).abs());
                }
            }
            for (i, ((a, b), z)) in ut.fine[d]
                .iter()
                .zip(&u0.fine[d])
                .zip(&uz.fine[d])
                .enumerate()
            {
                if mp.fine_class[d][i] != FaceClass::Closed {
                    change = change.max((a - b - z).abs());
                }
            }
        }
        assert!(st.converged);
        // Against the uniformly fine projection with the same cut (on the
        // fine faces of the patch, which coincide).
        let um = UniformMac::new(2 * n, Some(&sphere()), q_sph);
        let mut uu = um.sample(ustar_sph);
        let _ = um.project(&mut uu);
        let lo = mp.comp.lo;
        let off = [2 * lo[0], 2 * lo[1], 2 * lo[2]];
        let fg = mp.comp.fine;
        let (mut dmax, mut efine, mut ecomp) = (0.0f64, 0.0f64, 0.0f64);
        for d in 0..3 {
            let nx = fg.nx + usize::from(d == 0);
            let ny = fg.ny + usize::from(d == 1);
            for (i, v) in u.fine[d].iter().enumerate() {
                if mp.fine_class[d][i] != FaceClass::Fine {
                    continue;
                }
                let c = [i % nx, (i / nx) % ny, i / (nx * ny)];
                let gi = ufi(um.g, d, [c[0] + off[0], c[1] + off[1], c[2] + off[2]]);
                let x = mp.face_centre(true, d, i);
                let ex = udiv_sph(d, x);
                dmax = dmax.max((v - uu[d][gi]).abs());
                efine = efine.max((uu[d][gi] - ex).abs());
                ecomp = ecomp.max((v - ex).abs());
            }
        }
        eprintln!(
            "P3 n {n}: div after/before {:.2e}, cut cells {cut_cells}, translating: div {translate_div:.2e} change {change:.2e}; vs uniform fine: max diff {dmax:.3e} (fine err {efine:.3e}, composite err {ecomp:.3e})",
            r.div_rel
        );
        if let Some(f) = extra.as_mut() {
            writeln!(
                f,
                "{n},{:.3e},{cut_cells},{translate_div:.3e},{change:.3e},{dmax:.3e},{efine:.3e},{ecomp:.3e}",
                r.div_rel
            )
            .ok();
        }
        assert!(r.div_rel <= 1e-10, "P3 n {n}: div {:.2e}", r.div_rel);
        assert!(
            translate_div <= 1e-12,
            "P3 n {n}: translating div {translate_div:.2e}"
        );
        assert!(change <= 1e-10, "P3 n {n}: translating change {change:.2e}");
        rungs.push(r);
    }
    let orders = ladder_report("sphere", &rungs, &mut csv);
    let last = orders.last().expect("pair");
    let mut ok_all = true;
    for c in [0usize, 1] {
        let (o2, oi) = last[c];
        let ok = o2 >= 1.8;
        ok_all &= ok;
        eprintln!(
            "P3 GATE {} L2 {o2:.2} (Linf {oi:.2} reported) -> {}",
            CLASSES[c],
            if ok { "HELD" } else { "FAILED" }
        );
    }
    for c in [2usize, 3, 4, 5] {
        let (o2, oi) = last[c];
        eprintln!("P3 report {} L2 {o2:.2} Linf {oi:.2}", CLASSES[c]);
    }
    assert!(ok_all, "P3 interface order gate failed");
}

#[test]
#[ignore = "P4 cost (minutes)"]
fn projection_cost() {
    let mut csv = out_file("p4_cost.csv");
    if let Some(f) = csv.as_mut() {
        writeln!(
            f,
            "field,case,n,patch,unknowns,iters,build_s,solve_s,mac_s,total_s,l2_patch,linf_patch,div_after"
        )
        .ok();
    }
    let n = 64usize;
    let h = 1.0 / n as f64;
    type Case = (&'static str, Field3, Scalar);
    let fields: [Case; 2] = [("smooth", ustar_mms, zero), ("bump", ustar_bump, q_bump)];
    for ((field, ustar, dirichlet), (label, lo, hi)) in fields
        .into_iter()
        .flat_map(|f| [("half", 16usize, 48usize), ("quarter", 24, 40)].map(move |p| (f, p)))
    {
        let t0 = std::time::Instant::now();
        let g = Grid::cubic(n, n, n, h);
        let mut mp = projection(g, [lo; 3], [hi; 3], Interface::Quadratic, None, dirichlet);
        let build_s = t0.elapsed().as_secs_f64();
        let mut u = mp.sample(&ustar);
        let mut p = vec![0.0; mp.comp.unknowns()];
        let st = mp.project(&mut u, &mut p, None, 1e-10, 400);
        assert!(st.converged);
        let mut e = Err::default();
        let hf = mp.comp.fine.dx;
        for d in 0..3 {
            for (i, v) in u.fine[d].iter().enumerate() {
                let x = mp.face_centre(true, d, i);
                e.add(v - udiv(d, x), hf * hf * hf);
            }
        }
        let total = build_s + st.solve_s + st.mac_s;
        eprintln!(
            "P4 {field} composite {label}: {} unknowns, {} it, build+setup {build_s:.2} s, solve {:.2} s, mac {:.2} s, total {total:.2} s, patch L2 {:.3e} Linf {:.3e}, div {:.2e}",
            mp.comp.unknowns(),
            st.iterations,
            st.solve_s,
            st.mac_s,
            e.l2(),
            e.max,
            st.div_after
        );
        if let Some(f) = csv.as_mut() {
            writeln!(
                f,
                "{field},composite,{n},{lo}..{hi},{},{},{build_s:.3},{:.3},{:.3},{total:.3},{:.6e},{:.6e},{:.3e}",
                mp.comp.unknowns(),
                st.iterations,
                st.solve_s,
                st.mac_s,
                e.l2(),
                e.max,
                st.div_after
            )
            .ok();
        }
        for m in [2 * n, n] {
            let t0 = std::time::Instant::now();
            let um = UniformMac::new(m, None, dirichlet);
            let build_s = t0.elapsed().as_secs_f64();
            let mut uu = um.sample(ustar);
            let (setup_s, solve_s, mac_s, it, div) = um.project(&mut uu);
            let (x0, x1) = (lo as f64 * h, hi as f64 * h);
            let mut e = Err::default();
            let hm = um.g.dx;
            for d in 0..3 {
                let nx = m + usize::from(d == 0);
                let ny = m + usize::from(d == 1);
                for (i, v) in uu[d].iter().enumerate() {
                    let c = [i % nx, (i / nx) % ny, i / (nx * ny)];
                    let mut x = [
                        (c[0] as f64 + 0.5) * hm,
                        (c[1] as f64 + 0.5) * hm,
                        (c[2] as f64 + 0.5) * hm,
                    ];
                    x[d] = c[d] as f64 * hm;
                    if x.iter().all(|&v| v >= x0 - 1e-12 && v <= x1 + 1e-12) {
                        e.add(v - udiv(d, x), hm * hm * hm);
                    }
                }
            }
            let total = build_s + setup_s + solve_s + mac_s;
            let case = if m == 2 * n {
                "uniform-fine"
            } else {
                "uniform-coarse"
            };
            eprintln!(
                "P4 {field} {case} n {m} (region {label}): {} cells, {it} it, build+setup {:.2} s, solve {solve_s:.2} s, mac {mac_s:.2} s, total {total:.2} s, region L2 {:.3e} Linf {:.3e}, div {div:.2e}",
                m * m * m,
                build_s + setup_s,
                e.l2(),
                e.max
            );
            if let Some(f) = csv.as_mut() {
                writeln!(
                    f,
                    "{field},{case},{m},{lo}..{hi},{},{it},{:.3},{solve_s:.3},{mac_s:.3},{total:.3},{:.6e},{:.6e},{div:.3e}",
                    m * m * m,
                    build_s + setup_s,
                    e.l2(),
                    e.max
                )
                .ok();
            }
        }
    }
}