Files
CelerisLab/tests/validation/run_jet_circle_cv.py
T
ClaudeandClaude Opus 5.5 8064834c37 test: tighten historical-pin checks per review of 941eb28b
- test_refined_dfg::test_source_pin runs its mismatch negative controls
  first (now including window mode) so missing history cannot skip them.
- test_re100: renamed to test_manifest_matches_its_pinned_commit, plus an
  explicit record that today's src differs from the pin.
- Structural manifest checks run before any history-dependent check.
- Real tamper test: a flipped digest raises via verify_manifest, and a
  path absent at the commit raises ValueError.
- A blobless clone (tree present, blob missing) is reported as unavailable
  history, not as a missing path.
- All history skips go through history_unavailable(): fatal when
  CELERIS_REQUIRE_HISTORY=1.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-20 18:28:09 +00:00

677 lines
26 KiB
Python

#!/usr/bin/env python3
"""Bounded jet-on periodic cylinder: actual AMR Jet API vs blocked validation adapter."""
import argparse
import hashlib
import subprocess
import json
import os
from pathlib import Path
import sys
import time
import numpy as np
from source_provenance import checkout_provenance
PIN = "47b05bc8abfc957fe876b0c63bdb386e3c6cbaa5"
HERE = Path(__file__).resolve().parent
CX = np.array([0, 1, -1, 0, 0, 1, -1, 1, -1])
CY = np.array([0, 0, 0, 1, -1, 1, -1, -1, 1])
W = np.array([4 / 9] + [1 / 9] * 4 + [1 / 36] * 4)
def case(scale=1, re=100, steps=800):
u = 0.03 / scale
d = 24.084 * scale
return dict(
scale=scale,
re_d=re,
beta=0.1,
nx=192 * scale,
ny=128 * scale,
diameter=d,
center=[96 * scale, 64 * scale],
u=u,
jet_speed=0.1 * u,
nu=u * d / re,
tau_root=0.5 + 3 * u * d / re,
tau_fine=0.5 + 6 * u * d / re,
bounds=[48 * scale, 32 * scale, 144 * scale, 96 * scale],
contours=[
[72 * scale, 44 * scale, 120 * scale, 84 * scale],
[64 * scale, 40 * scale, 128 * scale, 88 * scale],
],
steps=steps * scale**2,
window=100 * scale**2,
dx_root=1 / scale,
dt_root=1 / scale**2,
fine_dx_root_units=0.5,
fine_dt_root_units=0.5,
mach=np.sqrt(3) * u,
final_t_star=steps * scale**2 * u / d,
)
def state(shape, u):
cu = CX[:, None, None] * u
return (
np.broadcast_to(
W[:, None, None] * (1 + 3 * cu + 4.5 * cu**2 - 1.5 * u * u), shape
)
.astype(np.float32)
.copy()
)
def moments(g):
f = np.asarray(g, np.float64)
rho = f.sum(axis=0)
momentum = np.stack([np.einsum("q,qyx->yx", v, f) for v in (CX, CY)])
return rho, momentum, momentum / rho
def cv(g, fluid, origin, rectangle, nu, h=0.5):
"""Cell-centred inventory and midpoint face quadrature, physical root units.
Cauchy stress from centred velocity gradients, NOT postcollision non-equilibrium.
Rectangle faces lie on cell edges; fluid mask excludes the staircase cylinder.
"""
return _cv_quantities(*moments(g), fluid, origin, rectangle, nu, h)
def cv_from_macros(rho, ux, uy, fluid, origin, rectangle, nu, h=0.5):
"""Same CV calculation from public macros, in one common length/time frame.
``origin`` is the lower cell edge, ``h`` the cell width, and ``nu`` uses
that frame. Returns [mass, Px, Py] inventory and outward total flux.
Crop fields/mask together and translate origin before calling if desired.
No boundary-exchange observation enters this independent estimator.
"""
rho = np.asarray(rho, np.float64)
if rho.ndim != 2 or any(np.shape(v) != rho.shape for v in (ux, uy, fluid)):
raise ValueError("CV macros and fluid mask must have matching 2D shapes")
vel = np.asarray([ux, uy], np.float64)
return _cv_quantities(rho, rho[None] * vel, vel, fluid, origin, rectangle, nu, h)
def _cv_quantities(rho, mom, vel, fluid, origin, rectangle, nu, h):
x0, y0, x1, y1 = np.rint((np.array(rectangle) - np.tile(origin, 2)) / h).astype(int)
np.testing.assert_allclose(
np.array([x0, y0, x1, y1]) * h + np.tile(origin, 2), rectangle, atol=1e-12
)
ny, nx = fluid.shape
if not (2 <= x0 < x1 <= nx - 2 and 2 <= y0 < y1 <= ny - 2):
raise ValueError("Contour needs two fluid-node stencil layers")
# Every face and its velocity-gradient stencil must be active fluid.
for ys, xs in [
(slice(y0 - 2, y1 + 2), slice(x0 - 2, x0 + 2)),
(slice(y0 - 2, y1 + 2), slice(x1 - 2, x1 + 2)),
(slice(y0 - 2, y0 + 2), slice(x0 - 2, x1 + 2)),
(slice(y1 - 2, y1 + 2), slice(x0 - 2, x1 + 2)),
]:
if not fluid[ys, xs].all():
raise ValueError("Contour stencil intersects solid/inactive region")
use = fluid[y0:y1, x0:x1]
inventory = (
np.array(
[
rho[y0:y1, x0:x1][use].sum(),
mom[0, y0:y1, x0:x1][use].sum(),
mom[1, y0:y1, x0:x1][use].sum(),
]
)
* h
* h
)
duy, dux = np.gradient(vel[0], h)
dvy, dvx = np.gradient(vel[1], h)
mu = rho * nu
# Subtract constant reference pressure, whose closed-contour integral is zero.
p = (rho - 1) / 3
tensor = np.array(
[
[
rho * vel[0] ** 2 + p - 2 * mu * dux,
rho * vel[0] * vel[1] - mu * (duy + dvx),
],
[
rho * vel[0] * vel[1] - mu * (duy + dvx),
rho * vel[1] ** 2 + p - 2 * mu * dvy,
],
]
)
transport = np.concatenate(
[mom[None], tensor], axis=0
) # [mass,Px,Py] x normal axis
flux = np.zeros(3)
for i, sign in [(x0, -1), (x1, 1)]:
flux += (
sign
* 0.5
* (transport[:, 0, y0:y1, i - 1] + transport[:, 0, y0:y1, i]).sum(axis=1)
* h
)
for j, sign in [(y0, -1), (y1, 1)]:
flux += (
sign
* 0.5
* (transport[:, 1, j - 1, x0:x1] + transport[:, 1, j, x0:x1]).sum(axis=1)
* h
)
if not np.isfinite(np.r_[inventory, flux]).all():
raise RuntimeError("Nonfinite CV quantities")
return inventory, flux
def aperture_probe(fluid, origin, center, radius, h=0.5):
"""Independent true-arc quadrature with nearest active-fluid sampling.
This is a first-order near-wall diagnostic, not certified wall reconstruction.
No target velocity or raw-population ledger is used to measure the flux.
"""
z, w = np.polynomial.legendre.leggauss(64)
result = []
for angle in (np.pi / 2, 3 * np.pi / 2):
theta = angle + np.pi / 12 * z
er = np.column_stack((np.cos(theta), np.sin(theta)))
points = np.asarray(center) + radius * er
selected = []
distance = []
for point in points:
ix, iy = np.rint((point - np.asarray(origin)) / h - 0.5).astype(int)
candidates = [
(y, x)
for y in range(iy - 3, iy + 4)
for x in range(ix - 3, ix + 4)
if 0 <= y < fluid.shape[0] and 0 <= x < fluid.shape[1] and fluid[y, x]
]
if not candidates:
raise ValueError("No active near-wall sampling node")
yy, xx = np.array(candidates).T
dist = np.linalg.norm(
np.column_stack((xx, yy)) * h + np.asarray(origin) + h / 2 - point,
axis=1,
)
k = np.argmin(dist)
selected.append(candidates[k])
distance.append(dist[k])
yy, xx = np.array(selected).T
result.append((yy, xx, er, w * radius * np.pi / 12, float(max(distance))))
return result
def measured_aperture_flux(g, probe):
_, momentum, _ = moments(g)
return np.array(
[
np.sum(np.einsum("in,ni->n", momentum[:, y, x], normal) * ds)
for y, x, normal, ds, _ in probe
]
)
def verify(repo):
manifest = json.loads(
(HERE / "data/jet_circle_cv/source_manifest.json").read_text()
)
if manifest["revision"] != PIN:
raise RuntimeError("Wrong pinned revision")
for path, sha in manifest["files"].items():
if hashlib.sha256((repo / path).read_bytes()).hexdigest() != sha:
raise RuntimeError("Source hash mismatch: " + path)
return manifest
def verify_historical(repo):
"""The pinned manifest describes the PIN tree (read with git, not the working tree)."""
from source_provenance import verify_manifest
manifest = json.loads((HERE / "data/jet_circle_cv/source_manifest.json").read_text())
if manifest["revision"] != PIN:
raise RuntimeError("Wrong pinned revision")
verify_manifest(repo, PIN, manifest["files"])
return manifest
def current_checkout(repo):
"""Provenance for a run on the current checkout; never labelled as the pinned source."""
files = subprocess.check_output(["git", "-C", str(repo), "ls-files", "src"],
text=True).splitlines()
return dict(source_mode="current_checkout",
source_capture="current checkout, not the frozen pinned qualification",
source_revision=checkout_provenance(repo)["checkout_revision"],
src_tree=subprocess.check_output(["git", "-C", str(repo), "rev-parse",
"HEAD:src"], text=True).strip(),
src_sha256={f: hashlib.sha256((repo / f).read_bytes()).hexdigest()
for f in files})
def run_case(c, ctx, output, provenance, continuation=False):
from CelerisLab.body.geometry.circle import CircleGeometry
from CelerisLab.body.objects import SimObject
from CelerisLab.body.jets import Jet, Sector, Radial
from CelerisLab.config import LBMConfig
from CelerisLab.lbm.blocks import BlockLayout2D
from CelerisLab.lbm.block_stepper import PeriodicBlockStepper
from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY
import pycuda.driver as cuda
np.testing.assert_array_equal(CX, D2Q9_EX)
np.testing.assert_array_equal(CY, D2Q9_EY)
nx, ny, u, nu = c["nx"], c["ny"], c["u"], c["nu"]
center = np.array(c["center"], float)
radius = c["diameter"] / 2
cfg = LBMConfig(
nx=nx,
ny=ny,
viscosity=nu,
velocity=0.0,
collision="SRT",
streaming="double_buffer",
store_precision="FP32",
ddf_shifting=False,
les_enabled=False,
)
fine_cfg = LBMConfig(
nx=2 * nx,
ny=2 * ny,
viscosity=2 * nu,
velocity=0.0,
collision="SRT",
streaming="double_buffer",
store_precision="FP32",
ddf_shifting=False,
les_enabled=False,
)
obj = SimObject(
0, CircleGeometry(*center, radius), center=tuple(center), radius=radius
)
fc = 2 * center - 0.5
fobj = SimObject(
0, CircleGeometry(*fc, 2 * radius), center=tuple(fc), radius=2 * radius
)
jets = [
Jet(
Sector(tuple(center), angle, np.pi / 12, radius + 0.1),
Radial(tuple(center), inward=False),
"uniform",
sign * c["jet_speed"],
)
for angle, sign in [(np.pi / 2, 1), (3 * np.pi / 2, -1)]
]
layout = BlockLayout2D.build(2 * ny, 2 * nx, 8, 8)
with (
PeriodicBlockStepper(fine_cfg, layout, context=ctx, objects=[fobj]) as uniform,
FixedPatchStepper(
cfg,
c["bounds"],
context=ctx,
base_block=(8, 8),
fine_block=(8, 8),
fine_objects=[obj],
jets=jets,
observe_raw_exchange=True,
) as amr,
):
b = amr._level_boundaries[1]
count = len(b.tables["cell"])
if np.any(b.tables["q"] <= 0.010001) or np.any(b.tables["fallback"]):
raise RuntimeError("Pilot requires unclamped q and no fallback")
links = fobj.geometry.build_cut_links(2 * nx, 2 * ny)
np.testing.assert_array_equal([v.dir for v in links], b.tables["dir"])
np.testing.assert_allclose(
[v.q for v in links], b.tables["q"], atol=2e-7, rtol=0
)
np.testing.assert_array_equal([v.fallback for v in links], b.tables["fallback"])
np.testing.assert_allclose(
np.array([[v.hit_x, v.hit_y] for v in links]) * 0.5 + 0.25,
b.host.curved_hits,
atol=2e-12,
rtol=0,
)
# Same pre-existing optional-input adapter as test_raw_exchange_gpu.py.
velocity = uniform._upload(b.tables["surface_velocity"])
rx = uniform._upload(b.tables["rx"])
ry = uniform._upload(b.tables["ry"])
raw = uniform._allocate(count * 16)
original = uniform._curved
samples = []
def boundary(*args, **kw):
original(*args[:7], velocity, rx, ry, raw, **kw)
uniform._stream.synchronize()
a = np.empty((count, 4), np.float32)
cuda.memcpy_dtoh(a, raw)
samples.append(a)
uniform._curved = boundary
def initialize():
uniform.initialize(state((9, 2 * ny, 2 * nx), u))
amr.initialize(*(state(l.shape, u) for l in amr._levels))
initialize()
# Compile/launch all stepping and diagnostic paths, then reset physical time.
amr.step()
uniform.step(2)
uniform.download()
amr.download()
amr.read_raw_exchange(0)
initialize()
samples.clear()
jet_id = b.jet_cache[1]
masks = {
"wall": jet_id < 0,
"aperture": jet_id >= 0,
"total": np.ones(count, bool),
"blowing": jet_id == 0,
"suction": jet_id == 1,
}
if not masks["blowing"].any() or not masks["suction"].any():
raise RuntimeError("Both apertures must resolve cut links")
fluid_u = ~uniform.boundary["solid"]
fluid_a = amr._leaf[1]
origin_a = np.array(amr.patch.fine_origin)
probes = {
"uniform_adapter": aperture_probe(fluid_u, (0.0, 0.0), center, radius),
"amr_api": aperture_probe(fluid_a, origin_a, center, radius),
}
recorder = None
if continuation:
from jet_cv_measurements import MeasurementRecorder, sample_steps
recorder = MeasurementRecorder(
c, output, fluid_u, amr._leaf, origin_a, CX, CY
)
scheduled = set(sample_steps(c).tolist())
else:
scheduled = set(range(c["steps"] + 1))
cumulative_raw = {
name: {k: np.zeros(4) for k in masks}
for name in ("uniform_adapter", "amr_api")
}
def capture(step):
ug = uniform.download()
levels = amr.download()
ag = levels[1]
data = {}
aperture = {}
for name, g, fluid, origin in [
("uniform_adapter", ug, fluid_u, (0.0, 0.0)),
("amr_api", ag, fluid_a, origin_a),
]:
data[name] = [cv(g, fluid, origin, r, nu) for r in c["contours"]]
aperture[name] = measured_aperture_flux(g, probes[name])
rho, _, _ = moments(g)
if not np.isfinite(rho[fluid]).all() or rho[fluid].min() <= 0:
raise RuntimeError("Nonfinite/nonpositive fluid density")
ox, oy = np.rint(2 * origin_a).astype(int)
ay, ax = fluid_a.shape
_, _, uv = moments(ug[:, oy : oy + ay, ox : ox + ax])
_, _, av = moments(ag)
error = float(
np.linalg.norm((av - uv)[:, fluid_a])
/ max(np.linalg.norm(uv[:, fluid_a]), 1e-30)
)
if recorder is not None:
recorder.capture(step, ug, levels, data, aperture, error)
recorder.records[-1]["raw_cumulative"] = np.array(
[
[cumulative_raw[name][k].copy() for k in masks]
for name in ("uniform_adapter", "amr_api")
]
)
return data, error, aperture
previous, _, previous_aperture = capture(0)
previous_sample_step = 0
initial = previous
integrated = {name: np.zeros((2, 3)) for name in previous}
integrated_aperture = {name: np.zeros(2) for name in previous}
accumulated = {name: {k: np.zeros(4) for k in masks} for name in previous}
rows = []
start = time.monotonic()
print(
json.dumps(
dict(
status="READY",
case=c,
links=count,
actual_amr_jet_api=True,
uniform_reference="private optional-input validation adapter",
q_range=[float(b.tables["q"].min()), float(b.tables["q"].max())],
fallback_count=int(np.count_nonzero(b.tables["fallback"])),
)
),
flush=True,
)
for step in range(1, c["steps"] + 1):
if time.monotonic() - start > c.get("post_ready_budget_seconds", 1800):
raise RuntimeError(
"Declared post-READY budget exceeded; partial report retained"
)
samples.clear()
amr.step()
uniform.step(2)
if len(samples) != 2:
raise RuntimeError("Uniform interval must contain two fine hops")
ur = sum(a.astype(float) * [0.25, 0.25, 0.125, 0.25] for a in samples)
ar = np.empty((count, 4), np.float64)
with ctx.activate():
cuda.memcpy_dtoh(ar, b.interval_integral)
if not np.isfinite(ar).all() or not np.isfinite(ur).all():
raise RuntimeError("Nonfinite raw exchange in completed base interval")
obs = amr.read_raw_exchange(0)
if obs["time_start"] != step - 1 or obs["time_end"] != step:
raise RuntimeError("Raw interval time mismatch")
np.testing.assert_allclose(
ar.sum(axis=0)[:2], obs["total"]["impulse"], atol=1e-13, rtol=1e-13
)
for name, rr in [("uniform_adapter", ur), ("amr_api", ar)]:
for k, mask in masks.items():
subtotal = rr[mask].sum(axis=0)
accumulated[name][k] += subtotal
cumulative_raw[name][k] += subtotal
if step not in scheduled:
continue
current, error, current_aperture = capture(step)
sample_dt = step - previous_sample_step
for name in current:
integrated_aperture[name] += (
0.5 * sample_dt * (previous_aperture[name] + current_aperture[name])
)
integrated[name] += (
0.5
* sample_dt
* (
np.array([x[1] for x in previous[name]])
+ np.array([x[1] for x in current[name]])
)
)
previous = current
previous_aperture = current_aperture
previous_sample_step = step
if step % c["window"] == 0 or step == c["steps"]:
end0 = rows[-1]["step_end"] if rows else 0
duration = step - end0
row = dict(
step_start=end0,
step_end=step,
t_star=step * u / c["diameter"],
fine_velocity_relative_l2=error,
elapsed_seconds=time.monotonic() - start,
measurements={},
)
for name in current:
delta = np.array([x[0] for x in current[name]]) - np.array(
[x[0] for x in initial[name]]
)
mass_cv = delta[:, 0] + integrated[name][:, 0]
reaction_cv = -delta[:, 1:] - integrated[name][:, 1:]
ledger = accumulated[name]
total = ledger["total"]
normal = 0.5 * u * u * c["diameter"]
row["measurements"][name] = dict(
raw_interval={k: v.tolist() for k, v in ledger.items()},
cv_mass_into_fluid=mass_cv.tolist(),
cv_reaction_impulse=reaction_cv.tolist(),
mass_balance_residual=(mass_cv - total[3]).tolist(),
reaction_balance_residual=(reaction_cv - total[:2]).tolist(),
reaction_residual_coefficient=(
(reaction_cv - total[:2]) / duration / normal
).tolist(),
raw_reaction_coefficient=(
total[:2] / duration / normal
).tolist(),
cv_reaction_coefficient=(
reaction_cv / duration / normal
).tolist(),
contour_reaction_spread_coefficient=(
(reaction_cv[1] - reaction_cv[0]) / duration / normal
).tolist(),
measured_arc_mass_rates=(
integrated_aperture[name] / duration
).tolist(),
arc_max_sample_distance=[v[-1] for v in probes[name]],
nominal_volume_fluxes=[
2 * np.pi / 12 * radius * c["jet_speed"],
-2 * np.pi / 12 * radius * c["jet_speed"],
],
aperture_mass_rates={
k: float(ledger[k][3] / duration)
for k in ("blowing", "suction", "wall", "total")
},
)
integrated[name].fill(0)
integrated_aperture[name].fill(0)
for v in accumulated[name].values():
v.fill(0)
rows.append(row)
initial = current
if recorder is not None:
recorder.flush()
payload = dict(
case=c,
**{"source_revision": PIN, **provenance},
status="PILOT_DIAGNOSTIC_ONLY",
windows=rows,
physics_validation_pass=False,
precision="FP32 evolution, FP64 host diagnostics",
)
(output / f"re{c['re_d']}_s{c['scale']}.json").write_text(
json.dumps(payload, indent=2) + "\n"
)
print(
json.dumps(
dict(
status="PROGRESS",
step=step,
steps=c["steps"],
elapsed_seconds=row["elapsed_seconds"],
fine_velocity_relative_l2=error,
)
),
flush=True,
)
return payload
def main():
p = argparse.ArgumentParser(description=__doc__)
p.add_argument("--repo", type=Path, required=True)
p.add_argument("--output", type=Path, required=True)
p.add_argument("--scale", type=int, choices=[1, 2, 4], default=1)
p.add_argument("--re", type=int, choices=[100, 400], default=100)
p.add_argument(
"--steps", type=int, default=800, help="Base s=1 steps; scales as s squared"
)
p.add_argument("--allow-re400", action="store_true")
p.add_argument(
"--continuation",
action="store_true",
help="Fixed measurement campaign: Re100, scale1/2, 3200 reference steps",
)
p.add_argument("--plan", action="store_true")
p.add_argument("--current-checkout", action="store_true",
help="Run on the current clean checkout and record its provenance "
"instead of requiring the frozen pinned source")
a = p.parse_args()
if a.steps <= 0 or a.steps > 3200:
p.error("steps must be 1..3200")
if a.re == 400 and not a.allow_re400:
p.error("Re400 requires explicit conditional-campaign opt-in")
if a.continuation:
if a.scale not in (1, 2) or a.re != 100 or a.steps != 3200:
p.error("--continuation requires --re 100 --scale 1/2 --steps 3200")
c = case(a.scale, a.re, a.steps)
if a.continuation:
from jet_cv_measurements import sample_steps, factors
c.update(
measurement_schema=1,
sample_count=len(sample_steps(c)),
startup_sample_dt_reference=1,
later_sample_dt_reference=4,
dense_startup_end_reference=100,
reference_conversion=factors(a.scale),
post_ready_budget_seconds=1200 if a.scale == 1 else 5400,
runtime_estimate_seconds=[210, 500] if a.scale == 1 else [900, 3200],
)
plan = dict(
case=c,
required_revision=None if a.current_checkout else PIN,
source_mode="current_checkout" if a.current_checkout else "pinned",
required_free_device_mib=1024,
estimated_peak_device_mib=128 * a.scale**2,
estimated_output_mib=64 if a.continuation else 2,
max_steps=c["steps"],
gpu_launch_owner="root; this runner does not select or reserve a GPU",
)
if a.plan:
print(json.dumps(plan, indent=2))
return
if (
not os.environ.get("CUDA_VISIBLE_DEVICES")
or "," in os.environ["CUDA_VISIBLE_DEVICES"]
):
p.error("Root must select one verified idle GPU via CUDA_VISIBLE_DEVICES")
repo = a.repo.resolve()
if a.current_checkout:
provenance = dict(**checkout_provenance(repo), **current_checkout(repo))
if provenance["checkout_dirty"]:
raise RuntimeError("Current-checkout runs require a clean checkout")
else:
source = verify(repo)
provenance = dict(
source_experiment=source["experiment"], source_mode="pinned",
**checkout_provenance(repo)
)
sys.path.insert(0, str(repo / "src"))
import CelerisLab
if Path(CelerisLab.__file__).resolve() != repo / "src/CelerisLab/__init__.py":
raise RuntimeError("Wrong source import")
from CelerisLab.cuda.context import CudaContext
import pycuda.driver as cuda
a.output.mkdir(parents=True, exist_ok=False)
(a.output / "plan.json").write_text(json.dumps(plan, indent=2) + "\n")
with CudaContext(0) as ctx:
with ctx.activate():
free, total = cuda.mem_get_info()
if free < max(1024, 2 * plan["estimated_peak_device_mib"]) * 1024**2:
raise RuntimeError("Insufficient free device memory")
run_case(c, ctx, a.output, provenance, continuation=a.continuation)
print(
json.dumps(
dict(status="COMPLETE_DIAGNOSTIC_ONLY", physics_validation_pass=False)
),
flush=True,
)
if __name__ == "__main__":
main()