Coverage for src / gwtransport / fronttracking / plot.py: 75%
337 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-27 06:32 +0000
1"""
2Visualization functions for front tracking.
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
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 matplotlib.pyplot as plt
13import numpy as np
14import pandas as pd
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
19from gwtransport.utils import step_plot_coords
22def plot_vt_diagram(
23 state: FrontTrackerState,
24 ax=None,
25 t_max: float | None = None,
26 figsize: tuple[float, float] = (14, 10),
27 *,
28 show_inactive: bool = False,
29 show_events: bool = False,
30):
31 """
32 Create V-t diagram showing all waves in space-time.
34 Plots characteristics (blue lines), shocks (red lines), and rarefactions
35 (green fans) in the (time, position) plane. This visualization shows how
36 waves propagate and interact throughout the simulation.
38 Parameters
39 ----------
40 state : FrontTrackerState
41 Complete simulation state containing all waves.
42 ax : matplotlib.axes.Axes, optional
43 Existing axes to plot into. If None, a new figure and axes are created
44 using ``figsize``.
45 t_max : float, optional
46 Maximum time to plot [days]. If None, uses final simulation time * 1.2.
47 figsize : tuple of float, optional
48 Figure size in inches (width, height). Default (14, 10).
49 show_inactive : bool, optional
50 Whether to show inactive waves (deactivated by interactions).
51 Default False.
52 show_events : bool, optional
53 Whether to show wave interaction events as markers.
54 Default False.
56 Returns
57 -------
58 fig : matplotlib.figure.Figure
59 Figure object containing the V-t diagram.
61 Notes
62 -----
63 - Characteristics appear as blue lines (constant velocity).
64 - Shocks appear as thick red lines (jump discontinuities).
65 - Rarefactions appear as green fans (smooth transition regions).
66 - Outlet position is shown as a horizontal dashed line.
67 - Only waves within domain [0, v_outlet] are plotted.
69 Examples
70 --------
71 ::
73 from gwtransport.fronttracking.solver import FrontTracker
75 tracker = FrontTracker(cin, flow, tedges, aquifer_pore_volume, sorption)
76 tracker.run()
77 fig = plot_vt_diagram(tracker.state)
78 fig.savefig("vt_diagram.png")
79 """
80 if t_max is None:
81 # Default to input data time range instead of simulation end time
82 # Convert tedges[-1] from Timestamp to days from tedges[0]
83 t_max = float((state.tedges[-1] - state.tedges[0]) / pd.Timedelta(days=1))
85 if ax is None:
86 _, ax = plt.subplots(figsize=figsize)
88 char_labeled = False
89 shock_labeled = False
90 raref_labeled = False
91 event_labeled = False
93 # Plot characteristics (blue lines)
94 for wave in state.waves:
95 if isinstance(wave, CharacteristicWave):
96 if not wave.is_active and not show_inactive:
97 continue
99 t_plot = np.linspace(wave.t_start, t_max, 100)
100 v_plot = []
101 t_plot_used = []
103 for t in t_plot:
104 # Compute position manually for inactive waves when show_inactive=True
105 v = wave.position_at_time(t) if wave.is_active else wave.v_start + wave.velocity() * (t - wave.t_start)
107 if v is not None and 0 <= v <= state.v_outlet:
108 v_plot.append(v)
109 t_plot_used.append(t)
110 elif v is not None and v > state.v_outlet:
111 # Wave crossed outlet - add exact intersection point
112 vel = wave.velocity()
113 if vel > 0:
114 t_outlet = wave.t_start + (state.v_outlet - wave.v_start) / vel
115 if wave.t_start <= t_outlet <= t_max:
116 v_plot.append(state.v_outlet)
117 t_plot_used.append(t_outlet)
118 break
119 else:
120 break
122 if len(v_plot) > 0:
123 alpha = 0.3 if not wave.is_active else 0.7
124 ax.plot(
125 t_plot_used,
126 v_plot,
127 "b-",
128 linewidth=0.5,
129 alpha=alpha,
130 label="Characteristic" if not char_labeled else "",
131 )
132 char_labeled = True
134 # Plot shocks (red lines)
135 for wave in state.waves:
136 if isinstance(wave, ShockWave):
137 if not wave.is_active and not show_inactive:
138 continue
140 t_plot = np.linspace(wave.t_start, t_max, 100)
141 v_plot = []
142 t_plot_used = []
144 for t in t_plot:
145 # Compute position manually for inactive waves when show_inactive=True
146 v = wave.position_at_time(t) if wave.is_active else wave.v_start + wave.velocity * (t - wave.t_start)
148 if v is not None and 0 <= v <= state.v_outlet:
149 v_plot.append(v)
150 t_plot_used.append(t)
151 elif v is not None and v > state.v_outlet:
152 # Wave crossed outlet - add exact intersection point
153 vel = wave.velocity
154 if vel is not None and vel > 0:
155 t_outlet = wave.t_start + (state.v_outlet - wave.v_start) / vel
156 if wave.t_start <= t_outlet <= t_max:
157 v_plot.append(state.v_outlet)
158 t_plot_used.append(t_outlet)
159 break
160 else:
161 break
163 if len(v_plot) > 0:
164 alpha = 0.5 if not wave.is_active else 1.0
165 ax.plot(
166 t_plot_used,
167 v_plot,
168 "r-",
169 linewidth=2,
170 alpha=alpha,
171 label="Shock" if not shock_labeled else "",
172 )
173 shock_labeled = True
175 # Plot rarefactions (green fans)
176 for wave in state.waves:
177 if isinstance(wave, RarefactionWave):
178 if not wave.is_active and not show_inactive:
179 continue
181 t_plot = np.linspace(wave.t_start, t_max, 100)
182 v_head_plot = []
183 v_tail_plot = []
184 t_plot_used = []
185 head_crossed = False
186 tail_crossed = False
188 for t in t_plot:
189 # Compute positions manually for inactive waves when show_inactive=True
190 if wave.is_active:
191 v_head = wave.head_position_at_time(t)
192 v_tail = wave.tail_position_at_time(t)
193 else:
194 # Manually compute positions for visualization
195 v_head = wave.v_start + wave.head_velocity() * (t - wave.t_start)
196 v_tail = wave.v_start + wave.tail_velocity() * (t - wave.t_start)
198 # Track time points
199 t_plot_used.append(t)
201 # Handle head
202 if v_head is not None and 0 <= v_head <= state.v_outlet:
203 v_head_plot.append(v_head)
204 elif v_head is not None and v_head > state.v_outlet and not head_crossed:
205 # Add exact outlet intersection for head
206 head_vel = wave.head_velocity()
207 if head_vel > 0:
208 t_outlet_head = wave.t_start + (state.v_outlet - wave.v_start) / head_vel
209 if wave.t_start <= t_outlet_head <= t_max:
210 # Insert the exact crossing point
211 v_head_plot.append(state.v_outlet)
212 head_crossed = True
213 v_head_plot.append(None)
214 else:
215 v_head_plot.append(None)
217 # Handle tail
218 if v_tail is not None and 0 <= v_tail <= state.v_outlet:
219 v_tail_plot.append(v_tail)
220 elif v_tail is not None and v_tail > state.v_outlet and not tail_crossed:
221 # Add exact outlet intersection for tail
222 tail_vel = wave.tail_velocity()
223 if tail_vel > 0:
224 t_outlet_tail = wave.t_start + (state.v_outlet - wave.v_start) / tail_vel
225 if wave.t_start <= t_outlet_tail <= t_max:
226 # Insert the exact crossing point
227 v_tail_plot.append(state.v_outlet)
228 tail_crossed = True
229 v_tail_plot.append(None)
230 else:
231 v_tail_plot.append(None)
233 # Plot head and tail boundaries
234 alpha = 0.5 if not wave.is_active else 0.8
235 label = "Rarefaction" if not raref_labeled else ""
237 # Plot head (faster boundary)
238 valid_head = [(t, v) for t, v in zip(t_plot_used, v_head_plot, strict=False) if v is not None]
239 if valid_head:
240 t_h, v_h = zip(*valid_head, strict=False)
241 ax.plot(t_h, v_h, "g-", linewidth=1.5, alpha=alpha, label=label)
242 raref_labeled = True
244 # Plot tail (slower boundary)
245 valid_tail = [(t, v) for t, v in zip(t_plot_used, v_tail_plot, strict=False) if v is not None]
246 if valid_tail:
247 t_t, v_t = zip(*valid_tail, strict=False)
248 ax.plot(t_t, v_t, "g--", linewidth=1.5, alpha=alpha)
250 # Fill between head and tail
251 if valid_head and valid_tail and len(valid_head) == len(valid_tail):
252 ax.fill_between(
253 t_h, # pyright: ignore[reportPossiblyUnboundVariable]
254 v_h, # pyright: ignore[reportPossiblyUnboundVariable]
255 v_t, # pyright: ignore[reportPossiblyUnboundVariable]
256 color="green",
257 alpha=0.1 if not wave.is_active else 0.2,
258 )
260 # Plot outlet position
261 ax.axhline(
262 state.v_outlet,
263 color="k",
264 linestyle="--",
265 linewidth=1,
266 alpha=0.5,
267 label=f"Outlet (V={state.v_outlet:.1f} m³)",
268 )
270 # Plot inlet position
271 ax.axhline(
272 0.0,
273 color="k",
274 linestyle=":",
275 linewidth=1,
276 alpha=0.5,
277 label="Inlet (V=0)",
278 )
280 # Plot wave interaction events as markers
281 if show_events and hasattr(state, "events") and state.events:
282 for event in state.events:
283 if "time" in event and "position" in event:
284 t_event = event["time"]
285 v_event = event["position"]
286 if 0 <= t_event <= t_max and 0 <= v_event <= state.v_outlet:
287 # Determine marker style based on event type
288 event_type = event.get("type", "unknown")
289 if "shock" in event_type.lower() or "collision" in event_type.lower():
290 marker = "X"
291 color = "red"
292 size = 100
293 elif "rarefaction" in event_type.lower():
294 marker = "o"
295 color = "green"
296 size = 80
297 elif "outlet" in event_type.lower():
298 marker = "s"
299 color = "black"
300 size = 80
301 else:
302 marker = "D"
303 color = "gray"
304 size = 60
306 ax.scatter(
307 t_event,
308 v_event,
309 marker=marker,
310 s=size,
311 color=color,
312 edgecolors="black",
313 linewidths=1.5,
314 alpha=0.8,
315 zorder=10,
316 label="Event" if not event_labeled else "",
317 )
318 event_labeled = True
320 ax.set_xlabel("Time [days]", fontsize=12)
321 ax.set_ylabel("Position (Pore Volume) [m³]", fontsize=12)
322 ax.set_title("V-t Diagram: Front Tracking Simulation", fontsize=14, fontweight="bold")
323 ax.grid(True, alpha=0.3)
324 ax.legend(loc="best")
325 ax.set_xlim(0, t_max)
326 ax.set_ylim(-state.v_outlet * 0.05, state.v_outlet * 1.05)
328 return ax
331def plot_breakthrough_curve(
332 state: FrontTrackerState,
333 ax=None,
334 t_max: float | None = None,
335 n_rarefaction_points: int = 50,
336 figsize: tuple[float, float] = (12, 6),
337 t_first_arrival: float | None = None,
338):
339 """
340 Plot exact analytical concentration breakthrough curve at outlet.
342 Uses wave segment information to plot the exact analytical solution
343 without discretization. Constant concentration regions are plotted
344 as horizontal lines, and rarefaction regions are plotted using their
345 exact self-similar solutions.
347 Parameters
348 ----------
349 state : FrontTrackerState
350 Complete simulation state containing all waves.
351 ax : matplotlib.axes.Axes, optional
352 Existing axes to plot into. If None, a new figure and axes are created
353 using ``figsize``.
354 t_max : float, optional
355 Maximum time to plot [days]. If None, uses final simulation time * 1.1.
356 n_rarefaction_points : int, optional
357 Number of points to use for plotting rarefaction segments (analytical
358 curves). Default 50 per rarefaction segment.
359 figsize : tuple of float, optional
360 Figure size in inches (width, height). Default (12, 6).
361 t_first_arrival : float, optional
362 First arrival time for marking spin-up period [days]. If None, spin-up
363 period is not plotted.
365 Returns
366 -------
367 fig : matplotlib.figure.Figure
368 Figure object containing the breakthrough curve.
370 Notes
371 -----
372 - Uses identify_outlet_segments to get exact analytical segment boundaries
373 - Constant concentration segments plotted as horizontal lines (no discretization)
374 - Rarefaction segments plotted using exact self-similar solution
375 - Shocks appear as instantaneous jumps at exact crossing times
376 - No bin averaging or discretization artifacts
378 Examples
379 --------
380 ::
382 from gwtransport.fronttracking.solver import FrontTracker
384 tracker = FrontTracker(cin, flow, tedges, aquifer_pore_volume, sorption)
385 tracker.run()
386 fig = plot_breakthrough_curve(tracker.state)
387 fig.savefig("exact_breakthrough.png")
388 """
389 if ax is None:
390 _, ax = plt.subplots(figsize=figsize)
392 if t_max is None:
393 # Default to input data time range instead of simulation end time
394 # Convert tedges[-1] from Timestamp to days from tedges[0]
395 t_max = float((state.tedges[-1] - state.tedges[0]) / pd.Timedelta(days=1))
397 # Use exact analytical segments
398 segments = identify_outlet_segments(0.0, t_max, state.v_outlet, state.waves, state.sorption)
400 for i, segment in enumerate(segments):
401 t_start = segment["t_start"]
402 t_end = segment["t_end"]
404 if segment["type"] == "constant":
405 # Constant concentration segment - plot as horizontal line
406 c_const = segment["concentration"]
407 ax.plot(
408 [t_start, t_end],
409 [c_const, c_const],
410 "b-",
411 linewidth=2,
412 label="Outlet concentration" if i == 0 else "",
413 )
414 elif segment["type"] == "rarefaction":
415 # Rarefaction segment - plot exact analytical curve
416 raref = segment["wave"]
417 t_raref = np.linspace(t_start, t_end, n_rarefaction_points)
418 c_raref = np.zeros_like(t_raref)
420 for j, t in enumerate(t_raref):
421 # Use the rarefaction wave's own concentration_at_point method
422 c_at_point = raref.concentration_at_point(state.v_outlet, t)
423 if c_at_point is not None:
424 c_raref[j] = c_at_point
425 else:
426 # Fallback to boundary values if not in fan
427 c_raref[j] = segment.get("c_start", raref.c_tail)
429 ax.plot(t_raref, c_raref, "b-", linewidth=2, label="Outlet concentration" if i == 0 else "")
431 # Mark first arrival time if provided
432 if t_first_arrival is not None and np.isfinite(t_first_arrival):
433 ax.axvline(
434 t_first_arrival,
435 color="r",
436 linestyle="--",
437 linewidth=1.5,
438 alpha=0.7,
439 label=f"First arrival (t={t_first_arrival:.2f} days)",
440 )
442 # Shade spin-up region
443 ax.axvspan(
444 0,
445 t_first_arrival,
446 alpha=0.1,
447 color="gray",
448 label="Spin-up period",
449 )
451 ax.set_xlabel("Time [days]", fontsize=12)
452 ax.set_ylabel("Concentration [mass/volume]", fontsize=12)
453 ax.set_title("Breakthrough Curve at Outlet (Exact Analytical)", fontsize=14, fontweight="bold")
454 ax.grid(True, alpha=0.3)
455 ax.legend(loc="best")
456 ax.set_xlim(0, t_max)
457 ax.set_ylim(bottom=0)
459 return ax
462def plot_wave_interactions(
463 state: FrontTrackerState,
464 figsize: tuple[float, float] = (14, 8),
465 ax=None,
466):
467 """
468 Plot event timeline showing wave interactions.
470 Creates a scatter plot showing when and where different types of wave
471 interactions occur during the simulation.
473 Parameters
474 ----------
475 state : FrontTrackerState
476 Complete simulation state containing all events.
477 ax : matplotlib.axes.Axes, optional
478 Existing axes to plot into. If None, a new figure and axes are created
479 using ``figsize``.
480 figsize : tuple of float, optional
481 Figure size in inches (width, height). Default (14, 8).
483 Returns
484 -------
485 fig : matplotlib.figure.Figure
486 Figure object containing the event timeline.
488 Notes
489 -----
490 - Each event type is shown with a different color and marker.
491 - Outlet crossings are shown separately from internal collisions.
492 - Event locations are plotted in the (time, position) plane.
494 Examples
495 --------
496 ::
498 from gwtransport.fronttracking.solver import FrontTracker
500 tracker = FrontTracker(cin, flow, tedges, aquifer_pore_volume, sorption)
501 tracker.run()
502 fig = plot_wave_interactions(tracker.state)
503 fig.savefig("wave_interactions.png")
504 """
505 if ax is None:
506 _, ax = plt.subplots(figsize=figsize)
508 # Group events by type
509 event_types = {}
510 for event_dict in state.events:
511 event_type = event_dict["type"]
512 if event_type not in event_types:
513 event_types[event_type] = {"times": [], "locations": []}
514 event_types[event_type]["times"].append(event_dict["time"])
515 event_types[event_type]["locations"].append(event_dict.get("location", 0.0))
517 # Define colors and markers for each event type
518 event_style = {
519 "CHAR_CHAR_COLLISION": {"color": "blue", "marker": "o", "label": "Char-Char"},
520 "SHOCK_SHOCK_COLLISION": {"color": "red", "marker": "s", "label": "Shock-Shock"},
521 "SHOCK_CHAR_COLLISION": {"color": "purple", "marker": "^", "label": "Shock-Char"},
522 "RAREF_CHAR_COLLISION": {"color": "green", "marker": "v", "label": "Raref-Char"},
523 "SHOCK_RAREF_COLLISION": {"color": "orange", "marker": "d", "label": "Shock-Raref"},
524 "RAREF_RAREF_COLLISION": {"color": "cyan", "marker": "p", "label": "Raref-Raref"},
525 "OUTLET_CROSSING": {"color": "black", "marker": "x", "label": "Outlet Crossing"},
526 "INLET_CHANGE": {"color": "gray", "marker": "+", "label": "Inlet Change"},
527 }
529 # Plot each event type
530 for event_type, data in event_types.items():
531 style = event_style.get(event_type, {"color": "gray", "marker": "o", "label": event_type})
532 ax.scatter(
533 data["times"],
534 data["locations"],
535 c=style["color"],
536 marker=style["marker"],
537 s=100,
538 alpha=0.7,
539 label=f"{style['label']} ({len(data['times'])})",
540 )
542 # Plot outlet line for reference
543 if state.events:
544 ax.axhline(
545 state.v_outlet,
546 color="k",
547 linestyle="--",
548 linewidth=1,
549 alpha=0.3,
550 label=f"Outlet (V={state.v_outlet:.1f} m³)",
551 )
553 ax.set_xlabel("Time [days]", fontsize=12)
554 ax.set_ylabel("Position (Pore Volume) [m³]", fontsize=12)
555 ax.set_title("Wave Interaction Events", fontsize=14, fontweight="bold")
556 ax.grid(True, alpha=0.3)
557 ax.legend(loc="best", ncol=2)
559 if state.events:
560 ax.set_xlim(left=0)
561 ax.set_ylim(-state.v_outlet * 0.05, state.v_outlet * 1.05)
563 return ax
566def plot_inlet_concentration(
567 tedges,
568 cin,
569 ax=None,
570 *,
571 t_first_arrival=None,
572 event_markers=None,
573 color="blue",
574 t_max=None,
575 xlabel="Time [days]",
576 ylabel="Concentration",
577 title="Inlet Concentration",
578 figsize=(8, 5),
579 **step_kwargs,
580):
581 """
582 Plot inlet concentration as step function with optional markers.
584 Parameters
585 ----------
586 tedges : pandas.DatetimeIndex
587 Time bin edges for inlet concentration.
588 Length = len(cin) + 1.
589 cin : array-like
590 Inlet concentration values.
591 Length = len(tedges) - 1.
592 ax : matplotlib.axes.Axes, optional
593 Existing axes to plot into. If None, creates new figure.
594 t_first_arrival : float, optional
595 First arrival time to mark with vertical line [days].
596 event_markers : list of dict, optional
597 Event markers to add. Each dict should have keys: 'time', 'label', 'color'.
598 color : str, optional
599 Color for inlet concentration line. Default 'blue'.
600 t_max : float, optional
601 Maximum time for x-axis [days]. If None, uses full range.
602 xlabel : str, optional
603 Label for x-axis. Default 'Time [days]'.
604 ylabel : str, optional
605 Label for y-axis. Default 'Concentration'.
606 title : str, optional
607 Plot title. Default 'Inlet Concentration'.
608 figsize : tuple, optional
609 Figure size if creating new figure. Default (8, 5).
610 **step_kwargs
611 Additional arguments passed to ax.plot().
613 Returns
614 -------
615 fig : matplotlib.figure.Figure
616 Figure object.
617 ax : matplotlib.axes.Axes
618 Axes object.
619 """
620 if ax is None:
621 _, ax = plt.subplots(figsize=figsize)
623 # Convert tedges to days from start
624 tedges_array = tedges.to_numpy() if hasattr(tedges, "to_numpy") else np.array(tedges)
626 t_days = (tedges_array - tedges_array[0]) / pd.Timedelta(days=1)
628 # Plot inlet concentration using step function
629 x_plot, y_plot = step_plot_coords(t_days, cin)
630 ax.plot(x_plot, y_plot, linewidth=2, color=color, label="Inlet", **step_kwargs)
632 # Add first arrival marker if provided
633 if t_first_arrival is not None and np.isfinite(t_first_arrival):
634 ax.axvline(
635 t_first_arrival,
636 color="green",
637 linestyle="--",
638 linewidth=1.5,
639 alpha=0.7,
640 label=f"First arrival ({t_first_arrival:.1f} days)",
641 )
643 # Add event markers if provided
644 if event_markers is not None:
645 for marker in event_markers:
646 t = marker.get("time")
647 label = marker.get("label", "")
648 marker_color = marker.get("color", "gray")
649 linestyle = marker.get("linestyle", "--")
651 if t is not None:
652 ax.axvline(
653 t,
654 color=marker_color,
655 linestyle=linestyle,
656 linewidth=1.5,
657 alpha=0.7,
658 label=label,
659 )
661 ax.set_xlabel(xlabel, fontsize=10)
662 ax.set_ylabel(ylabel, fontsize=10)
663 ax.set_title(title, fontsize=12, fontweight="bold")
664 ax.grid(True, alpha=0.3)
665 ax.legend(fontsize=8)
667 if t_max is not None:
668 ax.set_xlim(0, t_max)
669 else:
670 ax.set_xlim(0, t_days[-1])
672 return ax
675def plot_front_tracking_summary(
676 structure,
677 tedges,
678 cin,
679 cout_tedges,
680 cout,
681 *,
682 figsize=(16, 10),
683 show_exact=True,
684 show_bin_averaged=True,
685 show_events=True,
686 show_inactive=False,
687 t_max=None,
688 title=None,
689 inlet_color="blue",
690 outlet_exact_color="blue",
691 outlet_binned_color="red",
692 first_arrival_color="green",
693):
694 """
695 Create comprehensive 3-panel summary figure for front tracking simulation.
697 Creates a multi-panel visualization with:
698 - Top-left: V-t diagram showing wave propagation
699 - Top-right: Inlet concentration time series
700 - Bottom: Outlet concentration (exact and/or bin-averaged)
702 Parameters
703 ----------
704 structure : dict
705 Structure returned from infiltration_to_extraction_front_tracking_detailed.
706 Must contain keys: 'tracker_state', 't_first_arrival'.
707 tedges : pandas.DatetimeIndex
708 Time bin edges for inlet concentration.
709 Length = len(cin) + 1.
710 cin : array-like
711 Inlet concentration values.
712 Length = len(tedges) - 1.
713 cout_tedges : pandas.DatetimeIndex
714 Output time bin edges for bin-averaged concentration.
715 Length = len(cout) + 1.
716 cout : array-like
717 Bin-averaged output concentration values.
718 Length = len(cout_tedges) - 1.
719 figsize : tuple, optional
720 Figure size (width, height). Default (16, 10).
721 show_exact : bool, optional
722 Whether to show exact analytical breakthrough curve. Default True.
723 show_bin_averaged : bool, optional
724 Whether to show bin-averaged concentration. Default True.
725 show_events : bool, optional
726 Whether to show wave interaction events on V-t diagram. Default True.
727 show_inactive : bool, optional
728 Whether to show inactive waves on V-t diagram. Default False.
729 t_max : float, optional
730 Maximum time for plots [days]. If None, uses input data range.
731 title : str, optional
732 Overall figure title. If None, uses generic title.
733 inlet_color : str, optional
734 Color for inlet concentration. Default 'blue'.
735 outlet_exact_color : str, optional
736 Color for exact outlet curve. Default 'blue'.
737 outlet_binned_color : str, optional
738 Color for bin-averaged outlet. Default 'red'.
739 first_arrival_color : str, optional
740 Color for first arrival marker. Default 'green'.
742 Returns
743 -------
744 fig : matplotlib.figure.Figure
745 Figure object.
746 axes : dict
747 Dictionary with keys 'vt', 'inlet', 'outlet' containing axes objects.
748 """
749 # Create figure with 3-panel layout
750 fig = plt.figure(figsize=figsize)
751 gs = fig.add_gridspec(2, 2, hspace=0.3, wspace=0.3)
753 axes = {}
755 # Top left: V-t diagram
756 ax_vt = fig.add_subplot(gs[0, 0])
757 plot_vt_diagram(
758 structure["tracker_state"],
759 ax=ax_vt,
760 show_inactive=show_inactive,
761 show_events=show_events,
762 t_max=t_max,
763 )
764 ax_vt.set_title("V-t Diagram", fontsize=12, fontweight="bold")
765 axes["vt"] = ax_vt
767 # Top right: Inlet concentration
768 ax_inlet = fig.add_subplot(gs[0, 1])
769 plot_inlet_concentration(
770 tedges,
771 cin,
772 ax=ax_inlet,
773 t_first_arrival=structure["t_first_arrival"],
774 color=inlet_color,
775 t_max=t_max,
776 )
777 axes["inlet"] = ax_inlet
779 # Bottom: Outlet concentration (exact and bin-averaged)
780 ax_outlet = fig.add_subplot(gs[1, :])
782 # Compute time range
783 if t_max is None:
784 t_max = (tedges.to_numpy()[-1] - tedges.to_numpy()[0]) / pd.Timedelta(days=1)
786 # Exact breakthrough curve
787 if show_exact:
788 t_exact = np.linspace(0, t_max, 1000)
789 c_exact = [
790 concentration_at_point(
791 structure["tracker_state"].v_outlet,
792 t,
793 structure["tracker_state"].waves,
794 structure["tracker_state"].sorption,
795 )
796 for t in t_exact
797 ]
798 ax_outlet.plot(
799 t_exact,
800 c_exact,
801 color=outlet_exact_color,
802 linewidth=2.5,
803 label="Exact outlet concentration",
804 zorder=3,
805 )
807 # Bin-averaged outlet
808 if show_bin_averaged:
809 t_edges_days = ((cout_tedges - cout_tedges[0]) / pd.Timedelta(days=1)).values
810 xstep_cout, ystep_cout = step_plot_coords(t_edges_days, cout)
811 ax_outlet.plot(
812 xstep_cout,
813 ystep_cout,
814 color=outlet_binned_color,
815 linestyle="--",
816 linewidth=1.5,
817 alpha=0.7,
818 label="Bin-averaged outlet",
819 zorder=2,
820 )
822 # First arrival marker
823 t_first = structure["t_first_arrival"]
824 if np.isfinite(t_first):
825 ax_outlet.axvline(
826 t_first,
827 color=first_arrival_color,
828 linestyle="--",
829 linewidth=1.5,
830 alpha=0.7,
831 label=f"First arrival ({t_first:.1f} days)",
832 zorder=1,
833 )
835 ax_outlet.set_xlabel("Time [days]", fontsize=11)
836 ax_outlet.set_ylabel("Concentration", fontsize=11)
837 ax_outlet.set_title("Outlet Concentration: Exact vs Bin-Averaged", fontsize=12, fontweight="bold")
838 ax_outlet.grid(True, alpha=0.3)
839 ax_outlet.legend(fontsize=9)
840 ax_outlet.set_xlim(0, t_max)
841 axes["outlet"] = ax_outlet
843 # Overall title
844 if title is not None:
845 plt.suptitle(title, fontsize=14, fontweight="bold", y=0.995)
847 return fig
850def plot_sorption_comparison(
851 pulse_favorable_structure,
852 pulse_unfavorable_structure,
853 pulse_tedges,
854 pulse_cin,
855 dip_favorable_structure,
856 dip_unfavorable_structure,
857 dip_tedges,
858 dip_cin,
859 *,
860 figsize=(16, 12),
861 t_max_pulse=None,
862 t_max_dip=None,
863):
864 """
865 Compare how each inlet produces different outputs with n>1 vs n<1 sorption.
867 Creates a 2x3 grid:
868 - Row 1: Pulse inlet and its outputs with n>1 and n<1 sorption
869 - Row 2: Dip inlet and its outputs with n>1 and n<1 sorption
871 This demonstrates how the SAME inlet timeseries produces DIFFERENT breakthrough
872 curves depending on the sorption isotherm.
874 Parameters
875 ----------
876 pulse_favorable_structure : dict
877 Structure from pulse inlet with n>1 (higher C travels faster).
878 pulse_unfavorable_structure : dict
879 Structure from pulse inlet with n<1 (lower C travels faster).
880 pulse_tedges : pandas.DatetimeIndex
881 Time bin edges for pulse inlet.
882 Length = len(pulse_cin) + 1.
883 pulse_cin : array-like
884 Pulse inlet concentration (e.g., 0→10→0).
885 Length = len(pulse_tedges) - 1.
886 dip_favorable_structure : dict
887 Structure from dip inlet with n>1 (higher C travels faster).
888 dip_unfavorable_structure : dict
889 Structure from dip inlet with n<1 (lower C travels faster).
890 dip_tedges : pandas.DatetimeIndex
891 Time bin edges for dip inlet.
892 Length = len(dip_cin) + 1.
893 dip_cin : array-like
894 Dip inlet concentration (e.g., 10→2→10).
895 Length = len(dip_tedges) - 1.
896 figsize : tuple, optional
897 Figure size (width, height). Default (16, 12).
898 t_max_pulse : float, optional
899 Max time for pulse plots [days]. If None, auto-computed.
900 t_max_dip : float, optional
901 Max time for dip plots [days]. If None, auto-computed.
903 Returns
904 -------
905 fig : matplotlib.figure.Figure
906 Figure object.
907 axes : numpy.ndarray
908 2x3 array of axes objects.
909 """
910 fig, axes = plt.subplots(2, 3, figsize=figsize)
911 fig.suptitle(
912 "Sorption Comparison: How Each Inlet Responds to n>1 vs n<1 Sorption",
913 fontsize=15,
914 fontweight="bold",
915 y=0.995,
916 )
918 # Compute time ranges
919 if t_max_pulse is None:
920 t_max_pulse = (pulse_tedges.to_numpy()[-1] - pulse_tedges.to_numpy()[0]) / pd.Timedelta(days=1)
921 if t_max_dip is None:
922 t_max_dip = (dip_tedges.to_numpy()[-1] - dip_tedges.to_numpy()[0]) / pd.Timedelta(days=1)
924 # === ROW 1: Pulse inlet (0→10→0) ===
925 t_days_pulse = (pulse_tedges.to_numpy() - pulse_tedges.to_numpy()[0]) / pd.Timedelta(days=1)
927 # Column 1: Pulse inlet
928 ax_pulse_inlet = axes[0, 0]
929 x_pulse, y_pulse = step_plot_coords(t_days_pulse, pulse_cin)
930 ax_pulse_inlet.plot(x_pulse, y_pulse, linewidth=2.5, color="black")
931 ax_pulse_inlet.set_xlabel("Time [days]", fontsize=10)
932 ax_pulse_inlet.set_ylabel("Concentration", fontsize=10)
933 ax_pulse_inlet.set_title("Pulse Inlet\n(0→10→0)", fontsize=11, fontweight="bold")
934 ax_pulse_inlet.grid(True, alpha=0.3)
935 ax_pulse_inlet.set_xlim(0, t_max_pulse)
937 # Column 2: Pulse → n>1 outlet
938 ax_pulse_fav = axes[0, 1]
939 t_exact_pulse_fav = np.linspace(0, t_max_pulse, 1500)
940 c_exact_pulse_fav = [
941 concentration_at_point(
942 pulse_favorable_structure["tracker_state"].v_outlet,
943 t,
944 pulse_favorable_structure["tracker_state"].waves,
945 pulse_favorable_structure["tracker_state"].sorption,
946 )
947 for t in t_exact_pulse_fav
948 ]
949 ax_pulse_fav.plot(t_exact_pulse_fav, c_exact_pulse_fav, "b-", linewidth=2.5)
950 ax_pulse_fav.set_xlabel("Time [days]", fontsize=10)
951 ax_pulse_fav.set_ylabel("Concentration", fontsize=10)
952 ax_pulse_fav.set_title("n>1\nShock→Rarefaction", fontsize=11, fontweight="bold", color="darkblue")
953 ax_pulse_fav.grid(True, alpha=0.3)
954 ax_pulse_fav.set_xlim(0, t_max_pulse)
955 ax_pulse_fav.text(
956 0.05,
957 0.95,
958 "High C: FAST\nRise: Sharp\nFall: Smooth",
959 transform=ax_pulse_fav.transAxes,
960 verticalalignment="top",
961 bbox={"boxstyle": "round", "facecolor": "lightblue", "alpha": 0.7},
962 fontsize=8,
963 )
965 # Column 3: Pulse → n<1 outlet
966 ax_pulse_unfav = axes[0, 2]
967 t_exact_pulse_unfav = np.linspace(0, t_max_pulse, 1500)
968 c_exact_pulse_unfav = [
969 concentration_at_point(
970 pulse_unfavorable_structure["tracker_state"].v_outlet,
971 t,
972 pulse_unfavorable_structure["tracker_state"].waves,
973 pulse_unfavorable_structure["tracker_state"].sorption,
974 )
975 for t in t_exact_pulse_unfav
976 ]
977 ax_pulse_unfav.plot(t_exact_pulse_unfav, c_exact_pulse_unfav, "r-", linewidth=2.5)
978 ax_pulse_unfav.set_xlabel("Time [days]", fontsize=10)
979 ax_pulse_unfav.set_ylabel("Concentration", fontsize=10)
980 ax_pulse_unfav.set_title("n<1\nRarefaction→Shock", fontsize=11, fontweight="bold", color="darkred")
981 ax_pulse_unfav.grid(True, alpha=0.3)
982 ax_pulse_unfav.set_xlim(0, t_max_pulse)
983 ax_pulse_unfav.text(
984 0.05,
985 0.95,
986 "High C: SLOW\nRise: Smooth\nFall: Sharp",
987 transform=ax_pulse_unfav.transAxes,
988 verticalalignment="top",
989 bbox={"boxstyle": "round", "facecolor": "lightcoral", "alpha": 0.7},
990 fontsize=8,
991 )
993 # === ROW 2: Dip inlet (10→2→10) ===
994 t_days_dip = (dip_tedges.to_numpy() - dip_tedges.to_numpy()[0]) / pd.Timedelta(days=1)
996 # Column 1: Dip inlet
997 ax_dip_inlet = axes[1, 0]
998 x_dip, y_dip = step_plot_coords(t_days_dip, dip_cin)
999 ax_dip_inlet.plot(x_dip, y_dip, linewidth=2.5, color="black")
1000 ax_dip_inlet.set_xlabel("Time [days]", fontsize=10)
1001 ax_dip_inlet.set_ylabel("Concentration", fontsize=10)
1002 ax_dip_inlet.set_title("Dip Inlet\n(10→2→10)", fontsize=11, fontweight="bold")
1003 ax_dip_inlet.grid(True, alpha=0.3)
1004 ax_dip_inlet.set_xlim(0, t_max_dip)
1006 # Column 2: Dip → n>1 outlet
1007 ax_dip_fav = axes[1, 1]
1008 t_exact_dip_fav = np.linspace(0, t_max_dip, 1500)
1009 c_exact_dip_fav = [
1010 concentration_at_point(
1011 dip_favorable_structure["tracker_state"].v_outlet,
1012 t,
1013 dip_favorable_structure["tracker_state"].waves,
1014 dip_favorable_structure["tracker_state"].sorption,
1015 )
1016 for t in t_exact_dip_fav
1017 ]
1018 ax_dip_fav.plot(t_exact_dip_fav, c_exact_dip_fav, "b-", linewidth=2.5)
1019 ax_dip_fav.set_xlabel("Time [days]", fontsize=10)
1020 ax_dip_fav.set_ylabel("Concentration", fontsize=10)
1021 ax_dip_fav.set_title("n>1\nRarefaction→Shock", fontsize=11, fontweight="bold", color="darkblue")
1022 ax_dip_fav.grid(True, alpha=0.3)
1023 ax_dip_fav.set_xlim(0, t_max_dip)
1024 ax_dip_fav.text(
1025 0.05,
1026 0.95,
1027 "High C: FAST\nDrop: Smooth\nRise: Sharp",
1028 transform=ax_dip_fav.transAxes,
1029 verticalalignment="top",
1030 bbox={"boxstyle": "round", "facecolor": "lightblue", "alpha": 0.7},
1031 fontsize=8,
1032 )
1034 # Column 3: Dip → n<1 outlet
1035 ax_dip_unfav = axes[1, 2]
1036 t_exact_dip_unfav = np.linspace(0, t_max_dip, 1500)
1037 c_exact_dip_unfav = [
1038 concentration_at_point(
1039 dip_unfavorable_structure["tracker_state"].v_outlet,
1040 t,
1041 dip_unfavorable_structure["tracker_state"].waves,
1042 dip_unfavorable_structure["tracker_state"].sorption,
1043 )
1044 for t in t_exact_dip_unfav
1045 ]
1046 ax_dip_unfav.plot(t_exact_dip_unfav, c_exact_dip_unfav, "r-", linewidth=2.5)
1047 ax_dip_unfav.set_xlabel("Time [days]", fontsize=10)
1048 ax_dip_unfav.set_ylabel("Concentration", fontsize=10)
1049 ax_dip_unfav.set_title("n<1\nShock→Rarefaction", fontsize=11, fontweight="bold", color="darkred")
1050 ax_dip_unfav.grid(True, alpha=0.3)
1051 ax_dip_unfav.set_xlim(0, t_max_dip)
1052 ax_dip_unfav.text(
1053 0.05,
1054 0.95,
1055 "High C: SLOW\nDrop: Sharp\nRise: Smooth",
1056 transform=ax_dip_unfav.transAxes,
1057 verticalalignment="top",
1058 bbox={"boxstyle": "round", "facecolor": "lightcoral", "alpha": 0.7},
1059 fontsize=8,
1060 )
1062 return fig, axes