CelerisLab/tests/validation/run_dfg_re20.py

"""Steady DFG 2D-1 (Re 20): refined and uniform runs against the official reference intervals.

Not a full benchmark certification. Declared deviations from the official setup:
pressure-NEQ outlet (``neq_extrap``) instead of do-nothing, local Zou-He inlet.
Geometry, inflow profile, initial state, refined layout and observation reuse the
pinned Re100 runner (``run_refined_dfg``) and ``re100_helpers`` unchanged; only the
lattice viscosity is set for Re 20 (nu = U_lattice * D / 20 at the root level).

Units. The helpers use scaled DFG units with mean inflow speed 1, so
Cd = 2 Fx / (rho U^2 D) and Cl = 2 Fy / (rho U^2 D) are scale-free, and
``pressure_probes`` returns dp / (rho U^2). The official p_diff = p(0.15, 0.2) -
p(0.25, 0.2) uses U = 0.2 m/s, rho = 1, so p_diff = dp * 0.04. Both are reported;
only p_diff is compared with the official interval. Forces use the existing sign
convention: F is the force of the fluid on the cylinder, i.e. the refined raw
exchange impulse (window ledger, ``total.impulse``) or the uniform ``read_force``,
with x along the inflow and y up from the bottom wall. Cd must be positive, and the
official Cl is positive (centre 0.20 below the centreline 0.205), so a flipped sign
shows up as Cl outside its interval.

Sampling. One sample every SAMPLE_TSTAR (t* = U t / D; 0.25 s), taken at the end of
a completed root step with all levels synchronized; the step index is recorded. Cd
and Cl are means over the sample interval (window exchange ledger refined, per-batch
forces uniform), which filters acoustic and step-parity oscillation; p_diff is the
end-of-interval value at the exact surface points.

Steady rule (our engineering choice, frozen before any run; not official). Over the
trailing WINDOW_TSTAR (t* 50..60 inclusive, 21 samples), for each of Cd, Cl and
p_diff: range (max - min) <= 10 % and |least-squares slope * window| <= 2 % of the
width of its official interval. Absolute limits, so the small Cl is never a divisor;
the slope term bounds slow monotone drift that the range alone would pass. If the
rule fails at TSTAR_MAX the run reports "not steady" and gives no interval verdict;
limits are never relaxed after a run. The verdict uses the window mean, and the
distance of each mean to its nearest interval edge is reported (negative outside).
"""

import argparse
import json
import math
import os
from pathlib import Path
import sys
import time
from tempfile import TemporaryDirectory

import numpy as np

HERE = Path(__file__).resolve().parent
sys.path.insert(0, str(HERE))
import run_refined_dfg as rr  # noqa: E402

RE = 20.0
U_PHYS = 0.2                   # m/s mean inflow; rho = 1
TSTAR_MAX = 60.0               # 30 s physical
SAMPLE_TSTAR = 0.5             # 0.25 s
WINDOW_TSTAR = 10.0            # trailing 5 s, 21 samples
# Schaefer & Turek (1996), 2D-1 lower/upper bounds.
BANDS = dict(cd=(5.57, 5.59), cl=(0.0104, 0.0110), p_diff=(0.1172, 0.1176))
STEADY_RANGE = {k: (hi - lo) / 10 for k, (lo, hi) in BANDS.items()}
STEADY_SLOPE = {k: (hi - lo) / 50 for k, (lo, hi) in BANDS.items()}
DEVIATIONS = ("pressure-NEQ outlet (neq_extrap), not do-nothing", "local Zou-He inlet")

_re100_config = rr.config


def config(h, diameter, collision="SRT", jet_scenario="none"):
    """The pinned Re100 configuration with the lattice viscosity set for Re 20."""
    if jet_scenario != "none":
        raise ValueError("DFG 2D-1 runs have no Jets")
    cfg, body, p = _re100_config(h, diameter, collision, jet_scenario)
    nu = p["mean_speed"] * p["diameter"] / RE
    cfg["physics"]["viscosity"] = nu
    return cfg, body, dict(p, viscosity=nu, re=RE)


def lattice_parameters(diameter, speed=rr.SPEED):
    nu = speed * diameter / RE
    return dict(viscosity=nu, tau=0.5 + 3 * nu, steps=math.ceil(TSTAR_MAX * diameter / speed),
                interval=round(SAMPLE_TSTAR * diameter / speed))


def sample(impulse, n, diameter, dp):
    den = rr.SPEED**2 * diameter / 2
    return dict(cd=float(impulse[0] / n / den), cl=float(impulse[1] / n / den),
                dp_scaled=float(dp), p_diff=float(dp) * U_PHYS**2)


def summarize(rows):
    """Steady rule over the trailing window, then the interval verdict (or none)."""
    if not rows or rows[-1]["tstar"] < TSTAR_MAX - 1e-9:
        return dict(status="INCOMPLETE", steady=False, verdict=None)
    window = [r for r in rows if r["tstar"] >= TSTAR_MAX - WINDOW_TSTAR - 1e-9]
    t = np.array([r["tstar"] for r in window])
    ranges = {k: max(r[k] for r in window) - min(r[k] for r in window) for k in BANDS}
    drifts = {k: abs(np.polyfit(t, [r[k] for r in window], 1)[0]) * WINDOW_TSTAR for k in BANDS}
    steady = all(ranges[k] <= STEADY_RANGE[k] and drifts[k] <= STEADY_SLOPE[k] for k in BANDS)
    final = {k: float(np.mean([r[k] for r in window])) for k in BANDS}
    edge = {k: min(final[k] - BANDS[k][0], BANDS[k][1] - final[k]) for k in BANDS}
    result = dict(window_samples=len(window), ranges=ranges, drifts=drifts,
                  range_limits=STEADY_RANGE, slope_limits=STEADY_SLOPE, steady=steady,
                  window_mean=final, edge_distance=edge, bands=BANDS,
                  deviations=DEVIATIONS, cd_positive=final["cd"] > 0)
    if not steady:
        return dict(result, status="NOT_STEADY", verdict=None)
    inside = {k: BANDS[k][0] <= final[k] <= BANDS[k][1] for k in BANDS}
    return dict(result, status="WITHIN_INTERVALS" if all(inside.values()) else "OUTSIDE_INTERVALS",
                verdict=inside)


def healthy(row):
    return row["max_mach"] <= 0.2 and row["rho_min"] >= 0.5 and row["rho_max"] <= 1.5


def run_case(a, h, case):
    from contextlib import ExitStack
    from CelerisLab.cuda.context import CudaContext
    from CelerisLab.lbm.descriptors import FLUID
    from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
    from CelerisLab.simulation import Simulation

    refined = case == "refined"
    cfgdict, bodydict, p = config(h, 40 if refined else 80, a.collision)
    lp = lattice_parameters(p["diameter"])
    rows = []
    with ExitStack() as stack:
        if refined:
            cfg, pair, body, host = rr.host_geometry(a.source, h, a.collision)
            context = stack.enter_context(CudaContext(0))
            backend = stack.enter_context(FixedPatchStepper(
                cfg, context=context, coarse_boundary="channel", refinements=[rr.PLACEMENTS["wide"]],
                level_blocks=[rr.BASE_BLOCK, rr.FINE_BLOCK], level_objects=[(), (body,)],
                observe_raw_exchange=True, observe_window_exchange=True))
            from CelerisLab.lbm.channel import channel_masks
            from CelerisLab.lbm.refinement import ACTIVE
            from CelerisLab.lbm.descriptors import SOLID
            masks = [channel_masks(pair.coarse_shape)[0] & ~pair.covered,
                     (pair.fine_kind == ACTIVE)
                     & ((host.flags.reshape(pair.fine_shape) & SOLID) == 0)]
            backend.initialize(rr.initial_state(pair.coarse_shape, (0.0, 0.0), 1.0, root=True),
                               rr.initial_state(pair.fine_shape, pair.fine_origin, 0.5))
            probe_p = dict(p, center=[(80.5 - pair.fine_origin[0]) / 0.5 - 0.5,
                                      (81.0 - pair.fine_origin[1]) / 0.5 - 0.5], radius=40.0)
        else:
            tmp = Path(stack.enter_context(TemporaryDirectory()))
            (tmp / "lbm.json").write_text(json.dumps(cfgdict))
            (tmp / "body.json").write_text(json.dumps(bodydict))
            sim = stack.enter_context(Simulation(str(tmp / "lbm.json"), str(tmp / "body.json"),
                                                 device_id=0))
            sim.initialize()
            probe_p = p
            fluid = (sim.field.flag.reshape(p["ny"], p["nx"]) & FLUID) != 0
        begin = 0
        while begin < lp["steps"]:
            rr.check_deadline(a.deadline)
            n = min(lp["interval"], lp["steps"] - begin)
            if refined:
                impulse = rr.advance_refined_interval(backend, n, a.deadline, "window")
                fields = [rr.macros(g) for g in backend.download()]
                rho, ux, uy = fields[1]
                fluid = masks[1]
                vals = [tuple(v[m] for v in f) for m, f in zip(masks, fields)]
            else:
                impulse, done = np.zeros(2), 0
                while done < n:
                    batch = min(100, n - done)
                    sim.run(batch)
                    impulse += sim.read_force(0).astype(np.float64) * batch
                    done += batch
                m = sim.get_macroscopic()
                rho, ux, uy = [np.asarray(m[k], np.float64).reshape(p["ny"], p["nx"])
                               for k in ("rho", "ux", "uy")]
                vals = [(rho[fluid], ux[fluid], uy[fluid])]
            begin += n
            row = dict(step=begin, tstar=begin * rr.SPEED / p["diameter"],
                       **sample(impulse, n, p["diameter"],
                                h.pressure_probes(rho, fluid, probe_p)["dp"]),
                       rho_min=float(min(v[0].min() for v in vals)),
                       rho_max=float(max(v[0].max() for v in vals)),
                       max_mach=float(max(np.sqrt(3 * (v[1]**2 + v[2]**2)).max() for v in vals)))
            rr.ensure_finite(row)
            rows.append(row)
            if not healthy(row):
                break
    summary = summarize(rows) if healthy(rows[-1]) else dict(status="UNHEALTHY", steady=False,
                                                             verdict=None)
    return dict(case=case, collision=a.collision, lattice=lp, rows=rows, **summary)


def plan(collision="SRT"):
    return dict(re=RE, collision=collision, tstar_max=TSTAR_MAX, sample_tstar=SAMPLE_TSTAR,
                window_tstar=WINDOW_TSTAR, bands=BANDS, range_limits=STEADY_RANGE, slope_limits=STEADY_SLOPE,
                deviations=DEVIATIONS, pressure_scale=U_PHYS**2,
                refined_root=lattice_parameters(40), refined_fine_tau=0.5 + 3 * 2 * rr.SPEED * 40 / RE,
                uniform=lattice_parameters(80),
                label="compared against the official reference bands (non-canonical neq_extrap outlet, zou_he_local inlet); not a full benchmark certification")


def main():
    parser = argparse.ArgumentParser(description=__doc__.splitlines()[0])
    parser.add_argument("--case", choices=["refined", "uniform", "both"], default="refined")
    parser.add_argument("--collision", choices=["SRT", "MRT"], default="SRT")
    parser.add_argument("--source", type=Path, default=rr.DEFAULT_SOURCE)
    parser.add_argument("--output", type=Path, required=True)
    parser.add_argument("--plan", action="store_true", help="Write the frozen plan; no CUDA")
    parser.add_argument("--max-wall-seconds", type=float)
    a = parser.parse_args()
    a.source = a.source.resolve()
    a.output.mkdir(parents=True, exist_ok=False)
    (a.output / "plan.json").write_text(json.dumps(plan(a.collision), indent=2, allow_nan=False))
    if a.plan:
        return
    if not a.max_wall_seconds or a.max_wall_seconds <= 0 or not os.environ.get("CUDA_VISIBLE_DEVICES"):
        raise SystemExit("GPU run needs --max-wall-seconds > 0 and explicit CUDA_VISIBLE_DEVICES")
    h, _ = rr.helpers(a.source)
    rr.config = config  # host_geometry builds its root configuration through this name
    a.deadline = time.monotonic() + a.max_wall_seconds
    cases = ["refined", "uniform"] if a.case == "both" else [a.case]
    results = [run_case(a, h, case) for case in cases]
    (a.output / "result.json").write_text(json.dumps(results, indent=2, allow_nan=False))


if __name__ == "__main__":
    main()