Coverage for src/gwtransport/fronttracking/validation.py: 95%
105 statements
« prev ^ index » next coverage.py v7.15.3, created at 2026-08-04 21:13 +0000
1"""
2Physics validation utilities for front tracking in (V, θ) coordinates.
4This module provides functions to verify physical correctness of front-tracking
5simulations, including entropy conditions, concentration bounds, mass conservation,
6and event ordering. The solver runs in cumulative-flow coordinate
7``θ = ∫flow(t') dt'``; events on ``state.events`` carry ``"theta"`` (m³). Because
8``flow ≥ 0`` is enforced, θ is monotone non-decreasing in t, so θ-ordering and
9chronological ordering are equivalent.
11Available functions:
13- :func:`verify_physics` - Run seven checks on a completed front-tracking result and return a summary
14 dictionary with ``'all_passed'``, the check counts, the ``'failures'`` list, a per-check record list and a
15 one-line ``'summary'``: Lax entropy for every shock, no negative concentrations, output bounded by the
16 inlet maximum, finite first-arrival θ, no NaN after spin-up, θ-ordered events, and mass conservation
17 comparing an independent outlet integral plus the domain mass against the injected mass at ``θ_max``. The
18 mass-balance check uses ``max(rtol, _MASS_BALANCE_RTOL)`` because integrating a shock-bearing breakthrough
19 curve is only first-order accurate. Passing ``verbose=False`` suppresses printing but returns the same
20 dictionary.
22This file is part of gwtransport which is released under AGPL-3.0 license.
23See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details.
24"""
26import logging
28import numpy as np
29import numpy.typing as npt
30import pandas as pd
32from gwtransport._time import tedges_to_days
33from gwtransport.fronttracking.output import (
34 compute_breakthrough_curve,
35 compute_cumulative_inlet_mass,
36 compute_domain_mass,
37)
38from gwtransport.fronttracking.solver import FrontTrackerState
39from gwtransport.fronttracking.waves import ShockWave
41# Numerical tolerance constants
42EPSILON_CONCENTRATION_TOLERANCE = -1e-14 # Minimum allowed concentration (machine precision)
44# Mass-balance check (7) tolerance and grid.
45#
46# The independent outlet mass integrates the breakthrough curve with the trapezoid
47# rule (see ``_independent_outlet_mass``). For a shock-bearing, sharply-curved
48# breakthrough this is only first-order accurate: the measured relative error for the
49# canonical favorable-sorption pulse oscillates at ~5e-4 across a modest grid band
50# (2000-4000 points; it only falls reliably below 1e-4 above ~20000 points). The grid is
51# deliberately kept modest because a numerical DecayingShockWave makes
52# ``compute_breakthrough_curve`` slow (seconds to minutes for large grids), so we cannot
53# refine the integral to machine precision. ``_MASS_BALANCE_RTOL`` therefore bounds that
54# ~5e-4 grid noise with ~16x margin; a physical 30% inlet-mass error yields a relative
55# error of ~0.23 (= 1 - 1/1.3) to ~0.43 (= 1/0.7 - 1), i.e. ~20-40x this floor, so the
56# check still has strong teeth against a genuine conservation failure.
57_MASS_BALANCE_RTOL = 1e-2
58_MASS_BALANCE_GRID_POINTS = 3000
60logger = logging.getLogger(__name__)
63def _independent_outlet_mass(tracker_state: FrontTrackerState) -> float:
64 """Outlet-side mass total computed independently of the ``m_in - m_dom`` identity.
66 Integrated to θ_max (the last θ-bin edge). Sums the mass that has already left through
67 the outlet, ``∫₀^θ_max c_out(τ) dτ``, and the mass still in the domain, ``m_dom(θ_max)``.
68 The breakthrough integral uses :func:`compute_breakthrough_curve`, which dispatches
69 :func:`concentration_at_point` directly (pure wave evaluation), so this total never
70 references the conservation identity ``m_out = m_in − m_dom`` that the mass-balance
71 check is meant to test. Comparing it to :func:`compute_cumulative_inlet_mass` at θ_max
72 is therefore a genuine, non-tautological conservation check: for a pulse that has not
73 fully broken through by θ_max, the partial breakthrough integral plus the residual
74 domain mass still equals the cumulative inlet mass.
76 Parameters
77 ----------
78 tracker_state : FrontTrackerState
79 Solver state; must expose ``v_outlet``, ``sorption``, ``waves`` and
80 ``theta_edges``.
82 Returns
83 -------
84 float
85 Independent outlet-side mass total [mass].
86 """
87 v_outlet = tracker_state.v_outlet
88 sorption = tracker_state.sorption
89 waves = tracker_state.waves
90 theta_max = float(np.asarray(tracker_state.theta_edges, dtype=float)[-1])
92 if theta_max <= 0.0:
93 return compute_domain_mass(theta=theta_max, v_outlet=v_outlet, waves=waves, sorption=sorption)
95 theta_grid: npt.NDArray[np.floating] = np.linspace(0.0, theta_max, _MASS_BALANCE_GRID_POINTS)
96 breakthrough = compute_breakthrough_curve(theta_grid, v_outlet, waves, sorption)
97 mass_out = float(np.trapezoid(breakthrough, theta_grid))
98 mass_dom = compute_domain_mass(theta=theta_max, v_outlet=v_outlet, waves=waves, sorption=sorption)
99 return mass_out + mass_dom
102def verify_physics(
103 structure: dict,
104 cout: npt.ArrayLike,
105 cout_tedges: pd.DatetimeIndex,
106 cin: npt.ArrayLike,
107 *,
108 verbose: bool = True,
109 rtol: float = 1e-10,
110) -> dict:
111 """
112 Run comprehensive physics verification checks on front tracking results.
114 Performs the following checks:
116 1. Entropy condition for all shocks
117 2. No negative concentrations (within tolerance)
118 3. Output concentration <= input maximum
119 4. Finite first arrival θ
120 5. No NaN values after spin-up period
121 6. Events θ-ordered (equivalent to chronological under non-negative flow)
122 7. Mass conservation: independent outlet integral + domain mass == inlet mass at θ_max
124 Parameters
125 ----------
126 structure : dict
127 Structure returned from ``infiltration_to_extraction_nonlinear_sorption``.
128 Must contain keys: ``'waves'``, ``'theta_first_arrival'``, ``'events'``,
129 and optionally ``'tracker_state'``.
130 cout : array-like
131 Bin-averaged output concentrations.
132 cout_tedges : pandas.DatetimeIndex
133 Output time edges for bins (only used for the spin-up mask).
134 cin : array-like
135 Input concentrations.
136 verbose : bool, optional
137 If True, print detailed results. If False, only return summary. Default True.
138 rtol : float, optional
139 Relative tolerance for numerical checks. Default 1e-10. For the mass-balance
140 check (7) the effective tolerance is ``max(rtol, _MASS_BALANCE_RTOL)`` because
141 that check integrates a shock-bearing breakthrough curve and is only first-order
142 accurate (see ``_MASS_BALANCE_RTOL``).
144 Returns
145 -------
146 results : dict
147 Dictionary containing:
149 - ``'all_passed'``: bool - True if all checks passed
150 - ``'n_checks'``: int - Total number of checks performed
151 - ``'n_passed'``: int - Number of checks that passed
152 - ``'failures'``: list of str - Description of failed checks (empty if all passed)
153 - ``'checks'``: list of dict - Per-check result records; each has ``'name'``,
154 ``'passed'``, ``'message'`` keys.
155 - ``'summary'``: str - One-line summary
157 Examples
158 --------
159 .. disable_try_examples
161 ::
163 results = verify_physics(structure, cout, cout_tedges, cin, verbose=False)
164 print(results["summary"])
165 assert results["all_passed"]
166 """
167 cout = np.asarray(cout, dtype=float)
168 cin = np.asarray(cin, dtype=float)
169 failures: list[str] = []
170 checks: list[dict] = []
172 # Check 1: Entropy condition for all shocks (Lax in (V, θ): λ_θ(C_L) >= s >= λ_θ(C_R)).
173 shocks = [w for w in structure["waves"] if isinstance(w, ShockWave)]
174 entropy_violations = [s for s in shocks if not s.satisfies_entropy()]
175 check1_pass = len(entropy_violations) == 0
176 checks.append({
177 "name": "Shock entropy condition",
178 "passed": check1_pass,
179 "message": f"Entropy violations: {len(entropy_violations)}/{len(shocks)} shocks",
180 })
181 if not check1_pass:
182 failures.append(f"Entropy violations: {len(entropy_violations)} shocks violate entropy condition")
184 # Check 2: No negative concentrations (within tolerance)
185 valid_cout = cout[~np.isnan(cout)]
186 min_cout = np.min(valid_cout) if len(valid_cout) > 0 else 0.0
187 check2_pass = min_cout >= EPSILON_CONCENTRATION_TOLERANCE
188 checks.append({
189 "name": "Non-negative concentrations",
190 "passed": check2_pass,
191 "message": f"Minimum concentration: {min_cout:.2e}",
192 })
193 if not check2_pass:
194 failures.append(f"Negative concentrations found: min = {min_cout:.2e}")
196 # Check 3: Output doesn't exceed input (within tight tolerance)
197 max_cout = np.max(valid_cout) if len(valid_cout) > 0 else 0.0
198 max_cin = np.max(cin)
199 check3_pass = max_cout <= max_cin * (1.0 + rtol)
200 checks.append({
201 "name": "Output <= input maximum",
202 "passed": check3_pass,
203 "message": f"Max output: {max_cout:.2f}, Max input: {max_cin:.2f}",
204 })
205 if not check3_pass:
206 failures.append(f"Output exceeds input: {max_cout:.2f} > {max_cin:.2f}")
208 # Check 4: Finite first arrival θ
209 theta_first = structure["theta_first_arrival"]
210 check4_pass = np.isfinite(theta_first)
211 checks.append({
212 "name": "Finite first arrival θ",
213 "passed": check4_pass,
214 "message": f"First arrival: θ={theta_first:.2f}",
215 })
216 if not check4_pass:
217 failures.append(f"First arrival θ is not finite: {theta_first}")
219 # Check 5: No NaN values after spin-up
220 tracker_state = structure.get("tracker_state")
221 if tracker_state is not None and np.isfinite(theta_first):
222 # theta_at_t measures days from the input origin tracker_state.tedges[0]; the output
223 # grid must be referenced to that same origin (not its own first edge) or the mask
224 # shifts and NaNs after spin-up slip through.
225 t_days = tedges_to_days(cout_tedges, ref=tracker_state.tedges[0])[:-1]
226 theta_at_edge = tracker_state.theta_at_t_array(t_days)
227 mask_after_spinup = theta_at_edge >= theta_first
228 elif not np.isfinite(theta_first):
229 # No spin-up bound — every output row counts as "after spin-up".
230 mask_after_spinup = np.ones(len(cout), dtype=bool)
231 else:
232 # No tracker state to translate — nothing to check.
233 mask_after_spinup = np.zeros(len(cout), dtype=bool)
234 cout_after_spinup = cout[mask_after_spinup]
235 nan_count = np.sum(np.isnan(cout_after_spinup))
236 check5_pass = nan_count == 0
237 checks.append({
238 "name": "No NaN after spin-up",
239 "passed": check5_pass,
240 "message": f"NaN values after spin-up: {nan_count}/{len(cout_after_spinup)}",
241 })
242 if not check5_pass:
243 failures.append(f"Found {nan_count} NaN values after spin-up period")
245 # Check 6: Events θ-ordered. ``np.all(np.diff(...) >= 0)`` is vacuously True
246 # for an empty/singleton sequence, so the same expression covers the N/A case;
247 # only the message differs.
248 event_thetas = [e["theta"] for e in structure.get("events", [])]
249 check6_pass = bool(np.all(np.diff(event_thetas) >= 0))
250 checks.append({
251 "name": "Events θ-ordered",
252 "passed": check6_pass,
253 "message": f"{len(event_thetas)} events" if len(event_thetas) > 1 else f"{len(event_thetas)} events (N/A)",
254 })
255 if not check6_pass:
256 failures.append("Events are not θ-ordered")
258 # Check 7: Total integrated outlet mass vs total inlet mass (in θ-space).
259 #
260 # The outlet-side total is computed *independently* of the conservation identity
261 # ``m_out = m_in − m_dom``; using that identity on both sides would be an algebraic
262 # tautology that passes for any input. ``_independent_outlet_mass`` integrates the
263 # breakthrough curve and adds the spatial domain mass; both come from direct wave
264 # evaluation, so a mismatch with the cumulative inlet mass signals a real conservation
265 # failure. Integrated to θ_max (the last θ-bin edge); for pulses that have not fully
266 # broken through there, the partial breakthrough integral plus the residual domain mass
267 # still equals the cumulative inlet mass.
268 if tracker_state is not None and hasattr(tracker_state, "theta_edges"):
269 theta_edges_arr = np.asarray(tracker_state.theta_edges, dtype=float)
270 theta_integration_end = float(theta_edges_arr[-1])
272 total_mass_in = compute_cumulative_inlet_mass(theta=theta_integration_end, cin=cin, theta_edges=theta_edges_arr)
273 independent_mass_out = _independent_outlet_mass(tracker_state)
275 if total_mass_in > 0:
276 relative_error_total = abs(independent_mass_out - total_mass_in) / total_mass_in
277 else:
278 relative_error_total = abs(independent_mass_out - total_mass_in)
280 mass_balance_threshold = max(rtol, _MASS_BALANCE_RTOL)
281 check7_pass = relative_error_total <= mass_balance_threshold
282 checks.append({
283 "name": "Total integrated outlet mass",
284 "passed": check7_pass,
285 "message": (
286 f"Relative error: {relative_error_total:.2e} (independent outlet integral to "
287 f"θ={theta_integration_end:.1f}; threshold {mass_balance_threshold:.2e})"
288 ),
289 })
290 if not check7_pass:
291 failures.append(
292 f"Total outlet mass mismatch: relative_error={relative_error_total:.2e} > "
293 f"{mass_balance_threshold:.2e} (independent_mass_out={independent_mass_out:.6e}, "
294 f"total_mass_in={total_mass_in:.6e}, θ_integration_end={theta_integration_end:.1f})"
295 )
296 else:
297 check7_pass = True
298 checks.append({
299 "name": "Total integrated outlet mass",
300 "passed": True,
301 "message": "Skipped (tracker state not available)",
302 })
304 # Compile results
305 n_checks = len(checks)
306 n_passed = sum(c["passed"] for c in checks)
307 all_passed = len(failures) == 0
309 if all_passed:
310 summary = f"All {n_checks} physics checks passed"
311 else:
312 summary = f"{n_passed}/{n_checks} checks passed ({len(failures)} failures)"
314 results = {
315 "all_passed": all_passed,
316 "n_checks": n_checks,
317 "n_passed": n_passed,
318 "failures": failures,
319 "checks": checks,
320 "summary": summary,
321 }
323 if verbose:
324 logger.info("\nPhysics Verification:")
325 for i, check in enumerate(checks, 1):
326 status = "PASS" if check["passed"] else "FAIL"
327 logger.info(" %d. %s: %s %s", i, check["name"], status, check["message"])
329 if all_passed:
330 logger.info("\n%s", summary)
331 else:
332 logger.warning("\n%s", summary)
333 logger.warning("\nFailures:")
334 for i, failure in enumerate(failures, 1):
335 logger.warning(" %d. %s", i, failure)
337 return results