#!/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 historical_sha256 manifest = json.loads((HERE / "data/jet_circle_cv/source_manifest.json").read_text()) if manifest["revision"] != PIN: raise RuntimeError("Wrong pinned revision") if historical_sha256(repo, PIN, manifest["files"]) != manifest["files"]: raise RuntimeError("Historical manifest does not match its pinned revision") 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()