//! embedded3 gate 9a: the manufactured solution with an embedded sphere
//! (centre (0.6, 0.45, 0.5), r 0.2, off-centre so the exact force is not
//! zero by symmetry) carrying the exact field as its surface velocity, on
//! the binary ghost wall. The velocity error falls at the scheme's order,
//! every fluid cell is divergence-free, the compatibility correction
//! shrinks, and both load routes converge to the exact surface integral of
//! the manufactured stress (the control-volume route measures F − M with M
//! the momentum flux through the porous manufactured surface).

use rtx_cfd::solvers::incompressible::ConvectionScheme;
use rtx_cfd::solvers::incompressible::embedded3::{Body, Field, Fluid, Grid, Parameters, Solver};
use std::f64::consts::PI;

const RHO: f64 = 1.0;
const MU: f64 = 0.05;
const C: (f64, f64, f64) = (0.6, 0.45, 0.5);
const R: f64 = 0.2;

fn u3(x: f64, y: f64, z: f64) -> f64 {
    (PI * x).sin() * (PI * y).cos() * (PI * z).cos()
}
fn v3(x: f64, y: f64, z: f64) -> f64 {
    (PI * x).cos() * (PI * y).sin() * (PI * z).cos()
}
fn w3(x: f64, y: f64, z: f64) -> f64 {
    -2.0 * (PI * x).cos() * (PI * y).cos() * (PI * z).sin()
}
fn p3(x: f64, y: f64, z: f64) -> f64 {
    (PI * x).sin() * (PI * y).sin() * (PI * z).sin()
}
/// The velocity gradient ∂u_i/∂x_j and the pressure gradient.
fn grads(x: f64, y: f64, z: f64) -> ([[f64; 3]; 3], [f64; 3]) {
    let (sx, cx) = (PI * x).sin_cos();
    let (sy, cy) = (PI * y).sin_cos();
    let (sz, cz) = (PI * z).sin_cos();
    (
        [
            [PI * cx * cy * cz, -PI * sx * sy * cz, -PI * sx * cy * sz],
            [-PI * sx * sy * cz, PI * cx * cy * cz, -PI * cx * sy * sz],
            [
                2.0 * PI * sx * cy * sz,
                2.0 * PI * cx * sy * sz,
                -2.0 * PI * cx * cy * cz,
            ],
        ],
        [PI * cx * sy * sz, PI * sx * cy * sz, PI * sx * sy * cz],
    )
}
fn source3(x: f64, y: f64, z: f64) -> (f64, f64, f64) {
    let (g, gp) = grads(x, y, z);
    let u = [u3(x, y, z), v3(x, y, z), w3(x, y, z)];
    let lap = -3.0 * PI * PI;
    let conv = |i: usize| u[0] * g[i][0] + u[1] * g[i][1] + u[2] * g[i][2];
    (
        RHO * conv(0) + gp[0] - MU * lap * u[0],
        RHO * conv(1) + gp[1] - MU * lap * u[1],
        RHO * conv(2) + gp[2] - MU * lap * u[2],
    )
}
fn boundary3(x: f64, y: f64, z: f64) -> (f64, f64, f64) {
    let u = if x <= 0.0 || x >= 1.0 {
        0.0
    } else {
        u3(x, y, z)
    };
    let v = if y <= 0.0 || y >= 1.0 {
        0.0
    } else {
        v3(x, y, z)
    };
    let w = if z <= 0.0 || z >= 1.0 {
        0.0
    } else {
        w3(x, y, z)
    };
    (u, v, w)
}

/// Exact force `∮ (−p I + μ(∇u + ∇uᵀ)) n dA` and momentum flux `∮ ρ u (u·n) dA`
/// over the sphere by a fine Fibonacci quadrature.
fn exact_force_and_flux() -> ([f64; 3], [f64; 3]) {
    let n = 200_000;
    let golden = PI * (3.0 - 5.0_f64.sqrt());
    let (mut f, mut m) = ([0.0; 3], [0.0; 3]);
    let da = 4.0 * PI * R * R / n as f64;
    for k in 0..n {
        let zz = 1.0 - 2.0 * (k as f64 + 0.5) / n as f64;
        let rr = (1.0 - zz * zz).sqrt();
        let th = golden * k as f64;
        let nrm = [rr * th.cos(), rr * th.sin(), zz];
        let (x, y, z) = (C.0 + R * nrm[0], C.1 + R * nrm[1], C.2 + R * nrm[2]);
        let (g, _) = grads(x, y, z);
        let p = p3(x, y, z);
        let u = [u3(x, y, z), v3(x, y, z), w3(x, y, z)];
        let un = u[0] * nrm[0] + u[1] * nrm[1] + u[2] * nrm[2];
        for i in 0..3 {
            let mut t = -p * nrm[i];
            for j in 0..3 {
                t += MU * (g[i][j] + g[j][i]) * nrm[j];
            }
            f[i] += t * da;
            m[i] += RHO * u[i] * un * da;
        }
    }
    (f, m)
}

struct Measurement {
    l2_velocity: f64,
    max_div: f64,
    ghost_correction: f64,
    force_surface: [f64; 3],
    skipped: usize,
    force_cv: [f64; 3],
}

fn measure(n: usize) -> Measurement {
    let h = 1.0 / n as f64;
    let dt = 0.4 * (h * h / (4.0 * MU / RHO)).min(h);
    let mut solver = Solver::new(
        Fluid {
            density: RHO,
            viscosity: MU,
            reference_velocity: 1.0,
            reference_length: 1.0,
        },
        Parameters {
            corrector_steps: 2,
            tolerance: 1e-8,
            convection_scheme: ConvectionScheme::Upwind,
            ..Parameters::default()
        },
    );
    solver.set_momentum_source(|x, y, z, _t| source3(x, y, z));
    solver.set_boundary_velocity(|x, y, z, _t| boundary3(x, y, z));
    solver.set_body(
        Body::sphere(|_t| C, R)
            .with_surface_velocity(|x, y, z, _t| (u3(x, y, z), v3(x, y, z), w3(x, y, z))),
    );
    let g = Grid::cubic(n, n, n, h);
    let mut f = Field::new(g);
    solver.initialize(&mut f);
    for _ in 0..200_000 {
        let (bu, bv, bw) = (f.u.clone(), f.v.clone(), f.w.clone());
        solver.advance(&mut f, dt);
        let mut change = 0.0_f64;
        for (a, b) in
            f.u.iter()
                .zip(&bu)
                .chain(f.v.iter().zip(&bv))
                .chain(f.w.iter().zip(&bw))
        {
            change = change.max((a - b).abs());
        }
        if change / dt < 1e-6 {
            break;
        }
    }
    let mask = solver.mask().expect("mask");
    use rtx_cfd::solvers::incompressible::embedded3::FaceKind;
    let (mut sq, mut vol) = (0.0, 0.0);
    let dv = h * h * h;
    for k in 0..n {
        for j in 0..n {
            for i in 1..n {
                if mask.u_kind(g.uface(k, j, i)) == FaceKind::Fluid {
                    let e = f.u[g.uface(k, j, i)]
                        - u3(i as f64 * h, (j as f64 + 0.5) * h, (k as f64 + 0.5) * h);
                    sq += e * e * dv;
                    vol += dv;
                }
            }
        }
        for j in 1..n {
            for i in 0..n {
                if mask.v_kind(g.vface(k, j, i)) == FaceKind::Fluid {
                    let e = f.v[g.vface(k, j, i)]
                        - v3((i as f64 + 0.5) * h, j as f64 * h, (k as f64 + 0.5) * h);
                    sq += e * e * dv;
                    vol += dv;
                }
            }
        }
    }
    for k in 1..n {
        for j in 0..n {
            for i in 0..n {
                if mask.w_kind(g.wface(k, j, i)) == FaceKind::Fluid {
                    let e = f.w[g.wface(k, j, i)]
                        - w3((i as f64 + 0.5) * h, (j as f64 + 0.5) * h, k as f64 * h);
                    sq += e * e * dv;
                    vol += dv;
                }
            }
        }
    }
    let mut max_div = 0.0_f64;
    for k in 0..n {
        for j in 0..n {
            for i in 0..n {
                if mask.is_fluid_cell(g.cell(k, j, i)) {
                    let div = (f.u[g.uface(k, j, i + 1)] - f.u[g.uface(k, j, i)]) / h
                        + (f.v[g.vface(k, j + 1, i)] - f.v[g.vface(k, j, i)]) / h
                        + (f.w[g.wface(k + 1, j, i)] - f.w[g.wface(k, j, i)]) / h;
                    max_div = max_div.max(div.abs());
                }
            }
        }
    }
    let body = solver.body().expect("body");
    let surface = mask.surface_force(body, &f, MU, solver.time(), 0.5 * h);
    let (i0, i1) = (n / 8, n - n / 8);
    let src = |x: f64, y: f64, z: f64| source3(x, y, z);
    let force_cv = mask.control_volume_force(&f, dt, RHO, MU, Some(&src), (i0, i1, i0, i1, i0, i1));
    Measurement {
        l2_velocity: (sq / vol).sqrt(),
        max_div,
        ghost_correction: solver.ghost_correction().abs(),
        force_surface: surface.f,
        skipped: surface.skipped,
        force_cv,
    }
}

fn norm(a: [f64; 3]) -> f64 {
    (a[0] * a[0] + a[1] * a[1] + a[2] * a[2]).sqrt()
}

fn ladder(resolutions: &[usize]) {
    let (fe, m) = exact_force_and_flux();
    let f_scale = norm(fe);
    let fcv = [fe[0] - m[0], fe[1] - m[1], fe[2] - m[2]];
    println!(
        "  exact force {fe:.5?}; momentum flux {m:.5?}; the control-volume route measures {fcv:.5?}"
    );
    let ms: Vec<Measurement> = resolutions.iter().map(|&n| measure(n)).collect();
    let errors: Vec<f64> = ms.iter().map(|x| x.l2_velocity).collect();
    let mut se = Vec::new();
    let mut ce = Vec::new();
    for (k, (mm, &n)) in ms.iter().zip(resolutions).enumerate() {
        let rate = if k == 0 {
            "    -".to_string()
        } else {
            format!("{:5.2}", (errors[k - 1] / errors[k]).log2())
        };
        let s = norm([
            mm.force_surface[0] - fe[0],
            mm.force_surface[1] - fe[1],
            mm.force_surface[2] - fe[2],
        ]) / f_scale;
        let c = norm([
            mm.force_cv[0] - fcv[0],
            mm.force_cv[1] - fcv[1],
            mm.force_cv[2] - fcv[2],
        ]) / f_scale;
        println!(
            "  n = {n:3}  L2 u {:.4e} (order {rate})  max div {:.2e}  ghost corr {:.2e}  F_surface {:.4?} rel {s:.3e} (skipped {})  F_cv {:.4?} rel {c:.3e}",
            mm.l2_velocity,
            mm.max_div,
            mm.ghost_correction,
            mm.force_surface,
            mm.skipped,
            mm.force_cv
        );
        se.push(s);
        ce.push(c);
    }
    assert!(
        errors.windows(2).all(|w| w[1] < w[0]),
        "errors not monotone {errors:?}"
    );
    for w in errors.windows(2) {
        let rate = (w[0] / w[1]).log2();
        assert!(
            rate > 0.75 && rate < 2.3,
            "order {rate:.3} outside [0.75, 2.3]"
        );
    }
    for mm in &ms {
        assert!(mm.max_div < 1e-5, "max div {:.3e}", mm.max_div);
    }
    assert!(
        se.windows(2).all(|w| w[1] < w[0]),
        "surface-route error not falling {se:?}"
    );
    assert!(
        ce.windows(2).all(|w| w[1] < w[0]),
        "control-volume-route error not falling {ce:?}"
    );
}

#[test]
fn embedded_sphere_recovers_the_manufactured_solution() {
    ladder(&[12, 24]);
}

#[test]
#[ignore = "the three-rung ladder to n = 48 (minutes on the host)"]
fn embedded_sphere_three_rungs() {
    ladder(&[12, 24, 48]);
}