test: add bounded dense-root Jet CV temporal diagnostics

tests/validation/jet_cv_measurements.py

@@ -10,13 +10,63 @@ LINEAR = np.array([2.5, -1.5, 0.0])
 QUADRATIC = np.array([4.375, -5.25, 1.875])
 
 
+EVIDENCE_LIMIT = 32 * 1024**2
+
+
+def install_evidence(replacements, limit=None):
+    """Preflight retained bytes before replacing files; cap failure preserves old pair.
+
+    Staged compression bytes are temporary disk overhead, not retained evidence.
+    This guards budget failures, not a transactional crash-recovery filesystem.
+    """
+    temporary = set(replacements.values())
+    try:
+        directory = next(iter(replacements)).parent
+        retained = sum(
+            p.stat().st_size
+            for p in directory.iterdir()
+            if p.is_file() and p not in temporary
+        )
+        projected = retained - sum(p.stat().st_size for p in replacements if p.exists())
+        projected += sum(p.stat().st_size for p in temporary)
+        if limit is not None and projected > limit:
+            raise RuntimeError(
+                "Dense-root 32 MiB retained-evidence budget would be exceeded; last verified files preserved"
+            )
+        for target, staged in replacements.items():
+            staged.replace(target)
+    finally:
+        for p in temporary:
+            if p.exists():
+                p.unlink()
+
+
+def write_evidence_text(path, text, *, dense_root=False):
+    staged = path.with_name("." + path.name + ".pending")
+    staged.write_text(text)
+    install_evidence({path: staged}, EVIDENCE_LIMIT if dense_root else None)
+
+
 def sample_steps(c):
     """Same reference times at both scales; dense startup retained."""
+    if c.get("dense_root", False):
+        return np.arange(c["steps"] + 1, dtype=int)
+    return old_sample_steps(c)
+
+
+def old_sample_steps(c):
+    """Exact previous physical cadence, clipped to this campaign's endpoint."""
     s2 = c["scale"] ** 2
-    reference = sorted(set(range(0, 101)) | set(range(100, 3201, 4)))
+    end = c.get("steps", 3200 * s2) // s2
+    reference = sorted(set(range(0, min(100, end) + 1)) | set(range(100, end + 1, 4)))
     return np.array(reference, int) * s2
 
 
+def absolute_aperture_interval_mass(rows, blowing, suction):
+    """Sum absolute aperture NET masses for one root interval, not per-link flux."""
+    return float(abs(rows[blowing, 3].sum()) + abs(rows[suction, 3].sum()))
+
+
 def factors(scale):
     """Lref=one s=1 root length, Tref=one s=1 root time, rho_ref=1."""
     return dict(
@@ -105,7 +155,7 @@ def ray_geometry(fluid, origin, center, radius, h=0.5):
     )
 
 
-def reconstruct_aperture(rho, vel, geo):
+def reconstruct_aperture(rho, vel, geo, *, return_absolute=False):
     field = np.concatenate((rho[None], vel), axis=0)
     x, y, a, b = geo["ix"], geo["iy"], geo["a"], geo["b"]
     ray = (
@@ -115,12 +165,24 @@ def reconstruct_aperture(rho, vel, geo):
         + field[:, y + 1, x + 1] * a * b
     )
     estimates = []
+    absolute = []
     for coeff in (LINEAR, QUADRATIC):
         wall = np.einsum("cajk,k->caj", ray, coeff)
         un = np.einsum("caj,ajc->aj", wall[1:], geo["normal"])
         q = np.sum(un * geo["ds"], axis=1)
         mass = np.sum(wall[0] * un * geo["ds"], axis=1)
         estimates.append(np.stack((q, mass), axis=-1))
+        absolute.append(
+            np.stack(
+                (
+                    np.sum(abs(un) * geo["ds"], axis=1),
+                    np.sum(abs(wall[0] * un) * geo["ds"], axis=1),
+                ),
+                axis=-1,
+            )
+        )
+    if return_absolute:
+        return ray, np.array(estimates), np.array(absolute)
     return ray, np.array(estimates)
 
 
@@ -136,6 +198,8 @@ class MeasurementRecorder:
         self.cx = cx
         self.cy = cy
         self.records = []
+        self.sparse_rays = []
+        self.ray_steps = set(old_sample_steps(c).tolist())
         self.rays = [
             ray_geometry(fluid_u, (0.0, 0.0), np.array(c["center"]), c["diameter"] / 2),
             ray_geometry(leaves[1], origin_a, np.array(c["center"]), c["diameter"] / 2),
@@ -179,8 +243,16 @@ class MeasurementRecorder:
         record["cv_inventory"] = []
         record["cv_flux"] = []
         record["momentum_flux_terms"] = []
-        record["arc_ray_values"] = []
-        record["arc_estimates"] = []
+        save_ray = not c.get("dense_root", False) or step in self.ray_steps
+        ray_record = dict(
+            sample_index=len(self.records),
+            arc_ray_values=[],
+            arc_estimates=[],
+            arc_absolute_estimates=[],
+        )
+        if not c.get("dense_root", False):
+            record["arc_ray_values"] = []
+            record["arc_estimates"] = []
         record["nearest_mass_flux"] = [
             nearest_flux[k] for k in ("uniform_adapter", "amr_api")
         ]
@@ -203,9 +275,18 @@ class MeasurementRecorder:
                 rtol=2e-10,
             )
             record["momentum_flux_terms"].append(terms)
-            ray, estimates = reconstruct_aperture(rho, vel, self.rays[idx])
-            record["arc_ray_values"].append(ray)
-            record["arc_estimates"].append(estimates)
+            if save_ray:
+                ray, estimates, absolute = reconstruct_aperture(
+                    rho, vel, self.rays[idx], return_absolute=True
+                )
+                ray_record["arc_ray_values"].append(ray)
+                ray_record["arc_estimates"].append(estimates)
+                ray_record["arc_absolute_estimates"].append(absolute)
+                if not c.get("dense_root", False):
+                    record["arc_ray_values"].append(ray)
+                    record["arc_estimates"].append(estimates)
+        if c.get("dense_root", False) and save_ray:
+            self.sparse_rays.append(ray_record)
         self.records.append(record)
         if step // s**2 in (0, 100, 800, 1600, 3200) and step % s**2 == 0:
             fields = {}
@@ -223,7 +304,11 @@ class MeasurementRecorder:
                 fields[name + "_origin_reference"] = np.asarray(origin) / s
                 fields[name + "_h_reference"] = np.array(h / s)
             file = self.output / f"fields_s{s}_tref{step//s**2}.npz"
-            np.savez_compressed(file, **fields)
+            staged = file.with_name("." + file.name + ".pending.npz")
+            np.savez_compressed(staged, **fields)
+            install_evidence(
+                {file: staged}, EVIDENCE_LIMIT if c.get("dense_root", False) else None
+            )
             self.snapshots.append(
                 dict(
                     file=file.name,
@@ -236,19 +321,32 @@ class MeasurementRecorder:
     def flush(self):
         path = self.output / f"samples_s{self.c['scale']}.npz"
         arrays = {k: np.array([r[k] for r in self.records]) for k in self.records[0]}
+        if self.c.get("dense_root", False):
+            arrays["ray_sample_indices"] = np.array(
+                [r["sample_index"] for r in self.sparse_rays], dtype=int
+            )
+            for key in ("arc_ray_values", "arc_estimates", "arc_absolute_estimates"):
+                arrays[key] = np.array([r[key] for r in self.sparse_rays])
         for i, g in enumerate(self.rays):
             for k in ("normal", "ds", "positions"):
                 arrays[f"ray_{i}_{k}"] = g[k]
         arrays["ray_offsets_fine_cells"] = OFFSETS
         arrays["linear_coefficients"] = LINEAR
         arrays["quadratic_coefficients"] = QUADRATIC
-        np.savez_compressed(path, **arrays)
+        staged = path.with_name("." + path.name + ".pending.npz")
+        np.savez_compressed(staged, **arrays)
         metadata = dict(
-            schema=1,
+            schema=2 if self.c.get("dense_root", False) else 1,
             case=self.c,
+            dense_root=bool(self.c.get("dense_root", False)),
+            ray_sample_index_key=(
+                "ray_sample_indices" if self.c.get("dense_root", False) else None
+            ),
+            raw_absolute_throughput_definition="Cumulative sum of abs(blowing NET dm)+abs(suction NET dm) per root interval; unresolved within-aperture/substep reversals; not abs(window net).",
+            absolute_throughput_definition="Quadrature integral of abs(reconstructed u.n), abs(reconstructed rho*u.n); not abs(signed integral), not prescribed nominal, not exact wall truth.",
             reference_conversion=factors(self.c["scale"]),
             sample_file=path.name,
-            sample_sha256=hashlib.sha256(path.read_bytes()).hexdigest(),
+            sample_sha256=hashlib.sha256(staged.read_bytes()).hexdigest(),
             sample_count=len(self.records),
             snapshots=self.snapshots,
             level_stats_columns=[
@@ -280,6 +378,10 @@ class MeasurementRecorder:
             ],
             warning="Extrapolation and nearest-node values are estimator-dependent, not exact wall values. Fields/ray arrays use root-lattice units except explicitly named reference arrays.",
         )
-        (self.output / f"samples_s{self.c['scale']}.json").write_text(
-            json.dumps(metadata, indent=2) + "\n"
+        meta_path = self.output / f"samples_s{self.c['scale']}.json"
+        meta_staged = meta_path.with_name("." + meta_path.name + ".pending")
+        meta_staged.write_text(json.dumps(metadata, indent=2) + "\n")
+        install_evidence(
+            {path: staged, meta_path: meta_staged},
+            EVIDENCE_LIMIT if self.c.get("dense_root", False) else None,
         )