fix(mt3dms-p01): advection scheme was not behaving as expected

doc/sections/ex-gwt-mt3dms-p01.tex

@@ -10,15 +10,15 @@ \subsection{Example description}
 
 All four model scenarios have 101 columns, 1 row, and 1 layer. The first and last columns use constant-head boundaries to simulate steady flow in confined conditions. Because the analytical solution assumes an infinite 1-dimensional flow field, the last column is set far enough from the source to avoid interfering with the final solution after 2,000 days. Initially, the model domain is devoid of solute; however, the first column uses a constant concentration boundary condition to ensure that water entering the simulation has a unit concentration of 1. Additional model parameters are shown in table~\ref{tab:ex-gwt-mt3dms-p01-01}.
 
+For these simulations, all \mf and MT3DMS simulations use the upstream-weighted finite-difference scheme for representing solute advection.
+
 % add 2nd static parameter value table
 \input{../tables/ex-gwt-mt3dms-p01-01}
 
 % for examples without scenarios
 \subsection{Example Results}
 
-Currently no options are available with \mf for simulating solute transport using particle tracking methods [referred to as Method of Characteristics (MOC) in the MT3DMS manual \cite{zheng1999mt3dms}.  Thus, the \mf solution is compared to an MT3DMS solution that uses the third-order total variation diminishing (TVD) option for solving the advection-only problem rather than invoking one of the MOC options available within MT3DMS.  Owing to different approaches between the two codes, namely TVD scheme of  MT3DMS and the second-order approach of \mf, differences between the two solutions and reflected in figure~\ref{fig:ex-gwt-mt3dms-p01a} are expected.  However, the differences are within acceptable tolerances.
-
-The comparison of the MT3DMS and\mf solutions for problem 1a, an advection dominated problem, represents an end-member test as the migrating concentration front is sharp (i.e., discontinuous). In technical terms, the grid Peclet number is infinity for this problem ($P_e$ = $v\Delta x/D_{xx}$ = $\Delta x$/$\alpha_L$ = $\infty$).
+The comparison of the MT3DMS and \mf solutions for problem 1a, an advection dominated problem, represents an end-member test as the migrating concentration front is sharp. Results for both the \mf and MT3DMS results are shown in figure~\ref{fig:ex-gwt-mt3dms-p01a}.  The grid Peclet number is infinity for this problem ($P_e$ = $v\Delta x/D_{xx}$ = $\Delta x$/$\alpha_L$ = $\infty$).  Though not shown here, this advection-only problem can also be represented with the TVD schemes in both \mf and MT3DMS.  For this particular problem, the TVD scheme in MT3DMS (called the ULTIMATE scheme) is better able to minimize numerical dispersion than the TVD scheme in \mf.
 
 % a figure
 \begin{StandardFigure}

scripts/ex-gwe-barends.py

@@ -105,9 +105,9 @@
 
 
 # delr is a calculated value based on the number of desired columns
-assert (
-    L / ncol % 1 == 0
-), "reconsider specification of NCOL such that length of the simulation divided by NCOL results in an even number value"
+assert L / ncol % 1 == 0, (
+    "reconsider specification of NCOL such that length of the simulation divided by NCOL results in an even number value"
+)
 delr = L / ncol  # Width along the row ($m$)
 
 # Calculated values

scripts/ex-gwe-radial.py

@@ -381,9 +381,9 @@ def generate_control_volumes(basedata, verts, flat_list, idx, silent=True):
                         if xchk == id:
                             xcoll.append(subvertx[num])
                             ycoll.append(subverty[num])
-                assert (
-                    len(xcoll) == 2
-                ), "Should only be two points touching the well bore"
+                assert len(xcoll) == 2, (
+                    "Should only be two points touching the well bore"
+                )
 
                 # was getting divide by zero error in MF6 when putting the
                 # centroid right on the line connecting the two points that

scripts/ex-gwf-curvilinear-90.py

@@ -1276,8 +1276,7 @@ def get_grid(self, name="__main__"):
     def add_grid(self, name, grid):
         if name == "" or name == "__main__":
             raise RuntimeError(
-                "\nDisvGridMerger.add_grid:\n"
-                'name = "" or "__main__"\nis not allowed.'
+                '\nDisvGridMerger.add_grid:\nname = "" or "__main__"\nis not allowed.'
             )
         if isinstance(grid, DisvPropertyContainer):
             grid = grid.copy()

scripts/ex-gwf-curvilinear.py

@@ -1276,8 +1276,7 @@ def get_grid(self, name="__main__"):
     def add_grid(self, name, grid):
         if name == "" or name == "__main__":
             raise RuntimeError(
-                "\nDisvGridMerger.add_grid:\n"
-                'name = "" or "__main__"\nis not allowed.'
+                '\nDisvGridMerger.add_grid:\nname = "" or "__main__"\nis not allowed.'
             )
         if isinstance(grid, DisvPropertyContainer):
             grid = grid.copy()

scripts/ex-gwf-radial.py

@@ -632,7 +632,7 @@ def drawdown_times(
         for stp, ts in enumerate(ts_list):
             if show_progress:
                 print(
-                    f"Solving {stp+1:4d} of {nstp}; time = {self.ts2time(ts, r)}",
+                    f"Solving {stp + 1:4d} of {nstp}; time = {self.ts2time(ts, r)}",
                     end="",
                 )
 
@@ -727,7 +727,7 @@ def drawdown_times(
                 "up to:\n"
                 f"   {bessel_roots0} * 2^bessel_loop_limit\nwhere:\n"
                 "   bessel_loop_limit = {bessel_loop_limit}\n"
-                f"resulting in {1024*2**bessel_loop_limit} roots evaluated, "
+                f"resulting in {1024 * 2**bessel_loop_limit} roots evaluated, "
                 "with the last root being {root}\n"
                 f"(That is, the Neuman integral was solved form 0 to {root})\n\n"
                 "You can either ignore this warning\n"

scripts/ex-gwt-mt3dms-p01.py

@@ -127,7 +127,7 @@
 ibound[0, 0, -1] = -1
 
 # Set some static transport related model parameter values
-mixelm = 0  # TVD
+mixelm = 0  # upstream
 rhob = 0.25
 sp2 = 0.0  # red, but not used in this problem
 sconc = np.zeros((nlay, nrow, ncol), dtype=float)
@@ -168,7 +168,6 @@ def build_models(
     dispersivity=0.0,
     retardation=0.0,
     decay=0.0,
-    mixelm=0,
     silent=False,
 ):
     mt3d_ws = os.path.join(workspace, sim_name, "mt3d")