diff --git a/CHANGELOG.md b/CHANGELOG.md index de49dbe01..a0f0a5a6d 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/docs/design/SPECIFICATION_TECHNIQUE_FINALE_POPS_ARCHITECTURE.md b/docs/design/SPECIFICATION_TECHNIQUE_FINALE_POPS_ARCHITECTURE.md index e8bb1baf1..3e49a84cc 100644 --- a/docs/design/SPECIFICATION_TECHNIQUE_FINALE_POPS_ARCHITECTURE.md +++ b/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 : diff --git a/docs/design/hyqmom15-final-contract.md b/docs/design/hyqmom15-final-contract.md index 4a3f8ea26..1e13997b3 100644 --- a/docs/design/hyqmom15-final-contract.md +++ b/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. diff --git a/examples/final/EXEMPLE_SPEC_FINALE_15_MOMENTS_HYQMOM.py b/examples/final/EXEMPLE_SPEC_FINALE_15_MOMENTS_HYQMOM.py index caaec9bf2..0a952c87e 100644 --- a/examples/final/EXEMPLE_SPEC_FINALE_15_MOMENTS_HYQMOM.py +++ b/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,6 +487,7 @@ 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(), ) @@ -380,7 +495,8 @@ def _snapshot(simulation: Any) -> RuntimeSnapshot: 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), diff --git a/tests/python/examples/final/test_hyqmom15_final_example.py b/tests/python/examples/final/test_hyqmom15_final_example.py index 3090b9656..b0e1c1dad 100644 --- a/tests/python/examples/final/test_hyqmom15_final_example.py +++ b/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 diff --git a/tests/python/unit/moments/test_hyqmom15_final_contract.py b/tests/python/unit/moments/test_hyqmom15_final_contract.py index 022d8071a..58162afc3 100644 --- a/tests/python/unit/moments/test_hyqmom15_final_contract.py +++ b/tests/python/unit/moments/test_hyqmom15_final_contract.py @@ -20,11 +20,25 @@ from pops.domain import RectangleFrame from pops.frames import Cartesian2D from pops.physics import Model -from pops.time import ProjectAndRecheck, RejectAttempt +from pops.time import ALL_PROVISIONAL_STORES, ProjectAndRecheck, RejectAttempt ROOT = Path(__file__).resolve().parents[4] EXAMPLE = ROOT / "examples/final/EXEMPLE_SPEC_FINALE_15_MOMENTS_HYQMOM.py" +_STANDARDIZED_SAMPLE = { + "S03": -0.2, + "S04": 2.8, + "S11": 0.15, + "S12": -0.35, + "S13": 0.42, + "S20": 1.0, + "S21": 0.25, + "S22": 1.2, + "S30": 0.3, + "S31": -0.1, + "S40": 3.1, + "S02": 1.0, +} def test_moment_flux_generator_is_public_for_explicit_python_models() -> None: @@ -83,22 +97,7 @@ def wrong_order(_standardized): def test_hyqmom15_closure_matches_closure_s5_matlab_oracle() -> None: """Pin the six non-Gaussian polynomial relations used by closureS5.m.""" - standardized = { - "S03": -0.2, - "S04": 2.8, - "S11": 0.15, - "S12": -0.35, - "S13": 0.42, - "S20": 1.0, - "S21": 0.25, - "S22": 1.2, - "S30": 0.3, - "S31": -0.1, - "S40": 3.1, - "S02": 1.0, - } - - closed = HyQMOM15Closure()(standardized) + closed = HyQMOM15Closure()(_STANDARDIZED_SAMPLE) assert closed == pytest.approx({ "S50": 2.1345, @@ -148,12 +147,22 @@ def test_final_authoring_derives_field_storage_and_complete_generic_program() -> target = _load_example().build_authoring() assert type(target.model) is Model + assert type(target.closure) is LocalClosure + assert target.closure.contract_data() == { + "kind": "local_moment_closure", + "order": 4, + "name": "user_hyqmom15_closure", + } + assert target.closure(_STANDARDIZED_SAMPLE) == pytest.approx( + HyQMOM15Closure()(_STANDARDIZED_SAMPLE) + ) assert isinstance(target.model.frame, RectangleFrame) assert target.components == tuple(moment_names(4)) assert target.model.field_spaces()[target.field.local_id].components == ( "phi", "grad_x", "grad_y") assert target.field_provider == target.model.operators["fields"] assert target.program.transaction_plan() is not None + assert target.program.transaction_plan().stores == ALL_PROVISIONAL_STORES guards = target.program.transaction_plan().guards assert [guard.name for guard in guards] == [ "hyqmom15_realizability_density", @@ -169,6 +178,36 @@ def test_final_authoring_derives_field_storage_and_complete_generic_program() -> assert projection.kind == "projection" +def test_particle_number_diagnostic_integrates_m00_and_rejects_drift() -> None: + example = _load_example() + target = example.build_authoring() + state = example.build_initial_state(cells=4)["plasma"] + reference = example._particle_number(state) + + assert reference == pytest.approx(1.0) + diagnostics = example._require_physical_diagnostics( + state, + projection=target.realizability, + reference_particle_number=reference, + where="initial state", + ) + assert diagnostics.realizable is True + assert diagnostics.particle_number == pytest.approx(reference) + assert diagnostics.particle_number_relative_error == pytest.approx(0.0) + + drifted = state.copy() + drifted[HyQMOM15.components.index("M00")] += ( + 2.0 * example.PARTICLE_NUMBER_RELATIVE_TOLERANCE + ) + with pytest.raises(RuntimeError, match="changed particle number"): + example._require_physical_diagnostics( + drifted, + projection=target.realizability, + reference_particle_number=reference, + where="drifted state", + ) + + def test_hyqmom15_projection_checks_all_moments_and_refuses_to_manufacture_density() -> None: example = _load_example() target = example.build_authoring()