"""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.5 s), taken at the end of a completed root step with all levels synchronized; the step index is recorded. Cd and Cl are the mean over the sample interval (smooths step-parity/acoustic oscillation; equal to the instantaneous value at steady state); p_diff is the end-of-interval value at the exact surface points. Steady rule (frozen before any run). Over the trailing WINDOW_TSTAR (last 5 s, 11 samples), both the range (max - min) and |least-squares slope * window| of each of Cd, Cl and p_diff must not exceed 1 % of the width of its official interval. These are absolute limits, so the small Cl is never a divisor, and the slope term catches slow monotone drift. If the rule fails at TSTAR_MAX the run reports "not steady" and gives no interval verdict; the limits are never relaxed after a run. """ 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 = 1.0 # 0.5 s WINDOW_TSTAR = 10.0 # trailing 5 s # 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_LIMIT = {k: (hi - lo) / 100 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_LIMIT[k] and drifts[k] <= STEADY_LIMIT[k] for k in BANDS) final = {k: rows[-1][k] for k in BANDS} result = dict(window_samples=len(window), ranges=ranges, drifts=drifts, steady_limits=STEADY_LIMIT, steady=steady, final=final, 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, rr.PLACEMENTS["wide"], context=context, coarse_boundary="channel", base_block=rr.BASE_BLOCK, fine_block=rr.FINE_BLOCK, fine_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, steady_limits=STEADY_LIMIT, 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()