test: add optional public macro CV observations

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

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":