Coverage for src/gwtransport/fronttracking/plot.py: 0%
186 statements
« prev ^ index » next coverage.py v7.15.3, created at 2026-08-04 20:54 +0000
1"""
2Visualization functions for front tracking.
4This module provides plotting utilities for visualizing front-tracking simulations:
6- V-t diagrams showing wave propagation in space-time
7- Breakthrough curves showing concentration at outlet over time
9Internally the simulation uses cumulative-flow coordinates (V, θ). All plots
10remain in user-facing time t (days). Translation is done via the state's
11``t_at_theta`` / ``theta_at_t`` methods at the plotting boundary.
13Available functions:
15- :func:`plot_vt_diagram` - Draw every wave of a completed simulation in the (time, position) plane:
16 characteristics as thin blue lines, shocks as thick red lines, rarefactions as filled green fans, with the
17 inlet and the outlet marked as horizontal lines. Each straight-in-θ trajectory is sampled in θ and translated
18 to time t before plotting. Returns the axes; ``show_inactive`` adds the waves deactivated by interactions and
19 ``show_events`` marks the interaction events.
21- :func:`plot_inlet_concentration` - Draw ``cin`` on its ``tedges`` as a step function of time [days], optionally
22 marking a first-arrival time with a vertical line. Returns the axes.
24- :func:`plot_front_tracking_summary` - Three-panel figure for one simulation: the V-t diagram, the inlet
25 concentration, and across the bottom the outlet concentration as the exact analytical breakthrough curve
26 and/or the bin-averaged ``cout`` step trace, both referenced to the tracker's own time origin so the two
27 overlay without a shift. Returns the figure and a dict of axes keyed ``'vt'``, ``'inlet'`` and ``'outlet'``.
29- :func:`plot_sorption_comparison` - 2x3 figure contrasting a pulse inlet and a dip inlet with the exact outlet
30 response each produces under a favorable (``n>1``) and an unfavorable (``n<1``) isotherm: column 0 holds the
31 two inlets, columns 1-2 their responses. Returns the figure and the 2x3 array of axes.
33This file is part of gwtransport which is released under AGPL-3.0 license.
34See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details.
35"""
37import matplotlib.pyplot as plt
38import numpy as np
39import numpy.typing as npt
40import pandas as pd
41from matplotlib.axes import Axes
42from matplotlib.figure import Figure
44from gwtransport._time import tedges_to_days
45from gwtransport.fronttracking.output import compute_breakthrough_curve
46from gwtransport.fronttracking.solver import FrontTrackerState
47from gwtransport.fronttracking.waves import CharacteristicWave, RarefactionWave, ShockWave
48from gwtransport.utils import step_plot_coords
51def _wave_trajectory_in_t(
52 state: FrontTrackerState,
53 theta_start: float,
54 v_start: float,
55 speed: float,
56 t_max: float,
57 *,
58 n_points: int = 100,
59) -> tuple[list[float], list[float]]:
60 """Convert a straight-in-θ wave trajectory into (t, V) samples for plotting.
62 A wave with position ``V(θ) = v_start + speed * (θ - theta_start)`` is
63 sampled at ``n_points`` θ-values between ``theta_start`` and the θ that
64 corresponds to ``t_max`` (or the outlet, whichever comes first), then
65 each sample is translated back to user-facing time via ``state.t_at_theta``.
67 Parameters
68 ----------
69 state : FrontTrackerState
70 Simulation state providing θ↔t translation.
71 theta_start : float
72 θ at which the wave forms [m³].
73 v_start : float
74 V at which the wave forms [m³].
75 speed : float
76 Wave speed dV/dθ.
77 t_max : float
78 Maximum user-facing time [days].
79 n_points : int, optional
80 Number of θ-samples (before clipping to outlet). Default 100.
82 Returns
83 -------
84 t_samples : list of float
85 User-facing times [days], monotonic.
86 v_samples : list of float
87 Wave positions at those times [m³], clipped to ``[0, v_outlet]``.
88 """
89 theta_max = state.theta_at_t(t_max)
90 if theta_max <= theta_start:
91 return [], []
93 # If the wave will exit the domain before θ_max, clip θ to the outlet
94 # crossing θ so we plot exactly up to V = v_outlet.
95 if speed > 0:
96 theta_outlet = theta_start + (state.v_outlet - v_start) / speed
97 theta_end = min(theta_max, theta_outlet)
98 else:
99 theta_end = theta_max
101 if theta_end <= theta_start:
102 return [], []
104 thetas = np.linspace(theta_start, theta_end, n_points)
105 vs = v_start + speed * (thetas - theta_start)
107 mask = (vs >= 0) & (vs <= state.v_outlet)
108 t_arr = [state.t_at_theta(float(theta)) for theta in thetas[mask]]
109 # Callers test truthiness of the returned lists (``if v_head:``), so keep
110 # Python lists rather than arrays.
111 return t_arr, vs[mask].tolist()
114def plot_vt_diagram(
115 state: FrontTrackerState,
116 ax: Axes | None = None,
117 *,
118 t_max: float | None = None,
119 figsize: tuple[float, float] = (14, 10),
120 show_inactive: bool = False,
121 show_events: bool = False,
122) -> Axes:
123 """
124 Create V-t diagram showing all waves in space-time.
126 Plots characteristics (blue lines), shocks (red lines), and rarefactions
127 (green fans) in the (time, position) plane. This visualization shows how
128 waves propagate and interact throughout the simulation.
130 Internally the waves live in (V, θ); each wave's straight-line θ-trajectory
131 is converted back to user-facing time t via ``state.t_at_theta`` before
132 plotting.
134 Parameters
135 ----------
136 state : FrontTrackerState
137 Complete simulation state containing all waves.
138 ax : matplotlib.axes.Axes, optional
139 Existing axes to plot into. If None, a new figure and axes are created
140 using ``figsize``.
141 t_max : float, optional
142 Maximum time to plot [days]. If None, uses the input data time range.
143 figsize : tuple of float, optional
144 Figure size in inches (width, height). Default (14, 10).
145 show_inactive : bool, optional
146 Whether to show inactive waves (deactivated by interactions).
147 Default False.
148 show_events : bool, optional
149 Whether to show wave interaction events as markers.
150 Default False.
152 Returns
153 -------
154 ax : matplotlib.axes.Axes
155 Axes object containing the V-t diagram.
157 See Also
158 --------
159 plot_front_tracking_summary : Multi-panel summary combining these views.
160 gwtransport.advection.infiltration_to_extraction_nonlinear_sorption : Produces the tracker state.
162 Notes
163 -----
164 - Characteristics appear as blue lines (constant speed in θ).
165 - Shocks appear as thick red lines (jump discontinuities).
166 - Rarefactions appear as green fans (smooth transition regions).
167 - Outlet position is shown as a horizontal dashed line.
168 - Only waves within domain [0, v_outlet] are plotted.
170 Examples
171 --------
172 .. disable_try_examples
174 ::
176 from gwtransport.fronttracking.solver import FrontTracker
178 tracker = FrontTracker(cin, flow, tedges, aquifer_pore_volume, sorption)
179 tracker.run()
180 ax = plot_vt_diagram(tracker.state)
181 ax.figure.savefig("vt_diagram.png")
182 """
183 if t_max is None:
184 t_max = float((state.tedges[-1] - state.tedges[0]) / pd.Timedelta(days=1))
186 if ax is None:
187 _, ax = plt.subplots(figsize=figsize)
189 char_labeled = False
190 shock_labeled = False
191 raref_labeled = False
192 event_labeled = False
194 for wave in state.waves:
195 if isinstance(wave, CharacteristicWave):
196 if not wave.is_active and not show_inactive:
197 continue
199 t_plot, v_plot = _wave_trajectory_in_t(state, wave.theta_start, wave.v_start, wave.speed(), t_max)
201 if len(v_plot) > 0:
202 alpha = 0.3 if not wave.is_active else 0.7
203 ax.plot(
204 t_plot,
205 v_plot,
206 "b-",
207 linewidth=0.5,
208 alpha=alpha,
209 label="Characteristic" if not char_labeled else "",
210 )
211 char_labeled = True
213 for wave in state.waves:
214 if isinstance(wave, ShockWave):
215 if not wave.is_active and not show_inactive:
216 continue
218 t_plot, v_plot = _wave_trajectory_in_t(state, wave.theta_start, wave.v_start, wave.speed, t_max)
220 if len(v_plot) > 0:
221 alpha = 0.5 if not wave.is_active else 1.0
222 ax.plot(
223 t_plot,
224 v_plot,
225 "r-",
226 linewidth=2,
227 alpha=alpha,
228 label="Shock" if not shock_labeled else "",
229 )
230 shock_labeled = True
232 for wave in state.waves:
233 if isinstance(wave, RarefactionWave):
234 if not wave.is_active and not show_inactive:
235 continue
237 t_head, v_head = _wave_trajectory_in_t(state, wave.theta_start, wave.v_start, wave.head_speed(), t_max)
238 t_tail, v_tail = _wave_trajectory_in_t(state, wave.theta_start, wave.v_start, wave.tail_speed(), t_max)
240 alpha = 0.5 if not wave.is_active else 0.8
241 label = "Rarefaction" if not raref_labeled else ""
243 if v_head:
244 ax.plot(t_head, v_head, "g-", linewidth=1.5, alpha=alpha, label=label)
245 raref_labeled = True
247 if v_tail:
248 ax.plot(t_tail, v_tail, "g--", linewidth=1.5, alpha=alpha)
250 # Fill between head and tail. Both are sampled from the same set of
251 # θ values, so when neither is clipped at the outlet they correspond
252 # one-to-one in time; sample lengths can differ once one boundary
253 # hits the outlet earlier. Fill only the overlap region.
254 if v_head and v_tail:
255 n_fill = min(len(v_head), len(v_tail))
256 if n_fill > 1:
257 ax.fill_between(
258 t_head[:n_fill],
259 v_head[:n_fill],
260 v_tail[:n_fill],
261 color="green",
262 alpha=0.1 if not wave.is_active else 0.2,
263 )
265 ax.axhline(
266 state.v_outlet,
267 color="k",
268 linestyle="--",
269 linewidth=1,
270 alpha=0.5,
271 label=f"Outlet (V={state.v_outlet:.1f} m³)",
272 )
274 ax.axhline(
275 0.0,
276 color="k",
277 linestyle=":",
278 linewidth=1,
279 alpha=0.5,
280 label="Inlet (V=0)",
281 )
283 # Plot wave interaction events as markers. Event records carry ``"theta"``;
284 # translate to user-facing t for display via ``state.t_at_theta``.
285 if show_events and state.events:
286 for event in state.events:
287 if "theta" in event and "location" in event:
288 t_event = state.t_at_theta(event["theta"])
289 v_event = event["location"]
290 if 0 <= t_event <= t_max and 0 <= v_event <= state.v_outlet:
291 # Determine marker style based on event type
292 event_type = event.get("type", "unknown")
293 if "shock" in event_type.lower() or "collision" in event_type.lower():
294 marker = "X"
295 color = "red"
296 size = 100
297 elif "rarefaction" in event_type.lower():
298 marker = "o"
299 color = "green"
300 size = 80
301 elif "outlet" in event_type.lower():
302 marker = "s"
303 color = "black"
304 size = 80
305 else:
306 marker = "D"
307 color = "gray"
308 size = 60
310 ax.scatter(
311 t_event,
312 v_event,
313 marker=marker,
314 s=size,
315 color=color,
316 edgecolors="black",
317 linewidths=1.5,
318 alpha=0.8,
319 zorder=10,
320 label="Event" if not event_labeled else "",
321 )
322 event_labeled = True
324 ax.set_xlabel("Time [days]", fontsize=12)
325 ax.set_ylabel("Position (Pore Volume) [m³]", fontsize=12)
326 ax.set_title("V-t Diagram: Front Tracking Simulation", fontsize=14, fontweight="bold")
327 ax.grid(True, alpha=0.3)
328 ax.legend(loc="best")
329 ax.set_xlim(0, t_max)
330 ax.set_ylim(-state.v_outlet * 0.05, state.v_outlet * 1.05)
332 return ax
335def plot_inlet_concentration(
336 tedges: pd.DatetimeIndex,
337 cin: npt.ArrayLike,
338 ax: Axes | None = None,
339 *,
340 t_first_arrival: float | None = None,
341 color: str = "blue",
342 t_max: float | None = None,
343 figsize: tuple[float, float] = (8, 5),
344) -> Axes:
345 """
346 Plot inlet concentration as a step function.
348 Parameters
349 ----------
350 tedges : pandas.DatetimeIndex
351 Time bin edges for inlet concentration.
352 Length = len(cin) + 1.
353 cin : array-like
354 Inlet concentration values.
355 Length = len(tedges) - 1.
356 ax : matplotlib.axes.Axes, optional
357 Existing axes to plot into. If None, creates new figure.
358 t_first_arrival : float, optional
359 First arrival time to mark with vertical line [days].
360 color : str, optional
361 Color for inlet concentration line. Default 'blue'.
362 t_max : float, optional
363 Maximum time for x-axis [days]. If None, uses full range.
364 figsize : tuple of float, optional
365 Figure size if creating new figure. Default (8, 5).
367 Returns
368 -------
369 ax : matplotlib.axes.Axes
370 Axes object.
372 See Also
373 --------
374 plot_front_tracking_summary : Multi-panel summary that places this inlet panel.
375 """
376 if ax is None:
377 _, ax = plt.subplots(figsize=figsize)
379 t_days = tedges_to_days(tedges)
381 x_plot, y_plot = step_plot_coords(t_days, cin)
382 ax.plot(x_plot, y_plot, linewidth=2, color=color, label="Inlet")
384 if t_first_arrival is not None and np.isfinite(t_first_arrival):
385 ax.axvline(
386 t_first_arrival,
387 color="green",
388 linestyle="--",
389 linewidth=1.5,
390 alpha=0.7,
391 label=f"First arrival ({t_first_arrival:.1f} days)",
392 )
394 ax.set_xlabel("Time [days]", fontsize=10)
395 ax.set_ylabel("Concentration", fontsize=10)
396 ax.set_title("Inlet Concentration", fontsize=12, fontweight="bold")
397 ax.grid(True, alpha=0.3)
398 ax.legend(fontsize=8)
400 if t_max is not None:
401 ax.set_xlim(0, t_max)
402 else:
403 ax.set_xlim(0, t_days[-1])
405 return ax
408def _outlet_concentration_curve(
409 state: FrontTrackerState,
410 t_array: npt.NDArray[np.floating],
411) -> npt.NDArray[np.floating]:
412 """Sample the exact outlet concentration at the given user-facing times.
414 The ``t`` array is mapped to θ vectorially via ``state.theta_at_t_array``
415 and delegated to :func:`compute_breakthrough_curve` (the outlet body of
416 ``concentration_at_point`` over a θ-array).
418 Parameters
419 ----------
420 state : FrontTrackerState
421 Simulation state.
422 t_array : numpy.ndarray
423 User-facing time points [days].
425 Returns
426 -------
427 c_out : numpy.ndarray
428 Outlet concentrations matching ``t_array``.
429 """
430 theta_array = state.theta_at_t_array(t_array)
431 return compute_breakthrough_curve(theta_array, state.v_outlet, state.waves, state.sorption)
434def plot_front_tracking_summary(
435 structure: dict,
436 tedges: pd.DatetimeIndex,
437 cin: npt.ArrayLike,
438 cout_tedges: pd.DatetimeIndex,
439 cout: npt.ArrayLike,
440 *,
441 figsize: tuple[float, float] = (16, 10),
442 show_exact: bool = True,
443 show_bin_averaged: bool = True,
444 show_events: bool = True,
445 show_inactive: bool = False,
446 t_max: float | None = None,
447 title: str | None = None,
448) -> tuple[Figure, dict]:
449 """
450 Create comprehensive 3-panel summary figure for front tracking simulation.
452 Creates a multi-panel visualization with:
453 - Top-left: V-t diagram showing wave propagation
454 - Top-right: Inlet concentration time series
455 - Bottom: Outlet concentration (exact and/or bin-averaged)
457 Parameters
458 ----------
459 structure : dict
460 Structure returned from infiltration_to_extraction_nonlinear_sorption.
461 Must contain keys: 'tracker_state', 'theta_first_arrival'.
462 tedges : pandas.DatetimeIndex
463 Time bin edges for inlet concentration.
464 Length = len(cin) + 1.
465 cin : array-like
466 Inlet concentration values.
467 Length = len(tedges) - 1.
468 cout_tedges : pandas.DatetimeIndex
469 Output time bin edges for bin-averaged concentration.
470 Length = len(cout) + 1.
471 cout : array-like
472 Bin-averaged output concentration values.
473 Length = len(cout_tedges) - 1.
474 figsize : tuple of float, optional
475 Figure size (width, height). Default (16, 10).
476 show_exact : bool, optional
477 Whether to show exact analytical breakthrough curve. Default True.
478 show_bin_averaged : bool, optional
479 Whether to show bin-averaged concentration. Default True.
480 show_events : bool, optional
481 Whether to show wave interaction events on V-t diagram. Default True.
482 show_inactive : bool, optional
483 Whether to show inactive waves on V-t diagram. Default False.
484 t_max : float, optional
485 Maximum time for plots [days]. If None, uses input data range.
486 title : str, optional
487 Overall figure title. If None, uses generic title.
489 Returns
490 -------
491 fig : matplotlib.figure.Figure
492 Figure object.
493 axes : dict
494 Dictionary with keys 'vt', 'inlet', 'outlet' containing axes objects.
496 See Also
497 --------
498 plot_vt_diagram : The top-left sub-panel.
499 plot_inlet_concentration : The top-right sub-panel.
500 gwtransport.advection.infiltration_to_extraction_nonlinear_sorption : Produces ``structure``.
501 """
502 fig = plt.figure(figsize=figsize)
503 gs = fig.add_gridspec(2, 2, hspace=0.3, wspace=0.3)
505 axes: dict = {}
506 tracker_state: FrontTrackerState = structure["tracker_state"]
508 if t_max is None:
509 t_max = float((tedges[-1] - tedges[0]) / pd.Timedelta(days=1))
511 # Top left: V-t diagram
512 ax_vt = fig.add_subplot(gs[0, 0])
513 plot_vt_diagram(
514 tracker_state,
515 ax=ax_vt,
516 show_inactive=show_inactive,
517 show_events=show_events,
518 t_max=t_max,
519 )
520 ax_vt.set_title("V-t Diagram", fontsize=12, fontweight="bold")
521 axes["vt"] = ax_vt
523 # Top right: Inlet concentration
524 ax_inlet = fig.add_subplot(gs[0, 1])
525 plot_inlet_concentration(
526 tedges,
527 cin,
528 ax=ax_inlet,
529 t_first_arrival=tracker_state.t_at_theta(structure["theta_first_arrival"]),
530 t_max=t_max,
531 )
532 axes["inlet"] = ax_inlet
534 # Bottom: Outlet concentration (exact and bin-averaged)
535 ax_outlet = fig.add_subplot(gs[1, :])
537 if show_exact:
538 t_exact = np.linspace(0, t_max, 1000)
539 c_exact = _outlet_concentration_curve(tracker_state, t_exact)
540 ax_outlet.plot(
541 t_exact,
542 c_exact,
543 color="blue",
544 linewidth=2.5,
545 label="Exact outlet concentration",
546 zorder=3,
547 )
549 if show_bin_averaged:
550 # Share the exact curve's origin (tracker_state.tedges[0]); referencing the output
551 # grid to its own first edge shifts the overlay by (cout_tedges[0] - tedges[0]) days.
552 t_edges_days = tedges_to_days(cout_tedges, ref=tracker_state.tedges[0])
553 xstep_cout, ystep_cout = step_plot_coords(t_edges_days, cout)
554 ax_outlet.plot(
555 xstep_cout,
556 ystep_cout,
557 color="red",
558 linestyle="--",
559 linewidth=1.5,
560 alpha=0.7,
561 label="Bin-averaged outlet",
562 zorder=2,
563 )
565 t_first = tracker_state.t_at_theta(structure["theta_first_arrival"])
566 if np.isfinite(t_first):
567 ax_outlet.axvline(
568 t_first,
569 color="green",
570 linestyle="--",
571 linewidth=1.5,
572 alpha=0.7,
573 label=f"First arrival ({t_first:.1f} days)",
574 zorder=1,
575 )
577 ax_outlet.set_xlabel("Time [days]", fontsize=11)
578 ax_outlet.set_ylabel("Concentration", fontsize=11)
579 ax_outlet.set_title("Outlet Concentration: Exact vs Bin-Averaged", fontsize=12, fontweight="bold")
580 ax_outlet.grid(True, alpha=0.3)
581 ax_outlet.legend(fontsize=9)
582 ax_outlet.set_xlim(0, t_max)
583 axes["outlet"] = ax_outlet
585 if title is not None:
586 plt.suptitle(title, fontsize=14, fontweight="bold", y=0.995)
588 return fig, axes
591def plot_sorption_comparison(
592 pulse_favorable_structure: dict,
593 pulse_unfavorable_structure: dict,
594 pulse_tedges: pd.DatetimeIndex,
595 pulse_cin: npt.ArrayLike,
596 dip_favorable_structure: dict,
597 dip_unfavorable_structure: dict,
598 dip_tedges: pd.DatetimeIndex,
599 dip_cin: npt.ArrayLike,
600 *,
601 figsize: tuple[float, float] = (16, 12),
602 t_max_pulse: float | None = None,
603 t_max_dip: float | None = None,
604) -> tuple[Figure, npt.NDArray]:
605 """
606 Compare how each inlet produces different outputs with n>1 vs n<1 sorption.
608 Creates a 2x3 grid:
609 - Row 1: Pulse inlet and its outputs with n>1 and n<1 sorption
610 - Row 2: Dip inlet and its outputs with n>1 and n<1 sorption
612 This demonstrates how the SAME inlet timeseries produces DIFFERENT breakthrough
613 curves depending on the sorption isotherm.
615 Parameters
616 ----------
617 pulse_favorable_structure : dict
618 Structure from pulse inlet with n>1 (higher C travels faster).
619 pulse_unfavorable_structure : dict
620 Structure from pulse inlet with n<1 (lower C travels faster).
621 pulse_tedges : pandas.DatetimeIndex
622 Time bin edges for pulse inlet.
623 Length = len(pulse_cin) + 1.
624 pulse_cin : array-like
625 Pulse inlet concentration (e.g., 0->10->0).
626 Length = len(pulse_tedges) - 1.
627 dip_favorable_structure : dict
628 Structure from dip inlet with n>1 (higher C travels faster).
629 dip_unfavorable_structure : dict
630 Structure from dip inlet with n<1 (lower C travels faster).
631 dip_tedges : pandas.DatetimeIndex
632 Time bin edges for dip inlet.
633 Length = len(dip_cin) + 1.
634 dip_cin : array-like
635 Dip inlet concentration (e.g., 10->2->10).
636 Length = len(dip_tedges) - 1.
637 figsize : tuple of float, optional
638 Figure size (width, height). Default (16, 12).
639 t_max_pulse : float, optional
640 Max time for pulse plots [days]. If None, auto-computed.
641 t_max_dip : float, optional
642 Max time for dip plots [days]. If None, auto-computed.
644 Returns
645 -------
646 fig : matplotlib.figure.Figure
647 Figure object.
648 axes : numpy.ndarray
649 2x3 array of axes objects.
650 """
651 fig, axes = plt.subplots(2, 3, figsize=figsize)
652 fig.suptitle(
653 "Sorption Comparison: How Each Inlet Responds to n>1 vs n<1 Sorption",
654 fontsize=15,
655 fontweight="bold",
656 y=0.995,
657 )
659 if t_max_pulse is None:
660 t_max_pulse = float((pulse_tedges[-1] - pulse_tedges[0]) / pd.Timedelta(days=1))
661 if t_max_dip is None:
662 t_max_dip = float((dip_tedges[-1] - dip_tedges[0]) / pd.Timedelta(days=1))
664 # Column 0 renders each inlet; columns 1-2 render that inlet's outlet response under the
665 # favorable (n>1) and unfavorable (n<1) isotherm.
666 for ax, panel_tedges, panel_cin, panel_title, panel_t_max in (
667 (axes[0, 0], pulse_tedges, pulse_cin, "Pulse Inlet\n(0->10->0)", t_max_pulse),
668 (axes[1, 0], dip_tedges, dip_cin, "Dip Inlet\n(10->2->10)", t_max_dip),
669 ):
670 x_step, y_step = step_plot_coords(tedges_to_days(panel_tedges), panel_cin)
671 ax.plot(x_step, y_step, linewidth=2.5, color="black")
672 ax.set_xlabel("Time [days]", fontsize=10)
673 ax.set_ylabel("Concentration", fontsize=10)
674 ax.set_title(panel_title, fontsize=11, fontweight="bold")
675 ax.grid(True, alpha=0.3)
676 ax.set_xlim(0, panel_t_max)
678 for ax, structure, panel_t_max, favorable, wave_sequence, annotation in (
679 (
680 axes[0, 1],
681 pulse_favorable_structure,
682 t_max_pulse,
683 True,
684 "Shock->Rarefaction",
685 "High C: FAST\nRise: Sharp\nFall: Smooth",
686 ),
687 (
688 axes[0, 2],
689 pulse_unfavorable_structure,
690 t_max_pulse,
691 False,
692 "Rarefaction->Shock",
693 "High C: SLOW\nRise: Smooth\nFall: Sharp",
694 ),
695 (
696 axes[1, 1],
697 dip_favorable_structure,
698 t_max_dip,
699 True,
700 "Rarefaction->Shock",
701 "High C: FAST\nDrop: Smooth\nRise: Sharp",
702 ),
703 (
704 axes[1, 2],
705 dip_unfavorable_structure,
706 t_max_dip,
707 False,
708 "Shock->Rarefaction",
709 "High C: SLOW\nDrop: Sharp\nRise: Smooth",
710 ),
711 ):
712 line_style, title_color, box_color = (
713 ("b-", "darkblue", "lightblue") if favorable else ("r-", "darkred", "lightcoral")
714 )
715 t_exact = np.linspace(0, panel_t_max, 1500)
716 ax.plot(t_exact, _outlet_concentration_curve(structure["tracker_state"], t_exact), line_style, linewidth=2.5)
717 ax.set_xlabel("Time [days]", fontsize=10)
718 ax.set_ylabel("Concentration", fontsize=10)
719 ax.set_title(
720 f"{'n>1' if favorable else 'n<1'}\n{wave_sequence}", fontsize=11, fontweight="bold", color=title_color
721 )
722 ax.grid(True, alpha=0.3)
723 ax.set_xlim(0, panel_t_max)
724 ax.text(
725 0.05,
726 0.95,
727 annotation,
728 transform=ax.transAxes,
729 verticalalignment="top",
730 bbox={"boxstyle": "round", "facecolor": box_color, "alpha": 0.7},
731 fontsize=8,
732 )
734 return fig, axes