//! The embedded-sphere manufactured solution shared by the embedded3
//! wall gates (items 9–11): the fields, the source, the boundary data,
//! the exact surface integrals, and one steady march measured.

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

pub const RHO: f64 = 1.0;
pub const MU: f64 = 0.05;
/// The sphere's default centre (off-centre so the exact force is not zero by symmetry).
pub const C: (f64, f64, f64) = (0.6, 0.45, 0.5);
pub const R: f64 = 0.2;

pub fn u3(x: f64, y: f64, z: f64) -> f64 {
    (PI * x).sin() * (PI * y).cos() * (PI * z).cos()
}
pub fn v3(x: f64, y: f64, z: f64) -> f64 {
    (PI * x).cos() * (PI * y).sin() * (PI * z).cos()
}
pub fn w3(x: f64, y: f64, z: f64) -> f64 {
    -2.0 * (PI * x).cos() * (PI * y).cos() * (PI * z).sin()
}
pub 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.
pub 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],
    )
}
pub 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],
    )
}
pub 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.
pub fn exact_force_and_flux(c: (f64, f64, f64)) -> ([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)
}

pub struct Measurement {
    pub l2_velocity: f64,
    pub max_div: f64,
    pub ghost_correction: f64,
    /// The scheme's wall route (the traction sampler on the binary wall,
    /// the operator route on the cut wall).
    pub force_surface: [f64; 3],
    pub skipped: usize,
    pub force_cv: [f64; 3],
    /// The traction sampler on the cut wall (S2-1); the wall route again
    /// on the binary wall.
    pub force_sampler: [f64; 3],
}

/// March the manufactured solution with the sphere at `c` to steady state on grid `n`.
pub fn measure(n: usize, scheme: WallScheme, c: (f64, f64, f64)) -> 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,
            wall_scheme: scheme,
            momentum_volume_cell_mean: std::env::var("RTX_E3_CELL_MEAN").is_ok_and(|v| v == "1"),
            ..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(move |_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);
    let mut last = solver.advance(&mut f, dt);
    for _ in 0..200_000 {
        let (bu, bv, bw) = (f.u.clone(), f.v.clone(), f.w.clone());
        last = 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");
    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 body = solver.body().expect("body");
    let t = solver.time();
    // The apertured divergence per unit volume, the porous surface's flux
    // through the wall included (the plain divergence on the binary wall);
    // a virtually merged small cell's flux counts with its master's (only
    // the pair's continuity holds).
    let mut max_div = 0.0_f64;
    let mut at_vol = 1.0;
    let mut sum_flux = 0.0;
    let (wall_fluxes, _) = mask.wall_flux_table(body, t);
    let mut cell_flux = vec![0.0; g.cells()];
    for k in 0..n {
        for j in 0..n {
            for i in 0..n {
                let idx = g.cell(k, j, i);
                if mask.is_fluid_cell(idx) {
                    let flux = (mask.a_u(g.uface(k, j, i + 1)) * f.u[g.uface(k, j, i + 1)]
                        - mask.a_u(g.uface(k, j, i)) * f.u[g.uface(k, j, i)])
                        * h
                        * h
                        + (mask.a_v(g.vface(k, j + 1, i)) * f.v[g.vface(k, j + 1, i)]
                            - mask.a_v(g.vface(k, j, i)) * f.v[g.vface(k, j, i)])
                            * h
                            * h
                        + (mask.a_w(g.wface(k + 1, j, i)) * f.w[g.wface(k + 1, j, i)]
                            - mask.a_w(g.wface(k, j, i)) * f.w[g.wface(k, j, i)])
                            * h
                            * h
                        + wall_fluxes[idx];
                    cell_flux[mask.master(idx).unwrap_or(idx)] += flux;
                }
            }
        }
    }
    for (idx, &flux) in cell_flux.iter().enumerate() {
        sum_flux += flux.abs();
        if (flux / (h * h * h)).abs() > max_div {
            max_div = (flux / (h * h * h)).abs();
            at_vol = mask.vol(idx);
        }
    }
    println!(
        "    [{scheme:?} n {n}] max div {max_div:.2e} in a cell of fluid fraction {at_vol:.3e}; Σ|flux| {sum_flux:.2e}; last step residual {:.2e}",
        last.final_residual
    );
    let surface = match scheme {
        WallScheme::GhostBinary => mask.surface_force(body, &f, MU, t, 0.5 * h),
        WallScheme::CutCell => rtx_cfd::solvers::incompressible::embedded3::SurfaceForce {
            f: mask.cut_wall_force(body, &f, MU, t).expect("cut wall"),
            samples: 0,
            skipped: 0,
        },
    };
    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));
    let force_sampler = match scheme {
        WallScheme::GhostBinary => surface.f,
        WallScheme::CutCell => mask
            .cut_wall_force_reconstructed(body, &f, MU, t, None)
            .expect("reconstructed"),
    };
    Measurement {
        l2_velocity: (sq / vol).sqrt(),
        max_div,
        ghost_correction: solver.ghost_correction().abs(),
        force_surface: surface.f,
        skipped: surface.skipped,
        force_cv,
        force_sampler,
    }
}