feat(prt): add developer option for cycle detection (#2435)

Cyclic particle pathlines may be legitimate over multiple time steps, but within a given application of the tracking method i.e. in a single time step they so far reliably indicate a programming error.

Add a developer option DEV_CYCLE_DETECTION_WINDOW to PRT PRP enabling cycle detection. Disabled (set to 0) by default. When enabled, look at the last DEV_CYCLE_DETECTION_WINDOW events and terminate the simulation with error if a cycle is found. A cycle is defined as a duplicate cell exit (i.e. through the same face) within a single application of the tracking method.

For a speed improvement we could keep a HashTable alongside the list of recent events, where the key is the event string, then just a single lookup instead of iterating over the full window. But since this is meant for infrequent use as a developer debugging tool not a user option, and the window will be small, I will punt performance optimization for now.

Consequentially, this PR also makes it the responsibility of the event firer to deallocate the event when it's no longer needed. The event firing routines generally now deallocate the event after dispatch/handling. When cycle detection is off, the cell exit event firing routine does the same. When cycle detection is on, events cannot be deallocated until they are popped from the history, so they are cleaned up at that time by the list node's value destruction mechanism.

Author: wpbonelli

doc/ReleaseNotes/develop.toml

@@ -51,7 +51,7 @@ description = "Fixed a memory access exception that can occur when writing conve
 [[items]]
 section = "features"
 subsection = "solution"
-description = "Added an optional new COORDINATE\\_CHECK\\_METHOD option to the PRT model's Particle Release Package, accepting values 'EAGER' (default) or 'NONE'. 'NONE' disables release-time verification that release point coordinates and release point cell IDs match, which can significantly reduce runtime when releasing many particles."
+description = "Added an optional new COORDINATE\\_CHECK\\_METHOD option to the PRT model's Particle Release Package, accepting values EAGER (default) or NONE. NONE disables release-time verification that release point coordinates and release point cell IDs match, which can significantly reduce runtime when releasing many particles."
 
 [[items]]
 section = "features"
@@ -72,3 +72,8 @@ description = "Previously the Particle Tracking (PRT) Model's generalized DISV t
 section = "fixes"
 subsection = "basic"
 description = "Previously the NOCHECK option was used to directly set the internal isimcheck variable. When NOCHECK is specified in the simulation name file isimcheck should be set to 0, otherwise isimcheck should be 1. The code has been fixed to implement the intended behaviour."
+
+[[items]]
+section = "features"
+subsection = "solution"
+description = "Added a developer option DEV\\_CYCLE\\_DETECTION\\_WINDOW to the PRT Model's Particle Release (PRP) Package. When enabled, this feature maintains a history of recent cell transitions, the size of which is configurable. If a cycle is detected, the program will terminate with an error. A cycle is identified as a duplicate particle event within a single application of the tracking loop, respective to certain particle and cell properties, e.g., exiting through the same cell face twice. This option defaults to 0, disabling cycle detection by default."

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

@@ -293,6 +293,15 @@ longname coordinate checking method
 description approach for verifying that release point coordinates are in the cell with the specified ID. By default, release point coordinates are checked at release time.
 default eager
 
+block options
+name dev_cycle_detection_window
+type integer
+reader urword
+optional true
+mf6internal icycwin
+longname cycle detection window size
+description integer value defining the size of the window (number of consecutive exit events) used for cycle detection. Defaults to 0, disabling cycle detection. With detection enabled, particle pathlines with duplicate cell exit events (i.e., exiting the same cell through the same face twice) will cause the program to terminate with an error. A larger window size provides more robust cycle detection at the cost of more runtime operations per cell exit.
+
 # --------------------- prt prp dimensions ---------------------
 
 block dimensions

make/makefile

@@ -315,7 +315,6 @@ $(OBJDIR)/ParticleEvent.o \
 $(OBJDIR)/sort.o \
 $(OBJDIR)/FlowModelInterface.o \
 $(OBJDIR)/Cell.o \
-$(OBJDIR)/FeatExitEvent.o \
 $(OBJDIR)/WeakSinkEvent.o \
 $(OBJDIR)/UserTimeEvent.o \
 $(OBJDIR)/TimeStepEvent.o \
@@ -325,20 +324,21 @@ $(OBJDIR)/Subcell.o \
 $(OBJDIR)/ReleaseEvent.o \
 $(OBJDIR)/prt-fmi.o \
 $(OBJDIR)/ParticleEvents.o \
-$(OBJDIR)/SubcellExitEvent.o \
+$(OBJDIR)/FeatExitEvent.o \
 $(OBJDIR)/Method.o \
+$(OBJDIR)/CellExitEvent.o \
+$(OBJDIR)/SubcellExitEvent.o \
+$(OBJDIR)/MethodCell.o \
 $(OBJDIR)/TernarySolveUtils.o \
 $(OBJDIR)/SubcellTri.o \
 $(OBJDIR)/MethodSubcell.o \
 $(OBJDIR)/SubcellRect.o \
 $(OBJDIR)/MethodSubcellTernary.o \
 $(OBJDIR)/MethodSubcellPollock.o \
-$(OBJDIR)/CellExitEvent.o \
 $(OBJDIR)/TimeStepSelect.o \
 $(OBJDIR)/LinearAlgebraUtils.o \
 $(OBJDIR)/SfrCrossSectionUtils.o \
 $(OBJDIR)/MethodSubcellPool.o \
-$(OBJDIR)/MethodCell.o \
 $(OBJDIR)/CellPoly.o \
 $(OBJDIR)/CellRectQuad.o \
 $(OBJDIR)/CellRect.o \

src/Idm/prt-prpidm.f90

@@ -41,6 +41,7 @@ module PrtPrpInputModule
     logical :: rttol = .false.
     logical :: rtfreq = .false.
     logical :: ichkmeth = .false.
+    logical :: icycwin = .false.
     logical :: nreleasepts = .false.
     logical :: nreleasetimes = .false.
     logical :: irptno = .false.
@@ -589,6 +590,24 @@ module PrtPrpInputModule
     .false. & ! timeseries
     )
 
+  type(InputParamDefinitionType), parameter :: &
+    prtprp_icycwin = InputParamDefinitionType &
+    ( &
+    'PRT', & ! component
+    'PRP', & ! subcomponent
+    'OPTIONS', & ! block
+    'DEV_CYCLE_DETECTION_WINDOW', & ! tag name
+    'ICYCWIN', & ! fortran variable
+    'INTEGER', & ! type
+    '', & ! shape
+    'cycle detection window size', & ! longname
+    .false., & ! required
+    .false., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
   type(InputParamDefinitionType), parameter :: &
     prtprp_nreleasepts = InputParamDefinitionType &
     ( &
@@ -909,6 +928,7 @@ module PrtPrpInputModule
     prtprp_rttol, &
     prtprp_rtfreq, &
     prtprp_ichkmeth, &
+    prtprp_icycwin, &
     prtprp_nreleasepts, &
     prtprp_nreleasetimes, &
     prtprp_irptno, &

src/Model/ParticleTracking/prt-prp.f90

@@ -54,6 +54,7 @@ module PrtPrpModule
     integer(I4B), pointer :: ifrctrn => null() !< force ternary solution for quad grids
     integer(I4B), pointer :: iexmeth => null() !< method for iterative solution of particle exit location and time in generalized Pollock's method
     integer(I4B), pointer :: ichkmeth => null() !< method for checking particle release coordinates are in the specified cells, 0 = none, 1 = eager
+    integer(I4B), pointer :: icycwin => null() !< cycle detection window size
     real(DP), pointer :: extol => null() !< tolerance for iterative solution of particle exit location and time in generalized Pollock's method
     real(DP), pointer :: rttol => null() !< tolerance for coincident particle release times
     real(DP), pointer :: rtfreq => null() !< frequency for regularly spaced release times
@@ -168,6 +169,7 @@ subroutine prp_da(this)
     call mem_deallocate(this%ifrctrn)
     call mem_deallocate(this%iexmeth)
     call mem_deallocate(this%ichkmeth)
+    call mem_deallocate(this%icycwin)
     call mem_deallocate(this%extol)
     call mem_deallocate(this%rttol)
     call mem_deallocate(this%rtfreq)
@@ -259,6 +261,7 @@ subroutine prp_allocate_scalars(this)
     call mem_allocate(this%ifrctrn, 'IFRCTRN', this%memoryPath)
     call mem_allocate(this%iexmeth, 'IEXMETH', this%memoryPath)
     call mem_allocate(this%ichkmeth, 'ICHKMETH', this%memoryPath)
+    call mem_allocate(this%icycwin, 'ICYCWIN', this%memoryPath)
     call mem_allocate(this%extol, 'EXTOL', this%memoryPath)
     call mem_allocate(this%rttol, 'RTTOL', this%memoryPath)
     call mem_allocate(this%rtfreq, 'RTFREQ', this%memoryPath)
@@ -283,6 +286,7 @@ subroutine prp_allocate_scalars(this)
     this%ifrctrn = 0
     this%iexmeth = 0
     this%ichkmeth = 1
+    this%icycwin = 0
     this%extol = DEM5
     this%rttol = DSAME * DEP9
     this%rtfreq = DZERO
@@ -540,6 +544,7 @@ subroutine initialize_particle(this, particle, ip, trelease)
     particle%ifrctrn = this%ifrctrn
     particle%iexmeth = this%iexmeth
     particle%iextend = this%iextend
+    particle%icycwin = this%icycwin
     particle%extol = this%extol
   end subroutine initialize_particle
 
@@ -712,6 +717,7 @@ subroutine prp_options(this)
                        found%iexmeth)
     call mem_set_value(this%ichkmeth, 'ICHKMETH', this%input_mempath, &
                        coorcheck_method, found%ichkmeth)
+    call mem_set_value(this%icycwin, 'ICYCWIN', this%input_mempath, found%icycwin)
 
     ! update internal state and validate input
     if (found%idrymeth) then
@@ -757,6 +763,11 @@ subroutine prp_options(this)
       end if
     end if
 
+    if (found%icycwin) then
+      if (this%icycwin < 0) &
+        call store_error('CYCLE_DETECTION_WINDOW MUST BE NON-NEGATIVE')
+    end if
+
     ! fileout options
     if (found%trackfile) then
       this%itrkout = getunit()

src/Model/ParticleTracking/prt.f90

@@ -989,6 +989,7 @@ subroutine prt_solve(this)
         end do
       end select
     end do
+    call particle%destroy()
     deallocate (particle)
   end subroutine prt_solve
 

src/Solution/ParticleTracker/Method/Method.f90

@@ -185,6 +185,7 @@ subroutine release(this, particle)
 
     allocate (ReleaseEventType :: event)
     call this%events%dispatch(particle, event)
+    deallocate (event)
   end subroutine release
 
   !> @brief Particle terminates.
@@ -198,6 +199,7 @@ subroutine terminate(this, particle, status)
     if (present(status)) particle%istatus = status
     allocate (TerminationEventType :: event)
     call this%events%dispatch(particle, event)
+    deallocate (event)
   end subroutine terminate
 
   !> @brief Time step ends.
@@ -208,6 +210,7 @@ subroutine timestep(this, particle)
 
     allocate (TimeStepEventType :: event)
     call this%events%dispatch(particle, event)
+    deallocate (event)
   end subroutine timestep
 
   !> @brief Particle leaves a weak sink.
@@ -218,6 +221,7 @@ subroutine weaksink(this, particle)
 
     allocate (WeakSinkEventType :: event)
     call this%events%dispatch(particle, event)
+    deallocate (event)
   end subroutine weaksink
 
   !> @brief User-defined tracking time occurs.
@@ -228,6 +232,7 @@ subroutine usertime(this, particle)
 
     allocate (UserTimeEventType :: event)
     call this%events%dispatch(particle, event)
+    deallocate (event)
   end subroutine usertime
 
 end module MethodModule

src/Solution/ParticleTracker/Method/MethodCell.f90

@@ -1,12 +1,14 @@
 module MethodCellModule
 
   use KindModule, only: DP, I4B, LGP
+  use ErrorUtilModule, only: pstop
   use ConstantsModule, only: DONE, DZERO
   use MethodModule, only: MethodType
-  use ParticleModule, only: ParticleType
+  use ParticleModule, only: ParticleType, TERM_NO_EXITS, TERM_BOUNDARY
   use ParticleEventModule, only: ParticleEventType
   use CellExitEventModule, only: CellExitEventType
   use CellDefnModule, only: CellDefnType
+  use IteratorModule, only: IteratorType
   implicit none
 
   private
@@ -16,6 +18,8 @@ module MethodCellModule
   contains
     procedure, public :: assess
     procedure, public :: cellexit
+    procedure, public :: forms_cycle
+    procedure, public :: store_event
   end type MethodCellType
 
 contains
@@ -32,8 +36,7 @@ subroutine assess(this, particle, cell_defn, tmax)
     use TdisModule, only: endofsimulation, totimc, totim
     use ParticleModule, only: TERM_WEAKSINK, TERM_NO_EXITS, &
                               TERM_STOPZONE, TERM_INACTIVE
-    use ParticleEventModule, only: FEATEXIT, TERMINATE, &
-                                   TIMESTEP, WEAKSINK, USERTIME
+    use ParticleEventModule, only: TERMINATE, TIMESTEP, WEAKSINK, USERTIME
     ! dummy
     class(MethodCellType), intent(inout) :: this
     type(ParticleType), pointer, intent(inout) :: particle
@@ -172,13 +175,87 @@ subroutine cellexit(this, particle)
     class(MethodCellType), intent(inout) :: this
     type(ParticleType), pointer, intent(inout) :: particle
     class(ParticleEventType), pointer :: event
+    ! local
+    integer(I4B) :: i, nhist
+    class(*), pointer :: prev
 
     allocate (CellExitEventType :: event)
     select type (event)
     type is (CellExitEventType)
       event%exit_face = particle%iboundary(2)
     end select
     call this%events%dispatch(particle, event)
+    if (particle%icycwin == 0) then
+      deallocate (event)
+      return
+    end if
+    if (this%forms_cycle(particle, event)) then
+      ! print event history
+      print *, "Cyclic pathline detected"
+      nhist = particle%history%Count()
+      do i = 1, nhist
+        prev => particle%history%GetItem(i)
+        select type (prev)
+        class is (ParticleEventType)
+          print *, "Back ", nhist - i + 1, ": ", prev%get_text()
+        end select
+      end do
+      print *, "Current :", event%get_text()
+      call pstop(1, 'Cyclic pathline detected, aborting')
+    else
+      call this%store_event(particle, event)
+    end if
   end subroutine cellexit
 
+  !> @brief Check if the event forms a cycle in the particle path.
+  function forms_cycle(this, particle, event) result(found_cycle)
+    ! dummy
+    class(MethodCellType), intent(inout) :: this
+    type(ParticleType), pointer, intent(inout) :: particle
+    class(ParticleEventType), pointer, intent(in) :: event
+    ! local
+    class(IteratorType), allocatable :: itr
+    logical(LGP) :: found_cycle
+
+    found_cycle = .false.
+    select type (event)
+    type is (CellExitEventType)
+      itr = particle%history%Iterator()
+      do while (itr%has_next())
+        call itr%next()
+        select type (prev => itr%value())
+        class is (CellExitEventType)
+          if (event%icu == prev%icu .and. &
+              event%ilay == prev%ilay .and. &
+              event%izone == prev%izone .and. &
+              event%exit_face == prev%exit_face .and. &
+              event%exit_face /= 0) then
+            found_cycle = .true.
+            exit
+          end if
+        end select
+      end do
+    end select
+  end function forms_cycle
+
+  !> @brief Save the event in the particle's history.
+  !! Acts like a queue, the oldest event is removed
+  !! when the event count exceeds the maximum size.
+  subroutine store_event(this, particle, event)
+    ! dummy
+    class(MethodCellType), intent(inout) :: this
+    type(ParticleType), pointer, intent(inout) :: particle
+    class(ParticleEventType), pointer, intent(in) :: event
+    ! local
+    class(*), pointer :: p
+
+    select type (event)
+    type is (CellExitEventType)
+      p => event
+      call particle%history%Add(p)
+      if (particle%history%Count() > particle%icycwin) &
+        call particle%history%RemoveNode(1, .true.)
+    end select
+  end subroutine store_event
+
 end module MethodCellModule

src/Solution/ParticleTracker/Method/MethodCellPassToBot.f90

@@ -64,10 +64,11 @@ subroutine apply_ptb(this, particle, tmax)
     type is (DisvType)
       nlay = dis%nlay
     end select
-    if (particle%ilay == nlay) &
+    if (particle%ilay == nlay) then
       call this%terminate(particle, status=TERM_NO_EXITS)
-
-    call this%cellexit(particle)
+    else
+      call this%cellexit(particle)
+    end if
   end subroutine apply_ptb
 
 end module MethodCellPassToBotModule

src/Solution/ParticleTracker/Method/MethodSubcell.f90

@@ -1,15 +1,15 @@
 module MethodSubcellModule
   use KindModule, only: DP, I4B
-  use MethodModule, only: MethodType
+  use MethodCellModule, only: MethodCellType
   use ParticleModule, only: ParticleType
   use CellDefnModule, only: CellDefnType
   use ParticleEventModule, only: ParticleEventType
-  use SubcellExitEventModule, only: SubcellExitEventType
+  use SubCellExitEventModule, only: SubCellExitEventType
 
   private
   public :: MethodSubcellType
 
-  type, abstract, extends(MethodType) :: MethodSubcellType
+  type, abstract, extends(MethodCellType) :: MethodSubcellType
   contains
     procedure, public :: assess
     procedure, public :: subcellexit
@@ -32,9 +32,9 @@ subroutine subcellexit(this, particle)
     type(ParticleType), pointer, intent(inout) :: particle
     class(ParticleEventType), pointer :: event
 
-    allocate (SubcellExitEventType :: event)
+    allocate (SubCellExitEventType :: event)
     select type (event)
-    type is (SubcellExitEventType)
+    type is (SubCellExitEventType)
       event%isc = particle%idomain(3)
       event%exit_face = particle%iboundary(3)
     end select