Frank14f/CelerisLab
1
0 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

test: integrate pinned Re100 and jet control-volume validation

Browse Source
This commit is contained in:
TRON committed 2026-10-17 23:38:53 -07:00
1 parent d79e3f6f0a
commit 69e9b78de9
12 files changed
+2418
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

tests/validation/run_jet_circle_cv.py
+558
View File
@@ -0,0 +1,558 @@
#!/usr/bin/env python3
"""Bounded jet-on periodic cylinder: actual AMR Jet API vs blocked validation adapter."""
import argparse
import hashlib
import json
import os
from pathlib import Path
import sys
import time
import numpy as np
PIN = "d79e3f6f0a424f5bcc3a1a6d610313e49daf97a2"
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.
"""
rho, mom, vel = moments(g)
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 run_case(c, ctx, output):
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),
}
def capture():
ug = uniform.download()
ag = amr.download()[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)
)
return data, error, aperture
previous, _, previous_aperture = capture()
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 > 1800:
raise RuntimeError(
"1800-second post-READY pilot 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)
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
)
current, error, current_aperture = capture()
for name, rr in [("uniform_adapter", ur), ("amr_api", ar)]:
for k, mask in masks.items():
accumulated[name][k] += rr[mask].sum(axis=0)
integrated_aperture[name] += 0.5 * (
previous_aperture[name] + current_aperture[name]
)
integrated[name] += 0.5 * (
np.array([x[1] for x in previous[name]])
+ np.array([x[1] for x in current[name]])
)
previous = current
previous_aperture = current_aperture
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
payload = dict(
case=c,
source_revision=PIN,
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("--plan", action="store_true")
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")
c = case(a.scale, a.re, a.steps)
plan = dict(
case=c,
required_revision=PIN,
required_free_device_mib=1024,
estimated_peak_device_mib=128 * a.scale**2,
estimated_output_mib=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()
verify(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)
print(
json.dumps(
dict(status="COMPLETE_DIAGNOSTIC_ONLY", physics_validation_pass=False)
),
flush=True,
)
if __name__ == "__main__":
main()