Merge pull request #256 from MODFLOW-USGS/develop

Release

.github/workflows/ex-rtd.yml

@@ -157,7 +157,7 @@ jobs:
         run: python ci_create_examples_rst.py
 
       - name: Upload artifacts
-        uses: actions/upload-artifact@v3
+        uses: actions/upload-artifact@v4
         with:
           name: rtd-files-for-${{ needs.set_options.outputs.sha }}
           path: |

.github/workflows/ex-workflow.yml

@@ -107,7 +107,7 @@ jobs:
         shell: python
 
       - name: Upload zipfile
-        uses: actions/upload-artifact@v3
+        uses: actions/upload-artifact@v4
         with:
           name: ${{ env.DISTNAME }}.zip
           path: modflow6-examples/${{ env.DISTNAME }}.zip
@@ -201,7 +201,7 @@ jobs:
 
       - name: Upload PDF document
         if: matrix.python == '3.9'
-        uses: actions/upload-artifact@v3
+        uses: actions/upload-artifact@v4
         with:
           name: ${{ env.DISTNAME}}.pdf
           path: modflow6-examples/${{ env.DISTNAME }}.pdf
@@ -222,13 +222,13 @@ jobs:
           format: MM/DD/YYYY HH:mm
 
       - name: Download ${{ env.DISTNAME }}.pdf
-        uses: actions/download-artifact@v3
+        uses: actions/download-artifact@v4.1.7
         with:
           name: ${{ env.DISTNAME }}.pdf
           path: ${{ env.TAG }}
 
       - name: Download ${{ env.DISTNAME }}.zip
-        uses: actions/download-artifact@v3
+        uses: actions/download-artifact@v4.1.7
         with:
           name: ${{ env.DISTNAME }}.zip
           path: ${{ env.TAG }}

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

@@ -37,7 +37,7 @@
 # NOTE: In the current example, ncol is specified as 250.  If users prefer,
 #       this value may be increased to 500, 1,000, or more for a more refined
 #       solution.  Currently, ncol is kept low for faster run times.
-# -
+
 
 # +
 import os
@@ -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

@@ -183,7 +183,7 @@ def create_outer_ring_of_ctrl_vols(outer_v, verts, iverts, xc, yc, ivt, idx):
     # For each outer point, create a new ivert
     # Starts at 9 o'clock and proceeds counter-clockwise
     # As such, 'pt3' will be used by the subsequent iverts
-    # However, 'pt4' will be re-used by each subseqeuent ivert
+    # However, 'pt4' will be reused by each subseqeuent ivert
     pt3_id = pt3_rad = pt3_ang = None
     for ii in np.arange(len(outer_v) - 1):
         # Create all the points involved
@@ -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-csub-p02.py

@@ -199,8 +199,8 @@ def build_models(
     opth = f"{name}.csub.obs"
     csub_csv = opth + ".csv"
     obs = [
-        ("tcomp", "interbed-compaction", "ib1"),
-        ("sk", "sk", "ib1"),
+        ("tcomp", "interbed-compaction", (0,)),
+        ("sk", "sk", (0,)),
         ("qtop", "delay-flowtop", (0,)),
         ("qbot", "delay-flowbot", (0,)),
     ]

scripts/ex-gwf-csub-p03.py

@@ -507,10 +507,19 @@ def build_models(
     opth = f"{name}.csub.obs"
     csub_csv = opth + ".csv"
     obs = [
-        ("cunit", "interbed-compaction", "cunit"),
-        ("aquitard", "interbed-compaction", "aquitard"),
-        ("nodelay", "interbed-compaction", "nodelay"),
-        ("delay", "interbed-compaction", "delay"),
+        ("cunit1", "interbed-compaction", (0,)),
+        ("cunit2", "interbed-compaction", (1,)),
+        ("cunit3", "interbed-compaction", (2,)),
+        ("aquitard6", "interbed-compaction", (5,)),
+        ("aquitard7", "interbed-compaction", (6,)),
+        ("aquitard8", "interbed-compaction", (7,)),
+        ("nodelay4", "interbed-compaction", (3,)),
+        ("nodelay5", "interbed-compaction", (4,)),
+        ("nodelay9", "interbed-compaction", (8,)),
+        ("nodelay10", "interbed-compaction", (9,)),
+        ("delay11", "interbed-compaction", (10,)),
+        ("delay12", "interbed-compaction", (11,)),
+        ("delay13", "interbed-compaction", (12,)),
         ("es14", "estress-cell", (nlay - 1, 0, 0)),
     ]
     for k in (1, 2, 3, 4, 6, 7, 8, 9, 11, 13):
@@ -719,7 +728,7 @@ def get_obs_dataframe(file_name, hash):
         known_hash=f"md5:{hash}",
     )
     df = pd.read_csv(fpath, index_col=0)
-    df.index = pd.to_datetime(df.index.values)
+    df.index = pd.to_datetime(df.index.values, format="%m/%d/%y")
     df.rename({"mean": "observed"}, inplace=True, axis=1)
     return df
 
@@ -793,7 +802,9 @@ def process_csub_obs(fpth):
     # transfer data from temporary storage
     for key in pcomp:
         if key != "TOTAL" and key != "SKELETAL":
-            v[key] = tdata[key].copy()
+            for obs_key in tdata.dtype.names:
+                if key in obs_key:
+                    v[key] += tdata[obs_key]
 
     # calculate skeletal
     for key in tdata.dtype.names[1:]:
@@ -1144,6 +1155,8 @@ def plot_calibration(silent=True):
         name = list(parameters.keys())[1]
         fpath = os.path.join(workspace, name, f"{name}.csub.obs.csv")
         df_sim = get_sim_dataframe(fpath)
+        for key in pcomp:
+            df_sim[key] = df_sim[list(df_sim.filter(regex=key))].sum(axis=1)
         df_sim.rename({"TOTAL": "simulated"}, inplace=True, axis=1)
 
         fname = "boundary_heads.csv"
@@ -1198,7 +1211,7 @@ def plot_calibration(silent=True):
         ix0 = df_sim.index.get_loc("1992-10-01 12:00:00")
 
         # get initial simulated compaction
-        cstart = df_sim.simulated[ix0]
+        cstart = df_sim.simulated.iloc[ix0]
 
         # cut off initial portion of simulated compaction
         df_sim_pc = df_sim[ix0:].copy()
@@ -1372,7 +1385,7 @@ def plot_calibration(silent=True):
             else:
                 dxc = 0.001
                 dyc = 1
-            xc = df_iobs_pc.observed[ixc]
+            xc = df_iobs_pc.observed.iloc[ixc]
             yc = es_obs[ixc]
             styles.add_annotation(
                 ax=ax,

scripts/ex-gwf-csub-p04.py

@@ -142,12 +142,14 @@
 # Create interbed package data
 icsubno = 0
 csub_pakdata = []
+boundname_dict = {}
 for i in range(nrow):
     for j in range(ncol):
         if ib[i, j] < 1 or (i, j) in chd_locs:
             continue
         for k in range(nlay):
             boundname = f"{k + 1:02d}_{i + 1:02d}_{j + 1:02d}"
+            boundname_dict[boundname] = (icsubno,)
             ib_lst = [
                 icsubno,
                 (k, i, j),
@@ -262,14 +264,14 @@ def build_models():
     opth = f"{sim_name}.csub.obs"
     csub_csv = opth + ".csv"
     obs = [
-        ("w1l1", "interbed-compaction", "01_09_10"),
-        ("w1l2", "interbed-compaction", "02_09_10"),
-        ("w1l3", "interbed-compaction", "03_09_10"),
-        ("w1l4", "interbed-compaction", "04_09_10"),
-        ("w2l1", "interbed-compaction", "01_12_07"),
-        ("w2l2", "interbed-compaction", "02_12_07"),
-        ("w2l3", "interbed-compaction", "03_12_07"),
-        ("w2l4", "interbed-compaction", "04_12_07"),
+        ("w1l1", "interbed-compaction", boundname_dict["01_09_10"]),
+        ("w1l2", "interbed-compaction", boundname_dict["02_09_10"]),
+        ("w1l3", "interbed-compaction", boundname_dict["03_09_10"]),
+        ("w1l4", "interbed-compaction", boundname_dict["04_09_10"]),
+        ("w2l1", "interbed-compaction", boundname_dict["01_12_07"]),
+        ("w2l2", "interbed-compaction", boundname_dict["02_12_07"]),
+        ("w2l3", "interbed-compaction", boundname_dict["03_12_07"]),
+        ("w2l4", "interbed-compaction", boundname_dict["04_12_07"]),
         ("s1l1", "coarse-compaction", (0, 8, 9)),
         ("s1l2", "coarse-compaction", (1, 8, 9)),
         ("s1l3", "coarse-compaction", (2, 8, 9)),
@@ -302,13 +304,13 @@ def build_models():
         ("sk1l2", "ske-cell", (1, 8, 9)),
         ("sk2l4", "ske-cell", (3, 11, 6)),
         ("t1l2", "theta", (1, 8, 9)),
-        ("w1qie", "elastic-csub", "02_09_10"),
-        ("w1qii", "inelastic-csub", "02_09_10"),
+        ("w1qie", "elastic-csub", boundname_dict["02_09_10"]),
+        ("w1qii", "inelastic-csub", boundname_dict["02_09_10"]),
         ("w1qaq", "coarse-csub", (1, 8, 9)),
         ("w1qt", "csub-cell", (1, 8, 9)),
         ("w1wc", "wcomp-csub-cell", (1, 8, 9)),
-        ("w2qie", "elastic-csub", "04_12_07"),
-        ("w2qii", "inelastic-csub", "04_12_07"),
+        ("w2qie", "elastic-csub", boundname_dict["04_12_07"]),
+        ("w2qii", "inelastic-csub", boundname_dict["04_12_07"]),
         ("w2qaq", "coarse-csub", (3, 11, 6)),
         ("w2qt ", "csub-cell", (3, 11, 6)),
         ("w2wc", "wcomp-csub-cell", (3, 11, 6)),
@@ -821,7 +823,7 @@ def plot_results(sim, silent=True):
 
 
 # +
-def scenario(silent=True):
+def scenario(silent=False):
     sim = build_models()
     if write:
         write_models(sim, silent=silent)

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