339 lines
13 KiB
Python
339 lines
13 KiB
Python
"""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
|
|
|
|
|
|
EXPERIMENTS = {
|
|
"historical": "source_pins.json",
|
|
"dev8403956": "source_pins_dev8403956.json",
|
|
"a135252": "source_pins_a135252.json",
|
|
}
|
|
|
|
|
|
def source_guard(root, variant, experiment="historical"):
|
|
manifest_path = HERE / "data/les_flow" / EXPERIMENTS[experiment]
|
|
manifest = json.loads(manifest_path.read_text())
|
|
for name, sha in manifest["sources"].items():
|
|
expected = sha
|
|
if name == HEADER:
|
|
expected = (manifest["candidate_header_sha256"] if variant == "corrected"
|
|
else manifest.get("baseline_header_sha256", sha))
|
|
if hashlib.sha256((root / name).read_bytes()).hexdigest() != expected:
|
|
raise RuntimeError("Source pin mismatch: " + name)
|
|
identity = 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,
|
|
)
|
|
if experiment != "historical":
|
|
identity.update(
|
|
experiment=experiment,
|
|
source_manifest_sha256=hashlib.sha256(manifest_path.read_bytes()).hexdigest(),
|
|
)
|
|
return identity
|
|
|
|
|
|
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(
|
|
"--experiment", choices=tuple(EXPERIMENTS), default="historical",
|
|
help="Explicit source contract; historical pins remain the default",
|
|
)
|
|
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, args.experiment)
|
|
# 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())
|