Coverage for src / gwtransport / fronttracking / validation.py: 87%
102 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-04 18:51 +0000
1"""
2Physics validation utilities for front tracking.
4This module provides functions to verify physical correctness of
5front-tracking simulations, including entropy conditions, concentration
6bounds, and mass conservation.
8This file is part of gwtransport which is released under AGPL-3.0 license.
9See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details.
10"""
12import logging
14import numpy as np
15import pandas as pd
17from gwtransport.fronttracking.output import (
18 compute_cumulative_inlet_mass,
19 compute_total_outlet_mass,
20)
21from gwtransport.fronttracking.waves import RarefactionWave, ShockWave
23# Numerical tolerance constants
24EPSILON_CONCENTRATION_TOLERANCE = -1e-14 # Minimum allowed concentration (machine precision)
26logger = logging.getLogger(__name__)
29def verify_physics(structure, cout, cout_tedges, cin, *, verbose=True, rtol=1e-10):
30 """
31 Run comprehensive physics verification checks on front tracking results.
33 Performs the following checks:
34 1. Entropy condition for all shocks
35 2. No negative concentrations (within tolerance)
36 3. Output concentration ≤ input maximum
37 4. Finite first arrival time
38 5. No NaN values after spin-up period
39 6. Events chronologically ordered
40 7. Rarefaction concentration ordering (head vs tail)
41 8. Total integrated outlet mass (until all mass exits)
43 Parameters
44 ----------
45 structure : dict
46 Structure returned from infiltration_to_extraction_front_tracking_detailed.
47 Must contain keys: 'waves', 't_first_arrival', 'events', 'n_shocks',
48 'n_rarefactions', etc.
49 cout : array-like
50 Bin-averaged output concentrations.
51 cout_tedges : pandas.DatetimeIndex
52 Output time edges for bins.
53 cin : array-like
54 Input concentrations.
55 verbose : bool, optional
56 If True, print detailed results. If False, only return summary.
57 Default True.
58 rtol : float, optional
59 Relative tolerance for numerical checks. Default 1e-10.
61 Returns
62 -------
63 results : dict
64 Dictionary containing:
65 - 'all_passed': bool - True if all checks passed
66 - 'n_checks': int - Total number of checks performed
67 - 'n_passed': int - Number of checks that passed
68 - 'failures': list of str - Description of failed checks (empty if all passed)
69 - 'summary': str - One-line summary (e.g., "✓ All 8 checks passed")
71 Examples
72 --------
73 ::
75 results = verify_physics(structure, cout, cout_tedges, cin, verbose=False)
76 print(results["summary"])
77 # ✓ All 8 checks passed
78 assert results["all_passed"]
79 """
80 failures = []
81 checks = []
83 # Check 1: Entropy condition for all shocks
84 shocks = [w for w in structure["waves"] if isinstance(w, ShockWave)]
85 entropy_violations = [s for s in shocks if not s.satisfies_entropy()]
86 check1_pass = len(entropy_violations) == 0
87 checks.append({
88 "name": "Shock entropy condition",
89 "passed": check1_pass,
90 "message": f"Entropy violations: {len(entropy_violations)}/{len(shocks)} shocks",
91 })
92 if not check1_pass:
93 failures.append(f"Entropy violations: {len(entropy_violations)} shocks violate entropy condition")
95 # Check 2: No negative concentrations (within tolerance)
96 valid_cout = cout[~np.isnan(cout)]
97 min_cout = np.min(valid_cout) if len(valid_cout) > 0 else 0.0
98 check2_pass = min_cout >= EPSILON_CONCENTRATION_TOLERANCE
99 checks.append({
100 "name": "Non-negative concentrations",
101 "passed": check2_pass,
102 "message": f"Minimum concentration: {min_cout:.2e}",
103 })
104 if not check2_pass:
105 failures.append(f"Negative concentrations found: min = {min_cout:.2e}")
107 # Check 3: Output doesn't exceed input (within tight tolerance)
108 max_cout = np.max(valid_cout) if len(valid_cout) > 0 else 0.0
109 max_cin = np.max(cin)
110 check3_pass = max_cout <= max_cin * (1.0 + rtol)
111 checks.append({
112 "name": "Output ≤ input maximum",
113 "passed": check3_pass,
114 "message": f"Max output: {max_cout:.2f}, Max input: {max_cin:.2f}",
115 })
116 if not check3_pass:
117 failures.append(f"Output exceeds input: {max_cout:.2f} > {max_cin:.2f}")
119 # Check 4: Finite first arrival time
120 t_first = structure["t_first_arrival"]
121 check4_pass = np.isfinite(t_first)
122 checks.append({
123 "name": "Finite first arrival time",
124 "passed": check4_pass,
125 "message": f"First arrival: {t_first:.2f} days",
126 })
127 if not check4_pass:
128 failures.append(f"First arrival time is not finite: {t_first}")
130 # Check 5: No NaN values after spin-up period
131 cout_tedges_days = ((cout_tedges - cout_tedges[0]) / pd.Timedelta(days=1)).values
132 mask_after_spinup = cout_tedges_days[:-1] >= t_first
133 cout_after_spinup = cout[mask_after_spinup]
134 nan_count = np.sum(np.isnan(cout_after_spinup))
135 check5_pass = nan_count == 0
136 checks.append({
137 "name": "No NaN after spin-up",
138 "passed": check5_pass,
139 "message": f"NaN values after spin-up: {nan_count}/{len(cout_after_spinup)}",
140 })
141 if not check5_pass:
142 failures.append(f"Found {nan_count} NaN values after spin-up period")
144 # Check 6: Events chronologically ordered
145 event_times = [e["time"] for e in structure.get("events", [])]
146 if len(event_times) > 1:
147 is_ordered = all(event_times[i] <= event_times[i + 1] for i in range(len(event_times) - 1))
148 check6_pass = is_ordered
149 checks.append({
150 "name": "Events chronologically ordered",
151 "passed": check6_pass,
152 "message": f"{len(event_times)} events",
153 })
154 if not check6_pass:
155 failures.append("Events are not in chronological order")
156 else:
157 check6_pass = True
158 checks.append({
159 "name": "Events chronologically ordered",
160 "passed": True,
161 "message": f"{len(event_times)} events (N/A)",
162 })
164 # Check 7: Rarefaction concentration ordering
165 rarefactions = [w for w in structure["waves"] if isinstance(w, RarefactionWave)]
166 raref_ordering_violations = 0
167 for raref in rarefactions:
168 # For n>1 (higher C travels faster), head should be higher concentration (faster)
169 # For n<1 (lower C travels faster), head should be lower concentration (faster)
170 # We can check this via velocities: head_velocity should always be >= tail_velocity
171 if raref.head_velocity() < raref.tail_velocity() - 1e-10:
172 raref_ordering_violations += 1
174 check7_pass = raref_ordering_violations == 0
175 checks.append({
176 "name": "Rarefaction wave ordering",
177 "passed": check7_pass,
178 "message": f"Ordering violations: {raref_ordering_violations}/{len(rarefactions)} rarefactions",
179 })
180 if not check7_pass:
181 failures.append(f"{raref_ordering_violations} rarefactions have incorrect head/tail ordering")
183 # Check 8: Total integrated outlet mass (until all mass exits)
184 # This replaces the old mass balance check by computing total integrated outlet mass
185 # and comparing it to total inlet mass
186 tracker_state = structure.get("tracker_state")
187 if tracker_state is not None and hasattr(tracker_state, "flow"):
188 # Get simulation parameters
189 waves = structure["waves"]
190 v_outlet = tracker_state.v_outlet
191 sorption = tracker_state.sorption
192 flow_arr = tracker_state.flow
193 tedges_in = tracker_state.tedges
194 tedges_days = (tedges_in - tedges_in[0]) / pd.Timedelta(days=1)
196 # Total mass that entered domain
197 t_inlet_end = tedges_days.values[-1]
198 total_mass_in = compute_cumulative_inlet_mass(
199 t=t_inlet_end, cin=cin, flow=flow_arr, tedges_days=tedges_days.values
200 )
202 # Total mass that exits domain (integrating until all mass has exited)
203 total_mass_out, t_integration_end = compute_total_outlet_mass(
204 v_outlet=v_outlet, waves=waves, sorption=sorption, flow=flow_arr, tedges_days=tedges_days.values
205 )
207 # Check if inlet ends with explicit transition to C=0
208 # Without this, the front tracking solver assumes the inlet continues
209 # at the last concentration forever, leading to incorrect mass balance
210 epsilon_conc_zero = 1e-10 # Tolerance for checking if concentration is zero
211 if len(cin) > 0 and abs(cin[-1]) > epsilon_conc_zero and verbose:
212 msg = (
213 f"\n⚠️ WARNING: Inlet concentration ends at C={cin[-1]:.3f} (not zero).\n"
214 " For accurate mass balance, the inlet should explicitly end with C=0.\n"
215 " Consider appending cin[-1]=0.0 or adding a final time step with C=0.\n"
216 " Without this, the solver assumes inlet continues indefinitely,\n"
217 " causing mass balance errors in the validation check.\n"
218 )
219 print(msg) # noqa: T201
221 # Check if total outlet mass matches total inlet mass
222 if total_mass_in > 0:
223 relative_error_total = abs(total_mass_out - total_mass_in) / total_mass_in
224 else:
225 relative_error_total = abs(total_mass_out - total_mass_in)
227 # For integrated rarefaction mass, use slightly relaxed tolerance
228 # Analytical rarefaction integrals to infinity can have O(1e-7) precision
229 check8_pass = relative_error_total <= max(rtol, 1e-6)
230 checks.append({
231 "name": "Total integrated outlet mass",
232 "passed": check8_pass,
233 "message": f"Relative error: {relative_error_total:.2e} (integrated to t={t_integration_end:.1f} days)",
234 })
235 if not check8_pass:
236 failures.append(
237 f"Total outlet mass mismatch: relative_error={relative_error_total:.2e} > {rtol:.2e} "
238 f"(total_mass_out={total_mass_out:.6e}, total_mass_in={total_mass_in:.6e}, "
239 f"t_integration_end={t_integration_end:.1f} days)"
240 )
241 else:
242 # Skip if tracker state not available
243 check8_pass = True
244 checks.append({
245 "name": "Total integrated outlet mass",
246 "passed": True,
247 "message": "Skipped (tracker state not available)",
248 })
250 # Compile results
251 n_checks = len(checks)
252 n_passed = sum(c["passed"] for c in checks)
253 all_passed = len(failures) == 0
255 if all_passed:
256 summary = f"✓ All {n_checks} physics checks passed"
257 else:
258 summary = f"✗ {n_passed}/{n_checks} checks passed ({len(failures)} failures)"
260 results = {
261 "all_passed": all_passed,
262 "n_checks": n_checks,
263 "n_passed": n_passed,
264 "failures": failures,
265 "checks": checks,
266 "summary": summary,
267 }
269 # Log detailed output if verbose
270 if verbose:
271 logger.info("\nPhysics Verification:")
272 for i, check in enumerate(checks, 1):
273 status = "✓" if check["passed"] else "✗"
274 logger.info(" %d. %s: %s %s", i, check["name"], status, check["message"])
276 if all_passed:
277 logger.info("\n%s", summary)
278 else:
279 logger.warning("\n%s", summary)
280 logger.warning("\nFailures:")
281 for i, failure in enumerate(failures, 1):
282 logger.warning(" %d. %s", i, failure)
284 return results