Compare commits

..
12 changed files with 767 additions and 69 deletions
+35
View File
@@ -0,0 +1,35 @@
# Compatibility and API changes
Changes that can make a previously accepted setup or call behave differently.
Each entry names the commit, the affected entry points and the observable result.
## 87c3a31: cut-link ownership is validated for every body
`validate_cut_link_ownership` (`body/preprocess/flag_overlay.py`) now runs for
passive bodies as well as Jet/aperture bodies. Every cut link must have an
active fluid receiver (FLUID set, SOLID clear), and each `(receiver, direction)`
population may have only one writer.
Rejected with `ValueError`, before flags, populations or device buffers change:
1. Overlapping or touching passive solids whose cut links write the same
population ("Boundary requires unique cut-link writers").
2. A body whose cut-link receiver is a non-fluid cell, e.g. a channel wall row
or the inlet/outlet column ("Cut links must have active fluid receivers").
Both were previously accepted silently for passive bodies; the interpolated
boundary then had undefined results. The rules apply to the cut links of all
bodies together, so a passive body overlapping a Jet body is rejected as well.
The check runs in `Simulation.initialize`, runtime body sync (`sync_bodies`),
checkpoint save and checkpoint load (`common/checkpoint.py`). Load validates
before restoring anything, so a checkpoint whose saved geometry violates either
rule no longer restores; rebuild the geometry with separated bodies.
## 87c3a31: `Simulation.read_surface_exchange` accepts Jet-driven bodies
`read_surface_exchange(id)` previously required `SurfaceActuation` on the body.
It now also accepts bodies that own Jet aperture links, including Jets with zero
or negative amplitude (the aperture test is geometric). Bodies with neither
`SurfaceActuation` nor aperture links still raise
`ValueError("Body has no surface actuation")`. The value is unchanged: the latest
step's wall/aperture exchange, not total nozzle force.
+59
View File
@@ -56,3 +56,62 @@ Root must independently review and admit any execution. An admitted paired pilot
the existing watchdog, device selection, `--gpu-gates-approved`, positive
`--max-wall-seconds` and a new output directory, adding only `--exchange-window` to the
selected collision's pilot command. No GPU or remote action was performed for this patch.
## Public Jet diagnostic option
`--collision MRT --exchange-window --jet-scenario zero|fixed|reverse` selects
unshifted FP32, no LES, constant transfer on both existing public backends. Default
`none` and the historical manifests remain unchanged. `zero` registers both Jets at
zero amplitude and retains every existing DFG verdict, including pressure. `fixed`
uses outward radial parabolic sectors at pi/2 and 3pi/2, half-width pi/12, peaks
(+.001,-.001). Both use the existing D40/D80 geometry and U=.01. Sector origins are
(80.5,81) in root edge coordinates and (160,160.5) in uniform node coordinates;
radial extent is [0,1.005R]. `reverse` keeps those peaks through t*=300, updates both
public setters at that completed boundary, and continues to 305 with opposite signs.
Fixed diagnostic windows are 250–260, 260–270, 270–280, 280–290 and 290–300; reversal
adds 300–305. Histories retain exact interval endpoints and root/uniform lattice
impulses. The refined window uses total exchange; uniform uses public `read_force`
batch averages. Driven windows report mean/min/max/RMS Cd/Cl without applying DFG
reference or periodicity gates. Exit zero means healthy execution, not qualification;
`matched_case_pass` remains false. No independent CV or stress-load claim is made.
This opt-in records current source/helper hashes with the existing manifest format,
checkout revision/dirty state and runner hash. Pilot/full must match that snapshot
and the exact scenario. It does not relabel current source as historical evidence.
Use the existing bounded plan/pilot/full commands with the three option flags above;
full requires the matching completed t*=5 paired pilot. `step_seconds` includes force
readback; observation/overhead and total wall time are separate, not kernel throughput.
Add `--observe-cv` to record independent scalar control-volume samples from the
already downloaded public macros. The pilot exercises this path at t*=.02..5;
full records 501 endpoints at 290..300. Each `cv_history.jsonl` row contains actual
root time, preceding interval start/duration/impulse, and two inventories/outward
fluxes ordered [mass,Px,Py]. Exclude the first row's preceding impulse: that row
supplies the initial inventory. Uniform native impulse converts by .5² to root
units; refined ledger impulse is already in root units. Explicit native-aligned
contours differ by .25 root cells in x and are recorded in the plan. Cropping
retains the two-cell face-gradient stencils. No fields or per-link arrays are saved.
`cv_observation_seconds` is the CV portion of observation/overhead time. Helper hash
and the observation flag must match pilot/full. Time integration, temporal
coarsening and any physical assessment remain separate offline work.
## Bounded resolution and placement controls
Defaults remain `--case both --uniform-diameter 80 --placement wide`. Single-case
runs use `--case uniform` or `--case refined` and never report a matched-case pass.
Their own reference/periodicity or driven-diagnostic results remain available.
A single zero-Jet full run exits successfully only when its periodicity and the
unchanged complete DFG reference verdict pass; driven completion is not qualification.
The D160 option requires `--case uniform --uniform-diameter 160` and an explicit
public Jet scenario. U and Jet speeds stay .01 and ±.001; viscosity is .016,
grid 3521×658, center (320,320.5), radius80. CV uses root spacing/time factor .25,
origin (.375,.75), impulse factor .0625, root viscosity .004, and x-face shift .125.
`--case refined --placement compact` selects root bounds (48,48,128,128), fine
168×168 and origin (46,46), with the same D40/fineD80 body/Jet resolution. Its CV
contours are (52,52,112,112) and (50,50,120,120); complete face-gradient stencils
are clear of both solids and the transfer footprint. Compact is refined-only.
CV metadata gives each grid's spacing/origin/contours explicitly. Case, uniform
diameter and placement must match the newly completed pilot in every receipt.
+2 -5
View File
@@ -216,11 +216,8 @@ def surface_link_velocity(objects, cfg, flags, fluid_idx, dir, body_id, rx, ry,
enabled = [obj for obj in objects if obj.surface_actuation is not None]
if not enabled and not jets:
return None, None
from ..lbm.descriptors import FLUID, SOLID
if np.any((flags[fluid_idx] & FLUID) == 0) or np.any(flags[fluid_idx] & SOLID):
raise ValueError("Surface cut links must have active fluid receivers")
if len(np.unique(fluid_idx.astype(np.uint64)*cfg.nq+dir)) != len(dir):
raise ValueError("Surface boundary requires unique cut-link writers")
from .preprocess.flag_overlay import validate_cut_link_ownership
validate_cut_link_ownership(flags, fluid_idx, dir, cfg.nq)
surface_velocity = np.zeros((len(dir), 2), np.float32)
surface_aperture = np.zeros(len(dir), bool)
for obj in enabled:
+17 -12
View File
@@ -27,7 +27,7 @@ import pycuda.driver as cuda
from .objects import SimObject
from .registry import BodyRegistry
from .preprocess.flag_overlay import merge_flag_masks
from .preprocess.flag_overlay import merge_flag_masks, validate_cut_link_ownership
from ..cuda.compiler_v2 import obs_layout
from ..cuda.obs_layout import ObsLayout
from ..lbm.descriptors import (
@@ -369,9 +369,9 @@ class ObjectManager:
Args:
objects: List of SimObject instances to process.
domain_flags: Optional uint16 flags for donor-cell classification.
domain_flags: Optional uint16 flags for receiver/donor classification.
sensor_cell_counts_out: Int array to fill with per-object sensor
cell counts. If provided, ``fill(0)`` is called first.
cell counts after the full geometry has passed validation.
hits_out: List that receives one FP64 (n, 2) array of true hit
points row-aligned with the cut-link columns (for Jets).
@@ -394,8 +394,8 @@ class ObjectManager:
active_fluid = ((domain_flags & FLUID) != 0) & ((domain_flags & SOLID) == 0)
force_owners: set[int] = set()
if sensor_cell_counts_out is not None:
sensor_cell_counts_out.fill(0)
sensor_counts = (np.zeros_like(sensor_cell_counts_out)
if sensor_cell_counts_out is not None else None)
ez = None
if self.cfg and self.cfg.is_d3q19:
@@ -430,9 +430,9 @@ class ObjectManager:
if len(cells) > 0:
s_cells.append(cells)
s_ids.append(ids)
if (sensor_cell_counts_out is not None
and 0 <= obj.obj_id < sensor_cell_counts_out.size):
sensor_cell_counts_out[obj.obj_id] = int(len(cells))
if (sensor_counts is not None
and 0 <= obj.obj_id < sensor_counts.size):
sensor_counts[obj.obj_id] = int(len(cells))
if obj.is_force_region and hasattr(obj, 'get_sensor_list'):
cells, ids = obj.get_sensor_list(self.nx, self.ny, self.nz)
keep = active_fluid[cells]
@@ -461,8 +461,9 @@ class ObjectManager:
r_fr_ids = np.concatenate(fr_ids) if fr_ids else np.zeros(0, dtype=np.int32)
r_hit = np.concatenate(cl_hit) if cl_hit else np.zeros((0, 2))
if hits_out is not None:
hits_out.append(r_hit)
# Passive walls need the same ownership guarantees as driven walls.
# Overlapping solids are not a supported union-boundary operation.
validate_cut_link_ownership(domain_flags, r_idx, r_dir, self.nq)
if self.jets or any(obj.surface_actuation is not None for obj in objects):
# Reject unresolved apertures/ownership before sync or restart mutates state.
from ..lbm.curved_links import CurvedLinkSoA
@@ -470,6 +471,11 @@ class ObjectManager:
links.assign_host(r_idx, r_dir, r_q, r_bid, r_rx, r_ry, r_rz, r_fallback)
links.assign_surface(objects, self.cfg, domain_flags, jets=self.jets, hits=r_hit)
if sensor_cell_counts_out is not None:
sensor_cell_counts_out[:] = sensor_counts
if hits_out is not None:
hits_out.append(r_hit)
return (
r_idx, r_dir, r_q, r_bid,
r_rx, r_ry, r_rz, r_fallback,
@@ -487,8 +493,6 @@ class ObjectManager:
the caller has already prepared and uploaded the final flag field
for this geometry, such as during initialization.
"""
self._telemetry_field = field
next_flags = self.build_flags(field.build_channel_flags()) if rebuild_flags else field.flag
hits = []
@@ -499,6 +503,7 @@ class ObjectManager:
fr_cells, fr_obj_id,
) = self.build_compact_lists(domain_flags=next_flags, hits_out=hits)
self._telemetry_field = field
if rebuild_flags:
field.flag = next_flags
field.upload_flags()
+10 -1
View File
@@ -1,6 +1,6 @@
# CelerisLab/body/preprocess/flag_overlay.py
"""
Flag overlay: merge per-object flag masks onto a clean domain base.
Flag overlay and cut-link ownership checks against the merged domain.
"""
from __future__ import annotations
@@ -12,6 +12,15 @@ import numpy as np
from ...lbm.descriptors import FLUID, SOLID, SENSOR_FLAG, FRC_REGION
def validate_cut_link_ownership(flags, fluid_idx, directions, nq):
"""Reject solid receivers and multiple writers to one population."""
if np.any((flags[fluid_idx] & FLUID) == 0) or np.any(flags[fluid_idx] & SOLID):
raise ValueError("Cut links must have active fluid receivers")
writers = fluid_idx.astype(np.uint64) * nq + directions
if len(np.unique(writers)) != len(directions):
raise ValueError("Boundary requires unique cut-link writers")
def merge_flag_masks(
base: np.ndarray,
masks: List[np.ndarray],
+17 -8
View File
@@ -336,13 +336,19 @@ class Simulation:
def read_surface_exchange(self, id: int):
"""Latest step's wall/aperture (Fx,Fy,Tz) exchange, not total nozzle force.
Supports SurfaceActuation and Jet-driven bodies, including zero-amplitude
apertures. Passive bodies without an aperture are rejected.
Unlike read_body this is not accumulated/time-averaged. It resets after
geometry sync or checkpoint load; advance before interpreting it.
"""
if self.bodies.get(id).surface_actuation is None:
obj = self.bodies.get(id)
curved = self.field.curved
has_aperture = (curved.aperture is not None
and np.any(curved.aperture & (curved.body_id == id)))
if obj.surface_actuation is None and not has_aperture:
raise ValueError("Body has no surface actuation")
self.stream.synchronize()
return self.field.curved.surface_exchange(id)
return curved.surface_exchange(id)
@context_bound
def read_body(self, id: int, *,
@@ -451,18 +457,21 @@ class Simulation:
@context_bound
def initialize(self):
"""Initialize flow field and sync objects to GPU."""
# Recompile if objects were added after construction.
self._assert_allocation_contract()
# Validate the complete geometry before flags, populations or modules
# change. Compact-list construction publishes no outputs on failure.
next_flags = self.bodies.build_flags(self.field.build_channel_flags())
if self.bodies.jets and self.bodies.count == 0:
raise ValueError("Jets need bodies with cut links")
self.bodies.build_compact_lists_for(self.bodies.objects, next_flags)
# Recompile if objects were added after construction.
if self._read_generated_n_objects() != self.bodies.count:
self.recompile()
self.field.flag = self.field.build_channel_flags()
if self.bodies.count > 0:
self.field.flag = self.bodies.build_flags(self.field.flag)
self.field.flag = next_flags
self.field.upload_flags()
if self.bodies.jets and self.bodies.count == 0:
raise ValueError("Jets need bodies with cut links")
self.stepper.initialize()
if self.bodies.count > 0:
self.bodies.sync_to_gpu(self.field, rebuild_flags=False)
+34
View File
@@ -16,6 +16,7 @@ from run_jet_circle_cv import (
case,
state,
cv,
cv_from_macros,
aperture_probe,
measured_aperture_flux,
verify,
@@ -28,6 +29,39 @@ def eq(rho, u, v):
class Host(unittest.TestCase):
def test_public_macros_pressure_advection_and_cropped_frames(self):
for h, origin, shape, rect in (
(.5, (0., 0.), (48, 64), (5., 4., 25., 20.)),
(.5, (.25, .5), (48, 64), (5.25, 4., 25.25, 20.)),
(.25, (.375, .75), (96, 128), (5.125, 4., 25.125, 20.)),
):
y, x = np.indices(shape)
x = origin[0] + h * (x + .5)
k, u, v = 1e-5, .02, -.01
rho = 1 + 3 * k * x # p=(rho-1)/3=k*x
ux, uy, fluid = np.full(shape, u), np.full(shape, v), np.ones(shape, bool)
result = cv_from_macros(rho, ux, uy, fluid, origin, rect, .004, h)
expected = 320 * np.array([3*k*u, k + 3*k*u*u, 3*k*u*v])
np.testing.assert_allclose(result[1], expected, rtol=1e-11, atol=1e-13)
np.testing.assert_allclose(result, cv(eq(rho, u, v), fluid, origin, rect, .004, h), atol=1e-12)
x0, y0, x1, y1 = np.rint((np.array(rect)-np.tile(origin, 2))/h).astype(int)
selection = np.s_[y0-2:y1+2, x0-2:x1+2]
cropped_origin = np.asarray(origin) + h*np.array([x0-2, y0-2])
cropped = cv_from_macros(rho[selection], ux[selection], uy[selection],
fluid[selection], cropped_origin, rect, .004, h)
np.testing.assert_array_equal(cropped, result)
def test_public_macros_full_normal_viscous_stress(self):
y, x = np.indices((48, 64))
x = (x + .5) * .5
rho, ux, uy, fluid = np.ones_like(x), 1e-4*x*x, np.zeros_like(x), np.ones(x.shape, bool)
args = (rho, ux, uy, fluid, (0., 0.), (5., 4., 25., 20.))
_, viscous = cv_from_macros(*args, .004)
_, inviscid = cv_from_macros(*args, 0.)
np.testing.assert_allclose(viscous-inviscid, [0, -4*.004*1e-4*320, 0], atol=1e-13)
with self.assertRaisesRegex(ValueError, "matching 2D"):
cv_from_macros(rho, ux, uy[:-1], fluid, (0, 0), (5, 4, 25, 20), .004)
def test_descriptor_and_uniform_closed_flux(self):
from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY
+49
View File
@@ -188,6 +188,55 @@ def jet_case(device, jets):
return field, bodies
def test_public_surface_exchange_for_jet_preserves_aperture_at_zero_and_reverse(device):
from CelerisLab.simulation import Simulation
jet = Jet(Sector(CENTER, 0., .6, 3.), Radial(CENTER, inward=False), "quartic", 2e-3)
field, bodies = jet_case(device, [jet])
sim = SimpleNamespace(bodies=bodies, field=field, _initialized=True,
stream=SimpleNamespace(synchronize=device.synchronize))
cl = field.curved
assert bodies.get(0).surface_actuation is None
aperture, velocity = cl.aperture, cl.surface_velocity.copy()
assert aperture.any() and (~aperture).any()
# Synthetic per-link exchanges test readback/splitting, not kernel physics.
values = np.arange(cl.count * 3, dtype=np.float32).reshape(-1, 3)
device.htod(cl.link_exchange_gpu, values)
for amplitude in (2e-3, 0., -1e-3):
Simulation.set_jet_amplitude(sim, 0, amplitude)
result = Simulation.read_surface_exchange(sim, 0)
assert cl.aperture is aperture
np.testing.assert_array_equal(aperture, cl.jet_id >= 0)
np.testing.assert_allclose(cl.surface_velocity, velocity * (amplitude / 2e-3))
for name, mask in (("wall", ~aperture), ("aperture", aperture)):
np.testing.assert_array_equal(result[name], values[mask].sum(axis=0, dtype=np.float64))
def test_public_surface_exchange_keeps_legacy_and_passive_guards(device):
from CelerisLab.simulation import Simulation
field = make_field(device)
bodies = manager(field, [circle(SurfaceActuation(0., .6, .002), radius=2.)])
bodies.sync_to_gpu(field)
sim = SimpleNamespace(bodies=bodies, field=field,
stream=SimpleNamespace(synchronize=device.synchronize))
for value in Simulation.read_surface_exchange(sim, 0).values():
np.testing.assert_array_equal(value, 0.)
field.curved.free_gpu_columns()
with pytest.raises(RuntimeError, match="diagnostic is allocated"):
Simulation.read_surface_exchange(sim, 0)
# An exchange buffer elsewhere in the scene must not admit a passive body.
for jets in ([], [Jet(Sector(CENTER, 0., .6, 3.), Radial(CENTER))]):
field, bodies = jet_case(device, jets)
passive = bodies.add(circle(center=(2., 3.), radius=.6))
bodies.sync_to_gpu(field)
sim.bodies, sim.field = bodies, field
with pytest.raises(ValueError, match="Body has no surface actuation"):
Simulation.read_surface_exchange(sim, passive)
if not jets:
with pytest.raises(ValueError, match="Body has no surface actuation"):
Simulation.read_surface_exchange(sim, 0)
def test_manager_passes_true_hits_and_reuploads_only_velocity(device):
# CPU memory double only: checks host/device bookkeeping, not GPU numerics.
jet = Jet(Sector(CENTER, 0., .6, 3.), Radial(CENTER, inward=False), "quartic", 2e-3)
+174 -4
View File
@@ -146,23 +146,26 @@ def test_full_requires_completed_pilot():
r.verify_pilot(None, r.verify_source(r.DEFAULT_SOURCE, "MRT", "window"))
def pilot_fixture(path, observation_mode="window"):
def pilot_fixture(path, observation_mode="window", jet_scenario="none", case_selection="both", uniform_diameter=80, placement="wide"):
"""Synthetic contract fixture only; never a physical result or run receipt."""
manifest = (r.verify_source(r.DEFAULT_SOURCE, "MRT", "window")
manifest = (r.verify_source(r.DEFAULT_SOURCE, "MRT", "window", jet_scenario)
if observation_mode == "window" else
json.loads((r.HERE / "data/refined_dfg/source_manifest_mrt.json").read_text()))
selection = dict(case_selection=case_selection, uniform_diameter=uniform_diameter, placement=placement)
(path / "plan.json").write_text(
json.dumps(
dict(
source_manifest=manifest,
**selection,
collision="MRT",
jet_scenario=jet_scenario,
observation_mode=observation_mode,
bounds_root_edge=list(r.BOUNDS),
bounds_root_edge=list(r.PLACEMENTS[placement]),
runner_sha256=r.runner_hash(),
)
)
)
for case, diameter in (("uniform", 80), ("refined", 40)):
for case, diameter in r.selected_cases(case_selection, uniform_diameter):
dest = path / case
dest.mkdir()
interval = round(0.02 * diameter / r.SPEED)
@@ -189,7 +192,9 @@ def pilot_fixture(path, observation_mode="window"):
json.dumps(
dict(
case=case,
**selection,
collision="MRT",
jet_scenario=jet_scenario,
observation_mode=observation_mode,
mode="pilot",
status="pilot_complete",
@@ -206,7 +211,9 @@ def pilot_fixture(path, observation_mode="window"):
json.dumps(
dict(
case=case,
**selection,
collision="MRT",
jet_scenario=jet_scenario,
observation_mode=observation_mode,
mode="pilot",
step=interval,
@@ -451,6 +458,169 @@ def test_window_contract_rejection():
assert not backend.events
def test_public_jet_configuration_and_geometry():
for diameter, center, radius in ((40, (80.5, 81), 20), (80, (160, 160.5), 40)):
cfg, _, _ = r.config(h, diameter, "MRT", "fixed")
assert not cfg["method"]["ddf_shifting"] and not cfg["method"]["les"]["enabled"]
jets = r.public_jets(center, radius, "fixed")
assert [j.amplitude for j in jets] == [.001, -.001]
assert all(j.field.inward is False and j.profile == "parabolic" for j in jets)
assert all(j.region.r_max == radius * 1.005 for j in jets)
assert all(j.amplitude == 0 for j in r.public_jets(center, radius, "zero"))
cfg, _, _, host = r.host_geometry(r.DEFAULT_SOURCE, h, "MRT", "fixed")
assert not cfg.ddf_shifting and len(host.boundary_cells) > 0
assert r.public_jets((0, 0), 1, "none") == ()
with unittest.TestCase().assertRaises(ValueError):
r.config(h, 40, "SRT", "fixed")
def test_public_jet_reversal_boundary_and_windows():
from types import SimpleNamespace
calls = []
handle = SimpleNamespace(set_jet_amplitude=lambda i, a: calls.append((i, a)))
for diameter in (40, 80):
boundary = round(250 * diameter / r.SPEED)
assert r.jet_interval(handle, "fixed", boundary - 3, 80, diameter) == (3, None)
switch = round(300 * diameter / r.SPEED)
assert r.jet_interval(handle, "reverse", switch - 3, 80, diameter) == (3, None)
assert not calls
n, event = r.jet_interval(handle, "reverse", switch, 80, diameter)
assert n == 80 and event["step"] == switch and event["tstar"] == 300
assert calls == [(0, -.001), (1, .001)]
calls.clear()
assert r.jet_interval(handle, "reverse", switch + 80, 80, diameter) == (80, None)
assert r.jet_interval(handle, "fixed", switch, 80, diameter) == (80, None)
assert not calls
rows = [dict(tstar=t, tstar_force=t - .5, cd=2.0, cl=-.1) for t in range(251, 306)]
windows = r.jet_diagnostics(rows, "reverse")
assert len(windows) == 6 and all(w["complete"] for w in windows)
assert all(w["cd"]["mean"] == 2 and np.isclose(w["cl"]["mean"], -.1) for w in windows)
assert windows[-1]["covered_tstar"] == 5
with unittest.TestCase().assertRaisesRegex(ValueError, "crosses"):
r.jet_diagnostics([dict(tstar=250.5, tstar_force=250, cd=1, cl=0)], "fixed")
assert r.plan(h, "MRT", "window", "reverse")["modes"]["full_tstar"] == 305
def test_public_jet_source_and_pilot_separation():
from tempfile import TemporaryDirectory
with TemporaryDirectory() as tmp:
path = Path(tmp)
manifest = pilot_fixture(path, jet_scenario="fixed")
assert manifest["source_capture"].startswith("current checkout")
assert len(manifest["checkout_revision"]) == 40
assert r.verify_pilot(path, manifest, "MRT", "window", "fixed") == 3606.0
for scenario in ("none", "zero", "reverse"):
with unittest.TestCase().assertRaisesRegex(ValueError, "Jet scenario"):
r.verify_pilot(path, manifest, "MRT", "window", scenario)
for relative in ("uniform/summary.json", "refined/ready.json"):
receipt = path / relative
record = json.loads(receipt.read_text())
receipt.write_text(json.dumps(dict(record, jet_scenario="zero")))
with unittest.TestCase().assertRaises(ValueError):
r.verify_pilot(path, manifest, "MRT", "window", "fixed")
receipt.write_text(json.dumps(record))
with unittest.TestCase().assertRaises(ValueError):
r.verify_source(r.DEFAULT_SOURCE, "MRT", "last_base_step", "fixed")
def test_cv_snapshot_common_units_and_window_endpoints():
records = []
for refined, shape, origin, diameter, impulse in (
(True, (240, 568), (38., 22.), 40, [4., -8.]),
(False, (330, 1761), (.25, .5), 80, [16., -32.]),
(False, (658, 3521), (.375, .75), 160, [64., -128.]),
):
rho, ux, uy, fluid = np.ones(shape), np.full(shape, .01), np.zeros(shape), np.ones(shape, bool)
begin, n = round(290*diameter/r.SPEED), round(.02*diameter/r.SPEED)
record = r.cv_snapshot(rho, ux, uy, fluid, origin, refined, begin, n, diameter, impulse)
records.append(record)
np.testing.assert_allclose(record["outward_flux"], 0, atol=1e-12)
assert r.cv_selected(begin, diameter, "full") and not r.cv_selected(begin-n, diameter, "full")
end = round(300*diameter/r.SPEED)
assert r.cv_selected(end, diameter, "full") and not r.cv_selected(end+n, diameter, "full")
assert r.cv_selected(n, diameter, "pilot")
assert sum(r.cv_selected(step, diameter, "full") for step in range(begin-n, end+n+1, n)) == 501
for key in ("time_root", "interval_start_root", "interval_dt_root", "interval_impulse_root", "inventory"):
for record in records[1:]:
np.testing.assert_allclose(records[0][key], record[key], rtol=0, atol=1e-10)
assert records[0]["interval_dt_root"] == 80 and np.isclose(records[0]["tstar"], 290.02)
np.testing.assert_array_equal(records[0]["interval_impulse_root"], [4, -8])
def test_cv_pilot_mode_and_helper_provenance():
from tempfile import TemporaryDirectory
manifest = r.verify_source(r.DEFAULT_SOURCE, "MRT", "window", "fixed", True)
assert "tests/validation/run_jet_circle_cv.py" in manifest["sha256"]
with TemporaryDirectory() as tmp:
path = Path(tmp)
old = pilot_fixture(path, jet_scenario="fixed")
with unittest.TestCase().assertRaisesRegex(ValueError, "CV observation mode"):
r.verify_pilot(path, old, "MRT", "window", "fixed", True)
for relative in ("plan.json", "uniform/summary.json", "refined/summary.json",
"uniform/ready.json", "refined/ready.json"):
file = path / relative
record = json.loads(file.read_text())
record.update(observe_cv=True, cv_samples=250)
if relative == "plan.json":
record["source_manifest"] = manifest
file.write_text(json.dumps(record))
assert r.verify_pilot(path, manifest, "MRT", "window", "fixed", True) == 3606
def test_compact_geometry_and_cv_stencils():
from CelerisLab.lbm.refinement import ACTIVE
from CelerisLab.lbm.descriptors import SOLID
_, pair, _, host = r.host_geometry(r.DEFAULT_SOURCE, h, "MRT", "fixed", "compact")
assert pair.fine_shape == (168, 168) and tuple(pair.fine_origin) == (46, 46)
fluid = (pair.fine_kind == ACTIVE) & ((host.flags.reshape(pair.fine_shape) & SOLID) == 0)
rho, ux, uy = np.ones(pair.fine_shape), np.full(pair.fine_shape, .01), np.zeros(pair.fine_shape)
record = r.cv_snapshot(rho, ux, uy, fluid, pair.fine_origin, True, 1160000, 80, 40, [1, 2], "compact")
np.testing.assert_allclose(record["outward_flux"], 0, atol=1e-12)
stencil = np.zeros(pair.fine_shape, bool)
for rect in r.cv_contours(True, 40, "compact"):
x0, y0, x1, y1 = np.rint((rect - np.tile(pair.fine_origin, 2))/.5).astype(int)
stencil[y0-2:y1+2, x0-2:x0+2] = stencil[y0-2:y1+2, x1-2:x1+2] = True
stencil[y0-2:y0+2, x0-2:x1+2] = stencil[y1-2:y1+2, x0-2:x1+2] = True
assert fluid[stencil].all() and not stencil.ravel()[host.transfer_fine_cells].any()
def test_d160_config_and_single_case_receipts():
from tempfile import TemporaryDirectory
from unittest.mock import patch
cfg, _, p = r.config(h, 160, "MRT", "fixed")
assert (p["nx"], p["ny"], p["center"], p["radius"]) == (3521, 658, [320, 320.5], 80)
assert p["viscosity"] == .016 and p["mean_speed"] == .01 and not cfg["method"]["ddf_shifting"]
assert r.uniform_origin(160) == (.375, .75)
np.testing.assert_array_equal(r.cv_contours(False, 160), r.CV_CONTOURS + [.125, 0, .125, 0])
for scenario, expected_exit in (("fixed", 0), ("zero", 2)):
with TemporaryDirectory() as tmp:
path = Path(tmp)
manifest = pilot_fixture(path, jet_scenario=scenario, case_selection="uniform", uniform_diameter=160)
assert r.verify_pilot(path, manifest, "MRT", "window", scenario, False, "uniform", 160) == 1803
for selection, diameter, placement in (("both", 160, "wide"), ("uniform", 80, "wide"), ("uniform", 160, "compact")):
with unittest.TestCase().assertRaises(ValueError):
r.verify_pilot(path, manifest, "MRT", "window", scenario, False, selection, diameter, placement)
# Synthetic dispatch/result check: no CUDA imports or execution.
output = path / "full"
argv = ["runner", "--collision", "MRT", "--exchange-window", "--jet-scenario", scenario,
"--case", "uniform", "--uniform-diameter", "160", "--mode", "full",
"--pilot-results", str(path), "--output", str(output), "--max-wall-seconds", "10", "--gpu-gates-approved"]
result = dict(case="uniform", healthy=True, periodic=scenario == "zero", jet_diagnostic_windows=[],
status="driven_diagnostics_complete" if scenario == "fixed" else "periodic",
reference_comparison=dict(reference_agreement=False))
with patch.object(sys, "argv", argv), patch.dict(r.os.environ, CUDA_VISIBLE_DEVICES="mock-no-gpu"), \
patch.object(r, "run_case", return_value=result) as run:
assert r.main() == expected_exit and run.call_count == 1
comparison = json.loads((output / "comparison.json").read_text())
assert comparison["matched_case_pass"] is False and "refined_vs_uniform" not in comparison
assert comparison["status"].startswith("single_case_")
if __name__ == "__main__":
suite = unittest.TestSuite(
unittest.FunctionTestCase(fn)
+75
View File
@@ -238,6 +238,81 @@ def test_overlapping_force_regions_rejected(device):
bodies.build_compact_lists(bodies.build_flags(field.flag))
@pytest.mark.parametrize("second_x,error", [(14., "unique cut-link writers"),
(17., "active fluid receivers")])
def test_passive_overlap_rejected_before_publishing_lists(second_x, error):
cfg = LBMConfig(nx=32, ny=24)
cfg.validate()
bodies = ObjectManager(cfg.nx, cfg.ny, cfg.nz, cfg.nq, cfg)
for x in (14., second_x):
bodies.add(obj(x=x, y=12., r=4.3))
bodies.add(obj(x=4., y=4., r=2., is_sensor=True))
counts = np.full(bodies.count, 91, np.int32)
hits = []
flags = bodies.build_flags(np.full(cfg.nx * cfg.ny, FLUID, np.uint16))
for domain_flags in (flags, None):
with pytest.raises(ValueError, match=error):
bodies.build_compact_lists_for(bodies.objects, domain_flags,
sensor_cell_counts_out=counts, hits_out=hits)
np.testing.assert_array_equal(counts, 91)
assert hits == []
@pytest.mark.parametrize("centers", [((10., 12.), (22., 12.)),
((10.5, 12.), (19.5, 12.))])
def test_separated_passive_bodies_keep_valid_writers(centers):
# The second case has a subcell physical gap with a resolved fluid column.
cfg = LBMConfig(nx=32, ny=24)
cfg.validate()
bodies = ObjectManager(cfg.nx, cfg.ny, cfg.nz, cfg.nq, cfg)
for x, y in centers:
bodies.add(obj(x=x, y=y, r=4.3))
flags = bodies.build_flags(np.full(cfg.nx * cfg.ny, FLUID, np.uint16))
indices, directions, _, ids, *_ = bodies.build_compact_lists(flags)
assert set(ids) == {0, 1}
assert np.all(flags[indices] & FLUID) and not np.any(flags[indices] & SOLID)
assert len(np.unique(indices.astype(np.uint64) * cfg.nq + directions)) == len(indices)
@pytest.mark.parametrize("entry", ["sync", "initialize", "runtime_edit"])
@pytest.mark.parametrize("second_x", [5., 6.])
def test_passive_overlap_failure_keeps_live_state(device, entry, second_x):
from CelerisLab.simulation import Simulation
field = make_field(device)
bodies = manager(field, [obj(), obj(x=9., y=3., r=1., is_sensor=True)])
bodies.sync_to_gpu(field)
overlapping = obj(x=second_x)
if entry == "runtime_edit":
bodies.stage_add(overlapping)
else:
bodies.add(overlapping)
objects = bodies.objects
flags, curved, counts = field.flag.copy(), field.curved.fluid_idx.copy(), bodies.sensor_cell_counts.copy()
action = bodies.action.copy()
buffers = [(block, bytes(block.data), block.freed) for block in device.blocks]
writes = device.writes
sim = SimpleNamespace(bodies=bodies, field=field, lbm_cfg=field.cfg,
_initialized=entry == "runtime_edit",
_assert_allocation_contract=lambda: None,
stream=SimpleNamespace(synchronize=device.synchronize))
with pytest.raises(ValueError, match="active fluid receivers|unique cut-link writers"):
if entry == "sync":
bodies.sync_to_gpu(field)
elif entry == "initialize":
Simulation.initialize(sim)
else:
Simulation.sync_bodies(sim)
assert bodies.objects == objects
assert bodies.has_pending_edit() == (entry == "runtime_edit")
assert sim._initialized == (entry == "runtime_edit")
assert device.writes == writes and len(device.blocks) == len(buffers)
assert all(bytes(block.data) == data and block.freed == freed for block, data, freed in buffers)
np.testing.assert_array_equal(field.flag, flags)
np.testing.assert_array_equal(field.curved.fluid_idx, curved)
np.testing.assert_array_equal(bodies.sensor_cell_counts, counts)
np.testing.assert_array_equal(bodies.action, action)
def test_action_upload_and_scratch_are_pinned(device):
field = make_field(device, streaming="esopull")
bodies = manager(field, [obj(is_force_region=True)])
+19 -1
View File
@@ -74,7 +74,25 @@ def cv(g, fluid, origin, rectangle, nu, h=0.5):
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)
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
+276 -38
View File
@@ -17,8 +17,109 @@ 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) # 240 x 568 including four fine halo cells.
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):
@@ -39,26 +140,34 @@ def helpers(source):
return h, periodic
def config(h, diameter, collision="SRT"):
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"):
_, _, c = config(h, 40, collision)
_, _, f = config(h, 80, collision)
x0, y0, x1, y1 = BOUNDS
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),
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"],
@@ -78,19 +187,31 @@ def plan(h, collision="SRT", observation_mode="last_base_step"):
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)
/ (2 * f["nx"] * f["ny"]),
/ ((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="For origin (38,22), fine physical x=(j+75.5)*dx80 versus uniform i*dx80: half-fine-cell x offset. Fine physical y=(k+42.5)*dx80 equals uniform (l-.5)*dx80 at l=k+43: y samples align. Body, wall and source planes coincide.",
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": 300.0},
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 {}),
)
@@ -126,7 +247,7 @@ def initial_state(shape, origin, spacing, height=164.0, speed=SPEED, root=False)
return w * (1 + 3 * cu + 4.5 * cu * cu - 1.5 * velocity * velocity)
def host_geometry(source, h, collision="SRT"):
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
@@ -139,12 +260,12 @@ def host_geometry(source, h, collision="SRT"):
from CelerisLab.body.geometry.circle import CircleGeometry
from CelerisLab.body.objects import SimObject
cfgdict, _, _ = config(h, 40, collision)
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, BOUNDS)
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]
@@ -153,18 +274,41 @@ def host_geometry(source, h, collision="SRT"):
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,), _shifted_blocked=True)
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 verify_source(source, collision="SRT", observation_mode="last_base_step"):
def verify_source(source, collision="SRT", observation_mode="last_base_step", jet_scenario="none", observe_cv=False):
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")
if observe_cv and jet_scenario == "none":
raise ValueError("--observe-cv requires an explicit public Jet scenario")
name = ("source_manifest_window.json" if observation_mode == "window" else
"source_manifest.json" if collision == "SRT" else "source_manifest_mrt.json")
expected = json.loads((HERE / "data/refined_dfg" / name).read_text())
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
paths = set(expected["sha256"])
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():
if hashlib.sha256((source / name).read_bytes()).hexdigest() != digest:
raise RuntimeError("Source mismatch: " + name)
@@ -219,29 +363,33 @@ def run_case(a, h, periodic, case):
from contextlib import ExitStack
refined = case == "refined"
cfgdict, bodydict, p = config(h, 40 if refined else 80, a.collision)
cap = 5.0 if a.mode == "pilot" else 300.0
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)
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,
BOUNDS,
PLACEMENTS[a.placement],
context=context,
coarse_boundary="channel",
base_block=BASE_BLOCK,
fine_block=FINE_BLOCK,
fine_objects=(body,),
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 {}),
)
@@ -274,6 +422,9 @@ def run_case(a, h, periodic, case):
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
@@ -290,6 +441,10 @@ def run_case(a, h, periodic, case):
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:
@@ -359,6 +514,19 @@ def run_case(a, h, periodic, case):
/ 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)
@@ -402,14 +570,24 @@ def run_case(a, h, periodic, case):
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),
**(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(
@@ -425,7 +603,8 @@ def run_case(a, h, periodic, case):
else (
"pilot_complete"
if a.mode == "pilot"
else ("periodic" if result["periodic"] else "incomplete_periodicity")
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))
@@ -464,27 +643,39 @@ def check_deadline(deadline):
)
def verify_pilot(path, manifest, collision="SRT", observation_mode="last_base_step"):
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(
BOUNDS
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 (("uniform", 80), ("refined", 40)):
for case, diameter in selected_cases(case_selection, uniform_diameter):
summary = read_finite_json(path / case / "summary.json")
if (
summary.get("observation_mode") != observation_mode
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"
@@ -537,7 +728,10 @@ def verify_pilot(path, manifest, collision="SRT", observation_mode="last_base_st
raise ValueError("Unhealthy pilot observations")
ready = read_finite_json(path / case / "ready.json")
if (
ready.get("observation_mode") != observation_mode
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"
@@ -555,9 +749,9 @@ def verify_pilot(path, manifest, collision="SRT", observation_mode="last_base_st
or not 0 < stepping <= wall
):
raise ValueError("Invalid pilot calibration timings")
# Full physical duration is 60 times pilot; dense pilot sampling makes
# 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 + 60 * (wall - startup)
calibrated_seconds += startup + (61 if jet_scenario == "reverse" else 60) * (wall - startup)
return 3 * calibrated_seconds
@@ -566,6 +760,13 @@ def main():
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")
@@ -577,15 +778,26 @@ def main():
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)
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)
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 (
@@ -602,17 +814,31 @@ def main():
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)
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,
reference_gates=periodic.REFERENCE_GATES,
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,
)
host_geometry(a.source, h, a.collision)
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
@@ -621,8 +847,9 @@ def main():
"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 ("uniform", "refined"):
for case, _ in cases:
try:
result = run_case(a, h, periodic, case)
except TimeoutError as exc:
@@ -650,7 +877,12 @@ def main():
status="pilot_only" if a.mode == "pilot" else "incomplete",
matched_case_pass=False,
)
if len(results) == 2 and all(r["healthy"] and r["periodic"] for r in results):
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(
@@ -660,14 +892,20 @@ def main():
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) == 2
if len(results) == len(cases)
and all(r["healthy"] for r in results)
and (a.mode == "pilot" or comparison["matched_case_pass"])
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
)