//! A-P3, first task (`docs/overset_metal_campaign.md` §5.10): the wall
//! force on the curvilinear patch. On the phantom circle of the embedded
//! MMS ((0.6, 0.45), r = 0.2, the exact field as the wall velocity) the
//! patch's surface force must converge to the exact surface integral of
//! the manufactured stress, `F = ∮ (−p I + μ(∇u + ∇uᵀ)) n ds` — the
//! embedded solver's staircase reconstruction measured relative errors
//! 0.52 / 0.29 / 0.15 at n = 16 / 32 / 64 (first order). The patch runs
//! standalone with its acceptor ring stamped from the exact field (the S3
//! harness), so the wall force is the patch's own.

use rtx_cfd::mesh::PatchSide;
use rtx_cfd::mesh::patch_gen::annulus_skewed;
use rtx_cfd::solvers::incompressible::{
    CurvilinearParameters, CurvilinearPisoSolver, NormalDiffusion, PatchField,
};
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;

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 p_exact(x: f64, y: f64) -> f64 {
    (PI * x).sin() * (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(),
    )
}

/// `F = ∮ (−p I + μ(∇u + ∇uᵀ)) n ds` on the circle by fine quadrature
/// (`embedded_mms.rs`).
fn exact_force() -> (f64, f64) {
    let n = 20_000;
    let (mut fx, mut fy) = (0.0, 0.0);
    for k in 0..n {
        let theta = (k as f64 + 0.5) * 2.0 * PI / n as f64;
        let (s, c) = theta.sin_cos();
        let (x, y) = (CX + R0 * c, CY + R0 * s);
        let ux = PI * (PI * x).cos() * (PI * y).cos();
        let uy = -PI * (PI * x).sin() * (PI * y).sin();
        let vx = PI * (PI * x).sin() * (PI * y).sin();
        let vy = -PI * (PI * x).cos() * (PI * y).cos();
        let p = p_exact(x, y);
        let sxx = -p + 2.0 * MU * ux;
        let syy = -p + 2.0 * MU * vy;
        let sxy = MU * (uy + vx);
        let ds = 2.0 * PI * R0 / n as f64;
        fx += (sxx * c + sxy * s) * ds;
        fy += (sxy * c + syy * s) * ds;
    }
    (fx, fy)
}

async fn wall_force(n: usize) -> CfdResult<((f64, f64), usize)> {
    let mesh = annulus_skewed([CX, CY], R0, R1, 9 * n / 4, n / 4, 0.3, 3.0)?;
    let nu = MU / RHO;
    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 dt = 0.4 * (hs * hs / (4.0 * nu)).min(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 solver = CurvilinearPisoSolver::new(
        config,
        CurvilinearParameters {
            tolerance: 1e-5,
            normal_diffusion: NormalDiffusion::LineImplicit,
            ..CurvilinearParameters::default()
        },
        mesh,
    )?;
    solver.set_boundary_velocity(|x, y, _| (u_exact(x, y), v_exact(x, y)));
    solver.set_momentum_source(|x, y, _| source(x, y));
    solver.set_acceptor_ring(true);
    let (ns, nn) = (solver.mesh().ns(), solver.mesh().nn());
    let acc: Vec<(f64, f64, f64)> = (0..ns)
        .map(|i| {
            let xy = solver.mesh().centre(solver.mesh().cell(nn - 1, i));
            (
                u_exact(xy[0], xy[1]),
                v_exact(xy[0], xy[1]),
                p_exact(xy[0], xy[1]),
            )
        })
        .collect();
    let zeros = vec![0.0; ns];
    let mut field = PatchField::new(solver.mesh());
    solver.initialize(&mut field, |_, _| (0.0, 0.0));
    solver.stamp_acceptors(&mut field, &acc);
    solver.set_acceptor_correction(&zeros);
    let mut steady = f64::INFINITY;
    let mut steps = 0;
    for _ in 0..400_000 {
        let before = (field.u.clone(), field.v.clone());
        solver.advance(&mut field, dt).await?;
        solver.stamp_acceptors(&mut field, &acc);
        steps += 1;
        let change = field
            .u
            .iter()
            .zip(&before.0)
            .chain(field.v.iter().zip(&before.1))
            .map(|(a, b)| (a - b).abs())
            .fold(0.0, f64::max);
        steady = change / dt;
        if steady < 1e-6 {
            break;
        }
    }
    assert!(steady < 1e-6, "no steady state: {steady:.3e}");
    let load = solver.surface_force(&field, PatchSide::Inner, solver.time());
    let f = load.total();
    Ok(((f[0], f[1]), steps))
}

#[tokio::test]
async fn wall_force_on_the_phantom_circle_converges_to_the_exact_stress_integral() -> CfdResult<()>
{
    let (ex, ey) = exact_force();
    let scale = (ex * ex + ey * ey).sqrt();
    println!("  exact force ({ex:.6e}, {ey:.6e}), |F| {scale:.4e}");
    let mut errs = Vec::new();
    for n in [32usize, 64, 128] {
        let ((fx, fy), steps) = wall_force(n).await?;
        let rel = ((fx - ex).powi(2) + (fy - ey).powi(2)).sqrt() / scale;
        println!(
            "  patch n = {n} (ns {} nn {}): force ({fx:.6e}, {fy:.6e}), relative error {rel:.3e}, {steps} steps",
            9 * n / 4,
            n / 4
        );
        errs.push(rel);
    }
    let orders: Vec<f64> = errs.windows(2).map(|w| (w[0] / w[1]).log2()).collect();
    println!("  wall-force orders {orders:?} (embedded staircase: 0.52 / 0.29 / 0.15 at 16/32/64)");
    // Better than the staircase reconstruction at equal h, and converging.
    assert!(
        errs[0] < 0.29 && errs[1] < 0.15,
        "patch wall force worse than the staircase: {errs:?}"
    );
    assert!(
        orders.iter().all(|&o| o > 0.8),
        "wall-force orders {orders:?}"
    );
    Ok(())
}