Coverage for src / gwtransport / fronttracking / plot.py: 75%

1""" 

2Visualization functions for front tracking. 

3 

4This module provides plotting utilities for visualizing front-tracking simulations: 

5- V-t diagrams showing wave propagation in space-time 

6- Breakthrough curves showing concentration at outlet over time 

7 

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""" 

11 

12import matplotlib.pyplot as plt 

13import numpy as np 

14import pandas as pd 

15 

16from gwtransport.fronttracking.output import concentration_at_point, identify_outlet_segments 

17from gwtransport.fronttracking.solver import FrontTrackerState 

18from gwtransport.fronttracking.waves import CharacteristicWave, RarefactionWave, ShockWave 

19 

20 

21def plot_vt_diagram( 

22 state: FrontTrackerState, 

23 ax=None, 

24 t_max: float | None = None, 

25 figsize: tuple[float, float] = (14, 10), 

26 *, 

27 show_inactive: bool = False, 

28 show_events: bool = False, 

29): 

30 """ 

31 Create V-t diagram showing all waves in space-time. 

32 

33 Plots characteristics (blue lines), shocks (red lines), and rarefactions 

34 (green fans) in the (time, position) plane. This visualization shows how 

35 waves propagate and interact throughout the simulation. 

36 

37 Parameters 

38 ---------- 

39 state : FrontTrackerState 

40 Complete simulation state containing all waves. 

41 ax : matplotlib.axes.Axes, optional 

42 Existing axes to plot into. If None, a new figure and axes are created 

43 using ``figsize``. 

44 t_max : float, optional 

45 Maximum time to plot [days]. If None, uses final simulation time * 1.2. 

46 figsize : tuple of float, optional 

47 Figure size in inches (width, height). Default (14, 10). 

48 show_inactive : bool, optional 

49 Whether to show inactive waves (deactivated by interactions). 

50 Default False. 

51 show_events : bool, optional 

52 Whether to show wave interaction events as markers. 

53 Default False. 

54 

55 Returns 

56 ------- 

57 fig : matplotlib.figure.Figure 

58 Figure object containing the V-t diagram. 

59 

60 Notes 

61 ----- 

62 - Characteristics appear as blue lines (constant velocity). 

63 - Shocks appear as thick red lines (jump discontinuities). 

64 - Rarefactions appear as green fans (smooth transition regions). 

65 - Outlet position is shown as a horizontal dashed line. 

66 - Only waves within domain [0, v_outlet] are plotted. 

67 

68 Examples 

69 -------- 

70 :: 

71 

72 from gwtransport.fronttracking.solver import FrontTracker 

73 

74 tracker = FrontTracker(cin, flow, tedges, aquifer_pore_volume, sorption) 

75 tracker.run() 

76 fig = plot_vt_diagram(tracker.state) 

77 fig.savefig("vt_diagram.png") 

78 """ 

79 if t_max is None: 

80 # Default to input data time range instead of simulation end time 

81 # Convert tedges[-1] from Timestamp to days from tedges[0] 

82 t_max = (state.tedges[-1] - state.tedges[0]) / pd.Timedelta(days=1) 

83 

84 if ax is None: 

85 _, ax = plt.subplots(figsize=figsize) 

86 

87 # Plot characteristics (blue lines) 

88 for wave in state.waves: 

89 if isinstance(wave, CharacteristicWave): 

90 if not wave.is_active and not show_inactive: 

91 continue 

92 

93 t_plot = np.linspace(wave.t_start, t_max, 100) 

94 v_plot = [] 

95 t_plot_used = [] 

96 

97 for t in t_plot: 

98 # Compute position manually for inactive waves when show_inactive=True 

99 v = wave.position_at_time(t) if wave.is_active else wave.v_start + wave.velocity() * (t - wave.t_start) 

100 

101 if v is not None and 0 <= v <= state.v_outlet: 

102 v_plot.append(v) 

103 t_plot_used.append(t) 

104 elif v is not None and v > state.v_outlet: 

105 # Wave crossed outlet - add exact intersection point 

106 vel = wave.velocity() 

107 if vel > 0: 

108 t_outlet = wave.t_start + (state.v_outlet - wave.v_start) / vel 

109 if wave.t_start <= t_outlet <= t_max: 

110 v_plot.append(state.v_outlet) 

111 t_plot_used.append(t_outlet) 

112 break 

113 else: 

114 break 

115 

116 if len(v_plot) > 0: 

117 alpha = 0.3 if not wave.is_active else 0.7 

118 ax.plot( 

119 t_plot_used, 

120 v_plot, 

121 "b-", 

122 linewidth=0.5, 

123 alpha=alpha, 

124 label="Characteristic" if not hasattr(ax, "gw_char_labeled") else "", 

125 ) 

126 ax.gw_char_labeled = True # type: ignore[attr-defined] 

127 

128 # Plot shocks (red lines) 

129 for wave in state.waves: 

130 if isinstance(wave, ShockWave): 

131 if not wave.is_active and not show_inactive: 

132 continue 

133 

134 t_plot = np.linspace(wave.t_start, t_max, 100) 

135 v_plot = [] 

136 t_plot_used = [] 

137 

138 for t in t_plot: 

139 # Compute position manually for inactive waves when show_inactive=True 

140 v = wave.position_at_time(t) if wave.is_active else wave.v_start + wave.velocity * (t - wave.t_start) 

141 

142 if v is not None and 0 <= v <= state.v_outlet: 

143 v_plot.append(v) 

144 t_plot_used.append(t) 

145 elif v is not None and v > state.v_outlet: 

146 # Wave crossed outlet - add exact intersection point 

147 vel = wave.velocity 

148 if vel is not None and vel > 0: 

149 t_outlet = wave.t_start + (state.v_outlet - wave.v_start) / vel 

150 if wave.t_start <= t_outlet <= t_max: 

151 v_plot.append(state.v_outlet) 

152 t_plot_used.append(t_outlet) 

153 break 

154 else: 

155 break 

156 

157 if len(v_plot) > 0: 

158 alpha = 0.5 if not wave.is_active else 1.0 

159 ax.plot( 

160 t_plot_used, 

161 v_plot, 

162 "r-", 

163 linewidth=2, 

164 alpha=alpha, 

165 label="Shock" if not hasattr(ax, "gw_shock_labeled") else "", 

166 ) 

167 ax.gw_shock_labeled = True # type: ignore[attr-defined] 

168 

169 # Plot rarefactions (green fans) 

170 for wave in state.waves: 

171 if isinstance(wave, RarefactionWave): 

172 if not wave.is_active and not show_inactive: 

173 continue 

174 

175 t_plot = np.linspace(wave.t_start, t_max, 100) 

176 v_head_plot = [] 

177 v_tail_plot = [] 

178 t_plot_used = [] 

179 head_crossed = False 

180 tail_crossed = False 

181 

182 for t in t_plot: 

183 # Compute positions manually for inactive waves when show_inactive=True 

184 if wave.is_active: 

185 v_head = wave.head_position_at_time(t) 

186 v_tail = wave.tail_position_at_time(t) 

187 else: 

188 # Manually compute positions for visualization 

189 v_head = wave.v_start + wave.head_velocity() * (t - wave.t_start) 

190 v_tail = wave.v_start + wave.tail_velocity() * (t - wave.t_start) 

191 

192 # Track time points 

193 t_plot_used.append(t) 

194 

195 # Handle head 

196 if v_head is not None and 0 <= v_head <= state.v_outlet: 

197 v_head_plot.append(v_head) 

198 elif v_head is not None and v_head > state.v_outlet and not head_crossed: 

199 # Add exact outlet intersection for head 

200 head_vel = wave.head_velocity() 

201 if head_vel > 0: 

202 t_outlet_head = wave.t_start + (state.v_outlet - wave.v_start) / head_vel 

203 if wave.t_start <= t_outlet_head <= t_max: 

204 # Insert the exact crossing point 

205 v_head_plot.append(state.v_outlet) 

206 head_crossed = True 

207 v_head_plot.append(None) 

208 else: 

209 v_head_plot.append(None) 

210 

211 # Handle tail 

212 if v_tail is not None and 0 <= v_tail <= state.v_outlet: 

213 v_tail_plot.append(v_tail) 

214 elif v_tail is not None and v_tail > state.v_outlet and not tail_crossed: 

215 # Add exact outlet intersection for tail 

216 tail_vel = wave.tail_velocity() 

217 if tail_vel > 0: 

218 t_outlet_tail = wave.t_start + (state.v_outlet - wave.v_start) / tail_vel 

219 if wave.t_start <= t_outlet_tail <= t_max: 

220 # Insert the exact crossing point 

221 v_tail_plot.append(state.v_outlet) 

222 tail_crossed = True 

223 v_tail_plot.append(None) 

224 else: 

225 v_tail_plot.append(None) 

226 

227 # Plot head and tail boundaries 

228 alpha = 0.5 if not wave.is_active else 0.8 

229 label = "Rarefaction" if not hasattr(ax, "gw_raref_labeled") else "" 

230 

231 # Plot head (faster boundary) 

232 valid_head = [(t, v) for t, v in zip(t_plot_used, v_head_plot, strict=False) if v is not None] 

233 if valid_head: 

234 t_h, v_h = zip(*valid_head, strict=False) 

235 ax.plot(t_h, v_h, "g-", linewidth=1.5, alpha=alpha, label=label) 

236 ax.gw_raref_labeled = True # type: ignore[attr-defined] 

237 

238 # Plot tail (slower boundary) 

239 valid_tail = [(t, v) for t, v in zip(t_plot_used, v_tail_plot, strict=False) if v is not None] 

240 if valid_tail: 

241 t_t, v_t = zip(*valid_tail, strict=False) 

242 ax.plot(t_t, v_t, "g--", linewidth=1.5, alpha=alpha) 

243 

244 # Fill between head and tail 

245 if valid_head and valid_tail and len(valid_head) == len(valid_tail): 

246 ax.fill_between( 

247 t_h, 

248 v_h, 

249 v_t, 

250 color="green", 

251 alpha=0.1 if not wave.is_active else 0.2, 

252 ) 

253 

254 # Plot outlet position 

255 ax.axhline( 

256 state.v_outlet, 

257 color="k", 

258 linestyle="--", 

259 linewidth=1, 

260 alpha=0.5, 

261 label=f"Outlet (V={state.v_outlet:.1f} m³)", 

262 ) 

263 

264 # Plot inlet position 

265 ax.axhline( 

266 0.0, 

267 color="k", 

268 linestyle=":", 

269 linewidth=1, 

270 alpha=0.5, 

271 label="Inlet (V=0)", 

272 ) 

273 

274 # Plot wave interaction events as markers 

275 if show_events and hasattr(state, "events") and state.events: 

276 for event in state.events: 

277 if "time" in event and "position" in event: 

278 t_event = event["time"] 

279 v_event = event["position"] 

280 if 0 <= t_event <= t_max and 0 <= v_event <= state.v_outlet: 

281 # Determine marker style based on event type 

282 event_type = event.get("type", "unknown") 

283 if "shock" in event_type.lower() or "collision" in event_type.lower(): 

284 marker = "X" 

285 color = "red" 

286 size = 100 

287 elif "rarefaction" in event_type.lower(): 

288 marker = "o" 

289 color = "green" 

290 size = 80 

291 elif "outlet" in event_type.lower(): 

292 marker = "s" 

293 color = "black" 

294 size = 80 

295 else: 

296 marker = "D" 

297 color = "gray" 

298 size = 60 

299 

300 ax.scatter( 

301 t_event, 

302 v_event, 

303 marker=marker, 

304 s=size, 

305 color=color, 

306 edgecolors="black", 

307 linewidths=1.5, 

308 alpha=0.8, 

309 zorder=10, 

310 label="Event" if not hasattr(ax, "gw_event_labeled") else "", 

311 ) 

312 ax.gw_event_labeled = True # type: ignore[attr-defined] 

313 

314 ax.set_xlabel("Time [days]", fontsize=12) 

315 ax.set_ylabel("Position (Pore Volume) [m³]", fontsize=12) 

316 ax.set_title("V-t Diagram: Front Tracking Simulation", fontsize=14, fontweight="bold") 

317 ax.grid(True, alpha=0.3) 

318 ax.legend(loc="best") 

319 ax.set_xlim(0, t_max) 

320 ax.set_ylim(-state.v_outlet * 0.05, state.v_outlet * 1.05) 

321 

322 return ax 

323 

324 

325def plot_breakthrough_curve( 

326 state: FrontTrackerState, 

327 ax=None, 

328 t_max: float | None = None, 

329 n_rarefaction_points: int = 50, 

330 figsize: tuple[float, float] = (12, 6), 

331 t_first_arrival: float | None = None, 

332): 

333 """ 

334 Plot exact analytical concentration breakthrough curve at outlet. 

335 

336 Uses wave segment information to plot the exact analytical solution 

337 without discretization. Constant concentration regions are plotted 

338 as horizontal lines, and rarefaction regions are plotted using their 

339 exact self-similar solutions. 

340 

341 Parameters 

342 ---------- 

343 state : FrontTrackerState 

344 Complete simulation state containing all waves. 

345 ax : matplotlib.axes.Axes, optional 

346 Existing axes to plot into. If None, a new figure and axes are created 

347 using ``figsize``. 

348 t_max : float, optional 

349 Maximum time to plot [days]. If None, uses final simulation time * 1.1. 

350 n_rarefaction_points : int, optional 

351 Number of points to use for plotting rarefaction segments (analytical 

352 curves). Default 50 per rarefaction segment. 

353 figsize : tuple of float, optional 

354 Figure size in inches (width, height). Default (12, 6). 

355 t_first_arrival : float, optional 

356 First arrival time for marking spin-up period [days]. If None, spin-up 

357 period is not plotted. 

358 

359 Returns 

360 ------- 

361 fig : matplotlib.figure.Figure 

362 Figure object containing the breakthrough curve. 

363 

364 Notes 

365 ----- 

366 - Uses identify_outlet_segments to get exact analytical segment boundaries 

367 - Constant concentration segments plotted as horizontal lines (no discretization) 

368 - Rarefaction segments plotted using exact self-similar solution 

369 - Shocks appear as instantaneous jumps at exact crossing times 

370 - No bin averaging or discretization artifacts 

371 

372 Examples 

373 -------- 

374 :: 

375 

376 from gwtransport.fronttracking.solver import FrontTracker 

377 

378 tracker = FrontTracker(cin, flow, tedges, aquifer_pore_volume, sorption) 

379 tracker.run() 

380 fig = plot_breakthrough_curve(tracker.state) 

381 fig.savefig("exact_breakthrough.png") 

382 """ 

383 if ax is None: 

384 _, ax = plt.subplots(figsize=figsize) 

385 

386 if t_max is None: 

387 # Default to input data time range instead of simulation end time 

388 # Convert tedges[-1] from Timestamp to days from tedges[0] 

389 t_max = (state.tedges[-1] - state.tedges[0]) / pd.Timedelta(days=1) 

390 

391 # Use exact analytical segments 

392 segments = identify_outlet_segments(0.0, t_max, state.v_outlet, state.waves, state.sorption) 

393 

394 for i, segment in enumerate(segments): 

395 t_start = segment["t_start"] 

396 t_end = segment["t_end"] 

397 

398 if segment["type"] == "constant": 

399 # Constant concentration segment - plot as horizontal line 

400 c_const = segment["concentration"] 

401 ax.plot( 

402 [t_start, t_end], 

403 [c_const, c_const], 

404 "b-", 

405 linewidth=2, 

406 label="Outlet concentration" if i == 0 else "", 

407 ) 

408 elif segment["type"] == "rarefaction": 

409 # Rarefaction segment - plot exact analytical curve 

410 raref = segment["wave"] 

411 t_raref = np.linspace(t_start, t_end, n_rarefaction_points) 

412 c_raref = np.zeros_like(t_raref) 

413 

414 for j, t in enumerate(t_raref): 

415 # Use the rarefaction wave's own concentration_at_point method 

416 c_at_point = raref.concentration_at_point(state.v_outlet, t) 

417 if c_at_point is not None: 

418 c_raref[j] = c_at_point 

419 else: 

420 # Fallback to boundary values if not in fan 

421 c_raref[j] = segment.get("c_start", raref.c_tail) 

422 

423 ax.plot(t_raref, c_raref, "b-", linewidth=2, label="Outlet concentration" if i == 0 else "") 

424 

425 # Mark first arrival time if provided 

426 if t_first_arrival is not None and np.isfinite(t_first_arrival): 

427 ax.axvline( 

428 t_first_arrival, 

429 color="r", 

430 linestyle="--", 

431 linewidth=1.5, 

432 alpha=0.7, 

433 label=f"First arrival (t={t_first_arrival:.2f} days)", 

434 ) 

435 

436 # Shade spin-up region 

437 ax.axvspan( 

438 0, 

439 t_first_arrival, 

440 alpha=0.1, 

441 color="gray", 

442 label="Spin-up period", 

443 ) 

444 

445 ax.set_xlabel("Time [days]", fontsize=12) 

446 ax.set_ylabel("Concentration [mass/volume]", fontsize=12) 

447 ax.set_title("Breakthrough Curve at Outlet (Exact Analytical)", fontsize=14, fontweight="bold") 

448 ax.grid(True, alpha=0.3) 

449 ax.legend(loc="best") 

450 ax.set_xlim(0, t_max) 

451 ax.set_ylim(bottom=0) 

452 

453 return ax 

454 

455 

456def plot_wave_interactions( 

457 state: FrontTrackerState, 

458 figsize: tuple[float, float] = (14, 8), 

459 ax=None, 

460): 

461 """ 

462 Plot event timeline showing wave interactions. 

463 

464 Creates a scatter plot showing when and where different types of wave 

465 interactions occur during the simulation. 

466 

467 Parameters 

468 ---------- 

469 state : FrontTrackerState 

470 Complete simulation state containing all events. 

471 ax : matplotlib.axes.Axes, optional 

472 Existing axes to plot into. If None, a new figure and axes are created 

473 using ``figsize``. 

474 figsize : tuple of float, optional 

475 Figure size in inches (width, height). Default (14, 8). 

476 

477 Returns 

478 ------- 

479 fig : matplotlib.figure.Figure 

480 Figure object containing the event timeline. 

481 

482 Notes 

483 ----- 

484 - Each event type is shown with a different color and marker. 

485 - Outlet crossings are shown separately from internal collisions. 

486 - Event locations are plotted in the (time, position) plane. 

487 

488 Examples 

489 -------- 

490 :: 

491 

492 from gwtransport.fronttracking.solver import FrontTracker 

493 

494 tracker = FrontTracker(cin, flow, tedges, aquifer_pore_volume, sorption) 

495 tracker.run() 

496 fig = plot_wave_interactions(tracker.state) 

497 fig.savefig("wave_interactions.png") 

498 """ 

499 if ax is None: 

500 _, ax = plt.subplots(figsize=figsize) 

501 

502 # Group events by type 

503 event_types = {} 

504 for event_dict in state.events: 

505 event_type = event_dict["type"] 

506 if event_type not in event_types: 

507 event_types[event_type] = {"times": [], "locations": []} 

508 event_types[event_type]["times"].append(event_dict["time"]) 

509 event_types[event_type]["locations"].append(event_dict.get("location", 0.0)) 

510 

511 # Define colors and markers for each event type 

512 event_style = { 

513 "CHAR_CHAR_COLLISION": {"color": "blue", "marker": "o", "label": "Char-Char"}, 

514 "SHOCK_SHOCK_COLLISION": {"color": "red", "marker": "s", "label": "Shock-Shock"}, 

515 "SHOCK_CHAR_COLLISION": {"color": "purple", "marker": "^", "label": "Shock-Char"}, 

516 "RAREF_CHAR_COLLISION": {"color": "green", "marker": "v", "label": "Raref-Char"}, 

517 "SHOCK_RAREF_COLLISION": {"color": "orange", "marker": "d", "label": "Shock-Raref"}, 

518 "RAREF_RAREF_COLLISION": {"color": "cyan", "marker": "p", "label": "Raref-Raref"}, 

519 "OUTLET_CROSSING": {"color": "black", "marker": "x", "label": "Outlet Crossing"}, 

520 "INLET_CHANGE": {"color": "gray", "marker": "+", "label": "Inlet Change"}, 

521 } 

522 

523 # Plot each event type 

524 for event_type, data in event_types.items(): 

525 style = event_style.get(event_type, {"color": "gray", "marker": "o", "label": event_type}) 

526 ax.scatter( 

527 data["times"], 

528 data["locations"], 

529 c=style["color"], 

530 marker=style["marker"], 

531 s=100, 

532 alpha=0.7, 

533 label=f"{style['label']} ({len(data['times'])})", 

534 ) 

535 

536 # Plot outlet line for reference 

537 if state.events: 

538 ax.axhline( 

539 state.v_outlet, 

540 color="k", 

541 linestyle="--", 

542 linewidth=1, 

543 alpha=0.3, 

544 label=f"Outlet (V={state.v_outlet:.1f} m³)", 

545 ) 

546 

547 ax.set_xlabel("Time [days]", fontsize=12) 

548 ax.set_ylabel("Position (Pore Volume) [m³]", fontsize=12) 

549 ax.set_title("Wave Interaction Events", fontsize=14, fontweight="bold") 

550 ax.grid(True, alpha=0.3) 

551 ax.legend(loc="best", ncol=2) 

552 

553 if state.events: 

554 ax.set_xlim(left=0) 

555 ax.set_ylim(-state.v_outlet * 0.05, state.v_outlet * 1.05) 

556 

557 return ax 

558 

559 

560def plot_inlet_concentration( 

561 tedges, 

562 cin, 

563 ax=None, 

564 *, 

565 t_first_arrival=None, 

566 event_markers=None, 

567 color="blue", 

568 t_max=None, 

569 xlabel="Time [days]", 

570 ylabel="Concentration", 

571 title="Inlet Concentration", 

572 figsize=(8, 5), 

573 **step_kwargs, 

574): 

575 """ 

576 Plot inlet concentration as step function with optional markers. 

577 

578 Parameters 

579 ---------- 

580 tedges : pandas.DatetimeIndex 

581 Time bin edges for inlet concentration. 

582 Length = len(cin) + 1. 

583 cin : array-like 

584 Inlet concentration values. 

585 Length = len(tedges) - 1. 

586 ax : matplotlib.axes.Axes, optional 

587 Existing axes to plot into. If None, creates new figure. 

588 t_first_arrival : float, optional 

589 First arrival time to mark with vertical line [days]. 

590 event_markers : list of dict, optional 

591 Event markers to add. Each dict should have keys: 'time', 'label', 'color'. 

592 color : str, optional 

593 Color for inlet concentration line. Default 'blue'. 

594 t_max : float, optional 

595 Maximum time for x-axis [days]. If None, uses full range. 

596 xlabel : str, optional 

597 Label for x-axis. Default 'Time [days]'. 

598 ylabel : str, optional 

599 Label for y-axis. Default 'Concentration'. 

600 title : str, optional 

601 Plot title. Default 'Inlet Concentration'. 

602 figsize : tuple, optional 

603 Figure size if creating new figure. Default (8, 5). 

604 **step_kwargs 

605 Additional arguments passed to ax.plot(). 

606 

607 Returns 

608 ------- 

609 fig : matplotlib.figure.Figure 

610 Figure object. 

611 ax : matplotlib.axes.Axes 

612 Axes object. 

613 """ 

614 if ax is None: 

615 _, ax = plt.subplots(figsize=figsize) 

616 

617 # Convert tedges to days from start 

618 tedges_array = tedges.to_numpy() if hasattr(tedges, "to_numpy") else np.array(tedges) 

619 

620 t_days = (tedges_array - tedges_array[0]) / pd.Timedelta(days=1) 

621 

622 # Plot inlet concentration using repeat pattern for step function 

623 x_plot, y_plot = np.repeat(t_days, 2)[1:-1], np.repeat(cin, 2) 

624 ax.plot(x_plot, y_plot, linewidth=2, color=color, label="Inlet", **step_kwargs) 

625 

626 # Add first arrival marker if provided 

627 if t_first_arrival is not None and np.isfinite(t_first_arrival): 

628 ax.axvline( 

629 t_first_arrival, 

630 color="green", 

631 linestyle="--", 

632 linewidth=1.5, 

633 alpha=0.7, 

634 label=f"First arrival ({t_first_arrival:.1f} days)", 

635 ) 

636 

637 # Add event markers if provided 

638 if event_markers is not None: 

639 for marker in event_markers: 

640 t = marker.get("time") 

641 label = marker.get("label", "") 

642 marker_color = marker.get("color", "gray") 

643 linestyle = marker.get("linestyle", "--") 

644 

645 if t is not None: 

646 ax.axvline( 

647 t, 

648 color=marker_color, 

649 linestyle=linestyle, 

650 linewidth=1.5, 

651 alpha=0.7, 

652 label=label, 

653 ) 

654 

655 ax.set_xlabel(xlabel, fontsize=10) 

656 ax.set_ylabel(ylabel, fontsize=10) 

657 ax.set_title(title, fontsize=12, fontweight="bold") 

658 ax.grid(True, alpha=0.3) 

659 ax.legend(fontsize=8) 

660 

661 if t_max is not None: 

662 ax.set_xlim(0, t_max) 

663 else: 

664 ax.set_xlim(0, t_days[-1]) 

665 

666 return ax 

667 

668 

669def plot_front_tracking_summary( 

670 structure, 

671 tedges, 

672 cin, 

673 cout_tedges, 

674 cout, 

675 *, 

676 figsize=(16, 10), 

677 show_exact=True, 

678 show_bin_averaged=True, 

679 show_events=True, 

680 show_inactive=False, 

681 t_max=None, 

682 title=None, 

683 inlet_color="blue", 

684 outlet_exact_color="blue", 

685 outlet_binned_color="red", 

686 first_arrival_color="green", 

687): 

688 """ 

689 Create comprehensive 3-panel summary figure for front tracking simulation. 

690 

691 Creates a multi-panel visualization with: 

692 - Top-left: V-t diagram showing wave propagation 

693 - Top-right: Inlet concentration time series 

694 - Bottom: Outlet concentration (exact and/or bin-averaged) 

695 

696 Parameters 

697 ---------- 

698 structure : dict 

699 Structure returned from infiltration_to_extraction_front_tracking_detailed. 

700 Must contain keys: 'tracker_state', 't_first_arrival'. 

701 tedges : pandas.DatetimeIndex 

702 Time bin edges for inlet concentration. 

703 Length = len(cin) + 1. 

704 cin : array-like 

705 Inlet concentration values. 

706 Length = len(tedges) - 1. 

707 cout_tedges : pandas.DatetimeIndex 

708 Output time bin edges for bin-averaged concentration. 

709 Length = len(cout) + 1. 

710 cout : array-like 

711 Bin-averaged output concentration values. 

712 Length = len(cout_tedges) - 1. 

713 figsize : tuple, optional 

714 Figure size (width, height). Default (16, 10). 

715 show_exact : bool, optional 

716 Whether to show exact analytical breakthrough curve. Default True. 

717 show_bin_averaged : bool, optional 

718 Whether to show bin-averaged concentration. Default True. 

719 show_events : bool, optional 

720 Whether to show wave interaction events on V-t diagram. Default True. 

721 show_inactive : bool, optional 

722 Whether to show inactive waves on V-t diagram. Default False. 

723 t_max : float, optional 

724 Maximum time for plots [days]. If None, uses input data range. 

725 title : str, optional 

726 Overall figure title. If None, uses generic title. 

727 inlet_color : str, optional 

728 Color for inlet concentration. Default 'blue'. 

729 outlet_exact_color : str, optional 

730 Color for exact outlet curve. Default 'blue'. 

731 outlet_binned_color : str, optional 

732 Color for bin-averaged outlet. Default 'red'. 

733 first_arrival_color : str, optional 

734 Color for first arrival marker. Default 'green'. 

735 

736 Returns 

737 ------- 

738 fig : matplotlib.figure.Figure 

739 Figure object. 

740 axes : dict 

741 Dictionary with keys 'vt', 'inlet', 'outlet' containing axes objects. 

742 """ 

743 # Create figure with 3-panel layout 

744 fig = plt.figure(figsize=figsize) 

745 gs = fig.add_gridspec(2, 2, hspace=0.3, wspace=0.3) 

746 

747 axes = {} 

748 

749 # Top left: V-t diagram 

750 ax_vt = fig.add_subplot(gs[0, 0]) 

751 plot_vt_diagram( 

752 structure["tracker_state"], 

753 ax=ax_vt, 

754 show_inactive=show_inactive, 

755 show_events=show_events, 

756 t_max=t_max, 

757 ) 

758 ax_vt.set_title("V-t Diagram", fontsize=12, fontweight="bold") 

759 axes["vt"] = ax_vt 

760 

761 # Top right: Inlet concentration 

762 ax_inlet = fig.add_subplot(gs[0, 1]) 

763 plot_inlet_concentration( 

764 tedges, 

765 cin, 

766 ax=ax_inlet, 

767 t_first_arrival=structure["t_first_arrival"], 

768 color=inlet_color, 

769 t_max=t_max, 

770 ) 

771 axes["inlet"] = ax_inlet 

772 

773 # Bottom: Outlet concentration (exact and bin-averaged) 

774 ax_outlet = fig.add_subplot(gs[1, :]) 

775 

776 # Compute time range 

777 if t_max is None: 

778 t_max = (tedges.to_numpy()[-1] - tedges.to_numpy()[0]) / pd.Timedelta(days=1) 

779 

780 # Exact breakthrough curve 

781 if show_exact: 

782 t_exact = np.linspace(0, t_max, 1000) 

783 c_exact = [ 

784 concentration_at_point( 

785 structure["tracker_state"].v_outlet, 

786 t, 

787 structure["tracker_state"].waves, 

788 structure["tracker_state"].sorption, 

789 ) 

790 for t in t_exact 

791 ] 

792 ax_outlet.plot( 

793 t_exact, 

794 c_exact, 

795 color=outlet_exact_color, 

796 linewidth=2.5, 

797 label="Exact outlet concentration", 

798 zorder=3, 

799 ) 

800 

801 # Bin-averaged outlet 

802 if show_bin_averaged: 

803 t_edges_days = ((cout_tedges - cout_tedges[0]) / pd.Timedelta(days=1)).values 

804 xstep_cout, ystep_cout = np.repeat(t_edges_days, 2)[1:-1], np.repeat(cout, 2) 

805 ax_outlet.plot( 

806 xstep_cout, 

807 ystep_cout, 

808 color=outlet_binned_color, 

809 linestyle="--", 

810 linewidth=1.5, 

811 alpha=0.7, 

812 label="Bin-averaged outlet", 

813 zorder=2, 

814 ) 

815 

816 # First arrival marker 

817 t_first = structure["t_first_arrival"] 

818 if np.isfinite(t_first): 

819 ax_outlet.axvline( 

820 t_first, 

821 color=first_arrival_color, 

822 linestyle="--", 

823 linewidth=1.5, 

824 alpha=0.7, 

825 label=f"First arrival ({t_first:.1f} days)", 

826 zorder=1, 

827 ) 

828 

829 ax_outlet.set_xlabel("Time [days]", fontsize=11) 

830 ax_outlet.set_ylabel("Concentration", fontsize=11) 

831 ax_outlet.set_title("Outlet Concentration: Exact vs Bin-Averaged", fontsize=12, fontweight="bold") 

832 ax_outlet.grid(True, alpha=0.3) 

833 ax_outlet.legend(fontsize=9) 

834 ax_outlet.set_xlim(0, t_max) 

835 axes["outlet"] = ax_outlet 

836 

837 # Overall title 

838 if title is not None: 

839 plt.suptitle(title, fontsize=14, fontweight="bold", y=0.995) 

840 

841 return fig 

842 

843 

844def plot_sorption_comparison( 

845 pulse_favorable_structure, 

846 pulse_unfavorable_structure, 

847 pulse_tedges, 

848 pulse_cin, 

849 dip_favorable_structure, 

850 dip_unfavorable_structure, 

851 dip_tedges, 

852 dip_cin, 

853 *, 

854 figsize=(16, 12), 

855 t_max_pulse=None, 

856 t_max_dip=None, 

857): 

858 """ 

859 Compare how each inlet produces different outputs with n>1 vs n<1 sorption. 

860 

861 Creates a 2x3 grid: 

862 - Row 1: Pulse inlet and its outputs with n>1 and n<1 sorption 

863 - Row 2: Dip inlet and its outputs with n>1 and n<1 sorption 

864 

865 This demonstrates how the SAME inlet timeseries produces DIFFERENT breakthrough 

866 curves depending on the sorption isotherm. 

867 

868 Parameters 

869 ---------- 

870 pulse_favorable_structure : dict 

871 Structure from pulse inlet with n>1 (higher C travels faster). 

872 pulse_unfavorable_structure : dict 

873 Structure from pulse inlet with n<1 (lower C travels faster). 

874 pulse_tedges : pandas.DatetimeIndex