ADC-694: complete the generic HyQMOM15 closure and transaction proof

CHANGELOG.md

@@ -18,6 +18,12 @@ Format: [Keep a Changelog](https://keepachangelog.com/en/1.1.0/); versioning
 
 ### Changed
 
+- The final HyQMOM15 executable now checks realizability and the conserved `M00`
+  particle number for rejected, accepted, restored, and continued runtime snapshots. Its JSON
+  evidence reports the measured integral and maximum relative drift against the documented
+  `1e-10` acceptance threshold. It now authors the six fifth-order relations through a public
+  `@closure(4)` value and authenticates that every typed provisional store belongs to the rejected
+  Program transaction, without adding a HyQMOM-specific native route.
 - Strict AMR checkpoint payload v7 now persists the accepted shared-interface flux audit together
   with Program clocks, histories, tagging state, conservative ledger and synchronization report.
   Restart validates every fragment's topology epoch, level pair, exact clock window, resolved

docs/design/SPECIFICATION_TECHNIQUE_FINALE_POPS_ARCHITECTURE.md

@@ -1516,10 +1516,12 @@ Quatre scripts sont des tests d'acceptation, pas des esquisses :
    `AMRExecution.subcycled()`, regrid/reflux, HDF5/NPZ/ParaView, restart strict et continuation
    bit-identique ;
 4. `examples/final/EXEMPLE_SPEC_FINALE_15_MOMENTS_HYQMOM.py` : état 15 moments, layout Uniform,
-   `Program` IMEX explicite avec garde de réalisabilité dans sa transaction, champ de Poisson,
-   HDF5/ParaView et continuation bit-identique, sans branche de scénario dans le compilateur. Le
-   preset `pops.lib.time.IMEX` reste un constructeur d'un `Program` ordinaire ; il ne remplace pas
-   cette écriture explicite lorsqu'une garde scientifique spécifique doit être composée.
+   fermeture utilisateur `@closure(4)` abaissée dans le graphe de flux générique, `Program` IMEX
+   explicite avec garde de réalisabilité et ensemble complet des stores provisoires dans sa
+   transaction, champ de Poisson, conservation du nombre de particules, HDF5/ParaView et
+   continuation bit-identique, sans branche de scénario dans le compilateur. Le preset
+   `pops.lib.time.IMEX` reste un constructeur d'un `Program` ordinaire ; il ne remplace pas cette
+   écriture explicite lorsqu'une garde scientifique spécifique doit être composée.
 
 `scripts/final_release_contract.py` fixe cet ensemble exact : aucun cinquième script `.py` n'est admis
 dans `examples/final/`. Chaque script doit :

docs/design/hyqmom15-final-contract.md

@@ -16,12 +16,16 @@ gauge and multigrid solver remain separate `FieldDiscretization` choices on the
 
 ## Generic extension boundaries
 
-- `LocalClosure(order, name, evaluator)` is the closure extension interface. The evaluator executes
-  once on symbolic standardized moments during authoring and must return exactly the order `N + 1`
-  keys. It is absent from native execution.
-- `RealizabilityProjection` configures the smooth floors used by moment algebra. It does not pretend
-  to be a time-step acceptance guard. A future realizability rejection policy must implement the
-  ordinary typed `AcceptanceGuard` protocol and participate in the Program transaction explicitly.
+- `LocalClosure(order, name, evaluator)` is the closure extension interface. The final script writes
+  the six fifth-order HyQMOM relations under `@closure(4)` and passes that value to
+  `HyQMOM15.vlasov_lorentz(closure=...)`. The evaluator executes once on symbolic standardized
+  moments during authoring and must return exactly the order `N + 1` keys. Its arithmetic is folded
+  into the ordinary flux graph, so there is no Python callback or mutable closure state in native
+  execution; the installed Program hash authenticates the resulting graph across restart.
+- `RealizabilityProjection` configures the smooth floors and the complete 15-moment projection.
+  `guard_hyqmom15_candidate(...)` authors ordinary typed acceptance guards with
+  `ProjectAndRecheck(on_failure=RejectAttempt())` inside the `Program` transaction. Rejection and
+  rollback therefore use the shared runtime path rather than a HyQMOM-specific branch.
 - `Model.field_spaces()` derives solved storage from the generic field-output protocol. A scalar
   `FieldOutput` contributes one component; a Cartesian `GradientOutput` contributes two. This rule
   lets any provided or user model add a potential-plus-gradient solve without a model-specific
@@ -42,6 +46,12 @@ explicitly so its realizability guard is visibly inside the commit transaction.
 route, but it does not hide this model-specific scientific guard. The local solve is specialized from
 the resolved state manifest and therefore prepares exact 15 by 15 stack storage for the shared
 pivoted local provider, without an explicit inverse, eight-component fallback or family dispatch.
+The executable also requires the installed transaction plan to own every typed provisional store:
+states, fields, topology, flux ledgers, caches, solver warm starts, histories, clocks, schedules,
+consumers, diagnostics and external effects. Its forced non-realizable attempt compares the
+accepted state, solved fields, histories, Program identity and ConsumerGraph cursors before and
+after rejection and refuses any published artifact. This Uniform case has no non-empty AMR reflux
+ledger; non-empty multilevel ledger persistence remains the responsibility of the AMR final example.
 
 The example executes only:
 
@@ -53,4 +63,7 @@ One accepted step publishes authenticated HDF5, ParaView and scheduled checkpoin
 script reopens both scientific formats, creates a manual checkpoint, restores it into a fresh bind,
 compares the full 15-component state, solved field, clock, program identity and consumer cursors,
 then advances the uninterrupted and restarted instances one more step and requires exact equality.
-This is the final behavior, not a transition or compatibility example.
+Every retained state must remain realizable and conserve the integral of `M00` over the unit square
+within a relative tolerance of `1e-10`; the machine-readable report exposes the measured particle
+number and maximum relative error. This is the final behavior, not a transition or compatibility
+example.

examples/final/EXEMPLE_SPEC_FINALE_15_MOMENTS_HYQMOM.py

@@ -31,7 +31,7 @@
 from pops.lib.models.moments import HyQMOM15
 from pops.math import laplacian
 from pops.mesh import CartesianGrid, PeriodicAxes
-from pops.moments import RealizabilityProjection
+from pops.moments import RealizabilityProjection, closure
 from pops.numerics import DiscretizationPlan, reconstruction, riemann, variables
 from pops.numerics.reconstruction import limiters
 from pops.numerics.spatial import FiniteVolume
@@ -42,6 +42,7 @@
 from pops.solvers import DenseLU
 from pops.solvers.elliptic import GeometricMG
 from pops.time import (
+    ALL_PROVISIONAL_STORES,
     AdaptiveCFL,
     Dense,
     LocalLinear,
@@ -54,6 +55,45 @@
 
 DEFAULT_CELLS = 8
 DEFAULT_T_END = 1.0e-5
+PARTICLE_NUMBER_RELATIVE_TOLERANCE = 1.0e-10
+
+
+@closure(4)
+def user_hyqmom15_closure(standardized: Any) -> dict[str, Any]:
+    """Close the six fifth-order moments through the public local algebra contract."""
+
+    s03 = standardized["S03"]
+    s04 = standardized["S04"]
+    s11 = standardized["S11"]
+    s12 = standardized["S12"]
+    s13 = standardized["S13"]
+    s21 = standardized["S21"]
+    s22 = standardized["S22"]
+    s30 = standardized["S30"]
+    s31 = standardized["S31"]
+    s40 = standardized["S40"]
+    return {
+        "S50": 0.5 * s30 * (5.0 * s40 - 3.0 * s30 * s30 - 1.0),
+        "S41": (
+            -0.25 * s30 * (8.0 * s40 - 9.0 * s30 * s30 - 4.0) * s11
+            + 0.25 * (10.0 * s40 - 15.0 * s30 * s30 - 6.0) * s21
+            + 2.0 * s30 * s31
+        ),
+        "S32": (
+            0.5 * (2.0 * s40 - 3.0 * s30 * s30) * s12
+            + 0.5 * (3.0 * s22 - 1.0) * s30
+        ),
+        "S23": (
+            0.5 * (2.0 * s04 - 3.0 * s03 * s03) * s21
+            + 0.5 * (3.0 * s22 - 1.0) * s03
+        ),
+        "S14": (
+            -0.25 * s03 * (8.0 * s04 - 9.0 * s03 * s03 - 4.0) * s11
+            + 0.25 * (10.0 * s04 - 15.0 * s03 * s03 - 6.0) * s12
+            + 2.0 * s03 * s13
+        ),
+        "S05": 0.5 * s03 * (5.0 * s04 - 3.0 * s03 * s03 - 1.0),
+    }
 
 
 def _native_output_mode() -> ParallelMode:
@@ -93,6 +133,7 @@ class HyQMOM15Authoring:
     """All exact declarations retained across the public lifecycle."""
 
     model: Any
+    closure: Any
     case: Any
     state: Any
     state_instance: Any
@@ -116,10 +157,20 @@ class RuntimeSnapshot:
     fields: dict[str, np.ndarray]
     histories: dict[str, tuple[np.ndarray, ...]]
     program_hash: str
+    transaction_stores: tuple[str, ...]
     consumer_graph_identity: str
     consumer_cursors: dict[str, Any]
 
 
+@dataclass(frozen=True, slots=True)
+class PhysicalDiagnostics:
+    """Retained-state checks required by the HyQMOM15 specification."""
+
+    realizable: bool
+    particle_number: float
+    particle_number_relative_error: float
+
+
 @dataclass(frozen=True, slots=True)
 class ExecutionEvidence:
     """Scientific artifacts and exact states produced by one final execution."""
@@ -137,6 +188,8 @@ class ExecutionEvidence:
     restored: RuntimeSnapshot
     continuous: RuntimeSnapshot
     restarted: RuntimeSnapshot
+    reference_particle_number: float
+    physical_diagnostics: dict[str, PhysicalDiagnostics]
 
 
 def _guarded_imex_program(
@@ -220,6 +273,7 @@ def build_authoring(
         "unit_square", lower=(0.0, 0.0), upper=(1.0, 1.0),
     ).frame(Cartesian2D())
     model = HyQMOM15.vlasov_lorentz(
+        closure=user_hyqmom15_closure,
         q_over_m=ConstParam("q_over_m", -1.0),
         omega_c=ConstParam("omega_c", 0.5),
         projection=realizability,
@@ -295,6 +349,7 @@ def build_authoring(
     )))
     return HyQMOM15Authoring(
         model=model,
+        closure=user_hyqmom15_closure,
         case=case,
         state=state,
         state_instance=state_instance,
@@ -329,6 +384,55 @@ def build_initial_state(*, cells: int = DEFAULT_CELLS) -> dict[str, np.ndarray]:
     return {"plasma": state}
 
 
+def _particle_number(state: Any) -> float:
+    """Integrate ``M00`` over the unit square represented by cell averages."""
+
+    values = np.asarray(state, dtype=np.float64)
+    if values.ndim != 3 or values.shape[0] != len(HyQMOM15.components):
+        raise ValueError(
+            "HyQMOM15 diagnostics require a (15, ny, nx) cell-average state"
+        )
+    density = values[HyQMOM15.components.index("M00")]
+    if density.size == 0:
+        raise ValueError("HyQMOM15 diagnostics require at least one cell")
+    return float(np.sum(density, dtype=np.float64) / density.size)
+
+
+def _require_physical_diagnostics(
+    state: Any,
+    *,
+    projection: RealizabilityProjection,
+    reference_particle_number: float,
+    where: str,
+) -> PhysicalDiagnostics:
+    """Require finite, realizable moments and conservative particle number."""
+
+    values = np.asarray(state, dtype=np.float64)
+    if not np.isfinite(values).all():
+        raise RuntimeError("%s contains a non-finite moment" % where)
+    if (
+        not np.isfinite(reference_particle_number)
+        or reference_particle_number <= 0.0
+    ):
+        raise ValueError("reference particle number must be finite and positive")
+    realizable = bool(projection.is_hyqmom15_realizable(values))
+    if not realizable:
+        raise RuntimeError("%s is not HyQMOM15-realizable" % where)
+    particle_number = _particle_number(values)
+    scale = max(abs(reference_particle_number), np.finfo(np.float64).tiny)
+    relative_error = abs(particle_number - reference_particle_number) / scale
+    if relative_error > PARTICLE_NUMBER_RELATIVE_TOLERANCE:
+        raise RuntimeError(
+            "%s changed particle number by %.6e (limit %.6e)"
+            % (where, relative_error, PARTICLE_NUMBER_RELATIVE_TOLERANCE)
+        )
+    return PhysicalDiagnostics(
+        realizable=realizable,
+        particle_number=particle_number,
+        particle_number_relative_error=relative_error,
+    )
+
+
 def compile_final_case(
     *, cells: int = DEFAULT_CELLS, inject_nonrealizable: bool = False,
 ) -> tuple[HyQMOM15Authoring, Any, Any]:
@@ -355,6 +459,16 @@ def compile_final_case(
 
 
 def _snapshot(simulation: Any) -> RuntimeSnapshot:
+    program_report = simulation.program_report()
+    if not program_report.installed:
+        raise RuntimeError("HyQMOM15 runtime has no installed Program report")
+    transaction_stores = tuple(program_report.step_transaction.get("stores", ()))
+    expected_stores = tuple(store.value for store in ALL_PROVISIONAL_STORES)
+    if transaction_stores != expected_stores:
+        raise RuntimeError(
+            "HyQMOM15 transaction does not own every provisional store: %r"
+            % (transaction_stores,)
+        )
     fields = {
         slot: np.asarray(simulation.field_potential_global(slot), dtype=np.float64).copy()
         for slot in simulation.field_provider_slots()
@@ -373,14 +487,16 @@ def _snapshot(simulation: Any) -> RuntimeSnapshot:
         fields=fields,
         histories=histories,
         program_hash=str(simulation.installed_program_hash()),
+        transaction_stores=transaction_stores,
         consumer_graph_identity=simulation.consumer_graph.identity.token,
         consumer_cursors=simulation.consumer_cursors.to_data(),
     )
 
 
 def _require_same_snapshot(left: RuntimeSnapshot, right: RuntimeSnapshot, *, where: str) -> bool:
     for name in (
-        "time", "macro_step", "program_hash", "consumer_graph_identity", "consumer_cursors",
+        "time", "macro_step", "program_hash", "transaction_stores",
+        "consumer_graph_identity", "consumer_cursors",
     ):
         if getattr(left, name) != getattr(right, name):
             raise RuntimeError("%s changed %s across restart" % (where, name))
@@ -461,8 +577,9 @@ def run_and_restart(
     root.mkdir(parents=True, exist_ok=True)
     rejected_before, rejected_after, rejection_reason = \
         _run_rejected_nonrealizable_attempt(root, cells=cells)
-    _target, _resolved, artifact = compile_final_case(cells=cells)
+    target, _resolved, artifact = compile_final_case(cells=cells)
     initial = build_initial_state(cells=cells)
+    reference_particle_number = _particle_number(initial["plasma"])
     simulation = _bind_artifact(artifact, initial_state=initial)
     accepted_root = root / "accepted"
     run_report = pops.run(
@@ -499,6 +616,23 @@ def run_and_restart(
         resumed, t_end=final_time, max_steps=1, output_dir=root / "restarted")
     continuous, restarted = _snapshot(simulation), _snapshot(resumed)
     _require_same_snapshot(continuous, restarted, where="bit-identical continuation")
+    snapshots = {
+        "rejected_before": rejected_before,
+        "rejected_after": rejected_after,
+        "accepted": accepted,
+        "restored": restored,
+        "continuous": continuous,
+        "restarted": restarted,
+    }
+    physical_diagnostics = {
+        name: _require_physical_diagnostics(
+            snapshot.state,
+            projection=target.realizability,
+            reference_particle_number=reference_particle_number,
+            where=name.replace("_", " "),
+        )
+        for name, snapshot in snapshots.items()
+    }
 
     return ExecutionEvidence(
         hdf5_path=hdf5_path,
@@ -514,6 +648,8 @@ def run_and_restart(
         restored=restored,
         continuous=continuous,
         restarted=restarted,
+        reference_particle_number=reference_particle_number,
+        physical_diagnostics=physical_diagnostics,
     )
 
 
@@ -533,11 +669,20 @@ def main(argv: list[str] | None = None) -> None:
     print("checkpoint: %s" % evidence.manual_checkpoint_path)
     print("non-realizable rollback: %s" % rollback)
     print("bit-identical restart: True")
+    diagnostics = evidence.physical_diagnostics
+    restarted_diagnostics = diagnostics["restarted"]
     print("report: " + json.dumps({
         "finite": bool(np.isfinite(evidence.restarted.state).all()),
+        "realizable": all(value.realizable for value in diagnostics.values()),
         "n_moments": int(evidence.restarted.state.shape[0]),
+        "particle_number": restarted_diagnostics.particle_number,
+        "particle_number_reference": evidence.reference_particle_number,
+        "particle_number_relative_error": max(
+            value.particle_number_relative_error for value in diagnostics.values()),
+        "particle_number_relative_tolerance": PARTICLE_NUMBER_RELATIVE_TOLERANCE,
         "runtime_steps": evidence.restarted.macro_step,
         "runtime_time": evidence.restarted.time,
+        "rollback_stores": list(evidence.restarted.transaction_stores),
         "rejection_reason": evidence.rejection_reason,
         "nonrealizable_rollback": rollback,
         "scheduled_checkpoint": str(evidence.scheduled_checkpoint_path),

tests/python/examples/final/test_hyqmom15_final_example.py

@@ -7,6 +7,7 @@
 import sys
 
 import numpy as np
+import pytest
 
 
 ROOT = Path(__file__).resolve().parents[4]
@@ -32,11 +33,25 @@ def test_hyqmom15_example_runs_outputs_and_restarts_bit_identically(tmp_path) ->
     report_line = next(
         line for line in completed.stdout.splitlines() if line.startswith("report: "))
     report = json.loads(report_line.removeprefix("report: "))
+    from pops.time import ALL_PROVISIONAL_STORES
+
     assert report["finite"] is True
+    assert report["realizable"] is True
     assert report["n_moments"] == 15
+    assert report["particle_number"] == pytest.approx(
+        report["particle_number_reference"],
+        rel=report["particle_number_relative_tolerance"],
+    )
+    assert (
+        report["particle_number_relative_error"]
+        <= report["particle_number_relative_tolerance"]
+    )
     assert report["nonrealizable_rollback"] is True
     assert "hyqmom15_realizability_density" in report["rejection_reason"]
     assert report["runtime_steps"] == 2
+    assert report["rollback_stores"] == [
+        store.value for store in ALL_PROVISIONAL_STORES
+    ]
 
     from pops.output import HDF5, ParaView