From 5731c5171b42b4d0cbbbd479f5c48b491a6d780f Mon Sep 17 00:00:00 2001 From: Codex Date: Mon, 19 Oct 2026 07:38:03 -0700 Subject: [PATCH] test: add optional public macro CV observations --- docs/validation/refined_dfg_window.md | 13 +++++ tests/unit/test_jet_circle_cv.py | 34 ++++++++++++ tests/unit/test_refined_dfg.py | 43 +++++++++++++++ tests/validation/run_jet_circle_cv.py | 20 ++++++- tests/validation/run_refined_dfg.py | 76 ++++++++++++++++++++++++--- 5 files changed, 179 insertions(+), 7 deletions(-) diff --git a/docs/validation/refined_dfg_window.md b/docs/validation/refined_dfg_window.md index 0849c62..2e04fa6 100644 --- a/docs/validation/refined_dfg_window.md +++ b/docs/validation/refined_dfg_window.md @@ -82,3 +82,16 @@ and the exact scenario. It does not relabel current source as historical evidenc 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. diff --git a/tests/unit/test_jet_circle_cv.py b/tests/unit/test_jet_circle_cv.py index 38c616f..03d4377 100644 --- a/tests/unit/test_jet_circle_cv.py +++ b/tests/unit/test_jet_circle_cv.py @@ -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 diff --git a/tests/unit/test_refined_dfg.py b/tests/unit/test_refined_dfg.py index 8af3843..818e432 100644 --- a/tests/unit/test_refined_dfg.py +++ b/tests/unit/test_refined_dfg.py @@ -521,6 +521,49 @@ def test_public_jet_source_and_pilot_separation(): 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.]), + ): + 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"): + np.testing.assert_allclose(records[0][key], records[1][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 + + if __name__ == "__main__": suite = unittest.TestSuite( unittest.FunctionTestCase(fn) diff --git a/tests/validation/run_jet_circle_cv.py b/tests/validation/run_jet_circle_cv.py index bc71f33..b6f46d0 100644 --- a/tests/validation/run_jet_circle_cv.py +++ b/tests/validation/run_jet_circle_cv.py @@ -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 diff --git a/tests/validation/run_refined_dfg.py b/tests/validation/run_refined_dfg.py index 6f9bba7..ac8f59c 100644 --- a/tests/validation/run_refined_dfg.py +++ b/tests/validation/run_refined_dfg.py @@ -21,6 +21,32 @@ FINE_BLOCK = (8, 8) # 240 x 568 including four fine halo cells. 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.)]) + + +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): + """Cropped public macros only; all inventories, fluxes and impulses use root units.""" + from run_jet_circle_cv import cv_from_macros + + rectangles = CV_CONTOURS + (0 if refined else np.array([.25, 0., .25, 0.])) + origin = np.asarray(origin, dtype=float) + x0, y0, x1, y1 = np.rint((rectangles[-1] - np.tile(origin, 2)) / .5).astype(int) + selection = np.s_[y0-2:y1+2, x0-2:x1+2] + cropped_origin = origin + .5 * 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=.5) + for rect in rectangles] + h_native = 1.0 if refined else .5 + 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): @@ -233,11 +259,13 @@ def host_geometry(source, h, collision="SRT", jet_scenario="none"): return cfg, pair, body, host -def verify_source(source, collision="SRT", observation_mode="last_base_step", jet_scenario="none"): +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()) @@ -249,6 +277,8 @@ def verify_source(source, collision="SRT", observation_mode="last_base_step", je 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( @@ -317,12 +347,14 @@ def run_case(a, h, periodic, case): 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) context = stack.enter_context(CudaContext(0)) @@ -464,6 +496,15 @@ def run_case(a, h, periodic, case): 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 (.25, .5), + refined, begin, n, p["diameter"], impulse) + 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) @@ -522,6 +563,9 @@ def run_case(a, h, periodic, case): 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( @@ -577,12 +621,14 @@ def check_deadline(deadline): ) -def verify_pilot(path, manifest, collision="SRT", observation_mode="last_base_step", jet_scenario="none"): +def verify_pilot(path, manifest, collision="SRT", observation_mode="last_base_step", jet_scenario="none", observe_cv=False): if path is None: raise ValueError( "Full comparison requires --pilot-results from a completed matched pilot" ) previous = read_finite_json(path / "plan.json") + 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: @@ -599,7 +645,9 @@ def verify_pilot(path, manifest, collision="SRT", observation_mode="last_base_st for case, diameter in (("uniform", 80), ("refined", 40)): summary = read_finite_json(path / case / "summary.json") if ( - summary.get("jet_scenario", "none") != jet_scenario + 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 @@ -653,7 +701,8 @@ 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("jet_scenario", "none") != jet_scenario + 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 @@ -685,6 +734,8 @@ def main(): parser.add_argument("--collision", choices=["SRT", "MRT"], default="SRT") 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") @@ -698,14 +749,16 @@ def main(): a = parser.parse_args() 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) + 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") + if a.observe_cv: + a.provenance["observe_cv"] = True calibration = ( - verify_pilot(a.pilot_results, manifest, a.collision, a.observation_mode, a.jet_scenario) + verify_pilot(a.pilot_results, manifest, a.collision, a.observation_mode, a.jet_scenario, a.observe_cv) if a.mode == "full" else None ) if a.mode != "plan" and ( @@ -732,6 +785,17 @@ def main(): 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.tolist(), + uniform_contours_root=(CV_CONTOURS + [.25, 0., .25, 0.]).tolist(), + uniform_grid_origin_root=[.25, .5], refined_grid_origin_root=[38., 22.], + cell_width_root=.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", + ) host_geometry(a.source, h, a.collision, a.jet_scenario) (a.output / "plan.json").write_text(json.dumps(specification, indent=2)) if a.mode == "plan":