#!/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 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 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, context=ctx, refinements=[c["bounds"]], level_blocks=[(8, 8), (8, 8)], level_objects=[(), (obj,)], level_jets=[(), tuple(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") 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=PIN, 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() source = verify(repo) provenance = dict( source_experiment=source["experiment"], **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()