//! A-P4 (`docs/overset_metal_campaign.md` §2.2 P4, §5.11): Turek–Hron CFD1
//! (Re = 20, steady) on the OVERSET — the rigid harness's background
//! (`turek_hron_cfd.rs`: parabolic inflow, outlet, multigrid, upwind) with
//! the cylinder–flag O-grid as a static patch (no-slip wall, line-implicit
//! across). Loads by the patch's wall stress (`surface_force`) and by the
//! background's control-volume momentum balance (the two-route rule).
//! Reference (FEATFLOW level 6): drag 14.2929, lift 1.11905. The embedded
//! staircase measured drag 15.71 (surface) / 15.62 (CV) at ny = 41 (+10%).

use rtx_cfd::mesh::patch_gen::cylinder_flag_patch;
use rtx_cfd::mesh::PatchSide;
use rtx_cfd::solvers::incompressible::{
    AleBoundaries, CellClass, CurvilinearParameters, CurvilinearPisoSolver, EmbeddedParameters,
    EmbeddedPisoSolver, FlowField, MomentumResidual, NormalDiffusion, OversetField,
    OversetParameters, OversetPisoSolver, PatchConvection, PatchField, PoissonSolverKind,
    SideBoundary,
};
use rtx_cfd::{CfdConfig, CfdResult};

const L: f64 = 2.5;
const H: f64 = 0.41;
const RHO: f64 = 1000.0;
const NU: f64 = 1e-3;
const U_MEAN: f64 = 0.2;
const REF_DRAG: f64 = 14.2929;
const REF_LIFT: f64 = 1.11905;

/// Tag of a settled-field set: resolution, scheme, and the overlap depth
/// when it is not the default (`RTX_OVERSET_ROWS`).
fn field_tag(ny: usize) -> String {
    let rows = overlap_rows();
    format!(
        "ny{ny}_{}{}",
        if std::env::var("RTX_OVERSET_CFD1_TVD").is_ok() {
            "tvd"
        } else {
            "upwind"
        },
        if rows == OversetParameters::default().overlap_rows {
            String::new()
        } else {
            format!("_rows{rows}")
        }
    )
}

/// `RTX_OVERSET_ROWS`: patch rows kept non-hole below the acceptor row
/// (the overlap depth; default 4). The band probe of §5.11.
fn overlap_rows() -> usize {
    std::env::var("RTX_OVERSET_ROWS")
        .ok()
        .and_then(|v| v.parse().ok())
        .unwrap_or(OversetParameters::default().overlap_rows)
}

fn inflow(y: f64) -> f64 {
    1.5 * U_MEAN * y * (H - y) / (0.5 * H).powi(2)
}

struct Cfd1 {
    drag_surface: f64,
    lift_surface: f64,
    drag_cv: f64,
    lift_cv: f64,
    steps: usize,
    seconds: f64,
    rounds_mean: f64,
    dt: f64,
    residual: MomentumResidual,
}

async fn run_cfd1(ny: usize, max_steps: usize) -> CfdResult<Cfd1> {
    let h = H / ny as f64;
    let nx = (L / h).round() as usize;
    let mu = RHO * NU;
    let config = CfdConfig::new()
        .with_density(RHO)
        .with_viscosity(mu)
        .with_reference_velocity(U_MEAN)
        .with_reference_length(0.1);
    let mut background = EmbeddedPisoSolver::new(
        config.clone(),
        EmbeddedParameters {
            corrector_steps: 2,
            tolerance: 1e-7,
            boundaries: AleBoundaries {
                left: SideBoundary::Velocity,
                right: SideBoundary::PressureOutlet,
                bottom: SideBoundary::Velocity,
                top: SideBoundary::Velocity,
            },
            poisson_solver: PoissonSolverKind::Multigrid,
            ..EmbeddedParameters::default()
        },
    )?;
    background.set_boundary_velocity(|x, y, _| {
        if x <= 0.0 {
            (inflow(y), 0.0)
        } else {
            (0.0, 0.0)
        }
    });

    let (mesh, _) = cylinder_flag_patch(
        [0.2, 0.2],
        0.05,
        0.01,
        0.6,
        h,
        0.5 * 0.41 / 41.0,
        6.0 * h,
        12,
        4.0,
        500,
    )?;
    // The patch's explicit along-body diffusion limit (its wall row is
    // line-implicit); the harness's combined criterion for the background.
    let mut hs = f64::INFINITY;
    for c in 0..mesh.cell_count() {
        for (f, _) in mesh.cell_faces(c) {
            if mesh.is_sface(f) {
                let d = mesh.faces()[f].d;
                hs = hs.min((d[0] * d[0] + d[1] * d[1]).sqrt());
            }
        }
    }
    let u_peak = 1.5 * 1.5 * U_MEAN;
    let dt_bg = 0.25 / (2.0 * u_peak / h + 4.0 * NU / (h * h));
    let dt_patch = 0.4 * (hs * hs / (4.0 * NU)).min(hs / u_peak);
    let dt = dt_bg.min(dt_patch);

    // P4 step 2: `RTX_OVERSET_CFD1_TVD=1` puts the van Albada deferred
    // correction on the patch (the background stays upwind, as recorded).
    let convection = if std::env::var("RTX_OVERSET_CFD1_TVD").is_ok() {
        PatchConvection::TvdVanAlbada
    } else {
        PatchConvection::Upwind
    };
    let mut patch = CurvilinearPisoSolver::new(
        config,
        CurvilinearParameters {
            tolerance: 1e-5,
            convection,
            normal_diffusion: NormalDiffusion::LineImplicit,
            ..CurvilinearParameters::default()
        },
        mesh,
    )?;
    patch.set_side_velocity(PatchSide::Inner, |_, _, _| (0.0, 0.0));
    let mut patch_field = PatchField::new(patch.mesh());
    patch.initialize(&mut patch_field, |_, _| (0.0, 0.0));

    let mut bg_field = FlowField::new(nx, ny, h, h)?;
    for j in 0..ny {
        let u0 = inflow((j as f64 + 0.5) * h);
        for i in 0..=nx {
            bg_field.u[(j, i)] = u0;
        }
    }
    // A steady march: stop a corrector's rounds when two rounds make no
    // progress (the noise floor); measured at ny = 41 without it: 5.0 rounds
    // per corrector on average (the second corrector 8 every step), 2099 s.
    let params = OversetParameters {
        stall_rounds: std::env::var("RTX_OVERSET_STALL")
            .ok()
            .and_then(|v| v.parse().ok())
            .unwrap_or(2),
        // Cost question (P4): does the second corrector's ~9 rounds buy a
        // measurable load? `RTX_OVERSET_MAX_ROUNDS=3` caps every corrector.
        max_rounds: std::env::var("RTX_OVERSET_MAX_ROUNDS")
            .ok()
            .and_then(|v| v.parse().ok())
            .unwrap_or(OversetParameters::default().max_rounds),
        overlap_rows: overlap_rows(),
        ..OversetParameters::default()
    };
    let mut solver = OversetPisoSolver::new(background, patch, (nx, ny, h, h), params)?;
    let mut field = OversetField {
        background: bg_field,
        patch: patch_field,
    };
    solver.initialize(&mut field)?;

    let cv = (
        (0.10 / h).round() as usize,
        (0.75 / h).round() as usize,
        (0.05 / h).round() as usize,
        (0.36 / h).round() as usize,
    );
    let cv_force = |field: &OversetField, solver: &OversetPisoSolver| {
        solver
            .background()
            .mask()
            .expect("mask")
            .control_volume_force(
                &field.background.u,
                &field.background.v,
                &field.background.p,
                &field.background.u_old,
                &field.background.v_old,
                dt,
                RHO,
                mu,
                None,
                cv,
            )
    };

    // `RTX_OVERSET_CFD1_LOAD=dir`: settled fields saved by a previous run
    // (`RTX_OVERSET_CFD1_SAVE`) replace the march — the diagnostics below
    // run offline in seconds instead of the 20–50 min settle.
    let loaded = match std::env::var("RTX_OVERSET_CFD1_LOAD") {
        Ok(dir) => {
            let tag = field_tag(ny);
            let dir = std::path::Path::new(&dir);
            field.background = FlowField::load(&dir.join(format!("bg_{tag}.bin")))?;
            let read = |name: &str| -> Vec<f64> {
                let bytes = std::fs::read(dir.join(format!("patch_{tag}_{name}.bin")))
                    .unwrap_or_else(|e| panic!("load patch {name}: {e}"));
                bytes
                    .chunks_exact(8)
                    .map(|c| f64::from_le_bytes(c.try_into().expect("8 bytes")))
                    .collect()
            };
            field.patch.u = read("u");
            field.patch.v = read("v");
            field.patch.p = read("p");
            field.patch.flux = read("flux");
            assert_eq!(field.patch.u.len(), solver.patch().mesh().cell_count());
            assert_eq!(field.patch.flux.len(), solver.patch().mesh().faces().len());
            println!("    loaded settled fields {tag} from {}", dir.display());
            true
        }
        Err(_) => false,
    };
    let start = std::time::Instant::now();
    let flow_through = L / U_MEAN;
    let min_steps = (flow_through / dt).ceil() as usize;
    let mut history: Vec<f64> = Vec::new();
    let mut steps = 0;
    let mut rounds_total = 0usize;
    let mut correctors_total = 0usize;
    let trace_first = std::env::var("RTX_OVERSET_CFD1_TRACE").is_ok();
    if !loaded {
        loop {
            let r = solver.advance(&mut field, dt).await?;
            steps += 1;
            let every: usize = std::env::var("RTX_OVERSET_CFD1_TRACE_EVERY")
                .ok()
                .and_then(|v| v.parse().ok())
                .unwrap_or(0);
            if (trace_first && steps <= 6) || (every > 0 && steps % every == 0) {
                let pmax = field
                    .background
                    .p
                    .iter()
                    .fold(0.0_f64, |m, v| m.max(v.abs()));
                let ppmax = field.patch.p.iter().fold(0.0_f64, |m, v| m.max(v.abs()));
                let upmax = field.patch.u.iter().fold(0.0_f64, |m, v| m.max(v.abs()));
                println!(
                "      step {steps}: rounds {:?} converged {} stalled {} bg res {:.2e} patch div {:.2e} patch iters {} conv {} | max|p| bg {pmax:.3e} patch {ppmax:.3e} max|u| patch {upmax:.3e} | defect bg {:.2e} patch {:.2e}",
                r.rounds,
                r.schwarz_converged,
                r.schwarz_stalled,
                r.background_residual,
                r.patch_max_divergence,
                r.patch_poisson_iterations,
                r.patch_converged,
                r.background_mass_defect,
                r.patch_mass_defect
            );
            }
            if steps >= max_steps {
                break;
            }
            rounds_total += r.rounds.iter().sum::<usize>();
            correctors_total += r.rounds.len();
            if steps % 50 == 0 {
                let (fx, _) = cv_force(&field, &solver);
                history.push(fx);
                let load =
                    solver
                        .patch()
                        .surface_force(&field.patch, PatchSide::Inner, solver.time());
                let umax = field
                    .background
                    .u
                    .iter()
                    .fold(0.0_f64, |m, v| m.max(v.abs()));
                if steps % 500 == 0 || !umax.is_finite() {
                    println!(
                    "    ny = {ny}: step {steps} t = {:.2} s  drag_cv {fx:.4}  drag_wall {:.4} lift_wall {:.4}  max|u| {umax:.3}  rounds {:?} bg res {:.1e} patch div {:.1e}  [{:.0} s]",
                    solver.time(),
                    load.total()[0],
                    load.total()[1],
                    r.rounds,
                    r.background_residual,
                    r.patch_max_divergence,
                    start.elapsed().as_secs_f64()
                );
                }
                assert!(
                    umax.is_finite(),
                    "velocity became non-finite at step {steps}"
                );
                if steps >= min_steps && history.len() > 4 {
                    let now = history[history.len() - 1];
                    let then = history[history.len() - 5];
                    if ((now - then) / now).abs() < 1e-4 {
                        break;
                    }
                }
            }
            assert!(steps < 2_000_000, "CFD1 at ny = {ny} did not settle");
        }
    }
    let seconds = start.elapsed().as_secs_f64();
    let load = solver
        .patch()
        .surface_force(&field.patch, PatchSide::Inner, solver.time());
    let (drag_cv, lift_cv) = cv_force(&field, &solver);
    // P4 momentum-defect measurement: the force the background transmits
    // into the ring (fringe + hole), into the hole alone, and the patch's
    // wall force — consecutive differences are the active region's
    // residual, the fringe ring's momentum defect, and the patch region's.
    let ring = solver
        .overlap()
        .region_force(&field.background, RHO, mu, |c| c != CellClass::Active);
    let hole = solver
        .overlap()
        .region_force(&field.background, RHO, mu, |c| c == CellClass::Hole);
    let wall = load.total();
    println!(
        "    overlap rows {} (patch nn {})",
        overlap_rows(),
        solver.patch().mesh().nn()
    );
    println!(
        "    momentum routes ny = {ny}: CV box ({drag_cv:.4}, {lift_cv:.4}) | ring outer ({:.4}, {:.4}) | hole boundary ({:.4}, {:.4}) | wall ({:.4}, {:.4}); defects [% of wall drag]: active {:+.2} fringe ring {:+.2} patch region {:+.2}",
        ring.0,
        ring.1,
        hole.0,
        hole.1,
        wall[0],
        wall[1],
        100.0 * (drag_cv - ring.0) / wall[0],
        100.0 * (ring.0 - hole.0) / wall[0],
        100.0 * (hole.0 - wall[0]) / wall[0],
    );
    // P4 option B: the momentum residual of the solver's OWN staggered
    // upwind stencil on every background face at the settled state. Solved
    // faces read zero by construction (the pin below); the prescribed
    // faces' sum is the momentum the stamping injects, in the solver's
    // metric and without the staircase curves' face-formula error.
    let mr = solver.momentum_residual(&field, dt);
    let pct = |b: &rtx_cfd::solvers::incompressible::ResidualBucket| 100.0 * b.fx / wall[0];
    println!(
        "    momentum residual ny = {ny} [N/m, x / y; % of wall drag; faces evaluated/total]: solved far Σr ({:+.3e}, {:+.3e}) Σ|r| ({:.3e}, {:.3e}) {}/{} | solved near ring Σr ({:+.3e}, {:+.3e}) Σ|r| ({:.4}, {:.4}) max|r| ({:.3e}, {:.3e}) {}/{} | fringe–fringe ({:+.4}, {:+.4}) {:+.2}% {}/{} | fringe–hole ({:+.4}, {:+.4}) {:+.2}% {}/{} | hole–hole skipped {} (ghosts: {} cells, {} faces) | Σ|r| fringe–fringe ({:.4}, {:.4}) fringe–hole ({:.4}, {:.4}); ring total ({:+.4}, {:+.4}) {:+.2}% vs routes' ring defect {:+.4} ({:+.2}%)",
        mr.solved_far.fx, mr.solved_far.fy, mr.solved_far.abs_x, mr.solved_far.abs_y, mr.solved_far.evaluated, mr.solved_far.total,
        mr.solved_near.fx, mr.solved_near.fy, mr.solved_near.abs_x, mr.solved_near.abs_y, mr.solved_near.max_abs_x, mr.solved_near.max_abs_y, mr.solved_near.evaluated, mr.solved_near.total,
        mr.fringe_fringe.fx, mr.fringe_fringe.fy, pct(&mr.fringe_fringe), mr.fringe_fringe.evaluated, mr.fringe_fringe.total,
        mr.fringe_hole.fx, mr.fringe_hole.fy, pct(&mr.fringe_hole), mr.fringe_hole.evaluated, mr.fringe_hole.total,
        mr.hole_hole_skipped, mr.hole_ghosts, mr.ghost_faces,
        mr.fringe_fringe.abs_x, mr.fringe_fringe.abs_y, mr.fringe_hole.abs_x, mr.fringe_hole.abs_y,
        mr.fringe_fringe.fx + mr.fringe_hole.fx, mr.fringe_fringe.fy + mr.fringe_hole.fy,
        pct(&mr.fringe_fringe) + pct(&mr.fringe_hole),
        hole.0 - ring.0,
        100.0 * (hole.0 - ring.0) / wall[0],
    );
    // Where along the ring: the prescribed u faces' x-momentum residual in
    // x-bands (cylinder front, cylinder–flag junction, flag, trailing edge).
    let mut bands = [
        (0.0_f64, 0.20, 0.0_f64, 0usize),
        (0.20, 0.30, 0.0, 0),
        (0.30, 0.55, 0.0, 0),
        (0.55, 1.0, 0.0, 0),
    ];
    for f in mr.prescribed.iter().filter(|f| f.is_u && f.r.is_finite()) {
        let x = f.i as f64 * h;
        if let Some(b) = bands.iter_mut().find(|b| x >= b.0 && x < b.1) {
            b.2 += f.r;
            b.3 += 1;
        }
    }
    println!(
        "    ring x-momentum residual by x-band ny = {ny} (N/m, u faces): {}; level offset δ = {:.3e} Pa on {} + {} interface faces",
        bands
            .iter()
            .map(|b| format!("x {:.2}–{:.2}: {:+.4} ({} faces)", b.0, b.1, b.2, b.3))
            .collect::<Vec<_>>()
            .join(" | "),
        mr.level_offset(h),
        mr.interface_u,
        mr.interface_v
    );
    // `RTX_OVERSET_CFD1_FRONT=1`: every prescribed u face with x < 0.20
    // (the cylinder front, where the ring's source sits) with the pieces
    // of its residual, and the band's sums by piece.
    if std::env::var("RTX_OVERSET_CFD1_FRONT").is_ok() {
        let mut sums = [0.0_f64; 4];
        let mut consistency = 0.0_f64;
        for f in mr
            .prescribed
            .iter()
            .filter(|f| f.is_u && f.r.is_finite() && (f.i as f64) * h < 0.20)
        {
            let pc = f.pieces;
            let recon = pc[0] + pc[1] - pc[2] + pc[3];
            consistency = consistency.max((recon - f.r).abs());
            for k in 0..4 {
                sums[k] += pc[k];
            }
            let cls = |jj: usize, ii: usize| {
                format!("{:?}", solver.overlap().class(jj, ii))
                    .chars()
                    .next()
                    .unwrap()
            };
            println!(
                "      front u face j {:2} i {:2} (x {:.3} y {:.3}) {}|{}: r {:+.4e} = time {:+.2e} + conv {:+.2e} − diff {:+.2e} + pres {:+.2e}   u_old {:+.4} u {:+.4}",
                f.j,
                f.i,
                f.i as f64 * h,
                (f.j as f64 + 0.5) * h,
                cls(f.j, f.i - 1),
                cls(f.j, f.i),
                f.r,
                pc[0],
                pc[1],
                pc[2],
                pc[3],
                field.background.u_old[(f.j, f.i)],
                field.background.u[(f.j, f.i)]
            );
        }
        println!(
            "      front band sums: time {:+.4} conv {:+.4} −diff {:+.4} pres {:+.4} (pieces reconstruct r to {:.1e})",
            sums[0], sums[1], -sums[2], sums[3], consistency
        );
    }
    // `RTX_OVERSET_CFD1_RAW=1`: the same ring residual on the RAW
    // interpolated stamping — the prescribed faces re-stamped from the
    // patch without `balance_fringe_fluxes` (u and u_old alike) — against
    // the balanced one, and how far the balance moved the faces. If the
    // ring's source is the balance's, the raw residual is small.
    if std::env::var("RTX_OVERSET_CFD1_RAW").is_ok() {
        let mut raw = OversetField {
            background: field.background.clone(),
            patch: field.patch.clone(),
        };
        solver
            .overlap()
            .stamp_fringe_faces(&mut raw.background, &field.patch.u, &field.patch.v);
        let (mut moved, mut moved_max, mut scale) = (0.0_f64, 0.0_f64, 0.0_f64);
        for e in &solver.overlap().fringe_u {
            let d = raw.background.u[(e.j, e.i)] - field.background.u[(e.j, e.i)];
            moved += d.abs();
            moved_max = moved_max.max(d.abs());
            scale = scale.max(field.background.u[(e.j, e.i)].abs());
        }
        for e in &solver.overlap().fringe_v {
            let d = raw.background.v[(e.j, e.i)] - field.background.v[(e.j, e.i)];
            moved += d.abs();
            moved_max = moved_max.max(d.abs());
        }
        raw.background.u_old.copy_from(&raw.background.u);
        raw.background.v_old.copy_from(&raw.background.v);
        // Active faces keep u_old = u^n; only the prescribed faces changed.
        for j in 0..field.background.u.nrows() {
            for i in 0..field.background.u.ncols() {
                if solver.background().mask().expect("mask").u_kind(j, i)
                    == rtx_cfd::solvers::incompressible::FaceKind::Fluid
                {
                    raw.background.u_old[(j, i)] = field.background.u_old[(j, i)];
                    raw.background.u[(j, i)] = field.background.u[(j, i)];
                }
            }
        }
        for j in 0..field.background.v.nrows() {
            for i in 0..field.background.v.ncols() {
                if solver.background().mask().expect("mask").v_kind(j, i)
                    == rtx_cfd::solvers::incompressible::FaceKind::Fluid
                {
                    raw.background.v_old[(j, i)] = field.background.v_old[(j, i)];
                    raw.background.v[(j, i)] = field.background.v[(j, i)];
                }
            }
        }
        let mr_raw = solver.momentum_residual(&raw, dt);
        let mut front_raw = 0.0;
        for f in mr_raw
            .prescribed
            .iter()
            .filter(|f| f.is_u && f.r.is_finite() && (f.i as f64) * h < 0.20)
        {
            front_raw += f.r;
        }
        println!(
            "    raw-stamping ring residual ny = {ny}: fringe–fringe ({:+.4}, {:+.4}) fringe–hole ({:+.4}, {:+.4}) ring total x {:+.4} (balanced {:+.4}); front band x {:+.4} (balanced {:+.4}); the balance moved the prescribed faces by Σ|Δu| {:.3e} max {:.3e} (max |u| {:.3}); background mass defect raw {:.3e} balanced {:.3e}",
            mr_raw.fringe_fringe.fx,
            mr_raw.fringe_fringe.fy,
            mr_raw.fringe_hole.fx,
            mr_raw.fringe_hole.fy,
            mr_raw.fringe_fringe.fx + mr_raw.fringe_hole.fx,
            mr.fringe_fringe.fx + mr.fringe_hole.fx,
            front_raw,
            mr.prescribed
                .iter()
                .filter(|f| f.is_u && f.r.is_finite() && (f.i as f64) * h < 0.20)
                .map(|f| f.r)
                .sum::<f64>(),
            moved,
            moved_max,
            scale,
            solver.overlap().background_mass_defect(&raw.background),
            solver.overlap().background_mass_defect(&field.background),
        );
    }
    // The box force in the solver's own flux form on the five boxes (must
    // agree to rounding — the gate), and the bookkeeping it allows: the
    // numerical x-momentum source inside the box on the fluid is
    // wall − box; the ring's share is Σr; the rest is the patch region's
    // (hole boundary → wall) plus the two meshes' metric difference.
    let boxes = [
        (0.10, 0.75, 0.05, 0.36),
        (0.09, 0.70, 0.07, 0.34),
        (0.08, 1.00, 0.03, 0.38),
        (0.10, 1.50, 0.05, 0.36),
        (0.10, 0.75, 0.02, 0.39),
    ];
    let solver_forces: Vec<(f64, f64)> = boxes
        .iter()
        .map(|&(x0, x1, y0, y1)| {
            solver.solver_metric_force(
                &field,
                dt,
                (
                    (x0 / h).round() as usize,
                    (x1 / h).round() as usize,
                    (y0 / h).round() as usize,
                    (y1 / h).round() as usize,
                ),
            )
        })
        .collect();
    let sf = solver_forces[0];
    let spread = solver_forces.iter().fold(0.0_f64, |m, f| {
        m.max((f.0 - sf.0).abs()).max((f.1 - sf.1).abs())
    });
    let ring_sum = mr.fringe_fringe.fx + mr.fringe_hole.fx;
    println!(
        "    solver-metric box force ny = {ny}: ({:.4}, {:.4}) [5 boxes spread {:.2e}; CV formula {drag_cv:.4}]; x-momentum sources on the fluid inside the box [N/m, % of wall drag]: total wall − box {:+.4} ({:+.2}%) = ring Σr {:+.4} ({:+.2}%) + rest (patch region + metric) {:+.4} ({:+.2}%)",
        sf.0,
        sf.1,
        spread,
        wall[0] - sf.0,
        100.0 * (wall[0] - sf.0) / wall[0],
        ring_sum,
        100.0 * ring_sum / wall[0],
        wall[0] - sf.0 - ring_sum,
        100.0 * (wall[0] - sf.0 - ring_sum) / wall[0],
    );
    // The patch's own momentum balance on its solved cells (scheme
    // fluxes; the unsteady term is omitted — settled state): the balance
    // residual is the gate; flux_force − wall_force = the least-squares
    // pressure's non-conservation δP; the acceptor band's mismatch is
    // then (box + ring Σr) − flux_force, all in N/m.
    let pb = solver.patch().momentum_balance(&field.patch, solver.time());
    let ff = pb.flux_force();
    let fw = pb.wall_force();
    let dp = pb.pressure_defect();
    let bal = pb.balance();
    let hole_flux = sf.0 + ring_sum;
    println!(
        "    patch momentum balance ny = {ny} [N/m x / y; {} solved cells, {} interface faces, {} wall faces]: balance residual ({:+.3e}, {:+.3e}) | flux-form force through the interface ({:.4}, {:.4}) | wall force, scheme fluxes ({:.4}, {:.4}) | wall force, surface formula ({:.4}, {:.4}) | pressure defect δP = p_ls − p_face ({:+.4}, {:+.4}) [{:+.2}% of wall drag] | pieces: conv_acc ({:+.4}, {:+.4}) visc_acc ({:+.4}, {:+.4}) p_face_acc ({:+.4}, {:+.4}) visc_wall ({:+.4}, {:+.4}) p_face_wall ({:+.4}, {:+.4}) p_ls ({:+.4}, {:+.4}) | acceptor band: background hole flux {:.4} − patch interface {:.4} = {:+.4} ({:+.2}%)",
        pb.cells, pb.acc_faces, pb.wall_faces,
        bal[0], bal[1],
        ff[0], ff[1],
        fw[0], fw[1],
        wall[0], wall[1],
        dp[0], dp[1], 100.0 * dp[0] / wall[0],
        pb.conv_acc[0], pb.conv_acc[1], pb.visc_acc[0], pb.visc_acc[1], pb.p_face_acc[0], pb.p_face_acc[1],
        pb.visc_wall[0], pb.visc_wall[1], pb.p_face_wall[0], pb.p_face_wall[1], pb.p_ls[0], pb.p_ls[1],
        hole_flux, ff[0], hole_flux - ff[0], 100.0 * (hole_flux - ff[0]) / wall[0],
    );
    // `RTX_OVERSET_CFD1_SAVE=dir`: the settled fields, for offline
    // diagnostics without the march (background in `FlowField::save`'s
    // format; patch u, v, p as raw little-endian f64 vectors).
    if let Ok(dir) = std::env::var("RTX_OVERSET_CFD1_SAVE") {
        let tag = field_tag(ny);
        let dir = std::path::Path::new(&dir);
        std::fs::create_dir_all(dir).expect("save dir");
        field.background.save(&dir.join(format!("bg_{tag}.bin")))?;
        for (name, vals) in [
            ("u", &field.patch.u),
            ("v", &field.patch.v),
            ("p", &field.patch.p),
            ("flux", &field.patch.flux),
        ] {
            let bytes: Vec<u8> = vals.iter().flat_map(|x| x.to_le_bytes()).collect();
            std::fs::write(dir.join(format!("patch_{tag}_{name}.bin")), bytes).expect("save patch");
        }
        println!("    saved settled fields to {} as {tag}", dir.display());
    }
    // The wall load split and the fringe ring's extent (the tight box in
    // the sensitivity list must stay outside it).
    let (mut jmin, mut jmax, mut imin, mut imax) = (usize::MAX, 0, usize::MAX, 0);
    for e in &solver.overlap().fringe_cells {
        jmin = jmin.min(e.j);
        jmax = jmax.max(e.j);
        imin = imin.min(e.i);
        imax = imax.max(e.i);
    }
    println!(
        "    wall load split: pressure ({:.4}, {:.4}) viscous ({:.4}, {:.4}); fringe ring cells i {imin}–{imax} (x {:.3}–{:.3}) j {jmin}–{jmax} (y {:.3}–{:.3})",
        load.pressure[0],
        load.pressure[1],
        load.viscous[0],
        load.viscous[1],
        imin as f64 * h,
        (imax + 1) as f64 * h,
        jmin as f64 * h,
        (jmax + 1) as f64 * h
    );
    // Box sensitivity of the control-volume route: the same balance on
    // other rectangles, all outside the fringe ring (x 0.110–0.640, y
    // 0.110–0.290 at ny = 41 — a first list had a box at x0 = 0.12 cutting
    // through it and read −21 %). A route that moves with the box by more
    // than its own formula error cannot arbitrate the gap.
    for (x0, x1, y0, y1) in [
        (0.10, 0.75, 0.05, 0.36),
        (0.09, 0.70, 0.07, 0.34),
        (0.08, 1.00, 0.03, 0.38),
        (0.10, 1.50, 0.05, 0.36),
        (0.10, 0.75, 0.02, 0.39),
    ] {
        let boxc = (
            (x0 / h).round() as usize,
            (x1 / h).round() as usize,
            (y0 / h).round() as usize,
            (y1 / h).round() as usize,
        );
        let (bx, by) = solver
            .background()
            .mask()
            .expect("mask")
            .control_volume_force(
                &field.background.u,
                &field.background.v,
                &field.background.p,
                &field.background.u_old,
                &field.background.v_old,
                dt,
                RHO,
                mu,
                None,
                boxc,
            );
        println!(
            "    CV box x {x0:.2}–{x1:.2} y {y0:.2}–{y1:.2}: drag {bx:.4} ({:+.2}% vs wall) lift {by:.4}",
            100.0 * (bx - wall[0]) / wall[0]
        );
    }
    Ok(Cfd1 {
        drag_surface: load.total()[0],
        lift_surface: load.total()[1],
        drag_cv,
        lift_cv,
        steps,
        seconds,
        rounds_mean: rounds_total as f64 / correctors_total.max(1) as f64,
        dt,
        residual: mr,
    })
}

#[tokio::test]
async fn cfd1_on_the_overset_against_the_featflow_reference() -> CfdResult<()> {
    let resolutions: Vec<usize> = std::env::var("RTX_OVERSET_CFD1_NY").ok().map_or_else(
        || vec![41usize],
        |list| {
            list.split(',')
                .map(|t| t.trim().parse().expect("ny list"))
                .collect()
        },
    );
    for &ny in &resolutions {
        let max_steps: usize = std::env::var("RTX_OVERSET_CFD1_MAX_STEPS")
            .ok()
            .and_then(|v| v.parse().ok())
            .unwrap_or(2_000_000);
        let r = run_cfd1(ny, max_steps).await?;
        let rel = |a: f64, b: f64| 100.0 * (a - b) / b;
        println!(
            "  CFD1 overset ny = {ny} (h = {:.4}, dt = {:.2e}, patch {}): wall drag {:.4} ({:+.2}%) lift {:.4} ({:+.2}%); control volume drag {:.4} ({:+.2}%) lift {:.4}; routes differ {:.2}%; [{} steps, {:.0} s, Schwarz rounds mean {:.2}]  reference {REF_DRAG} / {REF_LIFT}; embedded staircase at ny=41: 15.71 / 15.62 (+10%)",
            H / ny as f64,
            r.dt,
            if std::env::var("RTX_OVERSET_CFD1_TVD").is_ok() {
                "tvd"
            } else {
                "upwind"
            },
            r.drag_surface,
            rel(r.drag_surface, REF_DRAG),
            r.lift_surface,
            rel(r.lift_surface, REF_LIFT),
            r.drag_cv,
            rel(r.drag_cv, REF_DRAG),
            r.lift_cv,
            100.0 * ((r.drag_surface - r.drag_cv) / r.drag_cv).abs(),
            r.steps,
            r.seconds,
            r.rounds_mean
        );
        assert!(r.drag_surface.is_finite() && r.drag_cv.is_finite());
    }
    Ok(())
}

/// The residual diagnostic is the solver's own operator: on every SOLVED
/// background face away from the ring the momentum residual (time term
/// included) is zero to rounding; on the solved faces NEXT to the ring it
/// is a constant per face that cancels in the sum — the composite's
/// pressure LEVEL offset between the active cells (whose `p'` had its mean
/// removed) and the fringe cells (re-stamped from the patch, which never
/// saw that shift); and the ring buckets are populated. This is what makes
/// the prescribed faces' sum readable as the stamping's momentum injection
/// in the solver's metric (§5.11, option B).
#[tokio::test]
async fn momentum_residual_vanishes_on_the_solved_faces() -> CfdResult<()> {
    let r = run_cfd1(41, 5).await?;
    let mr = &r.residual;
    let scale = r.drag_surface.abs().max(1.0);
    let far = &mr.solved_far;
    assert_eq!(
        far.evaluated, far.total,
        "every far solved face is evaluable"
    );
    assert!(
        far.abs_x <= 1e-9 * scale && far.abs_y <= 1e-9 * scale,
        "solved far: Σ|r| = ({:.3e}, {:.3e}) is not rounding against {scale:.3}",
        far.abs_x,
        far.abs_y
    );
    let near = &mr.solved_near;
    assert_eq!(
        near.evaluated, near.total,
        "every near solved face is evaluable"
    );
    assert!(
        near.fx.abs() <= 1e-9 * scale && near.fy.abs() <= 1e-9 * scale,
        "solved near: Σr = ({:.3e}, {:.3e}) does not cancel against {scale:.3}",
        near.fx,
        near.fy
    );
    // A pure level offset: every active–fringe INTERFACE face carries the
    // same |r| = δ·h and every other near face (one that only reads a
    // prescribed velocity) reads zero, so Σ|r| = N_interface · max|r| on
    // each lattice.
    assert!(
        (near.abs_x - mr.interface_u as f64 * near.max_abs_x).abs() <= 1e-6 * near.abs_x.max(1e-300)
            && (near.abs_y - mr.interface_v as f64 * near.max_abs_y).abs()
                <= 1e-6 * near.abs_y.max(1e-300),
        "solved near: not a uniform level offset on the interface — Σ|r| ({:.4e}, {:.4e}) vs N·max|r| ({:.4e}, {:.4e}) with N = ({}, {})",
        near.abs_x,
        near.abs_y,
        mr.interface_u as f64 * near.max_abs_x,
        mr.interface_v as f64 * near.max_abs_y,
        mr.interface_u,
        mr.interface_v
    );
    assert!(mr.fringe_fringe.evaluated > 0 && mr.fringe_hole.evaluated > 0);
    assert_eq!(
        mr.fringe_fringe.evaluated, mr.fringe_fringe.total,
        "every fringe–fringe face has a fully valid stencil"
    );
    assert_eq!(
        mr.fringe_hole.evaluated, mr.fringe_hole.total,
        "every fringe–hole face has a fully valid stencil with the ghost band"
    );
    Ok(())
}