"""DFG Re100 periodic reference comparison with declared pressure-NEQ outlet, not exact do-nothing.

CPU helpers import without CUDA. GPU execution verifies the explicitly pinned
AMR-channel experimental source snapshot. Historical d79 evidence remains tied to
its original commits. Default pilot is too short to certify periodic convergence.
"""

import argparse
import csv
import hashlib
import json
import math
import os
from pathlib import Path
import time
from tempfile import TemporaryDirectory
import numpy as np

ROOT = Path(__file__).resolve().parents[2]
PIN = "47b05bc8abfc957fe876b0c63bdb386e3c6cbaa5"


def geometry(n):
    if isinstance(n, bool) or n < 10 or n % 10:
        raise ValueError("Diameter must be a positive multiple of 10, at least 10")
    return dict(
        diameter=n,
        nx=22 * n + 1,
        ny=41 * n // 10 + 2,
        center=[2 * n, 2 * n + 0.5],
        radius=n / 2,
        dx=0.1 / n,
        bottom=0.5,
        top=41 * n / 10 + 0.5,
        inlet_source_x=0.0,
        outlet_source_x=22.0 * n,
    )


def parameters(n, speed):
    g = geometry(n)
    if not 0 < speed <= 0.04:
        raise ValueError("Require 0 < lattice mean speed <= .04")
    return dict(
        **g,
        mean_speed=speed,
        peak_speed=1.5 * speed,
        viscosity=speed * n / 100,
        dt=speed * g["dx"],
        peak_mach=1.5 * speed * math.sqrt(3),
        re=100.0,
    )


def make_config(n, speed, collision, outlet, inlet="regularized"):
    p = parameters(n, speed)
    cfg = json.loads((ROOT / "src/CelerisLab/configs/config_lbm.json").read_text())
    cfg["grid"].update(nx=p["nx"], ny=p["ny"], nz=1, lattice_model="D2Q9")
    cfg["physics"].update(velocity=speed, viscosity=p["viscosity"], rho=1.0)
    cfg["method"].update(
        collision=collision,
        streaming="double_buffer",
        store_precision="FP32",
        ddf_shifting=True,
        y_wall_bc="bounce_back",
    )
    cfg["method"]["les"]["enabled"] = False
    cfg["method"]["inlet"].update(profile="parabolic", scheme=inlet)
    cfg["method"]["outlet"].update(mode=outlet, backflow_clamp=False, blend_alpha=0.7)
    body = {
        "objects": [
            dict(type="cylinder", center=p["center"], radius=p["radius"], omega=0.0)
        ]
    }
    return cfg, body, p


def bilinear_fluid(a, fluid, x, y):
    ix, iy = math.floor(x), math.floor(y)
    tx, ty = x - ix, y - iy
    result = 0.0
    for xx, yy, w in (
        (ix, iy, (1 - tx) * (1 - ty)),
        (ix + 1, iy, tx * (1 - ty)),
        (ix, iy + 1, (1 - tx) * ty),
        (ix + 1, iy + 1, tx * ty),
    ):
        if w == 0:
            continue
        if not (0 <= yy < a.shape[0] and 0 <= xx < a.shape[1]) or not fluid[yy, xx]:
            raise ValueError("Pressure interpolation touched a non-fluid/outside node")
        result += w * float(a[yy, xx])
    return result


def pressure_probes(rho, fluid, p):
    # Linear extrapolation from distances h and 2h outside the cylinder to
    # each exact wall probe, after bilinear interpolation using only fluid.
    cx, cy = p["center"]
    r = p["radius"]
    values = []
    for sign in (-1, 1):
        wall = cx + sign * r
        near = bilinear_fluid(rho, fluid, wall + sign, cy)
        far = bilinear_fluid(rho, fluid, wall + 2 * sign, cy)
        values.append(2 * near - far)
    scale = (1.0 / p["mean_speed"]) ** 2 / 3
    return dict(
        rho_front=values[0], rho_back=values[1], dp=(values[0] - values[1]) * scale
    )


def outlet_diagnostic(rho, ux, uy, p):
    # Extrapolate fluid values to the declared x=NX-1 source plane. Diagnostic
    # excludes two rows beside each wall; not an exact discrete BC identity.
    rows = slice(3, -3)
    i = -2
    dudx = (5 * ux[rows, i] - 8 * ux[rows, i - 1] + 3 * ux[rows, i - 2]) / 2
    dvdx = (5 * uy[rows, i] - 8 * uy[rows, i - 1] + 3 * uy[rows, i - 2]) / 2
    rho_plane = 2 * rho[rows, i] - rho[rows, i - 1]
    pressure = (rho_plane - 1) / 3
    scale = p["mean_speed"] ** 2
    residual = np.stack((p["viscosity"] * dudx - pressure, p["viscosity"] * dvdx))
    return dict(
        rms_over_rho_U2=float(np.sqrt(np.mean(residual**2)) / scale),
        max_over_rho_U2=float(np.max(abs(residual)) / scale),
        definition="fluid-side one-sided estimate of nu*d_x(u)-p*n at declared outlet",
    )