fix(prt): distinguish cell exit from subcell exit events (#2438)

We currently define the PRT OC package's TRACK_EXIT keyword option as enabling reporting when particles exit any grid feature (model, cell or subcell), and report all such events with event code (i.e. reason code) 1. Thus subcell exits in the ternary method appear identical in output to cell transitions, despite the fact that subcells are an implementation detail of the ternary method and occur in the interior of user-defined grid cells.

Distinguish subcell exit events from cell exits. Introduce a new PRT OC option TRACK_SUBFEATURE_EXIT which applies to all grid features below the level of the user-defined discretization (currently only subcell exits in DISV grids, but this could apply later to e.g. subreaches in surface water models, etc). The old TRACK_EXIT option now only applies to grid-scale features and above, i.e. cells (and models once we implement particle exchanges).

While this changes the results reported by PRT under the default settings (see changed test snapshots for instance — fewer path points, now only on cell faces, which nicely matches mp7 when ternary is forced on a structured grid) we think that since subgrid features are not likely to be of interest to most users, this is reasonable to consider a fix, not a breaking change.

Users opting into subfeature transition reporting will still not be able to distinguish cell exits from subcell exits, but that may soon change pending a proposed new output column identifying the domain level at which an event occurs.

Also, move exit event firing routines to their respective method types, following on #2433.

Author: wpbonelli

autotest/__snapshots__/test_prt_disv_vert/test_mf6model[0-prtdisvvert].npy

@@ -311,7 +311,7 @@ def get_ireason_code(output_event):
         if output_event == "RELEASE"
         else (
             1
-            if output_event == "CELLEXIT"
+            if output_event == "EXITFEAT"
             else (
                 2
                 if output_event == "TIMESTEP"

autotest/__snapshots__/test_prt_ternary_methods/test_mf6model[0-prtterbr].npy

@@ -321,18 +321,20 @@ def compare_output(name, mf6_pls, mp7_pls, mp7_eps):
     del mp7_pls["node"]
 
     if "bprp" not in name:
-        # pollock's method should be (nearly) identical
+        # both methods should be (nearly) identical to mp7
         mf6_pls_plck = mf6_pls[mf6_pls.particlegroup == 1]
+        mf6_pls_tern = mf6_pls[mf6_pls.particlegroup == 2]
+
         assert mf6_pls_plck.shape == mp7_pls.shape
         assert np.allclose(mf6_pls_plck, mp7_pls, atol=1e-3)
 
-        # ternary method will have extra path points
-        # due to nudging, so just compare endpoints
-        mf6_eps_tern = mf6_eps[mf6_eps.particlegroup == 2]
-        mf6_eps_tern = mf6_eps_tern[["x", "y", "z", "time"]]
-        mp7_eps = mp7_eps[["x", "y", "z", "time"]]
-        assert mf6_eps_tern.shape == mp7_eps.shape
-        assert np.allclose(mf6_eps_tern, mp7_eps, atol=1e-3)
+        del mf6_pls_tern["particleid"]
+        del mf6_pls_tern["particlegroup"]
+        del mp7_pls["particleid"]
+        del mp7_pls["particlegroup"]
+
+        assert mf6_pls_tern.shape == mp7_pls.shape
+        assert np.allclose(mf6_pls_tern, mp7_pls, atol=1e-3)
 
 
 def check_output(idx, test):
@@ -498,9 +500,9 @@ def plot_output(idx, test):
     for ipl, ((iprp, irpt, trelease), pl) in enumerate(mf6_plines):
         pl.plot(
             title="MF6 pathlines",
-            linestyle="None" if "trst" in name else "--",
+            linestyle="--" if "trst" in name else None,
             marker="o",
-            markersize=2,
+            markersize=4,
             x="x",
             y="y",
             ax=ax[0],
@@ -523,6 +525,8 @@ def plot_output(idx, test):
             ax=ax[1],
             legend=False,
             color=cm.plasma(ipl / len(mp7_plines)),
+            marker="o",
+            markersize=4,
         )
 
     # view/save plot

autotest/__snapshots__/test_prt_ternary_methods/test_mf6model[1-prtterch].npy

@@ -244,7 +244,7 @@ def check_output(idx, test, snapshot):
     endpts = pls.sort_values("t").groupby(["imdl", "iprp", "irpt", "trelease"]).tail(1)
 
     # check pathline shape and endpoints
-    assert pls.shape == (116, 16)
+    assert pls.shape == (57, 16)
     assert endpts.shape == (2, 16)
     assert set(endpts.icell) == {111, 112}
 

autotest/prt_test_utils.py

@@ -299,7 +299,7 @@ def check_output(idx, test, snapshot):
     assert snapshot == endpts.drop("name", axis=1).round(3).to_records(index=False)
 
     if "r2l" in name:
-        assert pls.shape == (164, 16)
+        assert pls.shape == (84, 16)
         assert (pls.z == 0.5).all()
         rtol = 1e-6
         assert isclose(min(pls.x), 0, rel_tol=rtol)
@@ -310,7 +310,7 @@ def check_output(idx, test, snapshot):
         assert isclose(max(pls[pls.irpt == 2].y), 86, rel_tol=rtol)
         assert set(endpts.icell) == {130, 136}
     elif "diag" in name:
-        assert pls.shape == (232, 16)
+        assert pls.shape == (114, 16)
         assert endpts.shape == (4, 16)
         assert set(endpts.icell) == {111, 112}
 

autotest/test_prt_newton_disv1.py

@@ -61,4 +61,9 @@ description = "Add DEPENDENT\\_VARIABLE\\_SCALING option to scale the dependent
 [[items]]
 section = "fixes"
 subsection = "solution"
-description = "The Particle Tracking (PRT) Model previously reported cell exit events when a dry particle stranded in a partially saturated cell was dropped to the water table, as occurs with DRY\\_TRACKING\\_METHOD DROP (the default). An event is no longer reported at this time."
+description = "The Particle Tracking (PRT) Model previously reported cell exit events when a dry particle stranded in a partially saturated cell was dropped to the water table, as occurs with DRY\\_TRACKING\\_METHOD DROP (the default). An event is no longer reported at this time."
+
+[[items]]
+section = "fixes"
+subsection = "solution"
+description = "Previously the Particle Tracking (PRT) Model's generalized DISV tracking method did not distinguish subcell exit events from cell exits, reporting both with the same reason code (1). Since subcell exits occur at a scale below that defined by the model's grid discretization, they are an implementation detail of the generalized tracking method, and will therefore no longer be written by default to particle tracking output. The PRT Output Control (OC) package's TRACK\\_EXIT option will now apply exclusively to features at grid-scale or above. Writing particle transitions between sub-grid-scale features to output is now opt-in with a new keyword option TRACK\\_SUBFEATURE\\_EXIT."

autotest/test_prt_ternary_methods.py

@@ -156,8 +156,16 @@ name track_exit
 type keyword
 reader urword
 optional true
-longname track transitions
-description keyword to indicate that particle tracking output is to be written when a particle exits a feature (a model, cell, or subcell)
+longname track domain exits
+description keyword to indicate that particle tracking output is to be written when a particle exits a grid feature (e.g. model, cell)
+
+block options
+name track_subfeature_exit
+type keyword
+reader urword
+optional true
+longname track cell exits
+description keyword to indicate that particle tracking output is to be written when a particle exits a sub-grid-scale feature. Such features are considered implementation details of the tracking algorithm employed for particular grid types, and reporting for exit events from such features is disabled by default.
 
 block options
 name track_timestep

autotest/test_prt_triangle.py

@@ -121,7 +121,7 @@ \subsection{Particle Track Output}
 
 \begin{description} \itemsep0pt \parskip0pt \parsep0pt
 \item \texttt{0}: particle was released
-\item \texttt{1}: particle exited a cell
+\item \texttt{1}: particle exited a grid feature
 \item \texttt{2}: time step ended
 \item \texttt{3}: particle terminated
 \item \texttt{4}: particle entered a weak sink cell