refactor(prt): use enums for particle status and events (#2191)

Use enums for particle status and event codes, following #2188. Also remove an unnecessary/incorrect termination condition left over from #2185 in the base MethodType%check() routine — stationary particles get terminated at the end of the simulation by the tracking loop with timeout status, previously status 5 (no exit face) was getting reported.

src/Model/ModelUtilities/TrackFile.f90

@@ -1,7 +1,7 @@
 module TrackFileModule
   use KindModule, only: DP, I4B, LGP
   use ConstantsModule, only: DZERO, DPIO180
-  use ParticleModule, only: ParticleType
+  use ParticleModule, only: ParticleType, ACTIVE
   use GeomUtilModule, only: transform
 
   implicit none
@@ -18,51 +18,12 @@ module TrackFileModule
   !! state (e.g. tracking status, position) at a particular moment in time.
   !!
   !! Particles have no ID property. Particles can be uniquely identified
-  !! by composite key, i.e. combination of:
+  !! by composite key, i.e. combination of fields:
   !!
   !!   - imdl: originating model ID
   !!   - iprp: originating PRP ID
   !!   - irpt: particle release location ID
   !!   - trelease: particle release 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
-  !!     1: particle transitioned between cells
-  !!     2: current time step ended****
-  !!     3: particle terminated
-  !!     4: particle in weak sink
-  !!     5: user-specified tracking time
-  !!
-  !! Each record has an "istatus" property, which is the tracking status;
-  !! e.g., awaiting release, active, terminated. A particle may terminate
-  !! for several reasons. Status values greater than one imply termination.
-  !! Particle status strictly increases over time, starting at zero:
-  !!
-  !!     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: (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 due to stop time or end of simulation
-  !!
-  !! Comparison to MODPATH 7
-  !! -----------------------
-  !!
-  !! 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
@@ -99,7 +60,7 @@ subroutine save_record(iun, particle, kper, kstp, reason, csv)
 
     ! Set status
     if (particle%istatus .lt. 0) then
-      status = 1
+      status = ACTIVE
     else
       status = particle%istatus
     end if

src/Model/ParticleTracking/prt-prp.f90

@@ -13,7 +13,7 @@ module PrtPrpModule
   use BlockParserModule, only: BlockParserType
   use PrtFmiModule, only: PrtFmiType
   use ParticleModule, only: ParticleType, ParticleStoreType, &
-                            create_particle, allocate_particle_store
+                            create_particle, create_particle_store
   use SimModule, only: count_errors, store_error, store_error_unit, &
                        store_warning
   use SimVariablesModule, only: errmsg, warnmsg
@@ -181,8 +181,8 @@ subroutine prp_da(this)
     call mem_deallocate(this%rptm)
     call mem_deallocate(this%rptname, 'RPTNAME', this%memoryPath)
 
-    ! Deallocate particle store, release time and release step selections
-    call this%particles%deallocate(this%memoryPath)
+    ! Deallocate objects
+    call this%particles%destroy(this%memoryPath)
     call this%schedule%deallocate()
     deallocate (this%particles)
     deallocate (this%schedule)
@@ -211,7 +211,7 @@ subroutine prp_allocate_arrays(this, nodelist, auxvar)
 
     ! Allocate particle store, starting with the number
     ! of release points (arrays resized if/when needed)
-    call allocate_particle_store( &
+    call create_particle_store( &
       this%particles, &
       this%nreleasepoints, &
       this%memoryPath)
@@ -445,6 +445,7 @@ subroutine release(this, ip, trelease)
   end subroutine release
 
   subroutine initialize_particle(this, particle, ip, trelease)
+    use ParticleModule, only: TERM_UNRELEASED
     class(PrtPrpType), intent(inout) :: this !< this instance
     type(ParticleType), pointer, intent(inout) :: particle !< the particle
     integer(I4B), intent(in) :: ip !< particle index
@@ -480,7 +481,7 @@ subroutine initialize_particle(this, particle, ip, trelease)
     end select
     particle%ilay = ilay
     particle%izone = this%rptzone(ic)
-    particle%istatus = 0
+    particle%istatus = 0 ! status 0 until tracking starts
     ! If the cell is inactive, either drape the particle
     ! to the top-most active cell beneath it if drape is
     ! enabled, or else terminate permanently unreleased.
@@ -489,10 +490,13 @@ 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
+          ! negative unreleased status signals to the
+          ! tracking method that we haven't yet saved
+          ! a termination record, it needs to do so.
+          particle%istatus = -1 * TERM_UNRELEASED
         end if
       else
-        particle%istatus = -8
+        particle%istatus = -1 * TERM_UNRELEASED
       end if
     end if
 

src/Model/ParticleTracking/prt.f90

@@ -888,6 +888,7 @@ subroutine prt_solve(this)
     ! modules
     use TdisModule, only: kper, kstp, totimc, delt, endofsimulation
     use PrtPrpModule, only: PrtPrpType
+    use ParticleModule, only: ACTIVE, TERM_UNRELEASED, TERM_TIMEOUT
     ! dummy variables
     class(PrtModelType) :: this
     ! local variables
@@ -935,17 +936,17 @@ subroutine prt_solve(this)
           ! 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
+          if (particle%istatus == (-1 * TERM_UNRELEASED)) then
+            particle%istatus = TERM_UNRELEASED
             call this%method%save(particle, reason=3)
             call packobj%particles%put(particle, np)
           end if
 
           ! Skip terminated particles
-          if (particle%istatus > 1) cycle
+          if (particle%istatus > ACTIVE) cycle
 
           ! If particle was released this time step, record release
-          particle%istatus = 1
+          particle%istatus = ACTIVE
           if (particle%trelease >= totimc) &
             call this%method%save(particle, reason=0)
 
@@ -967,16 +968,17 @@ subroutine prt_solve(this)
           particle%izp = 0
 
           ! If the particle timed out, terminate it.
-          ! "Timeout" means
-          !   - the particle reached its stop time, or
+          ! "Timeout" means it remains active, but
+          !   - it 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. &
+          ! TODO maybe think about changing that?
+          if (particle%istatus <= ACTIVE .and. &
               (particle%ttrack == particle%tstop .or. &
                (endofsimulation .and. particle%iextend == 0))) then
-            particle%istatus = 10
+            particle%istatus = TERM_TIMEOUT
             call this%method%save(particle, reason=3)
           end if
 

src/Solution/ParticleTracker/Method.f90

@@ -194,6 +194,8 @@ end subroutine save
   subroutine check(this, particle, cell_defn, tmax)
     ! modules
     use TdisModule, only: endofsimulation, totimc, totim
+    use ParticleModule, only: TERM_WEAKSINK, TERM_NO_EXITS, &
+                              TERM_STOPZONE, TERM_INACTIVE
     ! dummy
     class(MethodType), intent(inout) :: this
     type(ParticleType), pointer, intent(inout) :: particle
@@ -213,22 +215,22 @@ subroutine check(this, particle, cell_defn, tmax)
     particle%izone = cell_defn%izone
     if (stop_zone) then
       particle%advancing = .false.
-      particle%istatus = 6
+      particle%istatus = TERM_STOPZONE
       call this%save(particle, reason=3)
       return
     end if
 
     if (no_exit_face .and. .not. dry_cell) then
       particle%advancing = .false.
-      particle%istatus = 5
+      particle%istatus = TERM_NO_EXITS
       call this%save(particle, reason=3)
       return
     end if
 
     if (weak_sink) then
       if (particle%istopweaksink > 0) then
         particle%advancing = .false.
-        particle%istatus = 3
+        particle%istatus = TERM_WEAKSINK
         call this%save(particle, reason=3)
         return
       else
@@ -244,7 +246,7 @@ subroutine check(this, particle, cell_defn, tmax)
       else if (particle%idrymeth == 1) then
         ! stop
         particle%advancing = .false.
-        particle%istatus = 7
+        particle%istatus = TERM_INACTIVE
         call this%save(particle, reason=3)
         return
       else if (particle%idrymeth == 2) then
@@ -280,7 +282,7 @@ subroutine check(this, particle, cell_defn, tmax)
 
         ! terminate if last period/step
         if (endofsimulation) then
-          particle%istatus = 5
+          particle%istatus = TERM_NO_EXITS
           particle%ttrack = ttrackmax
           call this%save(particle, reason=3)
           return
@@ -294,7 +296,7 @@ subroutine check(this, particle, cell_defn, tmax)
       else if (particle%idrymeth == 1) then
         ! stop
         particle%advancing = .false.
-        particle%istatus = 7
+        particle%istatus = TERM_INACTIVE
         call this%save(particle, reason=3)
         return
       else if (particle%idrymeth == 2) then
@@ -327,20 +329,12 @@ subroutine check(this, particle, cell_defn, tmax)
             if (t > ttrackmax) ttrackmax = t
           end do
         end if
-
-        ! terminate if last period/step
-        if (endofsimulation) then
-          particle%istatus = 5
-          particle%ttrack = ttrackmax
-          call this%save(particle, reason=3)
-        end if
-        return
       end if
     end if
 
     if (no_exit_face) then
       particle%advancing = .false.
-      particle%istatus = 5
+      particle%istatus = TERM_NO_EXITS
       call this%save(particle, reason=3)
       return
     end if

src/Solution/ParticleTracker/MethodCellPassToBot.f90

@@ -43,6 +43,7 @@ end subroutine deallocate
 
   !> @brief Pass particle vertically and instantaneously to the cell bottom
   subroutine apply_ptb(this, particle, tmax)
+    use ParticleModule, only: TERM_NO_EXITS
     ! dummy
     class(MethodCellPassToBotType), intent(inout) :: this
     type(ParticleType), pointer, intent(inout) :: particle
@@ -67,7 +68,7 @@ subroutine apply_ptb(this, particle, tmax)
     end select
     if (particle%ilay == nlay) then
       particle%advancing = .false.
-      particle%istatus = 5
+      particle%istatus = TERM_NO_EXITS
       call this%save(particle, reason=3)
     end if
 

src/Solution/ParticleTracker/MethodDis.f90

@@ -168,6 +168,7 @@ end subroutine load_dis
   !> @brief Load cell properties into the particle, including
   ! the z coordinate, entry face, and node and layer numbers.
   subroutine load_particle(this, cell, particle)
+    use ParticleModule, only: TERM_BOUNDARY
     ! dummy
     class(MethodDisType), intent(inout) :: this
     type(CellRectType), pointer, intent(inout) :: cell
@@ -210,7 +211,7 @@ subroutine load_particle(this, cell, particle)
       if (ic == particle%icp .and. inface == 7 .and. ilay < particle%ilay) then
         particle%advancing = .false.
         particle%idomain(2) = particle%icp
-        particle%istatus = 2
+        particle%istatus = TERM_BOUNDARY
         particle%izone = particle%izp
         call this%save(particle, reason=3)
         return
@@ -283,6 +284,7 @@ end subroutine update_flowja
 
   !> @brief Pass a particle to the next cell, if there is one
   subroutine pass_dis(this, particle)
+    use ParticleModule, only: TERM_BOUNDARY
     ! dummy
     class(MethodDisType), intent(inout) :: this
     type(ParticleType), pointer, intent(inout) :: particle
@@ -296,7 +298,7 @@ subroutine pass_dis(this, particle)
       ! boundary face, so terminate the particle.
       ! todo AMP: reconsider when multiple models supported
       if (cell%defn%facenbr(particle%iboundary(2)) .eq. 0) then
-        particle%istatus = 2
+        particle%istatus = TERM_BOUNDARY
         particle%advancing = .false.
         call this%save(particle, reason=3)
       else

src/Solution/ParticleTracker/MethodDisv.f90

@@ -141,6 +141,7 @@ end subroutine load_disv
   subroutine load_particle(this, cell, particle)
     ! modules
     use DisvModule, only: DisvType
+    use ParticleModule, only: TERM_BOUNDARY
     ! dummy
     class(MethodDisvType), intent(inout) :: this
     type(CellPolyType), pointer, intent(inout) :: cell
@@ -174,7 +175,7 @@ subroutine load_particle(this, cell, particle)
       if (ic == particle%icp .and. inface == 7 .and. ilay < particle%ilay) then
         particle%advancing = .false.
         particle%idomain(2) = particle%icp
-        particle%istatus = 2
+        particle%istatus = TERM_BOUNDARY
         particle%izone = particle%izp
         call this%save(particle, reason=3)
         return
@@ -223,6 +224,7 @@ end subroutine update_flowja
 
   !> @brief Pass a particle to the next cell, if there is one
   subroutine pass_disv(this, particle)
+    use ParticleModule, only: TERM_BOUNDARY
     ! dummy
     class(MethodDisvType), intent(inout) :: this
     type(ParticleType), pointer, intent(inout) :: particle
@@ -236,7 +238,7 @@ subroutine pass_disv(this, particle)
       ! boundary face, so terminate the particle.
       ! todo AMP: reconsider when multiple models supported
       if (cell%defn%facenbr(particle%iboundary(2)) .eq. 0) then
-        particle%istatus = 2
+        particle%istatus = TERM_BOUNDARY
         particle%advancing = .false.
         call this%save(particle, reason=3)
       else

src/Solution/ParticleTracker/MethodSubcellPollock.f90

@@ -84,6 +84,7 @@ end subroutine apply_msp
   !! this context and for any modifications or errors.
   !<
   subroutine track_subcell(this, subcell, particle, tmax)
+    use ParticleModule, only: ACTIVE, TERM_NO_EXITS_SUB
     ! dummy
     class(MethodSubcellPollockType), intent(inout) :: this
     class(SubcellRectType), intent(in) :: subcell
@@ -135,7 +136,7 @@ subroutine track_subcell(this, subcell, particle, tmax)
     !  Subcell has no exit face, terminate the particle
     !  todo: after initial release, consider ramifications
     if ((statusVX .eq. 3) .and. (statusVY .eq. 3) .and. (statusVZ .eq. 3)) then
-      particle%istatus = 9
+      particle%istatus = TERM_NO_EXITS_SUB
       particle%advancing = .false.
       call this%save(particle, reason=3)
       return
@@ -198,7 +199,7 @@ subroutine track_subcell(this, subcell, particle, tmax)
         particle%y = y * subcell%dy
         particle%z = z * subcell%dz
         particle%ttrack = t
-        particle%istatus = 1
+        particle%istatus = ACTIVE
         call this%save(particle, reason=5)
       end do
     end if
@@ -216,7 +217,7 @@ subroutine track_subcell(this, subcell, particle, tmax)
       z = new_x(vz, dvzdz, subcell%vz1, subcell%vz2, &
                 dt, initialZ, subcell%dz, statusVZ == 1)
       exitFace = 0
-      particle%istatus = 1
+      particle%istatus = ACTIVE
       particle%advancing = .false.
       reason = 2 ! timestep end
     else