fix(prt): prevent top exits from partially saturated cells (#2478)

Cyclic particle pathlines were possible in Newton models with DRY_TRACKING_METHOD DROP (the default) due to numerically insignificant upward flows in partially saturated cells. This occurred because PRT did not distinguish the water table from the cell top, so a particle at the water table in a partially saturated cell could obtain an exit solution through the cell's top face, jump to the dry cell above, then drop back down, and so on.

Prevent particles from exiting through the top face of partially saturated cells. If a particle obtains an exit solution upwards through the water table, its behavior now depends on whether the cell top is an assigned boundary face, and the DRY_TRACKING_METHOD setting. If the cell top is an assigned boundary, or if DRY_TRACKING_METHOD is STOP, the particle will terminate. Otherwise it will stay in that location. (Future work may expand on this behavior; this is a minimal bug fix.)

The issue is reproducible on the Keating example, so I added no new tests here. I added a snapshot comparison to Keating to check the new behavior.

This fix allows removing the old ad hoc trap for cycles between vertical faces. Also some other misc cleanup.

Author: wpbonelli

doc/ReleaseNotes/develop.toml

@@ -97,3 +97,8 @@ description = "Fixes a memory exception that has occurred when running parallel
 section = "fixes"
 subsection = "solution"
 description = "The PRT Model has previously been aware of boundary cell flows but not boundary faces. When an internal boundary face is assigned (an IFLOWFACE adjacent to another active cell), this can lead to cycles (infinite loops) in the particle's pathline, as the two adjacent cells fail to agree on the flow through their shared face. The PRT Model will now terminate particles at assigned boundary faces, whether they are internal or external, if the net boundary flow through the face is out of the cell."
+
+[[items]]
+section = "fixes"
+subsection = "solution"
+description = "Cyclic particle pathlines were possible in Newton models with DRY\\_TRACKING\\_METHOD DROP (the default) in the presence of numerically insignificant upward flows in partially saturated cells. This occurred because PRT did not distinguish the water table from the cell top, so a particle at the water table in a partially saturated cell could obtain an exit solution through the cell's top face, jump to the dry cell above, then drop back down, and so on. Particles are now prevented from exiting through the top face of partially saturated cells. If a particle reaches the water table, its behavior now depends on a combination of DRY\\_TRACKING\\_METHOD and whether the cell top is an assigned boundary face. If the cell top is an assigned boundary face, or if DRY\\_TRACKING\\_METHOD is STOP, the particle will terminate. Otherwise, it will stay there until tracking resumes the next time step."

src/Model/ParticleTracking/prt.f90

@@ -968,10 +968,6 @@ subroutine prt_solve(this)
           end if
           ! Apply the tracking method until the maximum time.
           call this%method%apply(particle, tmax)
-          ! Reset cell/zone one-backs, used for cycle detection.
-          ! TODO can remove when we have better cycle detection
-          particle%icp = 0
-          particle%izp = 0
           ! If the particle timed out, terminate it.
           ! "Timed out" means it's still active but
           !   - it reached its stop time, or

src/Solution/ParticleTracker/Method/ExitSolution.f90

@@ -8,7 +8,7 @@ module ExitSolutionModule
   !> @brief Exit status codes
   enum, bind(C)
     enumerator :: OK_EXIT = 0 !< exit found using velocity interpolation
-    ! below only used for linear interp
+    ! below only used for linear solution
     enumerator :: OK_EXIT_CONSTANT = 1 !< exit found, constant velocity
     enumerator :: NO_EXIT_STATIONARY = 2 !< no exit, zero velocity
     enumerator :: NO_EXIT_NO_OUTFLOW = 3 !< no exit, no outflow

src/Solution/ParticleTracker/Method/MethodCell.f90

@@ -29,23 +29,41 @@ module MethodCellModule
   subroutine try_pass(this, particle, nextlevel, advancing)
     class(MethodCellType), intent(inout) :: this
     type(ParticleType), pointer, intent(inout) :: particle
-    integer(I4B) :: nextlevel
-    logical(LGP) :: advancing
-
-    if (particle%advancing) then
-      ! if still advancing, pass to the next subdomain.
-      ! if that puts us on a boundary, then we're done.
-      ! raise a cell exit event.
-      call this%pass(particle)
-      if (particle%iboundary(nextlevel - 1) .ne. 0) then
-        advancing = .false.
-        call this%cellexit(particle)
-      end if
-    else
-      ! otherwise we're already done so
-      ! reset the domain boundary value.
+    integer(I4B) :: nextlevel, ic, iface
+    logical(LGP) :: advancing, on_face, on_top_face, partly_sat
+
+    if (.not. particle%advancing) then
       advancing = .false.
       particle%iboundary = 0
+      return
+    end if
+
+    call this%pass(particle)
+
+    ic = particle%itrdomain(LEVEL_FEATURE)
+    iface = particle%iboundary(LEVEL_FEATURE) - 1 ! cell is closed
+    on_face = iface >= 0
+    on_top_face = this%fmi%max_faces == iface
+    partly_sat = this%fmi%gwfsat(this%cell%defn%icell) < DONE
+
+    ! if at top and the cell is partially saturated,
+    ! we're at the water table, not the cell top, so
+    ! no exit event. if cell top is a boundary face,
+    ! terminate. if the dry tracking method is stop,
+    ! terminate. otherwise, leave the particle where
+    ! it is- keep tracking it on the next time step.
+    if (on_top_face .and. partly_sat) then
+      advancing = .false.
+      particle%advancing = .false.
+      if (this%fmi%is_net_out_boundary_face(ic, iface) .or. &
+          particle%idrymeth == 1) & ! dry_tracking_method stop
+        call this%terminate(particle, status=TERM_BOUNDARY)
+      return
+    end if
+
+    if (on_face) then
+      advancing = .false.
+      call this%cellexit(particle)
     end if
   end subroutine try_pass
 

src/Solution/ParticleTracker/Method/MethodCellPollock.f90

@@ -64,6 +64,7 @@ subroutine load_mcp(this, particle, next_level, submethod)
       call this%load_subcell(particle, subcell)
     end select
     call method_subcell_plck%init( &
+      fmi=this%fmi, &
       cell=this%cell, &
       subcell=this%subcell, &
       events=this%events, &
@@ -162,7 +163,8 @@ subroutine load_subcell(this, particle, subcell) !
 
     select type (cell => this%cell)
     type is (CellRectType)
-      ! Set subcell number to 1
+      ! Set cell/subcell numbers
+      subcell%icell = cell%defn%icell
       subcell%isubcell = 1
 
       ! Subcell calculations will be done in local subcell coordinates

src/Solution/ParticleTracker/Method/MethodDis.f90

@@ -204,18 +204,6 @@ subroutine load_particle(this, cell, particle)
       icu = dis%get_nodeuser(ic)
       call get_ijk(icu, dis%nrow, dis%ncol, dis%nlay, irow, icol, ilay)
 
-      ! TODO remove this cycle trap as soon as both cycle conditions are fixed!
-      if (ic == particle%icp .and. inface == 7 .and. ilay < particle%ilay) then
-        particle%itrdomain(LEVEL_FEATURE) = particle%icp
-        particle%izone = particle%izp
-        call this%terminate(particle, &
-                            status=TERM_BOUNDARY)
-        return
-      else
-        particle%icp = particle%itrdomain(LEVEL_FEATURE)
-        particle%izp = particle%izone
-      end if
-
       ! update node numbers and layer
       particle%itrdomain(LEVEL_FEATURE) = ic
       particle%icu = icu
@@ -288,7 +276,7 @@ subroutine pass_dis(this, particle)
     ! local
     type(CellRectType), pointer :: cell
     integer(I4B) :: iface
-    logical(LGP) :: no_neighbors, at_boundary
+    logical(LGP) :: at_boundary, no_neighbors
 
     select type (c => this%cell)
     type is (CellRectType)
@@ -304,7 +292,6 @@ subroutine pass_dis(this, particle)
       end if
 
       call this%load_particle(cell, particle)
-      if (.not. particle%advancing) return ! todo: remove after cycles fixed
       call this%update_flowja(cell, particle)
     end select
   end subroutine pass_dis

src/Solution/ParticleTracker/Method/MethodDisv.f90

@@ -166,17 +166,6 @@ subroutine load_particle(this, cell, particle)
       icu = dis%get_nodeuser(ic)
       call get_jk(icu, dis%ncpl, dis%nlay, icpl, ilay)
 
-      ! TODO remove this cycle trap as soon as both cycle conditions are fixed!
-      if (ic == particle%icp .and. inface == 7 .and. ilay < particle%ilay) then
-        particle%itrdomain(LEVEL_FEATURE) = particle%icp
-        particle%izone = particle%izp
-        call this%terminate(particle, status=TERM_BOUNDARY)
-        return
-      else
-        particle%icp = particle%itrdomain(LEVEL_FEATURE)
-        particle%izp = particle%izone
-      end if
-
       particle%itrdomain(LEVEL_FEATURE) = ic
       particle%icu = icu
       particle%ilay = ilay
@@ -225,7 +214,7 @@ subroutine pass_disv(this, particle)
     ! local
     type(CellPolyType), pointer :: cell
     integer(I4B) :: iface
-    logical(LGP) :: no_neighbors, at_boundary
+    logical(LGP) :: at_boundary, no_neighbors
 
     select type (c => this%cell)
     type is (CellPolyType)
@@ -241,7 +230,6 @@ subroutine pass_disv(this, particle)
       end if
 
       call this%load_particle(cell, particle)
-      if (.not. particle%advancing) return ! todo: remove after cycles fixed
       call this%update_flowja(cell, particle)
     end select
   end subroutine pass_disv

src/Solution/ParticleTracker/Method/MethodSubcellPollock.f90

@@ -115,15 +115,15 @@ subroutine track_subcell(this, subcell, particle, tmax)
     y0 = particle%y / subcell%dy
     z0 = particle%z / subcell%dz
 
-    ! Find exit solutions for each coordinate direction
+    ! Find exit solution in each direction
     call this%find_exits(particle, subcell)
 
     exit_x = this%exit_solutions(1)
     exit_y = this%exit_solutions(2)
     exit_z = this%exit_solutions(3)
 
     ! Subcell has no exit face, terminate the particle
-    ! todo: after initial release, consider ramifications
+    ! TODO: consider ramifications
     if (all([this%exit_solutions%status] == NO_EXIT_NO_OUTFLOW)) then
       call this%terminate(particle, status=TERM_NO_EXITS_SUB)
       return
@@ -260,19 +260,19 @@ function pick_exit(this, particle) result(exit_soln)
 
     exit_soln = 0
     dtmin = 1.0d+30
+
     if (this%exit_solutions(1)%status < 2) then
-      exit_soln = 1
+      exit_soln = 1 ! x
       dtmin = this%exit_solutions(1)%dt
     end if
     if (this%exit_solutions(2)%status < 2 .and. &
         this%exit_solutions(2)%dt < dtmin) then
-      exit_soln = 2
+      exit_soln = 2 ! y
       dtmin = this%exit_solutions(2)%dt
     end if
     if (this%exit_solutions(3)%status < 2 .and. &
         this%exit_solutions(3)%dt < dtmin) then
-      exit_soln = 3
-      dtmin = this%exit_solutions(3)%dt
+      exit_soln = 3 ! z
     end if
 
   end function pick_exit

src/Solution/ParticleTracker/Method/MethodSubcellTernary.f90

@@ -119,7 +119,7 @@ subroutine track_subcell(this, subcell, particle, tmax)
     exit_lateral = this%exit_solutions(2)
 
     ! If the subcell has no exit face, terminate the particle.
-    ! todo: after initial release, consider ramifications
+    ! TODO: consider ramifications
     if (exit_z%itopbotexit == 0 .and. &
         exit_lateral%itrifaceexit == 0) then
       call this%terminate(particle, status=TERM_NO_EXITS_SUB)
@@ -214,23 +214,14 @@ function pick_exit(this, particle) result(exit_soln)
     integer(I4B) :: exit_soln
 
     if (this%exit_solutions(1)%itopbotexit == 0) then
-      ! Exits through triangle face first
-      exit_soln = 2
-      this%exit_solutions(2)%iboundary = this%exit_solutions(2)%itrifaceexit
+      exit_soln = 2 ! lateral
     else if (this%exit_solutions(2)%itrifaceexit == 0 .or. &
              this%exit_solutions(1)%dt < this%exit_solutions(2)%dt) then
-      ! Exits through top/bottom first
-      exit_soln = 1
-      if (this%exit_solutions(1)%itopbotexit == -1) then
-        this%exit_solutions(1)%iboundary = 4
-      else
-        this%exit_solutions(1)%iboundary = 5
-      end if
+      exit_soln = 1 ! top/bottom
     else
-      ! Exits through triangle face first
-      exit_soln = 2
-      this%exit_solutions(2)%iboundary = this%exit_solutions(2)%itrifaceexit
+      exit_soln = 2 ! lateral
     end if
+
   end function pick_exit
 
   !> @brief Calculate exit solutions for each coordinate direction
@@ -264,7 +255,7 @@ subroutine find_exits(this, particle, domain)
     ntmax = 10000
     tol = particle%extol
 
-    ! Set solution method
+    ! Set lateral solution method
     if (particle%iexmeth == 0) then
       isolv = 1 ! default to Brent's
     else
@@ -310,7 +301,7 @@ subroutine find_exits(this, particle, domain)
       this%exit_solutions(1) = find_vertical_exit(subcell%vzbot, subcell%vztop, &
                                                   subcell%dz, zirel)
 
-      ! Calculate a lateral exit solution semi-analytically.
+      ! Calculate a semi-analytical lateral exit solution
       itrifaceenter = particle%iboundary(LEVEL_SUBFEATURE) - 1
       if (itrifaceenter == -1) itrifaceenter = 999
       this%exit_solutions(2) = find_lateral_exit(isolv, tol, &
@@ -324,6 +315,19 @@ subroutine find_exits(this, particle, domain)
       this%exit_solutions(2)%sxx = sxx
       this%exit_solutions(2)%sxy = sxy
       this%exit_solutions(2)%syy = syy
+
+      ! Set vertical solution exit face
+      if (this%exit_solutions(1)%itopbotexit /= 0) then
+        if (this%exit_solutions(1)%itopbotexit == -1) then
+          this%exit_solutions(1)%iboundary = 4
+        else
+          this%exit_solutions(1)%iboundary = 5
+        end if
+      end if
+
+      ! Set lateral solution exit face
+      if (this%exit_solutions(2)%itrifaceexit /= 0) &
+        this%exit_solutions(2)%iboundary = this%exit_solutions(2)%itrifaceexit
     end select
   end subroutine find_exits
 

src/Solution/ParticleTracker/Particle/Particle.f90

@@ -73,11 +73,9 @@ module ParticleModule
     ! state
     integer(I4B), public :: itrdomain(MAX_LEVEL) !< tracking domain indices
     integer(I4B), public :: iboundary(MAX_LEVEL) !< tracking domain boundary indices
-    integer(I4B), public :: icp !< previous cell number (reduced)
     integer(I4B), public :: icu !< user cell number
     integer(I4B), public :: ilay !< grid layer
     integer(I4B), public :: izone !< current zone number
-    integer(I4B), public :: izp !< previous zone number
     integer(I4B), public :: istatus !< tracking status
     real(DP), public :: x !< x coordinate
     real(DP), public :: y !< y coordinate
@@ -125,7 +123,6 @@ module ParticleModule
     integer(I4B), dimension(:), pointer, public, contiguous :: icu !< cell number (user)
     integer(I4B), dimension(:), pointer, public, contiguous :: ilay !< layer
     integer(I4B), dimension(:), pointer, public, contiguous :: izone !< current zone number
-    integer(I4B), dimension(:), pointer, public, contiguous :: izp !< previous zone number
     integer(I4B), dimension(:), pointer, public, contiguous :: istatus !< particle status
     real(DP), dimension(:), pointer, public, contiguous :: x !< model x coord of particle
     real(DP), dimension(:), pointer, public, contiguous :: y !< model y coord of particle
@@ -164,7 +161,6 @@ subroutine create_particle_store(store, np, mempath)
     call mem_allocate(store%icu, np, 'PLICU', mempath)
     call mem_allocate(store%ilay, np, 'PLILAY', mempath)
     call mem_allocate(store%izone, np, 'PLIZONE', mempath)
-    call mem_allocate(store%izp, np, 'PLIZP', mempath)
     call mem_allocate(store%istatus, np, 'PLISTATUS', mempath)
     call mem_allocate(store%x, np, 'PLX', mempath)
     call mem_allocate(store%y, np, 'PLY', mempath)
@@ -196,7 +192,6 @@ subroutine destroy(this, mempath)
     call mem_deallocate(this%icu, 'PLICU', mempath)
     call mem_deallocate(this%ilay, 'PLILAY', mempath)
     call mem_deallocate(this%izone, 'PLIZONE', mempath)
-    call mem_deallocate(this%izp, 'PLIZP', mempath)
     call mem_deallocate(this%istatus, 'PLISTATUS', mempath)
     call mem_deallocate(this%x, 'PLX', mempath)
     call mem_deallocate(this%y, 'PLY', mempath)
@@ -237,7 +232,6 @@ subroutine resize(this, np, mempath)
     call mem_reallocate(this%icu, np, 'PLICU', mempath)
     call mem_reallocate(this%ilay, np, 'PLILAY', mempath)
     call mem_reallocate(this%izone, np, 'PLIZONE', mempath)
-    call mem_reallocate(this%izp, np, 'PLIZP', mempath)
     call mem_reallocate(this%istatus, np, 'PLISTATUS', mempath)
     call mem_reallocate(this%x, np, 'PLX', mempath)
     call mem_reallocate(this%y, np, 'PLY', mempath)
@@ -279,11 +273,9 @@ subroutine get(this, particle, imdl, iprp, ip)
     particle%istopweaksink = this%istopweaksink(ip)
     particle%istopzone = this%istopzone(ip)
     particle%idrymeth = this%idrymeth(ip)
-    particle%icp = 0
     particle%icu = this%icu(ip)
     particle%ilay = this%ilay(ip)
     particle%izone = this%izone(ip)
-    particle%izp = this%izp(ip)
     particle%istatus = this%istatus(ip)
     particle%x = this%x(ip)
     particle%y = this%y(ip)
@@ -320,7 +312,6 @@ subroutine put(this, particle, ip)
     this%icu(ip) = particle%icu
     this%ilay(ip) = particle%ilay
     this%izone(ip) = particle%izone
-    this%izp(ip) = particle%izp
     this%istatus(ip) = particle%istatus
     this%x(ip) = particle%x
     this%y(ip) = particle%y