Release 6.7.0

doc/sections/ex-gwt-adv-schemes.tex

@@ -13,7 +13,7 @@ \subsection{Example description}
 
 \noindent where $C$ is the volumetric concentration ($M/L^3$), $S_w$ is water saturation (dimensionless), $\theta$ is porosity (dimensionless), and $\mathbf{q}$ is the specific discharge field ($L/T$). The problem is configured with no dispersion or diffusion terms, making it ideal for testing numerical scheme performance since an analytical solution exists.
 
-The model domain consists of a 100 cm $\times$ 100 cm square. A uniform inflow (flux) of magnitude 0.5 cm/s is applied to the inflow boundaries (left and bottom edges) resulting in a specific discharge field of $\vec{q} = (q_x, q_y) = (0.354, 0.354)$ cm/s, which is oriented at 45 degrees. Prescribed concentrations are applied on the inflow boundaries,while CHD boundary conditions are used on the outflow boundaries which extract any outflowing concentration.
+The model domain consists of a 100 cm $\times$ 100 cm square. A uniform inflow (flux) of magnitude 0.5 cm/s is applied to the inflow boundaries (left and bottom edges) resulting in a specific discharge field of $\vec{q} = (q_x, q_y) = (0.354, 0.354)$ cm/s, which is oriented at 45 degrees. Prescribed concentrations are applied on the inflow boundaries, while CHD boundary conditions are used on the outflow boundaries which extract any outflowing concentration.
 
 The spatial discretization uses three different grid types to test scheme performance across different grid geometries:
 \begin{itemize}
@@ -22,7 +22,7 @@ \subsection{Example description}
 \item \textbf{Voronoi grid}: Grid using Voronoi cells derived from triangular cells
 \end{itemize}
 
-\noindent The simulation time spans 300 seconds using adaptive time stepping with an initial time step of 5 seconds and a Courant number of 0.7. The Courant number is a dimensionless quantity representing the fraction of a cell that fluid travels in a single time step. A value less or equal to 1 is required for stability. A relatively high Courant number is chosen here specifically to challenge the central difference scheme and demonstrate its potential for oscillatory behavior on discontinuous functions. Four advection schemes are compared:
+\noindent The simulation time spans 300 seconds using adaptive time stepping with an initial time step of 5 seconds and a Courant number of 0.7. The Courant number is a dimensionless quantity representing the fraction of a cell that fluid travels in a single time step. A value less or equal to 1 is required for stability of the central-difference scheme. A relatively high Courant number is chosen here specifically to challenge the central difference scheme and demonstrate its potential for oscillatory behavior on discontinuous functions. Four advection schemes are compared:
 \begin{itemize}
 \item \textbf{Upstream}: First-order accurate, stable but diffusive scheme
 \item \textbf{Central difference}: Second-order accurate but prone to oscillations on discontinuities
@@ -52,23 +52,23 @@ \subsubsection{Flow Field}
 Figure~\ref{fig:ex-gwt-adv-schemes-flow} shows the computed head fields for all three grid types. The results verify that the flow field produces the expected uniform gradient with a 45° angle across all grid geometries, providing the necessary foundation for the pure advection transport simulations. Head contours are parallel straight lines as expected for uniform flow, confirming the correct implementation of the boundary conditions and hydraulic properties.
 
 \begin{StandardFigure}{
-    Computed head fields for the three grid types used in the advection scheme comparison. \textit{Left}: Structured rectangular grid. \textit{Center}: Triangular grid. \textit{Right}: Voronoi polygon grid. All grids produce the expected uniform head gradient at 45° angle.
+    Computed head fields for the three grid types used in the advection scheme comparison. \textit{Left}: Structured rectangular grid. \textit{Center}: Triangular grid. \textit{Right}: Voronoi polygon grid. All grids produce the expected uniform head gradient at a 45° angle.
 }{fig:ex-gwt-adv-schemes-flow}{../figures/ex-gwt-adv-schemes-flow.png}
 \end{StandardFigure}
 
 \subsubsection{Concentration Field Results}
 
 Figures~\ref{fig:ex-gwt-adv-schemes-sin²}, \ref{fig:ex-gwt-adv-schemes-step}, and \ref{fig:ex-gwt-adv-schemes-square} show the final concentration distributions at the end of the simulation, which are near the steady-state solution, for the three inflow concentration profiles across all scheme and grid combinations. These 2D concentration maps reveal the characteristic behavior of each numerical scheme, though the specific wave profiles are more clearly visualized in the cross-section analysis presented later. A mask has been applied to filter out all values below 0 and above 1. This has been done to make it easier to compare the results between different schemes.
 
-\paragraph{Sin² wave results} (Figure~\ref{fig:ex-gwt-adv-schemes-sin²}) reveal important differences in scheme performance for this smooth inflow concentration profile. While all schemes perform well on the structured grid, the upstream scheme is noticeable more diffusive. The central difference scheme becomes unstable on the triangular and Voronoi grids. On the structured grid, the central difference scheme is the least diffusive but exhibits undershoots resulting in negative concentrations due to the relatively high Courant number (0.7); lowering the Courant number eliminates this behavior. The UTVD scheme demonstrates robust performance across all grid types, maintaining both accuracy and stability.
+\paragraph{Sin² wave results} (Figure~\ref{fig:ex-gwt-adv-schemes-sin²}) reveal important differences in scheme performance for this smooth inflow concentration profile. While all schemes perform well on the structured grid, the upstream scheme is noticeably more diffusive. The central difference scheme becomes unstable on the triangular and Voronoi grids. On the structured grid, the central difference scheme is the least diffusive but exhibits undershoots resulting in negative concentrations due to the relatively high Courant number (0.7); lowering the Courant number eliminates this behavior. The UTVD scheme demonstrates robust performance across all grid types, maintaining both accuracy and stability.
 \begin{StandardFigure}{
     Final concentration distributions for the sin² wave inflow concentration profile. Rows represent the four advection schemes (upstream, central difference, TVD, UTVD) and columns represent the three grid types (structured, triangular, Voronoi).
-}{fig:ex-gwt-adv-schemes-sin²}{../figures/ex-gwt-adv-schemes-sin²-wave-conc.png}
+}{fig:ex-gwt-adv-schemes-sin²}{../figures/ex-gwt-adv-schemes-sin2-wave-conc.png}
 \end{StandardFigure}
 
-\paragraph{Step wave results} (Figure~\ref{fig:ex-gwt-adv-schemes-step}) show similar trends to the sin² wave results. On the structured grid, the central difference scheme outperforms UTVD in preserving the overall shape of the wave. However, while it only exhibited undershoot for the smooth sin² wave, the sharp discontinuity of the step wave causes both undershoot and overshoot. For the step wave, the central difference scheme becomes unstable on the triangular and Voronoi grids. The TVD scheme outperforms upstream on the structured grid but performs similarly on the triangular and Voronoi grids.The UTVD scheme emerges as the least dissipative scheme that maintains stability across all grid types.
+\paragraph{Step wave results} (Figure~\ref{fig:ex-gwt-adv-schemes-step}) show similar trends to the sin² wave results. On the structured grid, the central difference scheme outperforms UTVD in preserving the overall shape of the wave. However, while it only exhibited undershoot for the smooth sin² wave, the sharp discontinuity of the step wave causes both undershoot and overshoot. For the step wave, the central difference scheme becomes unstable on the triangular and Voronoi grids. The TVD scheme outperforms upstream on the structured grid but performs similarly on the triangular and Voronoi grids. The UTVD scheme emerges as the least dissipative scheme that maintains stability across all grid types.
 \begin{StandardFigure}{
-    Final concentration distributions for the step wave inflow concentration profile. The central difference scheme begins to show oscillatory behavior near the sharp transition, while TVD and UTVD schemes maintain stable solutions without spurious oscillations.
+    Final concentration distributions for the step wave inflow concentration profile. The central difference scheme begins to show oscillatory behavior near the sharp transition, while the TVD and UTVD schemes maintain stable solutions without spurious oscillations.
 }{fig:ex-gwt-adv-schemes-step}{../figures/ex-gwt-adv-schemes-step-wave-conc.png}
 \end{StandardFigure}
 
@@ -102,9 +102,9 @@ \subsubsection{Key Findings}
 \begin{enumerate}
 \item \textbf{Grid-dependent performance}: The central difference scheme performed reasonably on the structured grid but became unstable on the triangular and Voronoi grids, particularly for the discontinuous input concentration profiles. This demonstrates the critical importance of considering grid type when selecting numerical schemes.
 
-\item \textbf{Courant number sensitivity}: The central difference scheme's stability was highly sensitive to the Courant number. At Courant = 0.7, it exhibited undershoots leading to negative concentrations on smooth functions and oscillartory instability for the discontinuous inflow concentration profiles. Reducing the Courant number mitigated these issues for the smooth inflow concentration profile but at the cost of computational efficiency.
+\item \textbf{Courant number sensitivity}: The central difference scheme's stability was highly sensitive to the Courant number. At Courant = 0.7, it exhibited undershoots leading to negative concentrations on smooth functions and oscillartory instability for the discontinuous inflow concentration profiles. Reducing the Courant number by reducing the time step size mitigated these issues for the smooth inflow concentration profile but at the cost of computational efficiency.
 
-\item \textbf{UTVD scheme as overall best choice}: The UTVD scheme consistently outperformed all others across grid types and inflow concentration profiles .It was the least dissipative scheme that maintained stability, and is the only stable scheme that approached the maximum values of analytical solutions for sharp discontinuities.
+\item \textbf{UTVD scheme as overall best choice}: The UTVD scheme consistently outperformed all others across grid types and inflow concentration profiles. It was the least dissipative scheme that maintained stability, and is the only stable scheme that approached the maximum values of analytical solutions for sharp discontinuities.
 
 \item \textbf{TVD limitations on unstructured grids}: While TVD performed well on the structured grid and outperformed the upstream scheme for the discontinuous input concentration profiles, its advantage diminished on the unstructured grids, where it performed similarly to the upstream scheme.
 

scripts/ex-gwt-adv-schemes.py

@@ -8,7 +8,7 @@
 # where C is concentration [g/cm³] and **q** is the specific discharge field = (qx, qy) = (0.354, 0.354) cm/s at 45°. The problem is configured with no dispersion or diffusion terms, making it a perfect test case for numerical scheme performance since an analytical solution exists.
 #
 # **Problem Setup:**
-# - Domain: 100cm x 100cm square with uniform flow at 45° angle
+# - Domain: 100cm x 100cm square with uniform flow at a 45° angle
 # - Boundary conditions: Prescribed concentrations on inflow boundaries
 # - Time: 300 seconds with adaptive time stepping (initial dt = 5s)
 # - Physics: Pure advection without mixing processes (analytical solution available)
@@ -129,6 +129,11 @@
 schemes = ["upstream", "central", "tvd", "utvd"]  # 4 advection schemes
 wave_functions = ["sin²-wave", "step-wave", "square-wave"]  # 3 test functions
 
+# Abbreviations for model names (to fit MF6's 16-character limit)
+grid_abbrev = {"structured": "s", "triangle": "t", "voronoi": "v"}
+wave_abbrev = {"sin2-wave": "sin2", "step-wave": "step", "square-wave": "sqr"}
+scheme_abbrev = {"upstream": "up", "central": "cen", "tvd": "tvd", "utvd": "utvd"}
+
 # Compute discharge components
 angle = math.radians(angledeg)
 qx = specific_discharge * math.cos(angle)  # x-component of specific discharge (cm/s)
@@ -137,6 +142,7 @@
 AXES_FRACTION = "axes fraction"
 OFFSET_POINTS = "offset points"
 
+
 # %% [markdown]
 # # Analytical Solution
 #
@@ -151,8 +157,6 @@
 # - Rotated coordinate: y' = sin(-θ)·x + cos(-θ)·y
 # - Concentration: C(x,y,t) = inlet_signal(y' - v·t)
 # Note: At t=300s, the pattern has fully advected through the domain
-
-
 # %%
 def exact_solution_concentration(x, y, analytical):
     """Calculate exact concentration at any point in the domain.
@@ -378,7 +382,7 @@ def axis_aligned_segment_length(polygon, axis="y", value=0):
 
 # %%
 def build_mf6gwf(grid_type):
-    gwfname = f"flow_{grid_type}"
+    gwfname = f"gwf_{grid_type}"
     sim_ws = workspace / sim_name / Path(gwfname)
     sim = flopy.mf6.MFSimulation(sim_name=sim_name, sim_ws=sim_ws, exe_name="mf6")
 
@@ -509,10 +513,26 @@ def build_mf6gwf(grid_type):
     return sim
 
 
+def convert_superscript(text):
+    map = {
+        "¹": "1",
+        "²": "2",
+        "³": "3",
+    }
+    trans_table = str.maketrans(
+        "".join(map.keys()),
+        "".join(map.values()),
+    )
+    return text.translate(trans_table)
+
+
 def build_mf6gwt(grid_type, scheme, wave_func):
-    pathname = f"trans_{grid_type}_{wave_func}_{scheme}"
-    gwtname = "trans"
-    gwfname = f"flow_{grid_type}"
+    # Use abbreviated names for model name (16-char MF6 limit)
+    wave_converted = convert_superscript(wave_func)
+    gwtname = f"gwt_{grid_abbrev[grid_type]}_{wave_abbrev[wave_converted]}_{scheme_abbrev[scheme]}"
+    # Use descriptive name for workspace
+    pathname = f"gwt_{grid_type}_{wave_converted}_{scheme}"
+    gwfname = f"gwf_{grid_type}"
     sim_ws = workspace / sim_name / Path(pathname)
     sim = flopy.mf6.MFSimulation(sim_name=sim_name, sim_ws=sim_ws, exe_name="mf6")
 
@@ -644,7 +664,7 @@ def plot_flows(gwf_sims):
         fig, axs = plt.subplots(
             1, len(gwf_sims), dpi=300, figsize=figure_size, tight_layout=True
         )
-        fig.suptitle("Head - flow angle 45")
+        fig.suptitle("Head [cm] - flow angle 45 degrees")
 
         for idx, (grid, sim) in enumerate(gwf_sims.items()):
             plot_flow(sim, axs[idx])
@@ -699,7 +719,7 @@ def plot_concentrations(gwt_sims):
                 figsize=(7, 7 * len(schemes) / (len(grids) + 1)),
                 tight_layout=True,
             )
-            fig.suptitle(f"Concentration - {wave_func}")
+            fig.suptitle(f"Concentration [g/cm³] - {wave_func}")
 
             for idx_scheme, scheme in enumerate(schemes):
                 for idx_grid, grid in enumerate(grids):
@@ -736,7 +756,7 @@ def plot_concentrations(gwt_sims):
             if plot_show:
                 plt.show()
             if plot_save:
-                fname = f"{sim_name}-{wave_func}-conc.png"
+                fname = f"{sim_name}-{convert_superscript(wave_func)}-conc.png"
                 fpth = figs_path / fname
                 fig.savefig(fpth)