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.

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

autotest/test_prt_triangle.py

@@ -185,6 +185,7 @@ def build_prt_sim(idx, name, gwf_ws, prt_ws, targets):
             boundnames=True,
             stop_at_weak_sink=True,  # currently required for this problem
             exit_solve_tolerance=1e-5,
+            extend_tracking=True,
         )
     prt_track_file = f"{prtname}.trk"
     prt_track_csv_file = f"{prtname}.trk.csv"

autotest/test_prt_voronoi1.py

@@ -554,6 +554,7 @@ def build_prt_sim(idx, name, gwf_ws, prt_ws, targets, cell_ids):
         boundnames=True,
         stop_at_weak_sink=True,
         exit_solve_tolerance=1e-10,
+        extend_tracking=True,
     )
     prt_track_file = f"{prt_name}.trk"
     prt_track_csv_file = f"{prt_name}.trk.csv"

autotest/test_prt_voronoi2.py

@@ -173,6 +173,7 @@ def build_prt_sim(name, gwf_ws, prt_ws, targets, cell_ids):
         boundnames=True,
         stop_at_weak_sink=True,  # currently required for this problem
         exit_solve_tolerance=1e-10,
+        extend_tracking=True,
     )
     prt_track_file = f"{prt_name}.trk"
     prt_track_csv_file = f"{prt_name}.trk.csv"

autotest/test_prt_weak_sinks.py

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

doc/ReleaseNotes/develop.tex

@@ -24,6 +24,7 @@
 		\item The PRT model could hang upon particle termination due to no (sub)cell exit face. This occurred because a flag signaling particle advance was not set in the proper location.
 		\item Entangled with the previous issue, the ternary method could erroneously terminate a particle and report no exit face (before now this would hang) due to precision error in the exit time/position calculation. This could happen when two conditions are both met: the particle enters the subcell very close to one of its vertices, and flow very nearly parallels one of the subcell's faces. We have encountered similar situations before, solved by nudging the particle a small distance into the interior of the subcell before applying the tracking method. This particular case is resolved by increasing the padding distance from $\sisetup{input-digits = 0123456789\epsilon} \num{\epsilon e+2}$ to $\sisetup{input-digits = 0123456789\epsilon} \num{\epsilon e+5}$, where $\epsilon$ is machine precision.
 		\item For ASCII input files erroneously containing a mix of line endings, MODFLOW would sometimes proceed with unexpected results. The program was corrected to stop with an error message if an input line contained both carriage returns and line feeds.
+		\item 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\_TRACKING is provided. This is meant to provide a stronger guarantee that the program halts under default settings.
 	%	\item xxx
 	%	\item xxx
 	\end{itemize}

doc/mf6io/mf6ivar/dfn/prt-prp.dfn

@@ -42,6 +42,14 @@ optional true
 longname whether to use local z coordinates
 description indicates that ``zrpt'' defines the local z coordinate of the release point within the cell, with value of 0 at the bottom and 1 at the top of the cell.  If the cell is partially saturated at release time, the top of the cell is considered to be the water table elevation (the head in the cell) rather than the top defined by the user.
 
+block options
+name extend_tracking
+type keyword
+reader urword
+optional true
+longname whether to extend tracking beyond the end of the simulation
+description indicates that particles should be tracked beyond the end of the simulation's final time step (using that time step's flows) until particles terminate or reach a specified stop time.  By default, particles are terminated at the end of the simulation's final time step.
+
 block options
 name track_filerecord
 type record track fileout trackfile

doc/mf6io/mf6ivar/md/mf6ivar.md

@@ -1619,6 +1619,7 @@
 | PRT | PRP | OPTIONS | DEV_EXIT_SOLVE_METHOD | INTEGER | the method for iterative solution of particle exit location and time in the generalized Pollock's method.  1 Brent, 2 Chandrupatla.  The default is Brent. |
 | PRT | PRP | OPTIONS | EXIT_SOLVE_TOLERANCE | DOUBLE PRECISION | the convergence tolerance for iterative solution of particle exit location and time in the generalized Pollock's method.  A value of 0.00001 works well for many problems, but the value that strikes the best balance between accuracy and runtime is problem-dependent. |
 | PRT | PRP | OPTIONS | LOCAL_Z | KEYWORD | indicates that ``zrpt'' defines the local z coordinate of the release point within the cell, with value of 0 at the bottom and 1 at the top of the cell.  If the cell is partially saturated at release time, the top of the cell is considered to be the water table elevation (the head in the cell) rather than the top defined by the user. |
+| PRT | PRP | OPTIONS | EXTEND_TRACKING | KEYWORD | indicates that particles should be tracked beyond the end of the simulation's final time step (using that time step's flows) until particles terminate or reach a specified stop time.  By default, particles are terminated at the end of the simulation's final time step. |
 | PRT | PRP | OPTIONS | TRACK | KEYWORD | keyword to specify that record corresponds to a binary track output file.  Each PRP Package may have a distinct binary track output file. |
 | PRT | PRP | OPTIONS | FILEOUT | KEYWORD | keyword to specify that an output filename is expected next. |
 | PRT | PRP | OPTIONS | TRACKFILE | STRING | name of the binary output file to write tracking information. |

doc/mf6io/mf6ivar/tex/prt-prp-desc.tex

@@ -11,6 +11,8 @@
 
 \item \texttt{LOCAL\_Z}---indicates that ``zrpt'' defines the local z coordinate of the release point within the cell, with value of 0 at the bottom and 1 at the top of the cell.  If the cell is partially saturated at release time, the top of the cell is considered to be the water table elevation (the head in the cell) rather than the top defined by the user.
 
+\item \texttt{EXTEND\_TRACKING}---indicates that particles should be tracked beyond the end of the simulation's final time step (using that time step's flows) until particles terminate or reach a specified stop time.  By default, particles are terminated at the end of the simulation's final time step.
+
 \item \texttt{TRACK}---keyword to specify that record corresponds to a binary track output file.  Each PRP Package may have a distinct binary track output file.
 
 \item \texttt{FILEOUT}---keyword to specify that an output filename is expected next.

doc/mf6io/mf6ivar/tex/prt-prp-options.dat

@@ -3,6 +3,7 @@ BEGIN OPTIONS
   [PRINT_INPUT]
   EXIT_SOLVE_TOLERANCE <exit_solve_tolerance>
   [LOCAL_Z]
+  [EXTEND_TRACKING]
   [TRACK FILEOUT <trackfile>]
   [TRACKCSV FILEOUT <trackcsvfile>]
   [STOPTIME <stoptime>]

doc/mf6io/prt/prt.tex

@@ -12,7 +12,7 @@ \subsection{Units of Length and Time}
 The PRT Model formulates the particle trajectory equations without using prescribed length and time units. Any consistent units of length and time can be used when specifying the input data for a simulation. This capability gives a certain amount of freedom to the user, but care must be exercised to avoid mixing units.  The program cannot detect the use of inconsistent units.
 
 \subsection{Time Stepping}
-In \mf time step lengths are controlled by the user and specified in the Temporal Discretization (TDIS) input file.  When the flow model and particle-tracking model run in the same simulation, the time step length specified in TDIS is used for both models.  If the PRT Model runs in a separate simulation, the time discretization may differ.  Instructions for specifying time steps are described in the TDIS section of this user guide; additional information on GWF and PRT configurations are in the Flow Model Interface section.  
+In \mf time step lengths are controlled by the user and specified in the Temporal Discretization (TDIS) input file.  When the flow model and particle-tracking model run in the same simulation, the time step length specified in TDIS is used for both models.  If the PRT Model runs in a separate simulation, the time discretization may differ.  Instructions for specifying time steps are described in the TDIS section of this user guide; additional information on GWF and PRT configurations are in the Flow Model Interface section.  By default, the PRT model will terminate particles at the end of the simulation's final time step; PRT can be configured to extend tracking until particles terminate naturally (i.e. at boundary conditions, sinks, or specified stop times) via Particle Release Point (PRP) package settings.
 
 \subsection{Specifying Cell Face Flows using IFLOWFACE}
 

src/Model/ParticleTracking/prt-prp.f90

@@ -47,6 +47,7 @@ module PrtPrpModule
     integer(I4B), pointer :: itrkcsv => null() !< CSV track file
     integer(I4B), pointer :: irlstls => null() !< release time file
     logical(LGP), pointer :: localz => null() !< compute z coordinates local to the cell
+    integer(I4B), pointer :: extend => null() !< extend tracking beyond simulation's end
     logical(LGP), pointer :: rlsall => null() !< release in all time step
     logical(LGP), pointer :: rlsfirst => null() !< release in first time step
     logical(LGP), pointer :: rlstimelist => null() !< use global release time
@@ -146,6 +147,7 @@ subroutine prp_da(this)
     call mem_deallocate(this%rlsfirst)
     call mem_deallocate(this%rlstimelist)
     call mem_deallocate(this%localz)
+    call mem_deallocate(this%extend)
     call mem_deallocate(this%offset)
     call mem_deallocate(this%stoptime)
     call mem_deallocate(this%stoptraveltime)
@@ -246,6 +248,7 @@ subroutine prp_allocate_scalars(this)
     call mem_allocate(this%rlsfirst, 'RLSFIRST', this%memoryPath)
     call mem_allocate(this%rlstimelist, 'RELEASETIME', this%memoryPath)
     call mem_allocate(this%localz, 'LOCALZ', this%memoryPath)
+    call mem_allocate(this%extend, 'EXTEND', this%memoryPath)
     call mem_allocate(this%offset, 'OFFSET', this%memoryPath)
     call mem_allocate(this%stoptime, 'STOPTIME', this%memoryPath)
     call mem_allocate(this%stoptraveltime, 'STOPTRAVELTIME', this%memoryPath)
@@ -268,6 +271,7 @@ subroutine prp_allocate_scalars(this)
     this%rlsfirst = .false.
     this%rlstimelist = .false.
     this%localz = .false.
+    this%extend = 0
     this%offset = DZERO
     this%stoptime = huge(1d0)
     this%stoptraveltime = huge(1d0)
@@ -471,6 +475,7 @@ subroutine prp_ad(this)
         particle%ifrctrn = this%ifrctrn
         particle%iexmethod = this%iexmethod
         particle%extol = this%extol
+        particle%extend = this%extend
 
         ! -- Persist particle to particle store
         call this%particles%save_particle(particle, np)
@@ -828,6 +833,9 @@ subroutine prp_options(this, option, found)
     case ('LOCAL_Z')
       this%localz = .true.
       found = .true.
+    case ('EXTEND_TRACKING')
+      this%extend = 1
+      found = .true.
     case ('DEV_FORCETERNARY')
       call this%parser%DevOpt()
       this%ifrctrn = 1

src/Model/ParticleTracking/prt.f90

@@ -933,25 +933,30 @@ subroutine prt_solve(this)
           if (particle%istatus == 8) &
             call this%method%save(particle, reason=3) ! reason=3: termination
 
-          ! -- If particle is inactive or not yet to be released, cycle
+          ! If particle is inactive or not yet to be released, cycle
           if (particle%istatus > 1) cycle
 
-          ! -- If particle released this time step, record its initial state
+          ! If particle released this time step, record its initial state
           particle%istatus = 1
           if (particle%trelease >= totimc) &
             call this%method%save(particle, reason=0) ! reason=0: release
 
-          ! -- Maximum time is end of time step unless this is the last
-          !    time step in the simulation, which case it's the particle
-          !    stop time.
-          tmax = totimc + delt
-          if (nper == kper .and. nstp(kper) == kstp) &
+          ! Maximum time is the end of the time step or the particle
+          ! stop time, whichever comes first, unless it's the final
+          ! time step and the extend option is on, in which case
+          ! it's just the particle stop time.
+          if (nper == kper .and. &
+              nstp(kper) == kstp .and. &
+              particle%extend > 0) then
             tmax = particle%tstop
+          else
+            tmax = min(totimc + delt, particle%tstop)
+          end if
 
-          ! -- Get and apply the tracking method
+          ! Get and apply the tracking method
           call this%method%apply(particle, tmax)
 
-          ! -- Update particle storage
+          ! Update particle storage
           call packobj%particles%save_particle(particle, np)
         end do
       end select

src/Solution/ParticleTracker/MethodSubcellPollock.f90

@@ -181,15 +181,16 @@ subroutine track_subcell(this, subcell, particle, tmax)
 
     ! -- Select user tracking times to solve. If this is the first time step
     !    of the simulation, include all times before it begins; if it is the
-    !    last time step include all times after it ends, otherwise the times
-    !    within the current period and time step only.
+    !    last time step include all times after it ends only if the 'extend'
+    !    option is on, otherwise times in this period and time step only.
     call this%tracktimes%try_advance()
     tslice = this%tracktimes%selection
 
     if (all(tslice > 0)) then
       do i = tslice(1), tslice(2)
         t = this%tracktimes%times(i)
-        if (t < particle%ttrack .or. t >= texit .or. t >= tmax) cycle
+        if (t < particle%ttrack) cycle
+        if (t >= texit .or. t >= tmax) exit
         dt = t - t0
         x = new_x(vx, dvxdx, subcell%vx1, subcell%vx2, &
                   dt, initialX, subcell%dx, statusVX == 1)

src/Solution/ParticleTracker/MethodSubcellTernary.f90

@@ -226,16 +226,17 @@ subroutine track_subcell(this, subcell, particle, tmax)
     texit = particle%ttrack + dtexit
     t0 = particle%ttrack
 
-    ! Solve user-specified tracking times within the current stress period and
-    ! time step. If this is the last time step, solve all times after it ends,
-    ! as per MODPATH 7 with stop time option 'extend'.
-    ! todo: reconsider whether this should be default?
+    ! -- Select user tracking times to solve. If this is the first time step
+    !    of the simulation, include all times before it begins; if it is the
+    !    last time step include all times after it ends only if the 'extend'
+    !    option is on, otherwise times in this period and time step only.
     call this%tracktimes%try_advance()
     tslice = this%tracktimes%selection
     if (all(tslice > 0)) then
       do i = tslice(1), tslice(2)
         t = this%tracktimes%times(i)
-        if (t < particle%ttrack .or. t >= texit .or. t >= tmax) cycle
+        if (t < particle%ttrack) cycle
+        if (t >= texit .or. t >= tmax) exit
         dt = t - t0
         call calculate_xyz_position(dt, rxx, rxy, ryx, ryy, sxx, sxy, syy, &
                                     izstatus, x0, y0, az, vzi, vzbot, &