"""Frozen matched SRT/MRT DFG comparison. GPU execution requires explicit root admission.""" import argparse import csv import hashlib import json import math import os from pathlib import Path import sys import time from tempfile import TemporaryDirectory import numpy as np HERE = Path(__file__).resolve().parent DEFAULT_SOURCE = HERE.parents[1] BOUNDS = (40, 24, 320, 140) BASE_BLOCK = (2, 881) # Exact 166 x 881 domain; no geometric padding. FINE_BLOCK = (8, 8) # Wide 240 x 568 or compact 168 x 168, including fine halos. SPEED = 0.01 JET_PEAK = 0.001 JET_WINDOWS = [(t, t + 10) for t in range(250, 300, 10)] CV_CONTOURS = np.array([(48., 44., 160., 120.), (44., 36., 180., 128.)]) PLACEMENTS = {"wide": BOUNDS, "compact": (48, 48, 128, 128)} COMPACT_CV_CONTOURS = np.array([(52., 52., 112., 112.), (50., 50., 120., 120.)]) def selected_cases(selection="both", uniform_diameter=80): return [(case, diameter) for case, diameter in (("uniform", uniform_diameter), ("refined", 40)) if selection in ("both", case)] def uniform_origin(diameter): h = 40 / diameter return (.5 - .5*h, 1 - h) def cv_contours(refined, diameter, placement="wide"): contours = CV_CONTOURS if placement == "wide" else COMPACT_CV_CONTOURS return contours if refined else contours + [20/diameter, 0., 20/diameter, 0.] def cv_selected(step, diameter, mode): return mode == "pilot" or round(290 * diameter / SPEED) <= step <= round(300 * diameter / SPEED) def cv_snapshot(rho, ux, uy, fluid, origin, refined, begin, n, diameter, impulse, placement="wide"): """Cropped public macros only; all inventories, fluxes and impulses use root units.""" from run_jet_circle_cv import cv_from_macros rectangles = cv_contours(refined, diameter, placement) h_cell = .5 if refined else 40 / diameter origin = np.asarray(origin, dtype=float) x0, y0, x1, y1 = np.rint((rectangles[-1] - np.tile(origin, 2)) / h_cell).astype(int) selection = np.s_[y0-2:y1+2, x0-2:x1+2] cropped_origin = origin + h_cell * np.array([x0-2, y0-2]) fields = [v[selection] for v in (rho, ux, uy, fluid)] values = [cv_from_macros(*fields, cropped_origin, rect, .004, h=h_cell) for rect in rectangles] h_native = 1.0 if refined else 40 / diameter return dict(step=begin+n, tstar=(begin+n)*SPEED/diameter, time_root=(begin+n)*h_native, interval_start_root=begin*h_native, interval_dt_root=n*h_native, interval_impulse_root=(np.asarray(impulse)*h_native**2).tolist(), inventory=[inventory.tolist() for inventory, _ in values], outward_flux=[flux.tolist() for _, flux in values]) def public_jets(center, radius, scenario): """Same physical apertures; caller supplies its public geometry frame.""" if scenario == "none": return () from CelerisLab.body.jets import Jet, Sector, Radial peak = 0.0 if scenario == "zero" else JET_PEAK return tuple( Jet(Sector(center, angle, math.pi / 12, radius * 1.005), Radial(center, inward=False), "parabolic", amplitude) for angle, amplitude in ((math.pi / 2, peak), (3 * math.pi / 2, -peak)) ) def jet_interval(handle, scenario, begin, n, diameter): """Apply a reversal only at a completed interval boundary, before advancing.""" if scenario != "none": boundaries = [round(t * diameter / SPEED) for t in (250, 260, 270, 280, 290, 300, 305)] n = min([n] + [step - begin for step in boundaries if step > begin]) if scenario != "reverse": return n, None switch = round(300 * diameter / SPEED) if begin < switch: return min(n, switch - begin), None if begin == switch: handle.set_jet_amplitude(0, -JET_PEAK) handle.set_jet_amplitude(1, JET_PEAK) return n, dict(step=begin, tstar=300.0, amplitudes=[-JET_PEAK, JET_PEAK]) return n, None def jet_diagnostics(rows, scenario): """Raw interval-weighted diagnostics only; no driven-flow qualification gate.""" windows = JET_WINDOWS + ([(300, 305)] if scenario == "reverse" else []) result = [] for start, end in windows: selected = [r for r in rows if start <= r["tstar_force"] < end] if not selected: continue if any(2 * r["tstar_force"] - r["tstar"] < start - 1e-9 or r["tstar"] > end + 1e-9 for r in selected): raise ValueError("Jet diagnostic interval crosses a window boundary") weights = np.array([2 * (r["tstar"] - r["tstar_force"]) for r in selected]) record = dict(start_tstar=start, end_tstar=end, samples=len(selected), covered_tstar=float(weights.sum()), complete=bool(np.isclose(weights.sum(), end - start, rtol=0, atol=1e-9))) for key in ("cd", "cl"): values = np.array([r[key] for r in selected]) record[key] = dict(mean=float(np.average(values, weights=weights)), min=float(values.min()), max=float(values.max()), rms=float(np.sqrt(np.average(values**2, weights=weights)))) result.append(record) return result def helpers(source): sys.path.insert(0, str(source / "src")) sys.path.insert(0, str(source / "tests/validation")) import CelerisLab if ( Path(CelerisLab.__file__).resolve() != source.resolve() / "src/CelerisLab/__init__.py" ): raise RuntimeError("Wrong CelerisLab import") import re100_helpers as h import periodic if h.ROOT.resolve() != source.resolve(): raise RuntimeError("Helper module loaded from wrong source") return h, periodic def config(h, diameter, collision="SRT", jet_scenario="none"): cfg, body, p = h.make_config(diameter, SPEED, collision, "neq_extrap", "zou_he_local") cfg["method"]["outlet"]["srt_neq_damp"] = 0.5 if jet_scenario != "none": if collision != "MRT": raise ValueError("Public Jet scenarios require MRT") cfg["method"]["ddf_shifting"] = False return cfg, body, p def plan(h, collision="SRT", observation_mode="last_base_step", jet_scenario="none", uniform_diameter=80, placement="wide"): _, _, c = config(h, 40, collision, jet_scenario) _, _, f = config(h, uniform_diameter, collision, jet_scenario) bounds = PLACEMENTS[placement] x0, y0, x1, y1 = bounds fine_active = 4 * (x1 - x0) * (y1 - y0) coarse_uncovered = c["nx"] * c["ny"] - fine_active // 4 return dict( status="prepared_not_run", collision=collision, jet_scenario=jet_scenario, ddf_shifting=jet_scenario == "none", observation_mode=observation_mode, root_gpu_gates_required=True, root=c, uniform=f, bounds_root_edge=list(bounds), fine_origin_root_edge=[x0 - 2, y0 - 2], physical_patch_m=[ (x0 - 0.5) * c["dx"], (y0 - 1) * c["dx"], (x1 - 0.5) * c["dx"], (y1 - 1) * c["dx"], ], body_root_edge=[80.5, 81.0], radius_root=20.0, root_block=list(BASE_BLOCK), fine_block=list(FINE_BLOCK), root_wall_edges=[1.0, 165.0], root_inlet_outlet_edges=[0.5, 880.5], physical_domain_m=[2.2, 0.41], physical_center_m=[0.2, 0.2], physical_diameter_m=0.1, physical_mean_speed_m_s=1.0, physical_viscosity_m2_s=0.001, finest_dx_m=f["dx"], finest_dt_s=f["dt"], refined_finest_dx_m=c["dx"]/2, refined_finest_dt_s=c["dt"]/2, approximate_leaf_cell_ratio=(coarse_uncovered + fine_active) / (f["nx"] * f["ny"]), approximate_collision_work_ratio=(coarse_uncovered + 2 * fine_active) / ((uniform_diameter/40) * f["nx"] * f["ny"]), cost_caveat="Geometry-only estimates include ghosts/solids; exclude block halos, transfers, observation and launch overhead. Not measured speedup.", sampling_offset="Uniform nodes map to root edges X=.5+i*(40/D), Y=1+(j-.5)*(40/D). Refined centers are fine_origin_root_edge+(index+.5)*.5. Body, wall and source planes coincide; compare observables, not pointwise arrays.", outlet=( "SRT pressure NEQ, all nine fneq damped by .5, no backflow clamp; noncanonical DFG outlet" if collision == "SRT" else "MRT pressure NEQ, unknowns q2/6/8 extrapolated; no backflow clamp; noncanonical DFG outlet" ), modes={"pilot_tstar": 5.0, "full_tstar": 305.0 if jet_scenario == "reverse" else 300.0}, observation_dense_tstar=0.02, **(dict(jets=dict( uniform_origin_node=f["center"], refined_origin_root_edge=[80.5, 81.0], angles_rad=[math.pi / 2, 3 * math.pi / 2], half_width_rad=math.pi / 12, radial_extent_over_radius=[0.0, 1.005], field="Radial(outward)", profile="parabolic", condition="velocity", initial_amplitudes=[0.0, 0.0] if jet_scenario == "zero" else [JET_PEAK, -JET_PEAK], reversal_tstar=300.0 if jet_scenario == "reverse" else None, diagnostic_windows_tstar=JET_WINDOWS + ([(300, 305)] if jet_scenario == "reverse" else []), transfer="constant", driven_qualification="not assessed", )) if jet_scenario != "none" else {}), ) def macros(g): rho = np.sum(g, axis=0, dtype=np.float64) ux = (g[1] - g[2] + g[5] - g[6] + g[7] - g[8]) / rho uy = (g[3] - g[4] + g[5] - g[6] - g[7] + g[8]) / rho return rho, ux, uy def leaf_mass(arrays, dx, masks): """Finite public leaves only; covered coarse, solid and halos are excluded.""" if any(not np.isfinite(g[:, mask]).all() for g, mask in zip(arrays, masks)): raise ValueError("Nonfinite physical leaf population") return sum( float(np.sum(np.sum(g, axis=0, dtype=np.float64)[mask], dtype=np.float64)) * (dx / 2**level) ** 2 for level, (g, mask) in enumerate(zip(arrays, masks)) ) def initial_state(shape, origin, spacing, height=164.0, speed=SPEED, root=False): """Same physical parabolic initial condition as uniform InitTubeFlow_v2.""" y = origin[1] + spacing * (np.arange(shape[0]) + 0.5) u = 6 * speed * ((y - 1) / height) * (1 - (y - 1) / height) velocity = np.broadcast_to(u[:, None], shape).copy() if root: velocity[[0, -1], :] = 0 velocity[:, [0, -1]] = 0 cx = np.array([0, 1, -1, 0, 0, 1, -1, 1, -1])[:, None, None] w = np.array([4 / 9] + [1 / 9] * 4 + [1 / 36] * 4)[:, None, None] cu = cx * velocity return w * (1 + 3 * cu + 4.5 * cu * cu - 1.5 * velocity * velocity) def host_geometry(source, h, collision="SRT", jet_scenario="none", placement="wide"): """Production donor/transfer and complete cut-link checks, without CUDA.""" from CelerisLab.config import load_lbm_config from CelerisLab.lbm.refinement import FixedPatch2D, prepare_fine_boundary from CelerisLab.lbm.blocks import BlockLayout2D from CelerisLab.lbm.channel import ( validate_channel_config, validate_channel_transfers, ) from CelerisLab.lbm.refinement_tables import blocked_base_tables from CelerisLab.body.geometry.circle import CircleGeometry from CelerisLab.body.objects import SimObject cfgdict, _, _ = config(h, 40, collision, jet_scenario) with TemporaryDirectory() as tmp: p = Path(tmp) / "cfg.json" p.write_text(json.dumps(cfgdict)) cfg = load_lbm_config(str(p)) pair = FixedPatch2D.from_config(cfg, PLACEMENTS[placement]) layout = BlockLayout2D.build(cfg.ny, cfg.nx, *BASE_BLOCK) validate_channel_config(cfg, [BASE_BLOCK, FINE_BLOCK]) targets = blocked_base_tables(pair, layout, np.argsort(pair.coalesce_coarse))[1] validate_channel_transfers(pair, layout, targets) BlockLayout2D.build(*pair.fine_shape, *FINE_BLOCK, periodic=(False, False)) body = SimObject( 0, geometry=CircleGeometry(80.5, 81.0, 20.0), center=(80.5, 81.0), radius=20.0 ) host = prepare_fine_boundary( cfg, pair, (body,), jets=public_jets((80.5, 81.0), 20.0, jet_scenario), _shifted_blocked=jet_scenario == "none", ) return cfg, pair, body, host def source_manifest(collision="SRT", observation_mode="last_base_step"): """The frozen historical manifest for this collision/observation mode (no checkout read).""" if collision not in ("SRT", "MRT"): raise ValueError("Matched collision must be SRT or MRT") if observation_mode not in ("last_base_step", "window"): raise ValueError("Unknown refined observation mode") name = ("source_manifest_window.json" if observation_mode == "window" else "source_manifest.json" if collision == "SRT" else "source_manifest_mrt.json") return json.loads((HERE / "data/refined_dfg" / name).read_text()) def verify_source(source, collision="SRT", observation_mode="last_base_step", jet_scenario="none", observe_cv=False): expected = source_manifest(collision, observation_mode) if observe_cv and jet_scenario == "none": raise ValueError("--observe-cv requires an explicit public Jet scenario") if jet_scenario != "none": if collision != "MRT" or observation_mode != "window": raise ValueError("Public Jet scenarios require MRT and --exchange-window") import subprocess sys.path.insert(0, str(HERE)) from source_provenance import checkout_provenance # Harness files from the manifest; production files as this checkout lists them. paths = {name for name in expected["sha256"] if not name.startswith("src/")} if observe_cv: paths.add("tests/validation/run_jet_circle_cv.py") paths.update(subprocess.check_output( ["git", "-C", str(source), "ls-files", "src"], text=True).splitlines()) return dict( source_capture="current checkout, not historical DFG qualification", **checkout_provenance(source), observation_mode=observation_mode, sha256={name: hashlib.sha256((source / name).read_bytes()).hexdigest() for name in sorted(paths)}, ) for name, digest in expected["sha256"].items(): path = source / name if not path.is_file() or hashlib.sha256(path.read_bytes()).hexdigest() != digest: raise RuntimeError("Source mismatch: " + name) return expected def advance_refined_interval(backend, n, deadline, observation_mode): """Return the same root-unit interval impulse, never an already averaged force.""" impulse = np.zeros(2) if observation_mode == "window": begin = backend.step_count check_deadline(deadline) backend.start_exchange_window() backend.step(n) exchange = backend.read_exchange_window(0, level=1) if ( exchange["units"] != "base_lattice_per_unit_depth" or exchange["window"] != "window" or exchange["sample_count"] != 2 * n or exchange["time_start"] != begin or exchange["time_end"] != begin + n or exchange["dt"] != 0.5 ): raise RuntimeError("Unexpected window force observation contract") impulse += exchange["total"]["impulse"] check_deadline(deadline) return impulse if observation_mode != "last_base_step": raise ValueError("Unknown refined observation mode") for _ in range(n): check_deadline(deadline) backend.step(1) exchange = backend.read_raw_exchange( 0, level=1, window="last_base_step" ) if ( exchange["units"] != "base_lattice_per_unit_depth" or exchange["sample_count"] != 2 or exchange["time_end"] - exchange["time_start"] != 1 ): raise RuntimeError("Unexpected force observation contract") impulse += exchange["total"]["impulse"] check_deadline(deadline) return impulse def run_case(a, h, periodic, case): from CelerisLab.simulation import Simulation from CelerisLab.cuda.context import CudaContext from CelerisLab.lbm.refinement_stepper import FixedPatchStepper from CelerisLab.lbm.descriptors import FLUID from contextlib import ExitStack refined = case == "refined" cfgdict, bodydict, p = config(h, 40 if refined else a.uniform_diameter, a.collision, a.jet_scenario) cap = 5.0 if a.mode == "pilot" else (305.0 if a.jet_scenario == "reverse" else 300.0) max_steps = math.ceil(cap * p["diameter"] / SPEED) rows = [] healthy = True step_seconds = 0.0 cv_seconds, cv_samples = 0.0, 0 amplitude_changes = [] started = time.monotonic() check_deadline(a.deadline) dest = a.output / case dest.mkdir() with ExitStack() as stack: cv_out = stack.enter_context((dest / "cv_history.jsonl").open("w")) if a.observe_cv else None if refined: cfg, pair, body, host = host_geometry(a.source, h, a.collision, a.jet_scenario, a.placement) context = stack.enter_context(CudaContext(0)) backend = stack.enter_context( FixedPatchStepper( cfg, context=context, coarse_boundary="channel", refinements=[PLACEMENTS[a.placement]], level_blocks=[BASE_BLOCK, FINE_BLOCK], level_objects=[(), (body,)], level_jets=[(), public_jets((80.5, 81.0), 20.0, a.jet_scenario)], observe_raw_exchange=True, **({"observe_window_exchange": True} if a.observation_mode == "window" else {}), ) ) from CelerisLab.lbm.channel import channel_masks from CelerisLab.lbm.refinement import ACTIVE from CelerisLab.lbm.descriptors import SOLID masks = [ channel_masks(pair.coarse_shape)[0] & ~pair.covered, (pair.fine_kind == ACTIVE) & ((host.flags.reshape(pair.fine_shape) & SOLID) == 0), ] backend.initialize( initial_state(pair.coarse_shape, (0.0, 0.0), 1.0, root=True), initial_state(pair.fine_shape, pair.fine_origin, 0.5), ) probe_p = dict( p, center=[ (80.5 - pair.fine_origin[0]) / 0.5 - 0.5, (81.0 - pair.fine_origin[1]) / 0.5 - 0.5, ], radius=40.0, ) else: tmp = Path(stack.enter_context(TemporaryDirectory())) (tmp / "lbm.json").write_text(json.dumps(cfgdict)) (tmp / "body.json").write_text(json.dumps(bodydict)) sim = stack.enter_context( Simulation(str(tmp / "lbm.json"), str(tmp / "body.json"), device_id=0) ) for index, jet in enumerate(public_jets(p["center"], p["radius"], a.jet_scenario)): if sim.add_jet(jet) != index: raise RuntimeError("Unexpected public Jet registration index") sim.initialize() probe_p = p fluid = (sim.field.flag.reshape(p["ny"], p["nx"]) & FLUID) != 0 startup = time.monotonic() - started check_deadline(a.deadline) begin = 0 with (dest / "history.csv").open("w") as out: writer = None while begin < max_steps: check_deadline(a.deadline) width = ( 0.02 if a.mode == "pilot" or begin * SPEED / p["diameter"] >= 240 else 0.5 ) n = min(max(1, round(width * p["diameter"] / SPEED)), max_steps - begin) n, change = jet_interval(backend if refined else sim, a.jet_scenario, begin, n, p["diameter"]) if change is not None: amplitude_changes.append(change) impulse = np.zeros(2) tick = time.monotonic() if refined: impulse = advance_refined_interval( backend, n, a.deadline, a.observation_mode ) else: done = 0 while done < n: check_deadline(a.deadline) batch = min(100, n - done) sim.run(batch) impulse += sim.read_force(0).astype(np.float64) * batch done += batch check_deadline(a.deadline) step_seconds += time.monotonic() - tick if refined: arrays = backend.download() fields = [macros(g) for g in arrays] root_rho, root_ux, _ = fields[0] rho, ux, uy = fields[1] fluid = masks[1] if any( not np.isfinite(g[:, mask]).all() for g, mask in zip(arrays, masks) ): healthy = False break mass = leaf_mass(arrays, p["dx"], masks) vals = [ tuple(v[mask] for v in field) for mask, field in zip(masks, fields) ] else: m = sim.get_macroscopic() rho, ux, uy = [ np.asarray(m[k], np.float64).reshape(p["ny"], p["nx"]) for k in ("rho", "ux", "uy") ] root_rho, root_ux = rho, ux mass = float(rho[fluid].sum()) * p["dx"] ** 2 vals = [(rho[fluid], ux[fluid], uy[fluid])] healthy = ( all(np.isfinite(v).all() for vs in vals for v in vs) and np.isfinite(impulse).all() ) if not healthy: break den = SPEED**2 * p["diameter"] / 2 row = dict( step=begin + n, tstar=(begin + n) * SPEED / p["diameter"], tstar_force=(begin + n / 2) * SPEED / p["diameter"], cd=float(impulse[0] / n / den), cl=float(impulse[1] / n / den), mass=mass, rho_min=float(min(vs[0].min() for vs in vals)), rho_max=float(max(vs[0].max() for vs in vals)), max_mach=float( max(np.sqrt(3 * (vs[1] ** 2 + vs[2] ** 2)).max() for vs in vals) ), flux_in=float(np.sum(root_rho[1:-1, 2] * root_ux[1:-1, 2])) * p["dx"] / SPEED, flux_out=float(np.sum(root_rho[1:-1, -3] * root_ux[1:-1, -3])) * p["dx"] / SPEED, **h.pressure_probes(rho, fluid, probe_p), ) if a.jet_scenario != "none": row.update(step_start=begin, tstar_start=begin * SPEED / p["diameter"], interval_steps=n, impulse_x=float(impulse[0]), impulse_y=float(impulse[1])) if a.observe_cv and cv_selected(begin+n, p["diameter"], a.mode): tick = time.monotonic() record = cv_snapshot(rho, ux, uy, fluid, pair.fine_origin if refined else uniform_origin(p["diameter"]), refined, begin, n, p["diameter"], impulse, a.placement) cv_out.write(json.dumps(record, allow_nan=False) + "\n") cv_out.flush() cv_samples += 1 cv_seconds += time.monotonic() - tick check_deadline(a.deadline) ensure_finite(row) rows.append(row) if writer is None: writer = csv.DictWriter(out, fieldnames=row) writer.writeheader() writer.writerow(row) out.flush() healthy = ( row["max_mach"] <= 0.2 and row["rho_min"] >= 0.5 and row["rho_max"] <= 1.5 ) begin += n if not healthy: break if len(rows) == 1: ready = dict( case=case, collision=a.collision, observation_mode=a.observation_mode, **a.provenance, mode=a.mode, step=begin, tstar=row["tstar"], runner_sha256=runner_hash(), all_runtime_and_observation_paths_exercised=True, ) (dest / "ready.json").write_text( json.dumps(ready, indent=2, allow_nan=False) ) print("READY " + json.dumps(ready, allow_nan=False), flush=True) result = dict( case=case, collision=a.collision, observation_mode=a.observation_mode, **a.provenance, mode=a.mode, history=rows, healthy=bool(healthy), steps=begin, startup_seconds=startup, step_seconds=step_seconds, observation_and_overhead_seconds=time.monotonic() - started - startup - step_seconds, wall_seconds=time.monotonic() - started, force="raw stationary-body link momentum reaction, per-unit-depth; average over entire observation interval; normalized in each base lattice", mass_units="rho_reference * m^2 per unit depth", flux_units="rho_reference * m^2/s per unit depth", pressure="postcollision end-sample; physical pressure difference / (rho_reference * 1 m/s squared)", suite_validated=False, **(periodic.assess(rows) if a.jet_scenario in ("none", "zero") else dict(periodic=False, checks={}, cycles=[])), ) if a.jet_scenario != "none": result.update(jet_diagnostic_windows=jet_diagnostics(rows, a.jet_scenario), amplitude_changes=amplitude_changes, pressure_classification=("DFG reference verdict retained, including pressure" if a.jet_scenario == "zero" else "diagnostic only; no official DFG reference gate")) if a.observe_cv: result.update(cv_observation_seconds=cv_seconds, cv_samples=cv_samples, cv_qualification="not assessed; scalar observations only") if not healthy: result["periodic"] = False reference = json.loads( (a.source / "tests/validation/data/dfg_re100/reference_cycle.json").read_text() ) if result["cycles"]: result["reference_comparison"] = periodic.compare( result["cycles"][-1], reference ) result["status"] = ( "invalid_state" if not healthy else ( "pilot_complete" if a.mode == "pilot" else ("driven_diagnostics_complete" if a.jet_scenario in ("fixed", "reverse") else "periodic" if result["periodic"] else "incomplete_periodicity") ) ) (dest / "summary.json").write_text(json.dumps(result, indent=2, allow_nan=False)) return result def runner_hash(): return hashlib.sha256(Path(__file__).read_bytes()).hexdigest() def ensure_finite(value): if isinstance(value, dict): for item in value.values(): ensure_finite(item) elif isinstance(value, list): for item in value: ensure_finite(item) elif isinstance(value, (float, int)) and not isinstance(value, bool): if not math.isfinite(value): raise ValueError("Nonfinite JSON value") def read_finite_json(path): def reject(value): raise ValueError("Nonfinite JSON constant: " + value) value = json.loads(path.read_text(), parse_constant=reject) ensure_finite(value) return value def check_deadline(deadline): if time.monotonic() >= deadline: raise TimeoutError( "Campaign wall-time ceiling reached; no continuation or acceptance" ) def verify_pilot(path, manifest, collision="SRT", observation_mode="last_base_step", jet_scenario="none", observe_cv=False, case_selection="both", uniform_diameter=80, placement="wide"): if path is None: raise ValueError( "Full comparison requires --pilot-results from a completed matched pilot" ) previous = read_finite_json(path / "plan.json") options = dict(case_selection=(case_selection, "both"), uniform_diameter=(uniform_diameter, 80), placement=(placement, "wide")) if any(previous.get(k, default) != expected for k, (expected, default) in options.items()): raise ValueError("Pilot case, resolution or placement differs") if previous.get("observe_cv", False) != observe_cv: raise ValueError("Pilot CV observation mode differs") if previous.get("jet_scenario", "none") != jet_scenario: raise ValueError("Pilot Jet scenario differs") if previous.get("observation_mode") != observation_mode: raise ValueError("Pilot observation mode differs") if previous.get("collision") != collision: raise ValueError("Pilot collision differs") if previous["source_manifest"] != manifest or previous["bounds_root_edge"] != list( PLACEMENTS[placement] ): raise ValueError("Pilot source or frozen geometry differs") if previous.get("runner_sha256") != runner_hash(): raise ValueError("Pilot runner differs") calibrated_seconds = 0.0 for case, diameter in selected_cases(case_selection, uniform_diameter): summary = read_finite_json(path / case / "summary.json") if ( any(summary.get(k, default) != expected for k, (expected, default) in options.items()) or summary.get("observe_cv", False) != observe_cv or (observe_cv and summary.get("cv_samples") != 250) or summary.get("jet_scenario", "none") != jet_scenario or summary.get("observation_mode") != observation_mode or summary.get("collision") != collision or summary.get("case") != case or summary.get("mode") != "pilot" or summary.get("status") != "pilot_complete" or summary.get("healthy") is not True ): raise ValueError("Mislabeled or unhealthy pilot receipt") rows = summary.get("history") interval = round(0.02 * diameter / SPEED) total = round(5 * diameter / SPEED) if ( not isinstance(rows, list) or len(rows) != 250 or summary.get("steps") != total ): raise ValueError("Incomplete pilot history or step count") for i, row in enumerate(rows): step = (i + 1) * interval if type(row.get("step")) is not int or row["step"] != step: raise ValueError("Nonmonotone or incomplete pilot step history") for key, expected in ( ("tstar", step * SPEED / diameter), ("tstar_force", (step - interval / 2) * SPEED / diameter), ): value = row.get(key) if type(value) not in (float, int) or not math.isclose( value, expected, rel_tol=0, abs_tol=1e-12 ): raise ValueError("Incorrect pilot observation time") for key in ( "cd", "cl", "dp", "mass", "rho_min", "rho_max", "max_mach", "flux_in", "flux_out", "rho_front", "rho_back", ): if type(row.get(key)) not in (int, float): raise ValueError("Missing numeric pilot observation") if not ( 0 < row["mass"] and 0.5 <= row["rho_min"] <= row["rho_max"] <= 1.5 and 0 <= row["max_mach"] <= 0.2 ): raise ValueError("Unhealthy pilot observations") ready = read_finite_json(path / case / "ready.json") if ( any(ready.get(k, default) != expected for k, (expected, default) in options.items()) or ready.get("observe_cv", False) != observe_cv or ready.get("jet_scenario", "none") != jet_scenario or ready.get("observation_mode") != observation_mode or ready.get("collision") != collision or ready.get("case") != case or ready.get("mode") != "pilot" or ready.get("step") != interval or ready.get("runner_sha256") != runner_hash() or ready.get("all_runtime_and_observation_paths_exercised") is not True ): raise ValueError("Invalid pilot READY receipt") startup = summary.get("startup_seconds") wall = summary.get("wall_seconds") stepping = summary.get("step_seconds") if ( any(type(v) not in (float, int) for v in (startup, wall, stepping)) or not 0 <= startup < wall or not 0 < stepping <= wall ): raise ValueError("Invalid pilot calibration timings") # Full duration is 60 (or 61 with reversal) times pilot; dense sampling makes # this conservative for the sparser early full observations. Not a guarantee. calibrated_seconds += startup + (61 if jet_scenario == "reverse" else 60) * (wall - startup) return 3 * calibrated_seconds def main(): parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--exchange-window", action="store_true", help="Opt in to refined window-ledger observations on the pinned candidate") parser.add_argument("--collision", choices=["SRT", "MRT"], default="SRT") parser.add_argument("--case", choices=["both", "uniform", "refined"], default="both") parser.add_argument("--uniform-diameter", type=int, choices=[80, 160], default=80) parser.add_argument("--placement", choices=list(PLACEMENTS), default="wide") parser.add_argument("--jet-scenario", choices=["none", "zero", "fixed", "reverse"], default="none", help="Opt in to public unshifted MRT Jets; reverse switches at t*=300 and ends at 305") parser.add_argument("--observe-cv", action="store_true", help="Record independent cropped-macro CV samples during pilot or full t*=290..300") parser.add_argument("--source", type=Path, default=DEFAULT_SOURCE) parser.add_argument("--output", type=Path, required=True) parser.add_argument("--mode", choices=["plan", "pilot", "full"], default="plan") parser.add_argument("--gpu-gates-approved", action="store_true") parser.add_argument("--pilot-results", type=Path) parser.add_argument( "--max-wall-seconds", type=float, help="Required explicit campaign wall-time ceiling for GPU execution", ) a = parser.parse_args() if a.uniform_diameter == 160 and a.case != "uniform": raise ValueError("D160 requires --case uniform; no matched-resolution claim") if a.placement == "compact" and a.case != "refined": raise ValueError("Compact placement requires --case refined") if (a.uniform_diameter != 80 or a.placement != "wide") and a.jet_scenario == "none": raise ValueError("Resolution/placement variants require an explicit public Jet scenario") a.source = a.source.resolve() a.observation_mode = "window" if a.exchange_window else "last_base_step" manifest = verify_source(a.source, a.collision, a.observation_mode, a.jet_scenario, a.observe_cv) h, periodic = helpers(a.source) from source_provenance import checkout_provenance a.provenance = checkout_provenance(a.source) a.provenance.update(jet_scenario=a.jet_scenario, ddf_shifting=a.jet_scenario == "none") a.provenance.update(case_selection=a.case, uniform_diameter=a.uniform_diameter, placement=a.placement) if a.observe_cv: a.provenance["observe_cv"] = True calibration = ( verify_pilot(a.pilot_results, manifest, a.collision, a.observation_mode, a.jet_scenario, a.observe_cv, a.case, a.uniform_diameter, a.placement) if a.mode == "full" else None ) if a.mode != "plan" and ( a.max_wall_seconds is None or not math.isfinite(a.max_wall_seconds) or a.max_wall_seconds <= 0 ): raise ValueError("GPU execution requires finite positive --max-wall-seconds") budget = ( min(a.max_wall_seconds, calibration) if calibration is not None else a.max_wall_seconds ) if a.output.exists() and any(a.output.iterdir()): raise ValueError("Output must be empty/new") a.output.mkdir(parents=True, exist_ok=True) specification = plan(h, a.collision, a.observation_mode, a.jet_scenario, a.uniform_diameter, a.placement) specification.update( source_manifest=manifest, **a.provenance, gates=periodic.GATES if a.jet_scenario in ("none", "zero") else None, reference_gates=periodic.REFERENCE_GATES if a.jet_scenario in ("none", "zero") else None, runner_sha256=runner_hash(), wall_time_budget_seconds=budget, full_pilot_calibration_seconds=calibration, ) if a.observe_cv: specification["cv_observation"] = dict( full_tstar=[290, 300], pilot_tstar=[.02, 5], spacing_tstar=.02, refined_contours_root=cv_contours(True, 40, a.placement).tolist(), uniform_contours_root=cv_contours(False, a.uniform_diameter, a.placement).tolist(), uniform_grid_origin_root=uniform_origin(a.uniform_diameter), refined_grid_origin_root=specification["fine_origin_root_edge"], cell_width_root={"uniform": 40/a.uniform_diameter, "refined": .5}, viscosity_root=.004, components=["mass", "Px", "Py"], units="root lattice per unit depth; outward flux per root time", interval="Impulse belongs to the interval ending at each snapshot; first retained endpoint (pilot .02, full 290) supplies initial inventory only: exclude its preceding impulse", qualification="not assessed; trapezoidal integration and time-coarsening are offline", ) if a.case != "uniform": host_geometry(a.source, h, a.collision, a.jet_scenario, a.placement) (a.output / "plan.json").write_text(json.dumps(specification, indent=2)) if a.mode == "plan": return 0 if not a.gpu_gates_approved or not os.environ.get("CUDA_VISIBLE_DEVICES"): raise RuntimeError( "GPU execution blocked: root GPU gates approval and explicit device selection required" ) results = [] cases = selected_cases(a.case, a.uniform_diameter) a.deadline = time.monotonic() + budget for case, _ in cases: try: result = run_case(a, h, periodic, case) except TimeoutError as exc: failure = dict( status="wall_time_limit", collision=a.collision, observation_mode=a.observation_mode, **a.provenance, case=case, matched_case_pass=False, error=str(exc), ) (a.output / "comparison.json").write_text( json.dumps(failure, indent=2, allow_nan=False) ) print(json.dumps(failure, allow_nan=False), flush=True) return 3 results.append(result) if not result["healthy"]: break comparison = dict( collision=a.collision, observation_mode=a.observation_mode, **a.provenance, status="pilot_only" if a.mode == "pilot" else "incomplete", matched_case_pass=False, ) if a.jet_scenario in ("fixed", "reverse"): comparison.update(status=("pilot_only" if a.mode == "pilot" else "driven_diagnostics_complete") if len(results) == len(cases) and all(r["healthy"] for r in results) else "invalid_state", qualification="not assessed; execution completion is not physical acceptance", jet_diagnostic_windows={r["case"]: r["jet_diagnostic_windows"] for r in results}) if a.case == "both" and len(results) == 2 and all(r["healthy"] and r["periodic"] for r in results): comparison.update( status="periodic_comparison", refined_vs_uniform=periodic.compare( results[1]["cycles"][-1], results[0]["cycles"][-1] ), ) comparison["matched_case_pass"] = comparison["refined_vs_uniform"][ "reference_agreement" ] if a.case != "both": comparison.update(status="single_case_" + results[0]["status"], own_reference_comparison=results[0].get("reference_comparison")) (a.output / "comparison.json").write_text( json.dumps(comparison, indent=2, allow_nan=False) ) return ( 0 if len(results) == len(cases) and all(r["healthy"] for r in results) and (a.mode == "pilot" or comparison["matched_case_pass"] or a.jet_scenario in ("fixed", "reverse") or (a.case != "both" and results[0]["periodic"] and results[0].get("reference_comparison", {}).get("reference_agreement", False))) else 2 ) if __name__ == "__main__": raise SystemExit(main())