"""Bounded SRT/DB uniform-grid wall-bounded shear validation, Horsea only. Uses production Simulation/collision/streaming. Test-only raw population/flag adapter follows tests/integration/test_numerical_gpu.py; no periodic-y API. """ import argparse import hashlib import json import os from pathlib import Path import subprocess import sys from tempfile import TemporaryDirectory import time import numpy as np from les_flow_oracle import H, U, NU, DELTA, TIMES, solve, amplitude, relative HERE = Path(__file__).resolve().parent HEADER = "src/CelerisLab/lbm/kernels/operators/turbulence_smag.cuh" LIMITS = dict( linear_relative_l2=0.003, les_relative_l2=0.005, sgs_increment_relative_error=0.15, oracle_grid_relative=0.00008, mass_drift=0.00005, max_transverse=0.00001, x_invariance=0.000001, initial_relative_l2=0.00001, ) def decode_macros(data, w, cx, cy): """Reject invalid physical populations/density/macros before reporting norms.""" with np.errstate(over="ignore", invalid="ignore", divide="ignore"): f = data.astype(np.float64) + w[:, None, None].astype(np.float64) if not np.isfinite(f).all(): raise RuntimeError("Nonfinite populations") f = f[:, 1:-1, :] rho = f.sum(axis=0) if not np.isfinite(rho).all() or np.any(rho <= 0): raise RuntimeError("Nonpositive or nonfinite fluid density") vx = np.einsum("q,qyx->yx", cx, f) / rho vy = np.einsum("q,qyx->yx", cy, f) / rho profile = vx.mean(axis=1) if not ( np.isfinite(vx).all() and np.isfinite(vy).all() and np.isfinite(profile).all() ): raise RuntimeError("Nonfinite decoded velocity/profile") metrics = ( float(rho.mean()), float(np.max(np.abs(vy))), float(np.max(np.ptp(vx, axis=1))), ) if not np.isfinite(metrics).all(): raise RuntimeError("Nonfinite decoded diagnostics") return profile, *metrics def source_guard(root, variant): manifest = json.loads((HERE / "data/les_flow/source_pins.json").read_text()) for name, sha in manifest["sources"].items(): expected = ( manifest["candidate_header_sha256"] if name == HEADER and variant == "corrected" else sha ) if hashlib.sha256((root / name).read_bytes()).hexdigest() != expected: raise RuntimeError("Source pin mismatch: " + name) return dict( base_commit=manifest["base_commit"], production_origin_commit=manifest["production_origin_commit"], source_experiment=manifest["source_experiment"], checkout_commit=subprocess.check_output( ["git", "-C", str(root), "rev-parse", "HEAD"], text=True ).strip(), source_file_count=len(manifest["sources"]), header_sha256=hashlib.sha256((root / HEADER).read_bytes()).hexdigest(), variant=variant, ) def main(): ap = argparse.ArgumentParser(description=__doc__) ap.add_argument("--root", type=Path, required=True) ap.add_argument("--variant", choices=["baseline", "corrected"], required=True) ap.add_argument("--output", type=Path, required=True) ap.add_argument( "--device", type=int, default=0, help="Logical device in CUDA_VISIBLE_DEVICES" ) ap.add_argument("--run-gpu", action="store_true") args = ap.parse_args() if not args.run_gpu: ap.error("Explicit --run-gpu is required; authorized Horsea executor only") if not os.environ.get("CUDA_VISIBLE_DEVICES"): ap.error("Explicit CUDA_VISIBLE_DEVICES required") if args.output.exists(): ap.error("Output exists; choose a new evidence path") args.output.parent.mkdir(parents=True, exist_ok=True) root = args.root.resolve() pins = source_guard(root, args.variant) # Refuse installed-package shadowing; import only after explicit GPU opt-in. sys.path.insert(0, str(root / "src")) import CelerisLab if Path(CelerisLab.__file__).resolve() != root / "src/CelerisLab/__init__.py": raise RuntimeError("Wrong CelerisLab import") import pycuda.driver as cuda from CelerisLab.simulation import Simulation from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY, FLUID, SOLID, BC_WALL cx = np.asarray(D2Q9_EX) cy = np.asarray(D2Q9_EY) w = np.asarray([4 / 9] + [1 / 9] * 4 + [1 / 36] * 4, dtype=np.float32) h = int(H) nx = 8 ny = h + 2 y = np.arange(h) + 0.5 report = dict( source=pins, package=CelerisLab.__file__, cuda_visible_devices=os.environ["CUDA_VISIBLE_DEVICES"], logical_device=args.device, geometry=dict(nx=nx, ny=ny, halfway_wall_height=H), parameters=dict( nu=NU, U=U, delta_lattice=DELTA, tau0=0.5 + 3 * NU, re_height=U * H / NU, ma=U * np.sqrt(3), ), limits=LIMITS, criterion="standard |S|=sqrt(2 S:S); no baseline-failure masking", cases=[], ) def save(): args.output.write_text(json.dumps(report, indent=2, allow_nan=False) + "\n") started = time.monotonic() for cs in (0.0, 0.16, 1.0): cfg = json.loads((root / "src/CelerisLab/configs/config_lbm.json").read_text()) cfg["grid"].update(nx=nx, ny=ny, nz=1, lattice_model="D2Q9") cfg["physics"].update(viscosity=NU, velocity=0.0, rho=1.0) cfg["method"].update( collision="SRT", streaming="double_buffer", store_precision="FP32", ddf_shifting=True, y_wall_bc="bounce_back", ) cfg["method"]["les"].update(enabled=True, cs=cs, closed_form=True) with TemporaryDirectory(prefix="uniform-les-flow-") as temp: lbm = Path(temp) / "lbm.json" bodies = Path(temp) / "bodies.json" lbm.write_text(json.dumps(cfg)) bodies.write_text('{"objects":[]}') sim = Simulation(str(lbm), str(bodies), device_id=args.device) try: sim.initialize() # Exactly the existing integration-test spatial postcollision contract. # y-wall ghosts stay zero velocity; x edge fluid gives modulo x streaming. ux = np.zeros((ny, nx)) ux[1:-1] = U * np.sin(np.pi * y / H)[:, None] cu = cx[:, None, None] * ux[None, :, :] physical = w[:, None, None] * ( 1 + 3 * cu + 4.5 * cu * cu - 1.5 * ux[None, :, :] ** 2 ) shifted = np.ascontiguousarray( physical - w[:, None, None], dtype=np.float32 ) def install_initial(): with sim.ctx.activate(): cuda.memcpy_htod(sim.field.ddf_gpu, shifted) cuda.memcpy_htod(sim.field.temp_gpu, shifted) sim.stepper._step_count = ( 1 # count metadata; physical elapsed time is 0 ) sim.field.invalidate_host_ddf_cache() flags = sim.field.flag.reshape(ny, nx) flags[:] = FLUID flags[[0, -1], :] = SOLID | BC_WALL sim.field.upload_flags() sim.field.update_params(fx=0.0, fy=0.0, fz=0.0) def observe(): data = np.empty((9, ny, nx), np.float32) with sim.ctx.activate(): cuda.memcpy_dtoh(data, sim.field.ddf_gpu) return decode_macros(data, w, cx, cy) install_initial() observe() # reject invalid state before any advancing launch sim.run(1) # exercise the actual advancing kernel before claiming READY observe() # also require valid density/macros after the warm-up launch install_initial() # restore both (possibly swapped) DB buffers, flags and count initial, rho0, _, _ = observe() initial_error = relative(initial, U * np.sin(np.pi * y / H)) if ( not np.isfinite(initial_error) or initial_error > LIMITS["initial_relative_l2"] ): raise RuntimeError("Initial state adapter failed") with sim.ctx.activate(): device_name = cuda.Context.get_device().name() cuda.Context.synchronize() print( json.dumps( dict( status="READY", cs=cs, device=device_name, variant=args.variant, source=pins, warmup_steps=1, initial_state_restored=True, initial_relative_l2=initial_error, max_steps=TIMES[-1], mode="SRT double_buffer shifted FP32", ) ), flush=True, ) yf, ref = solve(cs, n=256) yc, coarse = solve(cs, n=128) _, wrong = solve(cs * 2 ** (-0.25), n=256) if cs else (yf, ref) case = dict( cs=cs, device=device_name, warmup_steps=1, initial_state_restored=True, configuration=cfg, initial_relative_l2=initial_error, observations=[], ) report["cases"].append(case) elapsed = 0 for target in TIMES: sim.run(target - elapsed) elapsed = target profile, mass, transverse, xspread = observe() standard = np.interp(y, yf, ref[target]) analytic = ( U * np.sin(np.pi * y / H) * np.exp(-NU * (np.pi / H) ** 2 * target) ) oracle_grid = relative(np.interp(yc, yf, ref[target]), coarse[target]) reference = standard if cs else analytic error = relative(profile, reference) obs = dict( elapsed_steps=elapsed, completed_count=sim.stepper.step_count, relative_l2=error, oracle_grid_relative=oracle_grid, mass_drift=abs(mass - rho0), max_transverse=transverse, x_invariance=xspread, amplitude=amplitude(y, profile), reference_amplitude=amplitude(y, reference), wrong_convention_relative_l2=relative( profile, np.interp(y, yf, wrong[target]) ), profile=profile.tolist(), reference_profile=reference.tolist(), ) obs["passed"] = bool( error <= LIMITS["les_relative_l2" if cs else "linear_relative_l2"] and oracle_grid <= LIMITS["oracle_grid_relative"] and abs(mass - rho0) <= LIMITS["mass_drift"] and transverse <= LIMITS["max_transverse"] and xspread <= LIMITS["x_invariance"] ) case["observations"].append(obs) save() print( json.dumps( {k: v for k, v in obs.items() if not k.endswith("profile")} ), flush=True, ) finally: sim.close() # Cancel the common molecular lattice error using the actual Cs=0 run. increments = [] for zero, les in zip( report["cases"][0]["observations"], report["cases"][-1]["observations"] ): measured = zero["amplitude"] - les["amplitude"] expected = zero["reference_amplitude"] - les["reference_amplitude"] increments.append( dict( elapsed_steps=zero["elapsed_steps"], measured=measured, expected=expected, relative_error=abs(measured - expected) / abs(expected), ) ) report.update(sgs_increments=increments, seconds=time.monotonic() - started) # Last two times are the declared signal gate; earliest is transient diagnostic. report["passed"] = bool( all(o["passed"] for c in report["cases"] for o in c["observations"]) and all( v["relative_error"] <= LIMITS["sgs_increment_relative_error"] for v in increments[1:] ) ) save() print( json.dumps( dict(status="DONE", passed=report["passed"], output=str(args.output)) ), flush=True, ) return 0 if report["passed"] else 1 if __name__ == "__main__": raise SystemExit(main())