"""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 (positive Cd along the inflow, positive Cl towards +y). Sampling. One sample every SAMPLE_TSTAR (t* = U t / D): Cd and Cl are the mean over that interval (window exchange ledger refined, per-batch forces uniform); p_diff is the end-of-interval value at the exact surface points. Steady rule (frozen before any run). Over the trailing WINDOW_TSTAR, the range (max - min) of each of Cd, Cl and p_diff must not exceed one tenth of the width of its official interval. Absolute ranges avoid dividing by the small Cl. If the rule fails at TSTAR_MAX the run reports "not steady" and gives no interval verdict. """ 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 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_RANGE = {k: (hi - lo) / 10 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] ranges = {k: max(r[k] for r in window) - min(r[k] for r in window) for k in BANDS} steady = all(ranges[k] <= STEADY_RANGE[k] for k in BANDS) final = {k: rows[-1][k] for k in BANDS} result = dict(window_samples=len(window), ranges=ranges, steady_limits=STEADY_RANGE, steady=steady, final=final, bands=BANDS, deviations=DEVIATIONS) 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_RANGE, 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 with official reference intervals; not a full benchmark certification") def main(): parser = argparse.ArgumentParser(description=__doc__.splitlines()[0]) parser.add_argument("--case", choices=["both", "refined", "uniform"], default="both") 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()