fix(prt): fix particle event reporting (MODFLOW-USGS#2185)

Reporting was incorrect for permanently unreleased particles. We want to report them once when they were scheduled for release. MODFLOW-USGS#2177 fixed an issue where they were reported every time step, but the solution was incorrect, only reporting them if the release occurred in the first time step.

Also fix reporting for timed-out particles. The new status code 10 introduced by MODFLOW-USGS#2177 covers not only the simulation ending but the particle reaching its stop time as well.

Update some comments and notes too.

doc/ReleaseNotes/develop.tex

@@ -24,7 +24,8 @@
 	\item GWT, GWE and PRT FMI Packages can now read a GWF Model's binary grid file via a new GWFGRID entry. This allows FMI Packages to perform the same grid equivalence checks as exchanges, which guarantees identical error-checking behavior whether a GWT, GWE or PRT Model is coupled to a GWF Model or running as a post-processor. The GWFGRID file entry is optional but recommended. A future version may make the GWFGRID entry mandatory.
 	\item A regression was recently introduced into the PRT model's generalized tracking method, in which a coordinate transformation was carried out prematurely. This could result in incorrect particle positions in and near quad-refined cells. This bug has been fixed.
 	\item The PRT model previously allowed particles to be released at any time. Release times falling outside the bounds of the simulation's time discretization could produce undefined behavior. Any release times occurring before the simulation begins (i.e. negative times) will now be skipped with a warning message. If EXTEND\_TRACKING is not enabled, release times occurring after the end of the simulation will now be skipped with a warning message as well. If EXTEND\_TRACKING is enabled, release times after the end of the simulation are allowed.
-	\item The PRT Model did not previously report all expected tracking events. In particular, time step end and termination events could go unreported with the DRY\_TRACKING\_METHOD options DROP and STAY (only relevant for Newton formulation models), and termination events were not always reported at the end of the simulation. Reporting has been corrected for the cases identified.
+	\item The PRT Model did not report time step end events correctly with the DRY\_TRACKING\_METHOD options DROP and STAY. This was only relevant for Newton formulation models.
+	\item The PRT Model did not report terminating events for particles still active at the end of the simulation. A corresponding termination event is now reported. This is according to the general expectation that all particle tracks will have exactly one terminating event. A new particle status code (10) has been introduced to indicate particle termination upon reaching a time boundary, either the particle's stop time or the simulation's end time.
 	\item In the generalized particle tracking method, a local Z coordinate could be calculated to fall slightly outside of the unit interval due to numerical imprecision. This could cause the vertical travel time calculation to fail with a floating point exception. Constrain the local Z coordinate to the unit interval to prevent this.
 	\item A profiling module is added for more enhanced performance diagnostics of the program. It can be activated through the PROFILE\_OPTION in the simulation name file.
 	\item Energy decay in the solid phase was added to the EST Package.  This capability is described in the Supplemental Technical Information document, but no actual support was provided in the source code.  Users can now specify distinct zeroth-order decay rates in either the aqueous or solid phases.

src/Model/ModelUtilities/TrackFile.f90

@@ -25,10 +25,10 @@ module TrackFileModule
   !!   - irpt: particle release location ID
   !!   - trelease: particle release time
   !!
-  !! Each record has an "ireason" property, which identifies the cause of
-  !! the record. The user selects 1+ conditions or events for recording.
-  !! Identical records (except "ireason") may be duplicated if multiple
-  !! reporting conditions apply to particles at the same moment in time.
+  !! Each record has an "ireason" property, which identifies the particle
+  !! event. The user selects among particle output events to be reported.
+  !! Records which are identical except for "ireason") may be written if
+  !! multiple reporting conditions apply to the particle at a single time.
   !! Each "ireason" value corresponds to an OC "trackevent" option value:
   !!
   !!     0: particle released
@@ -43,30 +43,26 @@ module TrackFileModule
   !! for several reasons. Status values greater than one imply termination.
   !! Particle status strictly increases over time, starting at zero:
   !!
-  !!     0: pending release*
+  !!     0: pending release (TODO is this necessary? will the user ever see it?)
   !!     1: active
   !!     2: terminated at boundary face
   !!     3: terminated in weak sink cell
-  !!     4: terminated in weak source cell
+  !!     4: (status code unused, see below)
   !!     5: terminated in cell with no exit face
   !!     6: terminated in cell with specified zone number
   !!     7: terminated in inactive cell
-  !!     8: permanently unreleased***
-  !!     9: terminated in subcell with no exit face*****
-  !!     10: terminated upon simulation's end
+  !!     8: permanently unreleased
+  !!     9: terminated in subcell with no exit face
+  !!     10: terminated due to stop time or end of simulation
   !!
-  !! PRT shares the same status enumeration as MODPATH 7. However, some
-  !! don't apply to PRT; for instance, MODPATH 7 distinguishes forwards
-  !! and backwards tracking, but status value 4 is not used by PRT.
+  !! Comparison to MODPATH 7
+  !! -----------------------
   !!
-  !! Notes
-  !! -----
-  !!
-  !!   * is this necessary?
-  !!   ** unnecessary since PRT makes no distinction between forwards/backwards tracking
-  !!   *** e.g., released into an inactive cell, a stop zone cell, or a termination zone
-  !!   **** this may coincide with termination, in which case two events are reported
-  !!   ***** PRT-specific status indicating a particle stopped within a cell subcell
+  !! PRT istatus codes 0-3 and 5-8 correspond directly to MODPATH 7 status codes.
+  !! Status code 4 does not apply to PRT because PRT does not distinguish forwards
+  !! from backwards tracking. Status code 9 provides more specific, subcell-
+  !! level information about a particle that terminated due to no exit face.
+  !! Status code 10 distinguishes particles which have terminated due to timeout.
   !<
   type :: TrackFileType
     private

src/Model/ParticleTracking/prt-prp.f90

@@ -489,10 +489,10 @@ subroutine initialize_particle(this, particle, ip, trelease)
       if (this%idrape > 0) then
         call this%dis%highest_active(ic, this%ibound)
         if (ic == ic_old .or. this%ibound(ic) == 0) then
-          particle%istatus = 8
+          particle%istatus = -8
         end if
       else
-        particle%istatus = 8
+        particle%istatus = -8
       end if
     end if
 

src/Model/ParticleTracking/prt.f90

@@ -928,18 +928,23 @@ subroutine prt_solve(this)
           ! Load particle from storage
           call packobj%particles%get(particle, this%id, iprp, np)
 
-          ! If particle is permanently unreleased, record
-          ! unreleased state if first time step and cycle
-          if (particle%istatus == 8) then
-            if (kper == 1 .and. kstp == 1) &
-              call this%method%save(particle, reason=3)
-            cycle
+          ! If particle is permanently unreleased, cycle.
+          ! Save a termination record if we haven't yet.
+          ! TODO: when we have generic dynamic vectors,
+          ! consider terminating permanently unreleased
+          ! in PRP instead of here. For now, status -8
+          ! indicates the permanently unreleased event
+          ! is not yet recorded, status 8 it has been.
+          if (particle%istatus == -8) then
+            particle%istatus = 8
+            call this%method%save(particle, reason=3)
+            call packobj%particles%put(particle, np)
           end if
 
-          ! If particle is inactive or not yet to be released, cycle
+          ! Skip terminated particles
           if (particle%istatus > 1) cycle
 
-          ! If particle released this time step, record its initial state
+          ! If particle was released this time step, record release
           particle%istatus = 1
           if (particle%trelease >= totimc) &
             call this%method%save(particle, reason=0)
@@ -961,10 +966,16 @@ subroutine prt_solve(this)
           particle%icp = 0
           particle%izp = 0
 
-          ! Terminate particles still active at end of simulation
-          if (endofsimulation .and. &
-              particle%iextend == 0 .and. &
-              particle%istatus < 2) then
+          ! If the particle timed out, terminate it.
+          ! "Timeout" means
+          !   - the particle reached its stop time, or
+          !   - the simulation is over and no extending.
+          ! We can't detect timeout within the tracking
+          ! method because the method just receives the
+          ! maximum time with no context on what it is.
+          if (particle%istatus < 2 .and. &
+              (particle%ttrack == particle%tstop .or. &
+               (endofsimulation .and. particle%iextend == 0))) then
             particle%istatus = 10
             call this%method%save(particle, reason=3)
           end if