test: integrate explicitly repinned uniform LES shear diagnostic
This commit is contained in:
TRON committed
1 parent
cd3a657f9d
commit
af984eab37
7 files changed
+755
No files matched your search
@@ -0,0 +1,321 @@
|
||||
"""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())
|
||||
Reference in new issue
Block a user