//! embedded3 item 12: Schäfer–Turek DFG 3D-2Z (the laminar cylinder in the
//! square channel, Re 20, steady) on the device cut-cell wall. Channel
//! 2.5 × 0.41 × 0.41, cylinder D 0.1 at (0.5, 0.2) across the width,
//! inflow U(y, z) = 16 U_m y z (H − y)(H − z)/H⁴ with U_m 0.45 (Ū = 0.2),
//! ρ 1, ν 1e-3. Coefficients `c = 2F/(ρ Ū² D H)`, Δp between the front
//! and back stagnation points. Reference (the DFG bar): c_D 6.05–6.25,
//! c_L 0.008–0.010, Δp 0.165–0.175.
//!
//! `RTX_E3_DFG_NY=62 RTX_CUDA_ARCH=sm_120 cargo test --release -p rtx-cfd --features cuda --test embedded3_dfg_2z -- --ignored --nocapture`
//! `RTX_E3_DFG_T` sets the flow time marched (default 10 s); `RTX_E3_DFG_VTK=<dir>`
//! writes the final instant; `RTX_E3_DFG_CSV=<path>` the load history.
#![cfg(feature = "cuda")]

use rtx_cfd::solvers::incompressible::ConvectionScheme;
use rtx_cfd::solvers::incompressible::embedded3::step::device::DeviceStep;
use rtx_cfd::solvers::incompressible::embedded3::{
    Body, Boundaries, Field, Fluid, Grid, Parameters, Side, Solver, WallScheme, write_vtk,
};
use std::io::Write as _;

const H: f64 = 0.41;
const L: f64 = 2.5;
const D: f64 = 0.1;
const CX: f64 = 0.5;
const CY: f64 = 0.2;
const U_M: f64 = 0.45;
const U_BAR: f64 = 4.0 / 9.0 * U_M;
const RHO: f64 = 1.0;
const NU: f64 = 1e-3;

fn env_f(name: &str, default: f64) -> f64 {
    std::env::var(name)
        .ok()
        .and_then(|v| v.parse().ok())
        .unwrap_or(default)
}

fn inflow(y: f64, z: f64) -> f64 {
    16.0 * U_M * y * z * (H - y) * (H - z) / (H * H * H * H)
}

#[test]
#[ignore = "item 12: the DFG rung on the device (minutes at ny 62, hours at ny 123)"]
fn dfg_3d_2z_on_the_device() {
    let ny = env_f("RTX_E3_DFG_NY", 62.0) as usize;
    let h = H / ny as f64;
    let nx = (L / h).round() as usize;
    let nz = ny;
    let t_end = env_f("RTX_E3_DFG_T", 10.0);
    // Explicit stability: CFL 0.3 on U_m and half the viscous limit.
    let dt = (0.3 * h / U_M).min(0.5 * h * h / (6.0 * NU));
    let mut solver = Solver::new(
        Fluid {
            density: RHO,
            viscosity: RHO * NU,
            reference_velocity: U_BAR,
            reference_length: D,
        },
        Parameters {
            corrector_steps: 2,
            tolerance: 1e-8,
            convection_scheme: ConvectionScheme::TvdVanAlbada,
            wall_scheme: WallScheme::CutCell,
            boundaries: Boundaries {
                x1: Side::PressureOutlet,
                ..Boundaries::default()
            },
            ..Parameters::default()
        },
    );
    solver.set_boundary_velocity(|x, y, z, _t| {
        if x <= 0.0 {
            (inflow(y, z), 0.0, 0.0)
        } else {
            (0.0, 0.0, 0.0)
        }
    });
    // The cylinder extruded across the width, with samples for the
    // traction route (S2-1).
    solver.set_body(Body::extruded(
        rtx_cfd::solvers::incompressible::EmbeddedBody::circle(CX, CY, 0.5 * D),
        H,
    ));
    let g = Grid::cubic(nx, ny, nz, h);
    let mut field = Field::new(g);
    // Start from the inflow profile everywhere (a faster approach to steady).
    for k in 0..nz {
        for j in 0..ny {
            let u0 = inflow((j as f64 + 0.5) * h, (k as f64 + 0.5) * h);
            for i in 0..=nx {
                field.u[g.uface(k, j, i)] = u0;
            }
        }
    }
    solver.initialize(&mut field);
    let mask_cells = solver.mask().map_or(0, |m| m.fluid_cells());
    println!(
        "  DFG 3D-2Z ny {ny}: {nx}×{ny}×{nz} = {} cells ({mask_cells} fluid), h {h:.4e}, dt {dt:.3e}, t_end {t_end}",
        g.cells()
    );
    let mut device = DeviceStep::new(solver, g);
    device.upload(&field);
    let steps = (t_end / dt).ceil() as usize;
    let coef = 2.0 / (RHO * U_BAR * U_BAR * D * H);
    let csv = std::env::var("RTX_E3_DFG_CSV").ok().map(|p| {
        let mut f = std::fs::File::create(p).expect("csv");
        writeln!(
            f,
            "t,cd_wall,cl_wall,cd_cv,cl_cv,dp,residual,cg,cd_sampler,cl_sampler"
        )
        .unwrap();
        f
    });
    let mut csv = csv;
    let sample_every = (steps / 100).max(1);
    let start = std::time::Instant::now();
    let mut last: Option<(f64, f64, f64, f64, f64)> = None;
    let mut last_sampler = (f64::NAN, f64::NAN);
    let mut settled = false;
    for step in 0..steps {
        let r = device.advance(dt);
        if (step + 1) % sample_every == 0 || step + 1 == steps {
            device.download(&mut field);
            let solver = &device.solver;
            let mask = solver.mask().expect("mask");
            let body = solver.body().expect("body");
            let t = solver.time();
            let fw = mask
                .cut_wall_force(body, &field, RHO * NU, t)
                .expect("wall");
            let margin = 3.0 * D;
            let ci = |x: f64| ((x / h).round() as usize).clamp(2, nx - 2);
            let cj = |y: f64| ((y / h).round() as usize).clamp(2, ny - 2);
            let bx = (
                ci(CX - margin),
                ci(CX + margin),
                cj(CY - 0.15),
                cj(CY + 0.15),
                0,
                nz,
            );
            let fcv =
                mask.control_volume_force_with_walls(&field, dt, RHO, RHO * NU, None, bx, true);
            // The reconstructed wall route (two probes on the cut polygons).
            let fr = mask
                .cut_wall_force_reconstructed(body, &field, RHO * NU, t, None)
                .expect("reconstructed");
            let fs = rtx_cfd::solvers::incompressible::embedded3::SurfaceForce {
                f: fr,
                samples: 0,
                skipped: 0,
            };
            let (cd_s, cl_s) = (coef * fs.f[0], coef * fs.f[1]);
            let zc = 0.5 * H;
            let p_front = mask
                .pressure_at(&field.p, CX - 0.5 * D, CY, zc)
                .unwrap_or(f64::NAN);
            let p_back = mask
                .pressure_at(&field.p, CX + 0.5 * D, CY, zc)
                .unwrap_or(f64::NAN);
            let dp = p_front - p_back;
            let (cd, cl, cd_cv, cl_cv) = (coef * fw[0], coef * fw[1], coef * fcv[0], coef * fcv[1]);
            println!(
                "  t {t:8.4}: c_D {cd:.4} (CV {cd_cv:.4}, reconstructed {cd_s:.4} skipped {})  c_L {cl:.5} (CV {cl_cv:.5}, reconstructed {cl_s:.5})  Δp {dp:.4}  residual {:.1e} CG {}  [{:.0} s]",
                fs.skipped,
                r.final_residual,
                r.poisson_iterations,
                start.elapsed().as_secs_f64()
            );
            if let Some(f) = csv.as_mut() {
                writeln!(
                    f,
                    "{t:.5},{cd:.6},{cl:.6},{cd_cv:.6},{cl_cv:.6},{dp:.6},{:.3e},{},{cd_s:.6},{cl_s:.6}",
                    r.final_residual, r.poisson_iterations
                )
                .unwrap();
            }
            last_sampler = (cd_s, cl_s);
            if let Some((pcd, pcl, _, _, pdp)) = last {
                let rel = ((cd - pcd) / cd)
                    .abs()
                    .max(((dp - pdp) / dp).abs())
                    .max((cl - pcl).abs() / 0.01);
                if rel < 1e-4 && t > 2.0 {
                    settled = true;
                }
            }
            last = Some((cd, cl, cd_cv, cl_cv, dp));
            if settled {
                println!("  settled (relative change < 1e-4 between samples)");
                break;
            }
        }
    }
    device.download(&mut field);
    let solver = &device.solver;
    let mask = solver.mask().expect("mask");
    if let Ok(dir) = std::env::var("RTX_E3_DFG_VTK") {
        let path = std::path::Path::new(&dir).join(format!("dfg_2z_ny{ny}.vtk"));
        write_vtk(&path, &field, Some(mask)).expect("vtk");
        println!("  instant written to {}", path.display());
    }
    let (cd, cl, cd_cv, cl_cv, dp) = last.expect("samples");
    {
        // S2-1 diagnosis: each wall route split into its pressure and shear parts.
        let body = solver.body().expect("body");
        let (po, so) = mask
            .cut_wall_force_parts(body, &field, RHO * NU, solver.time())
            .expect("parts");
        let (pr, sr) = mask
            .cut_wall_force_reconstructed_parts(body, &field, RHO * NU, solver.time(), None)
            .expect("parts");
        println!(
            "  SPLIT ny {ny}: operator c_D pressure {:.4} + shear {:.4}; reconstructed pressure {:.4} + shear {:.4}",
            coef * po[0],
            coef * so[0],
            coef * pr[0],
            coef * sr[0]
        );
    }
    println!(
        "  FINAL ny {ny}: c_D {cd:.4} (CV {cd_cv:.4}, routes {:.2e} apart; reconstructed {:.4}, {:.2e} from CV)  c_L {cl:.5} (CV {cl_cv:.5}, reconstructed {:.5})  Δp {dp:.4}  — reference c_D 6.05–6.25, c_L 0.008–0.010, Δp 0.165–0.175; {:.0} s",
        ((cd - cd_cv) / cd).abs(),
        last_sampler.0,
        ((last_sampler.0 - cd_cv) / cd_cv).abs(),
        last_sampler.1,
        start.elapsed().as_secs_f64()
    );
    if let Some(t) = device.timers() {
        println!("  timers: {t:?}");
    }
}