Coverage for src/gwtransport/fronttracking/events.py: 95%

1"""Event detection for front tracking in (V, θ) coordinates. 

2 

3All intersections are pure line/line geometry in the (V, θ) plane because 

4every wave speed dV/dθ is independent of flow. Functions return θ-coordinates 

5of intersections; the solver translates to user-facing t at the API boundary. 

6 

7Events include: 

8 

9- Characteristic-characteristic collisions 

10- Shock-shock collisions 

11- Shock-characteristic collisions 

12- Rarefaction boundary interactions 

13- Outlet crossings 

14 

15All calculations return exact floating-point results with machine precision. 

16""" 

17 

18from dataclasses import dataclass 

19from enum import Enum 

20 

21from gwtransport.fronttracking.math import characteristic_position, characteristic_speed 

22from gwtransport.fronttracking.waves import CharacteristicWave, DecayingShockWave, RarefactionWave, ShockWave 

23 

24# Numerical tolerance constants 

25EPSILON_SPEED = 1e-15 # Tolerance for checking if two speeds are equal (machine precision) 

26# A boundary state at/below the c_min retardation floor whose floored retardation 

27# exceeds this value is "pinned": for the n>1 dry-soil singularity R(c_min) is 

28# inflated to ~1e6, so the state advects orders of magnitude slower than any 

29# physical wave and its outlet crossing lands at a non-physical θ (~1e8) that only 

30# pollutes the diagnostic event record. n<1 clean water (R(c_min) ≈ 1, fast) stays 

31# well below this threshold and is NOT pinned, so its outlet crossing is kept. 

32OUTLET_PIN_RETARDATION = 1e4 

33 

34 

35def is_outlet_crossing_pinned(concentration: float, sorption) -> bool: 

36 """Whether a boundary state is pinned by the ``c_min`` retardation floor. 

37 

38 A crossing scheduled for such a state is a non-physical artifact (its speed is 

39 a floor artifact, not physics); the caller drops it so it does not pollute the 

40 solver's event record / ``theta_current``. 

41 

42 Parameters 

43 ---------- 

44 concentration : float 

45 Boundary-state concentration [mass/volume]. 

46 sorption : SorptionModel 

47 Sorption model (supplies ``c_min`` and ``retardation``). 

48 

49 Returns 

50 ------- 

51 bool 

52 ``True`` only when ``concentration`` is at/below ``c_min`` AND the floored 

53 retardation ``R(c_min)`` is inflated past :data:`OUTLET_PIN_RETARDATION`. 

54 """ 

55 c_min = getattr(sorption, "c_min", 0.0) 

56 if concentration > c_min: 

57 return False 

58 return float(sorption.retardation(c_min)) > OUTLET_PIN_RETARDATION 

59 

60 

61class EventType(Enum): 

62 """All possible event types in front tracking simulation.""" 

63 

64 CHAR_CHAR_COLLISION = "characteristic_collision" 

65 """Two characteristics intersect (will form shock).""" 

66 SHOCK_SHOCK_COLLISION = "shock_collision" 

67 """Two shocks collide (will merge).""" 

68 SHOCK_CHAR_COLLISION = "shock_characteristic_collision" 

69 """Shock catches or is caught by characteristic.""" 

70 RAREF_CHAR_COLLISION = "rarefaction_characteristic_collision" 

71 """Rarefaction boundary intersects with characteristic.""" 

72 SHOCK_RAREF_COLLISION = "shock_rarefaction_collision" 

73 """Shock intersects with rarefaction boundary.""" 

74 RAREF_RAREF_COLLISION = "rarefaction_rarefaction_collision" 

75 """Rarefaction boundary intersects with another rarefaction boundary.""" 

76 DSW_FAN_EXHAUSTED = "decaying_shock_fan_exhausted" 

77 """A decaying shock's fan is exhausted (c_decay reached c_fan_tail).""" 

78 OUTLET_CROSSING = "outlet_crossing" 

79 """Wave crosses outlet boundary.""" 

80 

81 

82@dataclass 

83class Event: 

84 """A single event in the simulation, ordered by cumulative flow θ. 

85 

86 The solver's priority queue orders ``(theta, counter, ...)`` tuples, not 

87 ``Event`` objects, so this dataclass intentionally defines no ordering. 

88 

89 Parameters 

90 ---------- 

91 theta : float 

92 Cumulative flow at which the event occurs [m³]. 

93 event_type : EventType 

94 Type of event. 

95 waves_involved : list 

96 List of wave objects involved in this event. 

97 location : float 

98 Volumetric position at which the event occurs [m³]. 

99 boundary_type : str or None 

100 Which rarefaction boundary collided: ``'head'`` or ``'tail'``. 

101 Set for rarefaction collision events. 

102 """ 

103 

104 theta: float 

105 event_type: EventType 

106 waves_involved: list # List[Wave] - can't type hint due to circular import 

107 location: float 

108 boundary_type: str | None = None 

109 

110 def __repr__(self): # noqa: D105 

111 return ( 

112 f"Event(θ={self.theta:.3f}, type={self.event_type.value}, " 

113 f"location={self.location:.3f}, n_waves={len(self.waves_involved)})" 

114 ) 

115 

116 

117def find_characteristic_intersection(char1, char2, theta_current: float) -> tuple[float, float] | None: 

118 """Find exact analytical intersection of two characteristics in (V, θ). 

119 

120 Returns (θ_intersect, V_intersect) if the intersection lies in the future 

121 (θ > θ_current) and both characteristics are active there; otherwise None. 

122 """ 

123 s1 = characteristic_speed(char1.concentration, char1.sorption) 

124 s2 = characteristic_speed(char2.concentration, char2.sorption) 

125 

126 if abs(s1 - s2) < EPSILON_SPEED: 

127 return None 

128 

129 theta_both_active = max(char1.theta_start, char2.theta_start, theta_current) 

130 

131 v1 = characteristic_position( 

132 char1.concentration, char1.sorption, char1.theta_start, char1.v_start, theta_both_active 

133 ) 

134 v2 = characteristic_position( 

135 char2.concentration, char2.sorption, char2.theta_start, char2.v_start, theta_both_active 

136 ) 

137 

138 if v1 is None or v2 is None: 

139 return None 

140 

141 # v1 + s1*dθ = v2 + s2*dθ 

142 dtheta = (v2 - v1) / (s1 - s2) 

143 

144 if dtheta <= 0: 

145 return None 

146 

147 theta_intersect = theta_both_active + dtheta 

148 v_intersect = v1 + s1 * dtheta 

149 

150 return (theta_intersect, v_intersect) 

151 

152 

153def find_shock_shock_intersection(shock1, shock2, theta_current: float) -> tuple[float, float] | None: 

154 """Find exact analytical intersection of two shocks in (V, θ).""" 

155 s1 = shock1.speed 

156 s2 = shock2.speed 

157 

158 if abs(s1 - s2) < EPSILON_SPEED: 

159 return None 

160 

161 theta_both_active = max(shock1.theta_start, shock2.theta_start, theta_current) 

162 

163 v1_ref = shock1.v_start + shock1.speed * (theta_both_active - shock1.theta_start) 

164 v2_ref = shock2.v_start + shock2.speed * (theta_both_active - shock2.theta_start) 

165 

166 dtheta = (v2_ref - v1_ref) / (s1 - s2) 

167 

168 if dtheta <= 0: 

169 return None 

170 

171 theta_intersect = theta_both_active + dtheta 

172 v_intersect = v1_ref + s1 * dtheta 

173 

174 return (theta_intersect, v_intersect) 

175 

176 

177def find_shock_characteristic_intersection(shock, char, theta_current: float) -> tuple[float, float] | None: 

178 """Find exact analytical intersection of a shock and a characteristic in (V, θ).""" 

179 s_shock = shock.speed 

180 s_char = characteristic_speed(char.concentration, char.sorption) 

181 

182 if abs(s_shock - s_char) < EPSILON_SPEED: 

183 return None 

184 

185 theta_both_active = max(shock.theta_start, char.theta_start, theta_current) 

186 

187 v_shock = shock.v_start + shock.speed * (theta_both_active - shock.theta_start) 

188 

189 v_char = characteristic_position( 

190 char.concentration, char.sorption, char.theta_start, char.v_start, theta_both_active 

191 ) 

192 

193 if v_char is None: 

194 return None 

195 

196 dtheta = (v_char - v_shock) / (s_shock - s_char) 

197 

198 if dtheta <= 0: 

199 return None 

200 

201 theta_intersect = theta_both_active + dtheta 

202 v_intersect = v_shock + s_shock * dtheta 

203 

204 return (theta_intersect, v_intersect) 

205 

206 

207def find_rarefaction_boundary_intersections(raref, other_wave, theta_current: float) -> list[tuple[float, float, str]]: 

208 """Intersections of a rarefaction's head/tail with another wave. 

209 

210 Both rarefaction boundaries propagate at characteristic speeds (head at 

211 ``1/R(c_head)``, tail at ``1/R(c_tail)``), so we synthesize temporary 

212 ``CharacteristicWave`` instances and reuse the analytical helpers. 

213 

214 Returns 

215 ------- 

216 list of tuple 

217 ``(θ_intersect, V_intersect, boundary_type)`` for each intersection, 

218 where boundary_type is ``'head'`` or ``'tail'``. 

219 """ 

220 intersections = [] 

221 

222 head_char = CharacteristicWave( 

223 theta_start=raref.theta_start, 

224 v_start=raref.v_start, 

225 concentration=raref.c_head, 

226 sorption=raref.sorption, 

227 is_active=raref.is_active, 

228 ) 

229 

230 tail_char = CharacteristicWave( 

231 theta_start=raref.theta_start, 

232 v_start=raref.v_start, 

233 concentration=raref.c_tail, 

234 sorption=raref.sorption, 

235 is_active=raref.is_active, 

236 ) 

237 

238 if isinstance(other_wave, CharacteristicWave): 

239 result = find_characteristic_intersection(head_char, other_wave, theta_current) 

240 if result: 

241 intersections.append((result[0], result[1], "head")) 

242 

243 result = find_characteristic_intersection(tail_char, other_wave, theta_current) 

244 if result: 

245 intersections.append((result[0], result[1], "tail")) 

246 

247 elif isinstance(other_wave, ShockWave): 

248 result = find_shock_characteristic_intersection(other_wave, head_char, theta_current) 

249 if result: 

250 intersections.append((result[0], result[1], "head")) 

251 

252 result = find_shock_characteristic_intersection(other_wave, tail_char, theta_current) 

253 if result: 

254 intersections.append((result[0], result[1], "tail")) 

255 

256 elif isinstance(other_wave, RarefactionWave): 

257 other_head_char = CharacteristicWave( 

258 theta_start=other_wave.theta_start, 

259 v_start=other_wave.v_start, 

260 concentration=other_wave.c_head, 

261 sorption=other_wave.sorption, 

262 is_active=other_wave.is_active, 

263 ) 

264 

265 other_tail_char = CharacteristicWave( 

266 theta_start=other_wave.theta_start, 

267 v_start=other_wave.v_start, 

268 concentration=other_wave.c_tail, 

269 sorption=other_wave.sorption, 

270 is_active=other_wave.is_active, 

271 ) 

272 

273 result = find_characteristic_intersection(head_char, other_head_char, theta_current) 

274 if result: 

275 intersections.append((result[0], result[1], "head")) 

276 

277 result = find_characteristic_intersection(head_char, other_tail_char, theta_current) 

278 if result: 

279 intersections.append((result[0], result[1], "head")) 

280 

281 result = find_characteristic_intersection(tail_char, other_head_char, theta_current) 

282 if result: 

283 intersections.append((result[0], result[1], "tail")) 

284 

285 result = find_characteristic_intersection(tail_char, other_tail_char, theta_current) 

286 if result: 

287 intersections.append((result[0], result[1], "tail")) 

288 

289 return intersections 

290 

291 

292def find_outlet_crossing(wave, v_outlet: float, theta_current: float) -> float | None: 

293 """Find the cumulative flow θ at which the wave crosses ``v_outlet``. 

294 

295 Handles ``CharacteristicWave``, ``ShockWave``, and ``DecayingShockWave``. 

296 Rarefaction outlet crossings are handled by the callers directly (the 

297 solver and ``output.py`` split them into head/tail boundary crossings), so 

298 a ``RarefactionWave`` never reaches this function and returns ``None``. 

299 

300 Assumes positive flow (waves always move toward larger V). Returns None if 

301 the wave has already passed the outlet, is not active, or moves backward. 

302 The "already past" check uses a relative tolerance so that a wave whose 

303 crossing event has just been processed (and is at v_outlet ± a few ULPs) 

304 does not re-emit a duplicate crossing one ULP later. 

305 """ 

306 if not wave.is_active: 

307 return None 

308 

309 # Suppress re-emission when v_current is within FP of v_outlet: the 

310 # crossing was already recorded on the prior iteration. 

311 tol = 1e-12 * max(abs(v_outlet), abs(wave.v_start), 1.0) 

312 

313 if isinstance(wave, CharacteristicWave): 

314 theta_eval = max(theta_current, wave.theta_start) 

315 v_current = characteristic_position( 

316 wave.concentration, wave.sorption, wave.theta_start, wave.v_start, theta_eval 

317 ) 

318 

319 if v_current is None or v_current >= v_outlet - tol: 

320 return None 

321 

322 speed = characteristic_speed(wave.concentration, wave.sorption) 

323 

324 # A c_min-floored (pinned) characteristic — R(c_min) inflated for n>1, 

325 # c→0 — advects too slowly to cross at any physical θ; suppress the 

326 # artifact crossing rather than scheduling it at θ~1e8. 

327 if speed <= 0 or is_outlet_crossing_pinned(wave.concentration, wave.sorption): 

328 return None 

329 

330 dtheta = (v_outlet - v_current) / speed 

331 return theta_eval + dtheta 

332 

333 if isinstance(wave, ShockWave): 

334 theta_eval = max(theta_current, wave.theta_start) 

335 v_current = wave.v_start + wave.speed * (theta_eval - wave.theta_start) 

336 

337 if v_current >= v_outlet - tol: 

338 return None 

339 

340 if wave.speed <= 0: 

341 return None 

342 

343 dtheta = (v_outlet - v_current) / wave.speed 

344 return theta_eval + dtheta 

345 

346 if isinstance(wave, DecayingShockWave): 

347 # Closed-form inverse V_s(theta) = v_outlet. 

348 theta_cross = wave.outlet_crossing_theta(v_outlet) 

349 if theta_cross is None: 

350 return None 

351 # Suppress re-emission within FP of the prior crossing (same convention 

352 # as the linear-shock branch above). 

353 if theta_cross <= theta_current + 1e-15 * max(abs(theta_current), 1.0): 

354 return None 

355 return theta_cross 

356 

357 return None