user-specified tracking times

autotest/prt/prt_test_utils.py

@@ -269,17 +269,17 @@ def plot_nodes_and_vertices(
 
     # plot vertices
     verts = mg.verts
-    ax.plot(verts[:, 0], verts[:, 1], "bo")
+    ax.plot(verts[:, 0], verts[:, 1], "bo", alpha=0.5)
     for i in range(ncpl):
         x, y = verts[i, 0], verts[i, 1]
-        ax.annotate(str(i + 1), verts[i, :], color="b")
+        ax.annotate(str(i + 1), verts[i, :], color="b", alpha=0.5)
 
     # plot nodes
     xc, yc = mg.get_xcellcenters_for_layer(0), mg.get_ycellcenters_for_layer(0)
     for i in range(ncpl):
         x, y = xc[i], yc[i]
-        ax.plot(x, y, "o", color="grey")
-        ax.annotate(str(i + 1), (x, y), color="grey")
+        ax.plot(x, y, "o", color="grey", alpha=0.5)
+        ax.annotate(str(i + 1), (x, y), color="grey", alpha=0.5)
 
     # create legend
     ax.legend(

autotest/prt/test_prt_disv.py → autotest/prt/test_prt_disv1.py

@@ -1,5 +1,6 @@
 """
-Tests particle tracking on a vertex (DISV) grid.
+Tests particle tracking on a vertex (DISV) grid
+that reduces to a regular grid.
 
 Two cases are provided, one with valid release
 position and cell correspondences, and another
@@ -30,8 +31,8 @@
 
 from framework import TestFramework
 
-simname = "prtfmi06"
-cases = [f"{simname}", f"{simname}bprp"]
+simname = "prtdisv1"
+cases = [f"{simname}", f"{simname}bprp", f"{simname}trts"]
 
 # model info
 nlay = 1
@@ -51,6 +52,7 @@
 nouter, ninner = 100, 300
 hclose, rclose, relax = 1e-9, 1e-3, 0.97
 porosity = 0.1
+tracktimes = list(np.linspace(0, 11, 1000))
 
 # vertex grid properties
 disvkwargs = get_disv_kwargs(
@@ -211,6 +213,7 @@ def build_prt_sim(idx, gwf_ws, prt_ws, mf6):
         pname="oc",
         track_filerecord=[prt_track_file],
         trackcsv_filerecord=[prt_track_csv_file],
+        tracktimerecord=tracktimes if "trts" in name else None,
     )
 
     # create the flow model interface
@@ -329,6 +332,8 @@ def check_output(idx, test):
 
     # load mf6 pathline results
     mf6_pls = pd.read_csv(prt_ws / prt_track_csv_file, na_filter=False)
+    if "trts" in name:
+        assert len(mf6_pls) == 9189
 
     # make sure pathline df has "name" (boundname) column and default values
     assert "name" in mf6_pls
@@ -359,52 +364,56 @@ def check_output(idx, test):
     spdis = bud.get_data(text="DATA-SPDIS")[0]
     qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(spdis, gwf)
 
-    # setup plot
-    fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(10, 10))
-    for a in ax:
-        a.set_aspect("equal")
-
-    # plot mf6 pathlines in map view
-    pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax[0])
-    pmv.plot_grid()
-    pmv.plot_array(hds[0], alpha=0.2)
-    pmv.plot_vector(qx, qy, normalize=True, color="white")
-    # set zoom area
-    # xmin, xmax = 2050, 4800
-    # ymin, ymax = 5200, 7550
-    plot_nodes_and_vertices(gwf, mg, None, mg.ncpl, ax[0])
-    mf6_plines = mf6_pls.groupby(["iprp", "irpt", "trelease"])
-    for ipl, ((iprp, irpt, trelease), pl) in enumerate(mf6_plines):
-        pl.plot(
-            title="MF6 pathlines",
-            kind="line",
-            x="x",
-            y="y",
-            ax=ax[0],
-            legend=False,
-            color=cm.plasma(ipl / len(mf6_plines)),
-        )
-
-    # plot mp7 pathlines in map view
-    pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax[1])
-    pmv.plot_grid()
-    pmv.plot_array(hds[0], alpha=0.2)
-    pmv.plot_vector(qx, qy, normalize=True, color="white")
-    mp7_plines = mp7_pls.groupby(["particleid"])
-    for ipl, (pid, pl) in enumerate(mp7_plines):
-        pl.plot(
-            title="MP7 pathlines",
-            kind="line",
-            x="x",
-            y="y",
-            ax=ax[1],
-            legend=False,
-            color=cm.plasma(ipl / len(mp7_plines)),
-        )
-
-    # view/save plot
-    # plt.show()
-    plt.savefig(gwf_ws / f"test_{simname}.png")
+    if "bprp" not in name:
+        # setup plot
+        fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(10, 10))
+        for a in ax:
+            a.set_aspect("equal")
+
+        # plot mf6 pathlines in map view
+        pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax[0])
+        pmv.plot_grid()
+        pmv.plot_array(hds[0], alpha=0.2)
+        pmv.plot_vector(qx, qy, normalize=True, color="white")
+        # set zoom area
+        # xmin, xmax = 2050, 4800
+        # ymin, ymax = 5200, 7550
+        # plot labeled nodes and vertices
+        plot_nodes_and_vertices(gwf, mg, None, mg.ncpl, ax[0])
+        mf6_plines = mf6_pls.groupby(["iprp", "irpt", "trelease"])
+        for ipl, ((iprp, irpt, trelease), pl) in enumerate(mf6_plines):
+            pl.plot(
+                title="MF6 pathlines",
+                linestyle="None" if "trst" in name else "--",
+                marker="o",
+                markersize=2,
+                x="x",
+                y="y",
+                ax=ax[0],
+                legend=False,
+                color=cm.plasma(ipl / len(mf6_plines)),
+            )
+
+        # plot mp7 pathlines in map view
+        pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax[1])
+        pmv.plot_grid()
+        pmv.plot_array(hds[0], alpha=0.2)
+        pmv.plot_vector(qx, qy, normalize=True, color="white")
+        mp7_plines = mp7_pls.groupby(["particleid"])
+        for ipl, (pid, pl) in enumerate(mp7_plines):
+            pl.plot(
+                title="MP7 pathlines",
+                kind="line",
+                x="x",
+                y="y",
+                ax=ax[1],
+                legend=False,
+                color=cm.plasma(ipl / len(mp7_plines)),
+            )
+
+        # view/save plot
+        # plt.show()
+        plt.savefig(gwf_ws / f"test_{simname}.png")
 
     # convert mf6 pathlines to mp7 format
     mf6_pls = to_mp7_pathlines(mf6_pls)
@@ -428,8 +437,9 @@ def check_output(idx, test):
     del mp7_pls["node"]
 
     # compare mf6 / mp7 pathline data
-    assert mf6_pls.shape == mp7_pls.shape
-    assert np.allclose(mf6_pls, mp7_pls, atol=1e-3)
+    if "trts" not in name:
+        assert mf6_pls.shape == mp7_pls.shape
+        assert np.allclose(mf6_pls, mp7_pls, atol=1e-3)
 
 
 @pytest.mark.parametrize("idx, name", enumerate(cases))

autotest/prt/test_prt_track_events.py

@@ -44,7 +44,7 @@
 
 from framework import TestFramework
 
-simname = "prtfmi02"
+simname = "prtevnt"
 cases = [
     f"{simname}all",
     f"{simname}rel",
@@ -53,6 +53,7 @@
     f"{simname}term",
     f"{simname}wksk",
     f"{simname}mult",
+    f"{simname}trts",
 ]
 releasepts_prt = {
     "a": [
@@ -78,6 +79,7 @@
         for i in range(5)
     ],
 }
+tracktimes = list(np.linspace(0, 50, 1000))
 
 
 def get_output_events(case_name) -> List[str]:
@@ -96,6 +98,8 @@ def get_output_events(case_name) -> List[str]:
         if "term" in case_name
         else ["RELEASE", "TERMINATE"]
         if "mult" in case_name
+        else ["TRACKTIME"]
+        if "trts" in case_name
         else ["ALL"]  # default
     )
 
@@ -170,6 +174,7 @@ def build_prt_sim(name, gwf_ws, prt_ws, mf6):
         track_filerecord=[prt_track_file],
         trackcsv_filerecord=[prt_track_csv_file],
         trackeventrecord=events,
+        tracktimerecord=tracktimes if "trts" in name else None,
     )
 
     # create the flow model interface
@@ -314,12 +319,11 @@ def check_output(idx, test):
     if "TIMESTEP" in events:
         # todo: timestep test case
         pass
-    if "TIMESERIES" in events:
-        # todo: timeseries test case
-        pass
     if "WEAKSINK" in events:
         # todo: weak sink test case
         pass
+    if "TRACKTIME" in events:
+        expected_len += 5324  # hardcoded
     assert len(mf6_pls) == expected_len
 
     # make sure mf6 pathline data have correct
@@ -373,7 +377,9 @@ def all_equal(col, val):
     for ipl, ((iprp, irpt, trelease), pl) in enumerate(mf6_plines):
         pl.plot(
             title="MF6 pathlines",
-            kind="line",
+            marker="o",
+            markersize=2,
+            linestyle="None",
             x="x",
             y="y",
             ax=ax[0],

autotest/prt/test_prt_voronoi01.py → autotest/prt/test_prt_voronoi1.py

@@ -33,8 +33,10 @@
 
 from framework import TestFramework
 
-simname = "prtvor01"
+simname = "prtvor1"
 cases = [f"{simname}l2r", f"{simname}welp", f"{simname}weli"]
+times = [True, False, False]
+tracktimes = list(np.linspace(0, 40000, 100))
 xmin = 0.0
 xmax = 2000.0
 ymin = 0.0
@@ -165,7 +167,7 @@ def build_gwf_sim(name, ws, targets):
     return sim
 
 
-def build_prt_sim(name, gwf_ws, prt_ws, targets):
+def build_prt_sim(idx, name, gwf_ws, prt_ws, targets):
     prt_ws = Path(prt_ws)
     gwf_name = get_model_name(name, "gwf")
     prt_name = get_model_name(name, "prt")
@@ -209,15 +211,15 @@ def build_prt_sim(name, gwf_ws, prt_ws, targets):
     prpdata = [
         # index, (layer, cell index), x, y, z
         (i, (0, vgrid.intersect(p[0], p[1])), p[0], p[1], p[2])
-        for i, p in enumerate(rpts[1:2])  # first release point crashes
+        for i, p in enumerate(rpts[1:])  # first release point crashes
     ]
     prp_track_file = f"{prt_name}.prp.trk"
     prp_track_csv_file = f"{prt_name}.prp.trk.csv"
     flopy.mf6.ModflowPrtprp(
         prt,
         pname="prp1",
         filename=f"{prt_name}_1.prp",
-        nreleasepts=1,  # len(prpdata),
+        nreleasepts=len(prpdata),
         packagedata=prpdata,
         perioddata={0: ["FIRST"]},
         track_filerecord=[prp_track_file],
@@ -232,6 +234,8 @@ def build_prt_sim(name, gwf_ws, prt_ws, targets):
         pname="oc",
         track_filerecord=[prt_track_file],
         trackcsv_filerecord=[prt_track_csv_file],
+        trackeventrecord=["TRACKTIME"] if times[idx] else ["ALL"],
+        tracktimerecord=tracktimes if times[idx] else None,
     )
     gwf_budget_file = gwf_ws / f"{gwf_name}.bud"
     gwf_head_file = gwf_ws / f"{gwf_name}.hds"
@@ -254,7 +258,7 @@ def build_prt_sim(name, gwf_ws, prt_ws, targets):
 def build_models(idx, test):
     gwf_sim = build_gwf_sim(test.name, test.workspace, test.targets)
     prt_sim = build_prt_sim(
-        test.name, test.workspace, test.workspace / "prt", test.targets
+        idx, test.name, test.workspace, test.workspace / "prt", test.targets
     )
     return gwf_sim, prt_sim
 
@@ -282,17 +286,15 @@ def check_output(idx, test):
     )
 
     if "l2r" in name:
-        assert pls.shape == (212, 16)
+        # assert pls.shape == (212, 16)
         assert (pls.z == 0.5).all()  # no z change
         # path should be horizontal from left to right
         assert isclose(min(pls.x), 20, rel_tol=1e-4)
-        assert isclose(max(pls.x), 1979, rel_tol=1e-4)
-        y = 21
-        assert isclose(min(pls.y), y, rel_tol=1e-4)
-        assert isclose(max(pls.y), y, rel_tol=1e-4)
-        assert set(endpts.icell) == {921}
+        assert isclose(max(pls.x), 1980.571, rel_tol=1e-4)
+        assert isclose(min(pls.y), 21, rel_tol=1e-4)
+        assert isclose(max(pls.y), 981, rel_tol=1e-4)
 
-    plot_2d = True
+    plot_2d = False
     if plot_2d:
         # plot in 2d with mpl
         fig = plt.figure(figsize=(16, 10))
@@ -347,13 +349,16 @@ def check_output(idx, test):
             pl.plot(
                 title=title,
                 kind="line",
+                linestyle="--",
+                marker="o",
+                markersize=2,
                 x="x",
                 y="y",
                 ax=ax,
                 legend=False,
                 color="black",
             )
-        # plt.show()
+        plt.show()
         plt.savefig(prt_ws / f"{name}.png")
 
     plot_3d = False

autotest/prt/test_prt_voronoi02.py → autotest/prt/test_prt_voronoi2.py

@@ -32,7 +32,7 @@
 
 from framework import TestFramework
 
-simname = "prtvor02"
+simname = "prtvor2"
 cases = [simname]
 xmin = 0.0
 xmax = 2000.0

autotest/test_gwf_ats02.py

@@ -202,7 +202,7 @@ def make_plot(test):
             label=f"Layer {ilay + 1}",
         )
     plt.legend()
-    plt.show()
+    # plt.show()
 
 
 def check_output(idx, test):