//! Shared harness of the overset A-P2 tests (`overset_mms.rs`,
//! `overset_moving.rs`): the P2 manufactured problem, the composite
//! builder, the two-mesh error measure and the steady march.
#![allow(dead_code)]

use rtx_cfd::mesh::PatchMesh;
use rtx_cfd::mesh::PatchSide;
use rtx_cfd::mesh::patch_gen::{annulus_skewed, stadium};
use rtx_cfd::solvers::incompressible::{
    CellClass, CurvilinearParameters, CurvilinearPisoSolver, EmbeddedParameters,
    EmbeddedPisoSolver, FaceKind, FlowField, NormalDiffusion, OversetField, OversetParameters,
    OversetPisoSolver, PatchField, PoissonSolverKind,
};
use rtx_cfd::{CfdConfig, CfdResult};
use std::f64::consts::PI;

const RHO: f64 = 1.0;
const MU: f64 = 0.05;
const CX: f64 = 0.6;
const CY: f64 = 0.45;
const R0: f64 = 0.2;
const R1: f64 = 0.354;

pub fn p_exact(x: f64, y: f64) -> f64 {
    (PI * x).sin() * (PI * y).sin()
}
fn u_exact(x: f64, y: f64) -> f64 {
    (PI * x).sin() * (PI * y).cos()
}
fn v_exact(x: f64, y: f64) -> f64 {
    -(PI * x).cos() * (PI * y).sin()
}
fn source(x: f64, y: f64) -> (f64, f64) {
    let conv = RHO * 0.5 * PI;
    (
        conv * (2.0 * PI * x).sin()
            + 2.0 * PI * PI * MU * u_exact(x, y)
            + PI * (PI * x).cos() * (PI * y).sin(),
        conv * (2.0 * PI * y).sin()
            + 2.0 * PI * PI * MU * v_exact(x, y)
            + PI * (PI * x).sin() * (PI * y).cos(),
    )
}
fn boundary_exact(x: f64, y: f64) -> (f64, f64) {
    let u = if x <= 0.0 || x >= 1.0 {
        0.0
    } else {
        u_exact(x, y)
    };
    let v = if y <= 0.0 || y >= 1.0 {
        0.0
    } else {
        v_exact(x, y)
    };
    (u, v)
}

pub fn p2_patch(n: usize) -> CfdResult<PatchMesh> {
    annulus_skewed([CX, CY], R0, R1, 9 * n / 4, n / 4, 0.3, 3.0)
}

pub struct Measurement {
    pub l2_background: f64,
    pub l2_patch: f64,
    /// Mean-shifted L2 pressure errors (background active cells; patch
    /// interior), and the pressure mismatch between the meshes: the largest
    /// |p_bg − p_exact-shift| on active cells adjacent to a fringe cell and
    /// the largest |p_fringe(stamped from the patch) − neighbours' own p|,
    /// both relative to the exact pressure range (2).
    pub l2_p_background: f64,
    pub l2_p_patch: f64,
    pub p_offset_between_meshes: f64,
    pub p_fringe_jump_rel: f64,
    pub steps: usize,
    pub mean_rounds: f64,
    pub max_rounds: usize,
    pub worst_bg_residual: f64,
    pub worst_patch_div_rel: f64,
    pub worst_defect_bg_rel: f64,
    pub worst_defect_patch_rel: f64,
    /// The defects at the steady state (last step), relative to the overlap flux.
    pub final_defect_bg_rel: f64,
    pub final_defect_patch_rel: f64,
    pub schwarz_failures: usize,
}

/// Build the composite for background resolution `n` on `patch_mesh`.
pub fn build(
    n: usize,
    patch_mesh: PatchMesh,
    schwarz_tol: f64,
) -> CfdResult<(OversetPisoSolver, OversetField)> {
    let h = 1.0 / n as f64;
    let config = CfdConfig::new()
        .with_density(RHO)
        .with_viscosity(MU)
        .with_reference_velocity(1.0)
        .with_reference_length(1.0);
    let mut background = EmbeddedPisoSolver::new(
        config.clone(),
        EmbeddedParameters {
            corrector_steps: 2,
            tolerance: 1e-8,
            poisson_solver: PoissonSolverKind::Multigrid,
            ..EmbeddedParameters::default()
        },
    )?;
    background.set_momentum_source(|x, y, _| source(x, y));
    background.set_boundary_velocity(|x, y, _| boundary_exact(x, y));

    let mut patch = CurvilinearPisoSolver::new(
        config,
        CurvilinearParameters {
            tolerance: 1e-5,
            // Line-implicit across the patch on request: the wall cells are
            // h/3 thick and the explicit limit at n = 128 (dt ≈ 3e-6) would
            // need > 1e6 steps to a steady state.
            normal_diffusion: if std::env::var("RTX_OVERSET_LINE").is_ok() {
                NormalDiffusion::LineImplicit
            } else {
                NormalDiffusion::Explicit
            },
            ..CurvilinearParameters::default()
        },
        patch_mesh,
    )?;
    patch.set_side_velocity(PatchSide::Inner, |x, y, _| (u_exact(x, y), v_exact(x, y)));
    patch.set_momentum_source(|x, y, _| source(x, y));
    let mut patch_field = PatchField::new(patch.mesh());
    patch.initialize(&mut patch_field, |_, _| (0.0, 0.0));

    let mut bg_field = FlowField::new(n, n, h, h)?;
    for j in 0..n {
        let y = (j as f64 + 0.5) * h;
        bg_field.u[(j, 0)] = boundary_exact(0.0, y).0;
        bg_field.u[(j, n)] = boundary_exact(1.0, y).0;
    }
    for i in 0..n {
        let x = (i as f64 + 0.5) * h;
        bg_field.v[(0, i)] = boundary_exact(x, 0.0).1;
        bg_field.v[(n, i)] = boundary_exact(x, 1.0).1;
    }
    let params = OversetParameters {
        schwarz_tolerance: schwarz_tol,
        // A steady march: stop rounds that make no progress over two
        // rounds (the step's p' is at the noise floor near the steady
        // state). Never on a transient — see `OversetParameters::stall_rounds`.
        stall_rounds: 2,
        ..OversetParameters::default()
    };
    let mut solver = OversetPisoSolver::new(background, patch, (n, n, h, h), params)?;
    let mut field = OversetField {
        background: bg_field,
        patch: patch_field,
    };
    solver.initialize(&mut field)?;
    Ok((solver, field))
}

/// L2 velocity error on the background's fluid faces and on the patch's
/// interior (non-acceptor) cells.
pub fn errors(solver: &OversetPisoSolver, field: &OversetField, n: usize) -> (f64, f64) {
    let h = 1.0 / n as f64;
    let mask = solver.background().mask().expect("mask");
    let (mut sq, mut area) = (0.0, 0.0);
    for j in 0..n {
        for i in 1..n {
            if mask.u_kind(j, i) == FaceKind::Fluid {
                let e = field.background.u[(j, i)] - u_exact(i as f64 * h, (j as f64 + 0.5) * h);
                sq += e * e * h * h;
                area += h * h;
            }
        }
    }
    for j in 1..n {
        for i in 0..n {
            if mask.v_kind(j, i) == FaceKind::Fluid {
                let e = field.background.v[(j, i)] - v_exact((i as f64 + 0.5) * h, j as f64 * h);
                sq += e * e * h * h;
                area += h * h;
            }
        }
    }
    let l2_bg = (sq / area).sqrt();
    let mesh = solver.patch().mesh();
    let (mut sq, mut vol) = (0.0, 0.0);
    for c in 0..mesh.cell_count() {
        if solver.patch().is_acceptor(c) {
            continue;
        }
        let xy = mesh.centre(c);
        let eu = field.patch.u[c] - u_exact(xy[0], xy[1]);
        let ev = field.patch.v[c] - v_exact(xy[0], xy[1]);
        sq += (eu * eu + ev * ev) * mesh.area(c);
        vol += mesh.area(c);
    }
    (l2_bg, (sq / vol).sqrt())
}

/// The coupled march's stationarity floor: the Schwarz stop (1e-3 of the
/// step's p') and the inner solvers' absolute stops leave |du/dt| noise
/// of ≈ 1e-5 at n = 64 (measured: hovering 5e-6–3e-5 for 70k steps while
/// both L2 errors held to four digits). 3e-5, or the embedded solver's
/// own floor `1e-6 (n/64)²` if larger; `RTX_OVERSET_STEADY` overrides.
pub fn steady_tolerance(n: usize) -> f64 {
    std::env::var("RTX_OVERSET_STEADY")
        .ok()
        .and_then(|v| v.parse().ok())
        .unwrap_or_else(|| (1e-6 * (n as f64 / 64.0).powi(2)).max(3e-5))
}

pub async fn march(n: usize, schwarz_tol: f64) -> CfdResult<Measurement> {
    let steady_tol = steady_tolerance(n);
    let h = 1.0 / n as f64;
    let nu = MU / RHO;
    let (mut solver, mut field) = build(n, p2_patch(n)?, schwarz_tol)?;
    // The patch's wall cells are h/3 thick; the explicit limit is theirs.
    // Line-implicit across the patch, only the along-body spacing counts.
    let line = std::env::var("RTX_OVERSET_LINE").is_ok();
    let mesh = solver.patch().mesh();
    let mut hp = f64::INFINITY;
    for c in 0..mesh.cell_count() {
        for (f, _) in mesh.cell_faces(c) {
            if line && !mesh.is_sface(f) {
                continue;
            }
            let d = mesh.faces()[f].d;
            hp = hp.min((d[0] * d[0] + d[1] * d[1]).sqrt());
        }
    }
    let dt = 0.4 * (hp * hp / (4.0 * nu)).min(h);

    let mut m = Measurement {
        l2_background: 0.0,
        l2_patch: 0.0,
        l2_p_background: 0.0,
        l2_p_patch: 0.0,
        p_offset_between_meshes: 0.0,
        p_fringe_jump_rel: 0.0,
        steps: 0,
        mean_rounds: 0.0,
        max_rounds: 0,
        worst_bg_residual: 0.0,
        worst_patch_div_rel: 0.0,
        worst_defect_bg_rel: 0.0,
        worst_defect_patch_rel: 0.0,
        final_defect_bg_rel: 0.0,
        final_defect_patch_rel: 0.0,
        schwarz_failures: 0,
    };
    let mut total_rounds = 0usize;
    let mut total_correctors = 0usize;
    let trace = std::env::var("RTX_OVERSET_TRACE").is_ok();
    let max_steps: usize = std::env::var("RTX_OVERSET_MAX_STEPS")
        .ok()
        .and_then(|v| v.parse().ok())
        .unwrap_or(400_000);
    let mut steady = f64::INFINITY;
    for step in 0..max_steps {
        let before = (
            field.background.u.clone(),
            field.background.v.clone(),
            field.patch.u.clone(),
            field.patch.v.clone(),
        );
        let r = solver.advance(&mut field, dt).await?;
        assert!(
            r.patch_converged,
            "step {step}: patch pressure solve did not converge: {r:?}"
        );
        m.steps = step + 1;
        total_rounds += r.rounds.iter().sum::<usize>();
        total_correctors += r.rounds.len();
        m.max_rounds = m
            .max_rounds
            .max(r.rounds.iter().copied().max().unwrap_or(0));
        if !r.schwarz_converged {
            m.schwarz_failures += 1;
        }
        m.worst_bg_residual = m.worst_bg_residual.max(r.background_residual);
        let flux_scale: f64 = field
            .patch
            .flux
            .iter()
            .map(|f| f.abs())
            .sum::<f64>()
            .max(1e-300);
        m.worst_patch_div_rel = m
            .worst_patch_div_rel
            .max(r.patch_max_divergence / flux_scale);
        let scale = r.overlap_flux_scale.max(1e-300);
        m.worst_defect_bg_rel = m.worst_defect_bg_rel.max(r.background_mass_defect / scale);
        m.worst_defect_patch_rel = m.worst_defect_patch_rel.max(r.patch_mass_defect / scale);
        m.final_defect_bg_rel = r.background_mass_defect / scale;
        m.final_defect_patch_rel = r.patch_mass_defect / scale;
        let change = (&field.background.u - &before.0)
            .abs()
            .max()
            .max((&field.background.v - &before.1).abs().max())
            .max(
                field
                    .patch
                    .u
                    .iter()
                    .zip(&before.2)
                    .chain(field.patch.v.iter().zip(&before.3))
                    .map(|(a, b)| (a - b).abs())
                    .fold(0.0, f64::max),
            );
        steady = change / dt;
        if trace && step % 500 == 0 {
            let (eb, ep) = errors(&solver, &field, n);
            println!(
                "    step {step} t={:.3}: |du/dt| {steady:.3e} rounds {:?} bg res {:.2e} defect bg {:.2e} patch {:.2e} L2 bg {eb:.3e} patch {ep:.3e}",
                solver.time(),
                r.rounds,
                r.background_residual,
                r.background_mass_defect / scale,
                r.patch_mass_defect / scale
            );
        }
        if steady < steady_tol {
            break;
        }
    }
    assert!(
        steady < steady_tol,
        "no steady state: |du/dt| = {steady:.3e}"
    );
    m.mean_rounds = total_rounds as f64 / total_correctors.max(1) as f64;
    let (eb, ep) = errors(&solver, &field, n);
    m.l2_background = eb;
    m.l2_patch = ep;
    let pe = pressure_errors(&solver, &field, n);
    m.l2_p_background = pe.0;
    m.l2_p_patch = pe.1;
    m.p_offset_between_meshes = pe.2;
    m.p_fringe_jump_rel = pe.3;
    Ok(m)
}

/// Pressure diagnostics: mean-shifted L2 errors on each mesh against the
/// exact pressure, the difference of the two meshes' mean shifts (a level
/// offset between them), and the largest jump between a fringe cell's
/// stamped pressure and the mean of its active neighbours' own pressure,
/// relative to the exact range (2).
pub fn pressure_errors(
    solver: &OversetPisoSolver,
    field: &OversetField,
    n: usize,
) -> (f64, f64, f64, f64) {
    let h = 1.0 / n as f64;
    let map = solver.overlap();
    let (mut sum_d, mut cnt) = (0.0, 0usize);
    for j in 0..n {
        for i in 0..n {
            if map.class(j, i) == CellClass::Active {
                sum_d += field.background.p[(j, i)]
                    - p_exact((i as f64 + 0.5) * h, (j as f64 + 0.5) * h);
                cnt += 1;
            }
        }
    }
    let shift_bg = sum_d / cnt as f64;
    let mut sq = 0.0;
    for j in 0..n {
        for i in 0..n {
            if map.class(j, i) == CellClass::Active {
                let e = field.background.p[(j, i)]
                    - shift_bg
                    - p_exact((i as f64 + 0.5) * h, (j as f64 + 0.5) * h);
                sq += e * e;
            }
        }
    }
    let l2_bg = (sq / cnt as f64).sqrt();
    let mesh = solver.patch().mesh();
    let (mut sum_d, mut vol) = (0.0, 0.0);
    for c in 0..mesh.cell_count() {
        if !solver.patch().is_acceptor(c) {
            let xy = mesh.centre(c);
            sum_d += (field.patch.p[c] - p_exact(xy[0], xy[1])) * mesh.area(c);
            vol += mesh.area(c);
        }
    }
    let shift_patch = sum_d / vol;
    let mut sq = 0.0;
    for c in 0..mesh.cell_count() {
        if !solver.patch().is_acceptor(c) {
            let xy = mesh.centre(c);
            let e = field.patch.p[c] - shift_patch - p_exact(xy[0], xy[1]);
            sq += e * e * mesh.area(c);
        }
    }
    let l2_patch = (sq / vol).sqrt();
    let mut jump = 0.0_f64;
    for e in &map.fringe_cells {
        let (j, i) = (e.j, e.i);
        let (mut ps, mut pc) = (0.0, 0usize);
        for (jj, ii) in [
            (j, i + 1),
            (j, i.wrapping_sub(1)),
            (j + 1, i),
            (j.wrapping_sub(1), i),
        ] {
            if jj < n && ii < n && map.class(jj, ii) == CellClass::Active {
                ps += field.background.p[(jj, ii)];
                pc += 1;
            }
        }
        if pc > 0 {
            jump = jump.max((field.background.p[(j, i)] - ps / pc as f64).abs());
        }
    }
    (l2_bg, l2_patch, shift_bg - shift_patch, jump / 2.0)
}

pub fn orders(errs: &[f64]) -> Vec<f64> {
    errs.windows(2).map(|p| (p[0] / p[1]).log2()).collect()
}

/// The falsifier plate (0.35 × 0.02 m) as a stadium O-grid at background
/// spacing `h`: 16 cells per end arc, straights graded 0.30 h → h at 1.15,
/// offset 6 h, 12 rows stretched 4× (§5.10).
pub fn plate_patch(centre: [f64; 2], h: f64) -> CfdResult<PatchMesh> {
    stadium(centre, 0.175, 0.01, 6.0 * h, 16, h, 1.15, 12, 4.0)
}