fix(prt): make extended tracking opt-in (#1888)

Previously the PRT model's default behavior was to track particles until termination, as with MODPATH 7 stop time option 2 (extend). Under extended tracking, the program may not halt if particles enter a cycle in the flow system. This PR changes the default to the equivalent of MP7 stop time option 1 (final), terminating at simulation end unless a new Particle Release Point (PRP) package keyword option EXTEND is provided. This is meant to provide a stronger guarantee that the program halts under default settings.

The extend option is stored as an integer rather than logical in case we need to enumerate more alternatives in future.

This PR also avoids unnecessary iterations in MethodSubcellPollock/Ternary, and changes the default visibility of ParticleStore members to private, as per style conventions.

Author: wpbonelli

autotest/test_prt_budget.py

@@ -89,6 +89,7 @@ def build_prt_sim(name, gwf_ws, prt_ws, mf6):
         stop_at_weak_sink=False,
         boundnames=True,
         exit_solve_tolerance=DEFAULT_EXIT_SOLVE_TOL,
+        extend_tracking=True,
     )
 
     # create output control package

autotest/test_prt_disv1.py

@@ -236,6 +236,7 @@ def build_prt_sim(idx, gwf_ws, prt_ws, mf6):
             print_input=True,
             dev_forceternary=i == 1,
             exit_solve_tolerance=1e-10,
+            extend_tracking=True,
         )
 
     # create output control package

autotest/test_prt_drape.py

@@ -190,6 +190,7 @@ def build_prt_sim(name, gwf_ws, prt_ws, mf6):
         trackcsv_filerecord=[prp_track_csv_file],
         drape="drp" in name,
         exit_solve_tolerance=DEFAULT_EXIT_SOLVE_TOL,
+        extend_tracking=True,
     )
 
     # create output control package

autotest/test_prt_exg.py

@@ -79,6 +79,7 @@ def build_mf6_sim(idx, test):
         perioddata={0: ["FIRST"]},
         boundnames="bnms" in name,
         exit_solve_tolerance=DEFAULT_EXIT_SOLVE_TOL,
+        extend_tracking=True,
     )
 
     # create output control package

autotest/test_prt_fmi.py

@@ -48,8 +48,8 @@
     DEFAULT_EXIT_SOLVE_TOL,
 )
 
-simname = "prtfmi01"
-cases = [simname, f"{simname}saws", f"{simname}bprp"]
+simname = "prtfmi"
+cases = [simname, f"{simname}saws", f"{simname}bprp", f"{simname}noext"]
 
 
 def build_prt_sim(name, gwf_ws, prt_ws, mf6):
@@ -123,6 +123,7 @@ def build_prt_sim(name, gwf_ws, prt_ws, mf6):
         stop_at_weak_sink="saws" in prt_name,
         boundnames=True,
         exit_solve_tolerance=DEFAULT_EXIT_SOLVE_TOL,
+        extend_tracking="noext" not in prt_name,
     )
 
     # create output control package
@@ -346,28 +347,32 @@ def check_output(idx, test):
         plt.show()
         plt.savefig(gwf_ws / f"test_{simname}.png")
 
-    # convert mf6 pathlines to mp7 format
-    mf6_pls = to_mp7_pathlines(mf6_pls)
-
-    # sort both dataframes by particleid and time
-    mf6_pls.sort_values(by=["particleid", "time"], inplace=True)
-    mp7_pls.sort_values(by=["particleid", "time"], inplace=True)
-
-    # drop columns for which there is no direct correspondence between mf6 and mp7
-    del mf6_pls["sequencenumber"]
-    del mf6_pls["particleidloc"]
-    del mf6_pls["xloc"]
-    del mf6_pls["yloc"]
-    del mf6_pls["zloc"]
-    del mp7_pls["sequencenumber"]
-    del mp7_pls["particleidloc"]
-    del mp7_pls["xloc"]
-    del mp7_pls["yloc"]
-    del mp7_pls["zloc"]
-
-    # compare mf6 / mp7 pathline data
-    assert mf6_pls.shape == mp7_pls.shape
-    assert np.allclose(mf6_pls, mp7_pls, atol=1e-3)
+    if "noext" in name:
+        # maximum tracking time should be simulation stop time
+        assert mf6_pls.t.max() == FlopyReadmeCase.nper * FlopyReadmeCase.perlen
+    else:
+        # convert mf6 pathlines to mp7 format
+        mf6_pls = to_mp7_pathlines(mf6_pls)
+
+        # sort both dataframes by particleid and time
+        mf6_pls.sort_values(by=["particleid", "time"], inplace=True)
+        mp7_pls.sort_values(by=["particleid", "time"], inplace=True)
+
+        # drop columns for which there is no direct correspondence between mf6 and mp7
+        del mf6_pls["sequencenumber"]
+        del mf6_pls["particleidloc"]
+        del mf6_pls["xloc"]
+        del mf6_pls["yloc"]
+        del mf6_pls["zloc"]
+        del mp7_pls["sequencenumber"]
+        del mp7_pls["particleidloc"]
+        del mp7_pls["xloc"]
+        del mp7_pls["yloc"]
+        del mp7_pls["zloc"]
+
+        # compare mf6 / mp7 pathline data
+        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/test_prt_quad_refinement.py

@@ -154,6 +154,7 @@ def build_mf6_sim(idx, test, **kwargs):
         perioddata={0: ["FIRST"]},
         exit_solve_tolerance=DEFAULT_EXIT_SOLVE_TOL,
         dev_forceternary=tracking_method == "ternary",
+        extend_tracking=True,
     )
     prt_track_file = f"{prt_name}.trk"
     prt_track_csv_file = f"{prt_name}.trk.csv"

autotest/test_prt_release_timing.py

@@ -177,6 +177,7 @@ def build_prt_sim(name, gwf_ws, prt_ws, mf6, fraction=None):
         ),
         print_input=True,
         exit_solve_tolerance=DEFAULT_EXIT_SOLVE_TOL,
+        extend_tracking=True,
     )
 
     # create output control package

autotest/test_prt_stop_zones.py

@@ -131,6 +131,7 @@ def build_prt_sim(name, gwf_ws, prt_ws, mf6):
         perioddata={0: ["FIRST"]},
         istopzone=1,
         exit_solve_tolerance=DEFAULT_EXIT_SOLVE_TOL,
+        extend_tracking=True,
     )
 
     # create output control package

autotest/test_prt_ternary_methods.py

@@ -100,6 +100,7 @@ def build_prt_sim(
         stop_at_weak_sink=True,  # currently required for this problem
         dev_exit_solve_method=methods[idx],
         exit_solve_tolerance=exit_solve_tolerance,
+        extend_tracking=True,
     )
     prt_track_file = f"{prtname}.trk"
     prt_track_csv_file = f"{prtname}.trk.csv"

autotest/test_prt_track_events.py

@@ -143,6 +143,7 @@ def build_prt_sim(name, gwf_ws, prt_ws, mf6):
             packagedata=releasepts_prt[grp],
             perioddata={0: ["FIRST"]},
             exit_solve_tolerance=DEFAULT_EXIT_SOLVE_TOL,
+            extend_tracking=True,
         )
         for grp in ["a", "b"]
     ]