//! The multigrid projection solves the SOR projection's system.
//!
//! `PisoSolver` and `EmbeddedPisoSolver` can run their pressure-correction
//! projection either by point SOR (the historical path) or by the
//! multigrid-preconditioned CG solver of `solvers/incompressible/poisson.rs`
//! (`PoissonSolverKind::Multigrid`). The two are handed the same five-point
//! coefficients, the same mass-imbalance source, the same anchor and the
//! same true-residual stop, so every benchmark of this suite must land on
//! the same discrete answer whichever is selected — to the inner tolerance,
//! which for the steady states below is far below the discretisation error
//! being measured. A multigrid branch that assembled a different system
//! (a wrong outlet coefficient, a dropped anchor, a coefficient across a
//! prescribed face) would move the steady solution by far more than the
//! tolerances here, and would not be caught by the solver's own unit tests,
//! which only see the `PoissonProblem` it is given.
//!
//! Four settings:
//! 1. the fixed-grid PISO manufactured steady solution (`tests/mms_piso.rs`
//!    at n = 32): same L2 velocity error to `1e-6` relative, divergence-free;
//! 2. the decaying Taylor–Green vortex (`tests/taylor_green.rs` at n = 32):
//!    divergence-free on every step, the energy-decay error no worse;
//! 3. the embedded-circle manufactured solution (`tests/embedded_mms.rs` at
//!    n = 32): same L2 velocity error to `1e-5` relative, and the no-body
//!    degeneracy (embedded solver == PISO to the bit) holds with multigrid
//!    on both;
//! 4. a channel with a pressure outlet and an embedded circle — the
//!    Turek–Hron configuration in miniature, which exercises the outlet's
//!    Dirichlet column and the level-free (un-anchored) system: the steady
//!    fields agree to `1e-6` relative.

use rtx_cfd::solvers::incompressible::{
    AleBoundaries, BoundaryConditions, EmbeddedBody, EmbeddedParameters, EmbeddedPisoSolver,
    FaceKind, FlowField, IncompressibleSolver, MgPrecision, PisoParameters, PisoSolver,
    PoissonSolverKind, SideBoundary,
};
use rtx_cfd::{CfdConfig, CfdResult};
use std::f64::consts::PI;

// ---------------------------------------------------------------------------
// Shared manufactured field (tests/mms_piso.rs, tests/embedded_mms.rs).
// ---------------------------------------------------------------------------

const RHO: f64 = 1.0;
const MU: f64 = 0.05;

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 fx = RHO * 0.5 * PI * (2.0 * PI * x).sin()
        + 2.0 * PI * PI * MU * u_exact(x, y)
        + PI * (PI * x).cos() * (PI * y).sin();
    let fy = RHO * 0.5 * PI * (2.0 * PI * y).sin()
        + 2.0 * PI * PI * MU * v_exact(x, y)
        + PI * (PI * x).sin() * (PI * y).cos();
    (fx, fy)
}

fn mms_config() -> CfdConfig {
    CfdConfig::new()
        .with_density(RHO)
        .with_viscosity(MU)
        .with_reference_velocity(1.0)
        .with_reference_length(1.0)
}

fn mms_time_step(n: usize) -> f64 {
    let dx = 1.0 / n as f64;
    let nu = MU / RHO;
    0.4 * (dx * dx / (4.0 * nu)).min(dx)
}

fn max_divergence(field: &FlowField, fluid: impl Fn(usize, usize) -> bool) -> f64 {
    let (nx, ny, dx, dy) = field.grid_info();
    let mut max_div: f64 = 0.0;
    for j in 0..ny {
        for i in 0..nx {
            if !fluid(j, i) {
                continue;
            }
            let div = (field.u[(j, i + 1)] - field.u[(j, i)]) / dx
                + (field.v[(j + 1, i)] - field.v[(j, i)]) / dy;
            max_div = max_div.max(div.abs());
        }
    }
    max_div
}

fn max_change(a: &FlowField, b: &FlowField) -> f64 {
    let mut m: f64 = 0.0;
    for (x, y) in a.u.iter().zip(b.u.iter()) {
        m = m.max((x - y).abs());
    }
    for (x, y) in a.v.iter().zip(b.v.iter()) {
        m = m.max((x - y).abs());
    }
    m
}

fn rel(a: f64, b: f64) -> f64 {
    ((a - b) / b).abs()
}

// ---------------------------------------------------------------------------
// 1. Fixed-grid PISO manufactured steady solution.
// ---------------------------------------------------------------------------

struct SteadyMeasurement {
    l2_velocity: f64,
    max_div: f64,
    steps: usize,
    seconds: f64,
}

/// `tests/mms_piso.rs::measure`, with the inner solver selectable.
async fn piso_mms(n: usize, kind: PoissonSolverKind) -> CfdResult<SteadyMeasurement> {
    let dx = 1.0 / n as f64;
    let dt = mms_time_step(n);
    let mut solver = PisoSolver::new(
        mms_config(),
        PisoParameters {
            corrector_steps: 2,
            time_step: dt,
            tolerance: 1e-8,
            poisson_solver: kind,
            poisson_precision: rtx_cfd::solvers::incompressible::MgPrecision::F64,
        },
    )?;
    solver.set_momentum_source(source);
    solver.set_wall_velocity(|x, y| (u_exact(x, y), v_exact(x, y)));

    let mut field = FlowField::new(n, n, dx, dx)?;
    for j in 0..n {
        let y = (j as f64 + 0.5) * dx;
        field.u[(j, 0)] = u_exact(0.0, y);
        field.u[(j, n)] = u_exact(1.0, y);
    }
    for i in 0..n {
        let x = (i as f64 + 0.5) * dx;
        field.v[(0, i)] = v_exact(x, 0.0);
        field.v[(n, i)] = v_exact(x, 1.0);
    }

    let empty = BoundaryConditions::new();
    let start = std::time::Instant::now();
    let mut steady_residual = f64::INFINITY;
    let mut steps = 0;
    for _ in 0..200_000 {
        let before = field.clone();
        solver.solve_time_step(&mut field, &empty, dt).await?;
        steps += 1;
        steady_residual = max_change(&field, &before) / dt;
        if steady_residual < 1e-6 {
            break;
        }
    }
    assert!(
        steady_residual < 1e-6,
        "PISO ({kind:?}) did not reach a steady state: |du/dt| = {steady_residual:.3e}"
    );

    let mut squared = 0.0;
    let mut volume = 0.0;
    for j in 0..n {
        for i in 1..n {
            let e = field.u[(j, i)] - u_exact(i as f64 * dx, (j as f64 + 0.5) * dx);
            squared += e * e * dx * dx;
            volume += dx * dx;
        }
    }
    for j in 1..n {
        for i in 0..n {
            let e = field.v[(j, i)] - v_exact((i as f64 + 0.5) * dx, j as f64 * dx);
            squared += e * e * dx * dx;
            volume += dx * dx;
        }
    }
    Ok(SteadyMeasurement {
        l2_velocity: (squared / volume).sqrt(),
        max_div: max_divergence(&field, |_, _| true),
        steps,
        seconds: start.elapsed().as_secs_f64(),
    })
}

/// Both inner solvers march the manufactured problem to the same discrete
/// steady state: L2 velocity errors equal to `1e-6` relative (the errors
/// themselves are `~2e-2`, so this is agreement to four orders below the
/// discretisation error), and every cell divergence-free under either.
#[tokio::test]
async fn piso_manufactured_steady_state_is_solver_independent() -> CfdResult<()> {
    let n = 32;
    let sor = piso_mms(n, PoissonSolverKind::Sor).await?;
    let mg = piso_mms(n, PoissonSolverKind::Multigrid).await?;
    println!(
        "  PISO MMS n = {n}: SOR L2 {:.6e} (div {:.2e}, {} steps, {:.1} s)   MG L2 {:.6e} \
         (div {:.2e}, {} steps, {:.1} s)   relative difference {:.2e}",
        sor.l2_velocity,
        sor.max_div,
        sor.steps,
        sor.seconds,
        mg.l2_velocity,
        mg.max_div,
        mg.steps,
        mg.seconds,
        rel(mg.l2_velocity, sor.l2_velocity)
    );
    assert!(
        rel(mg.l2_velocity, sor.l2_velocity) < 1e-6,
        "L2 velocity error differs between inner solvers: SOR {:.8e}, multigrid {:.8e}",
        sor.l2_velocity,
        mg.l2_velocity
    );
    for (name, m) in [("SOR", &sor), ("multigrid", &mg)] {
        assert!(
            m.max_div < 1e-5,
            "{name}: max |div u| = {:.3e}, the projection is not removing the divergence",
            m.max_div
        );
    }
    Ok(())
}

// ---------------------------------------------------------------------------
// 2. Taylor–Green.
// ---------------------------------------------------------------------------

const TG_NU: f64 = 0.02;
const TG_T_END: f64 = 0.25;

fn tg_amplitude(t: f64) -> f64 {
    (-2.0 * TG_NU * PI * PI * t).exp()
}

fn tg_u(x: f64, y: f64, t: f64) -> f64 {
    tg_amplitude(t) * (PI * x).sin() * (PI * y).cos()
}

fn tg_v(x: f64, y: f64, t: f64) -> f64 {
    -tg_amplitude(t) * (PI * x).cos() * (PI * y).sin()
}

fn tg_p(x: f64, y: f64, t: f64) -> f64 {
    let a = tg_amplitude(t);
    -RHO * a * a / 4.0 * ((2.0 * PI * x).cos() + (2.0 * PI * y).cos())
}

fn kinetic_energy(field: &FlowField, n: usize, dx: f64) -> f64 {
    let mut energy = 0.0;
    for j in 0..n {
        for i in 1..n {
            energy += 0.5 * RHO * field.u[(j, i)] * field.u[(j, i)] * dx * dx;
        }
    }
    for j in 1..n {
        for i in 0..n {
            energy += 0.5 * RHO * field.v[(j, i)] * field.v[(j, i)] * dx * dx;
        }
    }
    energy
}

struct TaylorGreen {
    energy_ratio: f64,
    /// Largest |div u| seen after any step.
    max_div_any_step: f64,
    steps: usize,
}

/// `tests/taylor_green.rs::measure`, with the inner solver selectable and
/// the divergence checked after every step rather than at the end only.
async fn taylor_green(n: usize, kind: PoissonSolverKind) -> CfdResult<TaylorGreen> {
    let dx = 1.0 / n as f64;
    let dt = 0.4 * dx * dx / (4.0 * TG_NU);
    let steps = (TG_T_END / dt).ceil() as usize;
    let dt = TG_T_END / steps as f64;

    let config = CfdConfig::new()
        .with_density(RHO)
        .with_viscosity(RHO * TG_NU)
        .with_reference_velocity(1.0)
        .with_reference_length(1.0);
    let mut solver = PisoSolver::new(
        config,
        PisoParameters {
            corrector_steps: 60,
            time_step: dt,
            tolerance: 1e-9,
            poisson_solver: kind,
            poisson_precision: rtx_cfd::solvers::incompressible::MgPrecision::F64,
        },
    )?;

    let mut field = FlowField::new(n, n, dx, dx)?;
    for j in 0..n {
        let y = (j as f64 + 0.5) * dx;
        for i in 0..=n {
            field.u[(j, i)] = tg_u(i as f64 * dx, y, 0.0);
        }
    }
    for j in 0..=n {
        let y = j as f64 * dx;
        for i in 0..n {
            field.v[(j, i)] = tg_v((i as f64 + 0.5) * dx, y, 0.0);
        }
    }
    for j in 0..n {
        for i in 0..n {
            field.p[(j, i)] = tg_p((i as f64 + 0.5) * dx, (j as f64 + 0.5) * dx, 0.0);
        }
    }
    let initial_energy = kinetic_energy(&field, n, dx);

    let empty = BoundaryConditions::new();
    let mut max_div_any_step: f64 = 0.0;
    for step in 0..steps {
        let t = step as f64 * dt;
        solver.set_wall_velocity(move |x, y| (tg_u(x, y, t), tg_v(x, y, t)));
        let result = solver.solve_time_step(&mut field, &empty, dt).await?;
        assert!(
            result.solver_result.converged,
            "{kind:?} step {step}: projection left mass residual {:.3e}",
            result.solver_result.final_residual
        );
        max_div_any_step = max_div_any_step.max(max_divergence(&field, |_, _| true));
    }
    Ok(TaylorGreen {
        energy_ratio: kinetic_energy(&field, n, dx) / initial_energy,
        max_div_any_step,
        steps,
    })
}

/// With the multigrid projection every step of the decaying vortex is
/// divergence-free and the kinetic-energy decay is no further from the
/// closed form `e^(-4 nu pi^2 T)` than with SOR — and equal to it to `1e-6`
/// relative, since both solve the same projection to the same stop.
#[tokio::test]
async fn taylor_green_is_divergence_free_every_step_with_multigrid() -> CfdResult<()> {
    let n = 32;
    let exact_ratio = (-4.0 * TG_NU * PI * PI * TG_T_END).exp();
    let sor = taylor_green(n, PoissonSolverKind::Sor).await?;
    let mg = taylor_green(n, PoissonSolverKind::Multigrid).await?;
    let deficit_sor = exact_ratio - sor.energy_ratio;
    let deficit_mg = exact_ratio - mg.energy_ratio;
    println!(
        "  Taylor-Green n = {n} ({} steps): E(T)/E(0) SOR {:.8} MG {:.8} (exact {exact_ratio:.8}); \
         deficits SOR {:.3e} MG {:.3e}; max |div u| over all steps SOR {:.2e} MG {:.2e}",
        mg.steps,
        sor.energy_ratio,
        mg.energy_ratio,
        deficit_sor,
        deficit_mg,
        sor.max_div_any_step,
        mg.max_div_any_step
    );
    assert!(
        mg.max_div_any_step < 1e-5,
        "multigrid: a step left max |div u| = {:.3e}",
        mg.max_div_any_step
    );
    assert!(
        deficit_mg.abs() <= deficit_sor.abs() * (1.0 + 1e-6) + 1e-12,
        "multigrid energy-decay error {deficit_mg:.6e} is worse than SOR's {deficit_sor:.6e}"
    );
    assert!(
        rel(mg.energy_ratio, sor.energy_ratio) < 1e-6,
        "energy ratios differ between inner solvers: SOR {:.10}, multigrid {:.10}",
        sor.energy_ratio,
        mg.energy_ratio
    );
    Ok(())
}

// ---------------------------------------------------------------------------
// 3. Embedded circle.
// ---------------------------------------------------------------------------

const CX: f64 = 0.6;
const CY: f64 = 0.45;
const R: f64 = 0.2;

/// The manufactured field on the box boundary with the normal components
/// snapped to their exact analytic zero (see `tests/embedded_mms.rs`).
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)
}

fn mms_initial_field(n: usize) -> CfdResult<FlowField> {
    let dx = 1.0 / n as f64;
    let mut field = FlowField::new(n, n, dx, dx)?;
    for j in 0..n {
        let y = (j as f64 + 0.5) * dx;
        field.u[(j, 0)] = boundary_exact(0.0, y).0;
        field.u[(j, n)] = boundary_exact(1.0, y).0;
    }
    for i in 0..n {
        let x = (i as f64 + 0.5) * dx;
        field.v[(0, i)] = boundary_exact(x, 0.0).1;
        field.v[(n, i)] = boundary_exact(x, 1.0).1;
    }
    Ok(field)
}

/// `tests/embedded_mms.rs::measure`, velocity error and divergence only,
/// with the inner solver selectable.
async fn embedded_mms(n: usize, kind: PoissonSolverKind) -> CfdResult<SteadyMeasurement> {
    embedded_mms_with(n, kind, MgPrecision::F64).await
}

/// `embedded_mms` with the multigrid V-cycle precision selectable (the M1
/// precision probe of `overset_metal_campaign.md` §3.2).
async fn embedded_mms_with(
    n: usize,
    kind: PoissonSolverKind,
    precision: MgPrecision,
) -> CfdResult<SteadyMeasurement> {
    let dx = 1.0 / n as f64;
    let dt = mms_time_step(n);
    let mut solver = EmbeddedPisoSolver::new(
        mms_config(),
        EmbeddedParameters {
            corrector_steps: 2,
            tolerance: 1e-8,
            poisson_solver: kind,
            poisson_precision: precision,
            ..EmbeddedParameters::default()
        },
    )?;
    solver.set_momentum_source(|x, y, _| source(x, y));
    solver.set_boundary_velocity(|x, y, _| boundary_exact(x, y));
    solver.set_body(
        EmbeddedBody::circle(CX, CY, R)
            .with_surface_velocity(|x, y, _| (u_exact(x, y), v_exact(x, y))),
    );
    let mut field = mms_initial_field(n)?;
    solver.initialize(&mut field)?;

    let start = std::time::Instant::now();
    let mut steady_residual = f64::INFINITY;
    let mut steps = 0;
    for _ in 0..200_000 {
        let before = field.clone();
        solver.advance(&mut field, dt).await?;
        steps += 1;
        steady_residual = max_change(&field, &before) / dt;
        if steady_residual < 1e-6 {
            break;
        }
    }
    assert!(
        steady_residual < 1e-6,
        "embedded PISO ({kind:?}) did not reach a steady state at n = {n}: |du/dt| = \
         {steady_residual:.3e}"
    );
    let mask = solver.mask().expect("mask built");
    let mut squared = 0.0;
    let mut volume = 0.0;
    for j in 0..n {
        for i in 1..n {
            if mask.u_kind(j, i) == FaceKind::Fluid {
                let e = field.u[(j, i)] - u_exact(i as f64 * dx, (j as f64 + 0.5) * dx);
                squared += e * e * dx * dx;
                volume += dx * dx;
            }
        }
    }
    for j in 1..n {
        for i in 0..n {
            if mask.v_kind(j, i) == FaceKind::Fluid {
                let e = field.v[(j, i)] - v_exact((i as f64 + 0.5) * dx, j as f64 * dx);
                squared += e * e * dx * dx;
                volume += dx * dx;
            }
        }
    }
    Ok(SteadyMeasurement {
        l2_velocity: (squared / volume).sqrt(),
        max_div: max_divergence(&field, |j, i| mask.is_fluid_cell(j, i)),
        steps,
        seconds: start.elapsed().as_secs_f64(),
    })
}

/// The embedded-circle manufactured steady state (`tests/embedded_mms.rs`
/// records L2 velocity 8.489e-3 at n = 32 with SOR) is the same under
/// multigrid to `1e-5` relative, and divergence-free on every fluid cell.
#[tokio::test]
async fn embedded_circle_steady_state_is_solver_independent() -> CfdResult<()> {
    let n = 32;
    let sor = embedded_mms(n, PoissonSolverKind::Sor).await?;
    let mg = embedded_mms(n, PoissonSolverKind::Multigrid).await?;
    println!(
        "  embedded MMS n = {n}: SOR L2 {:.6e} (div {:.2e}, {} steps, {:.1} s)   MG L2 {:.6e} \
         (div {:.2e}, {} steps, {:.1} s)   relative difference {:.2e}   (recorded SOR value 8.4892e-3)",
        sor.l2_velocity,
        sor.max_div,
        sor.steps,
        sor.seconds,
        mg.l2_velocity,
        mg.max_div,
        mg.steps,
        mg.seconds,
        rel(mg.l2_velocity, sor.l2_velocity)
    );
    assert!(
        rel(mg.l2_velocity, sor.l2_velocity) < 1e-5,
        "L2 velocity error differs between inner solvers: SOR {:.8e}, multigrid {:.8e}",
        sor.l2_velocity,
        mg.l2_velocity
    );
    for (name, m) in [("SOR", &sor), ("multigrid", &mg)] {
        assert!(
            m.max_div < 1e-5,
            "{name}: a fluid cell is not divergence-free, max |div u| = {:.3e}",
            m.max_div
        );
    }
    Ok(())
}

/// `tests/embedded_mms.rs::without_a_body_the_embedded_solver_is_piso_to_the_bit`
/// with the multigrid projection on both solvers: the no-body embedded
/// solver assembles exactly the fixed-grid PISO's `PoissonProblem` (all
/// cells active, anchor `(1, 1)`, no outlet), so the fields must still agree
/// to the bit.
/// M1 precision probe (`overset_metal_campaign.md` §3.2 M1, §5.2): with
/// the V-cycle in single precision inside the f64 conjugate gradient, the
/// embedded-circle steady state must be the same field to the projection
/// stop — the f64 CG owns the true residual, so the preconditioner's
/// precision may cost iterations, never accuracy. The measured gap is
/// printed; the assert is at the tolerance scale, not at f32 eps.
#[tokio::test]
async fn embedded_circle_steady_state_survives_an_f32_vcycle() -> CfdResult<()> {
    let n = 32;
    let f64_arm = embedded_mms_with(n, PoissonSolverKind::Multigrid, MgPrecision::F64).await?;
    let f32_arm = embedded_mms_with(n, PoissonSolverKind::Multigrid, MgPrecision::F32).await?;
    println!(
        "f32 V-cycle vs f64: L2 error {:.6e} vs {:.6e} (rel {:.2e}), max div {:.2e} vs {:.2e}, \
         steps {} vs {}, wall {:.2} s vs {:.2} s",
        f32_arm.l2_velocity,
        f64_arm.l2_velocity,
        rel(f32_arm.l2_velocity, f64_arm.l2_velocity),
        f32_arm.max_div,
        f64_arm.max_div,
        f32_arm.steps,
        f64_arm.steps,
        f32_arm.seconds,
        f64_arm.seconds
    );
    assert!(
        rel(f32_arm.l2_velocity, f64_arm.l2_velocity) < 1e-5,
        "f32 V-cycle changed the steady error: {:.6e} vs {:.6e}",
        f32_arm.l2_velocity,
        f64_arm.l2_velocity
    );
    assert!(
        f32_arm.max_div < 1e-5,
        "divergence with the f32 V-cycle {:.3e}",
        f32_arm.max_div
    );
    Ok(())
}

#[tokio::test]
async fn without_a_body_the_embedded_solver_is_piso_to_the_bit_with_multigrid() -> CfdResult<()> {
    let n = 16;
    let dt = mms_time_step(n);
    let mut piso = PisoSolver::new(
        mms_config(),
        PisoParameters {
            corrector_steps: 2,
            time_step: dt,
            tolerance: 1e-8,
            poisson_solver: PoissonSolverKind::Multigrid,
            poisson_precision: rtx_cfd::solvers::incompressible::MgPrecision::F64,
        },
    )?;
    piso.set_momentum_source(source);
    piso.set_wall_velocity(boundary_exact);
    let mut embedded = EmbeddedPisoSolver::new(
        mms_config(),
        EmbeddedParameters {
            corrector_steps: 2,
            tolerance: 1e-8,
            poisson_solver: PoissonSolverKind::Multigrid,
            poisson_precision: rtx_cfd::solvers::incompressible::MgPrecision::F64,
            ..EmbeddedParameters::default()
        },
    )?;
    embedded.set_momentum_source(|x, y, _| source(x, y));
    embedded.set_boundary_velocity(|x, y, _| boundary_exact(x, y));

    let mut a = mms_initial_field(n)?;
    let mut b = mms_initial_field(n)?;
    let empty = BoundaryConditions::new();
    for _ in 0..200 {
        piso.solve_time_step(&mut a, &empty, dt).await?;
        embedded.advance(&mut b, dt).await?;
    }
    let mut max_diff: f64 = max_change(&a, &b);
    for (x, y) in a.p.iter().zip(b.p.iter()) {
        max_diff = max_diff.max((x - y).abs());
    }
    // Both must also have actually moved the field — a pair of solvers that
    // both do nothing agree to the bit too.
    let moved = max_change(&a, &mms_initial_field(n)?);
    assert!(
        moved > 1e-3,
        "the solvers did not advance the field ({moved:.3e})"
    );
    assert!(
        max_diff == 0.0,
        "embedded solver without a body differs from PISO by {max_diff:.3e} under multigrid"
    );
    Ok(())
}

// ---------------------------------------------------------------------------
// 4. Channel with a pressure outlet and an embedded circle.
// ---------------------------------------------------------------------------

/// Steady channel flow past a circle with a pressure outlet on the right —
/// the outlet's Dirichlet column makes the system non-singular (no anchor),
/// the circle masks cells: every arm of the embedded assembly is exercised.
/// Returns the steady `u`, `v`, `p` fields.
async fn channel_with_circle(kind: PoissonSolverKind) -> CfdResult<(FlowField, usize, f64)> {
    let (length, height) = (2.0, 0.5);
    let ny = 20;
    let h = height / ny as f64;
    let nx = (length / h).round() as usize;
    let (rho, nu, u_mean) = (1.0, 0.01, 1.0);
    let u_peak = 1.5 * 1.5 * u_mean;
    let dt = 0.25 / (2.0 * u_peak / h + 4.0 * nu / (h * h));
    let config = CfdConfig::new()
        .with_density(rho)
        .with_viscosity(rho * nu)
        .with_reference_velocity(u_mean)
        .with_reference_length(height);
    let mut solver = EmbeddedPisoSolver::new(
        config,
        EmbeddedParameters {
            corrector_steps: 2,
            tolerance: 1e-8,
            boundaries: AleBoundaries {
                left: SideBoundary::Velocity,
                right: SideBoundary::PressureOutlet,
                bottom: SideBoundary::Velocity,
                top: SideBoundary::Velocity,
            },
            poisson_solver: kind,
            poisson_precision: rtx_cfd::solvers::incompressible::MgPrecision::F64,
            ..EmbeddedParameters::default()
        },
    )?;
    let inflow = move |y: f64| 1.5 * u_mean * y * (height - y) / (0.5 * height).powi(2);
    solver.set_boundary_velocity(move |x, y, _| {
        if x <= 0.0 {
            (inflow(y), 0.0)
        } else {
            (0.0, 0.0)
        }
    });
    solver.set_body(EmbeddedBody::circle(0.5, 0.27, 0.1));
    let mut field = FlowField::new(nx, ny, h, h)?;
    for j in 0..ny {
        let u0 = inflow((j as f64 + 0.5) * h);
        for i in 0..=nx {
            field.u[(j, i)] = u0;
        }
    }
    solver.initialize(&mut field)?;

    let start = std::time::Instant::now();
    let mut steps = 0;
    let mut steady_residual = f64::INFINITY;
    for _ in 0..200_000 {
        let before = field.clone();
        solver.advance(&mut field, dt).await?;
        steps += 1;
        steady_residual = max_change(&field, &before) / dt;
        if steady_residual < 1e-6 * u_mean {
            break;
        }
    }
    assert!(
        steady_residual < 1e-6 * u_mean,
        "channel ({kind:?}) did not reach a steady state: |du/dt| = {steady_residual:.3e}"
    );
    let mask = solver.mask().expect("mask built");
    let max_div = max_divergence(&field, |j, i| mask.is_fluid_cell(j, i));
    assert!(
        max_div < 1e-5 * u_mean / h,
        "channel ({kind:?}): a fluid cell is not divergence-free, max |div u| = {max_div:.3e}"
    );
    Ok((field, steps, start.elapsed().as_secs_f64()))
}

/// The steady channel-with-circle fields agree between the two inner
/// solvers to `1e-6` relative (RMS of the difference over the RMS of the
/// field), for velocity and for pressure — the outlet column and the
/// un-anchored system are assembled as the SOR loop forms them.
#[tokio::test]
async fn outlet_channel_with_circle_steady_state_is_solver_independent() -> CfdResult<()> {
    let (sor, sor_steps, sor_seconds) = channel_with_circle(PoissonSolverKind::Sor).await?;
    let (mg, mg_steps, mg_seconds) = channel_with_circle(PoissonSolverKind::Multigrid).await?;
    let sums = |a: &nalgebra::DMatrix<f64>, b: &nalgebra::DMatrix<f64>| {
        let diff: f64 = a.iter().zip(b.iter()).map(|(x, y)| (x - y) * (x - y)).sum();
        let scale: f64 = b.iter().map(|y| y * y).sum();
        (diff, scale)
    };
    // Velocity: both components against the velocity scale (v alone is
    // small in a channel); pressure against its own RMS (the outlet fixes
    // the level, so the RMS is a scale and not an arbitrary offset).
    let (du, su) = sums(&mg.u, &sor.u);
    let (dv, sv) = sums(&mg.v, &sor.v);
    let (dp, sp) = sums(&mg.p, &sor.p);
    let velocity = ((du + dv) / (su + sv)).sqrt();
    let pressure = (dp / sp).sqrt();
    println!(
        "  outlet channel with circle: SOR {sor_steps} steps {sor_seconds:.1} s, MG {mg_steps} steps \
         {mg_seconds:.1} s; relative RMS differences velocity {velocity:.2e} pressure {pressure:.2e}"
    );
    assert!(
        velocity < 1e-6 && pressure < 1e-6,
        "steady fields differ between inner solvers: velocity {velocity:.3e}, pressure {pressure:.3e}"
    );
    Ok(())
}