fix(prt): record user tracking times for stationary particles

Author: wpbonelli

autotest/__snapshots__/test_prt_dry/test_mf6model[4-prtdry_stay].npy

@@ -70,6 +70,8 @@
 for i in range(nper):
     period_data.append((perlen[i], nstp[i], tsmult[i]))
 
+user_time = 100.0
+
 
 def build_gwf_sim(name, gwf_ws, mf6, newton=False):
     sim = flopy.mf6.MFSimulation(
@@ -319,6 +321,8 @@ def build_prt_sim(name, gwf, prt_ws, mf6, drape=False, dry_tracking_method=False
         trackcsv_filerecord=[trackcsvfile],
         saverecord=[("BUDGET", "ALL")],
         pname="oc",
+        ntracktimes=1 if "stay" in name else 0,
+        tracktimes=[(user_time,)] if "stay" in name else None,
     )
 
     rel_prt_folder = os.path.relpath(gwf_ws, start=prt_ws)
@@ -373,6 +377,8 @@ def check_output(idx, test, snapshot):
     # ignore particle 4, it terminates early with optimization=2 when built with ifort
     if "drop" in name:
         actual = actual.drop(actual[actual.irpt == 4].index)
+    if "stay" in name:
+        assert len(actual[actual.t == user_time]) == 5
     assert snapshot == actual.to_records(index=False)
 
     plot_pathlines = False

autotest/test_prt_dry.py

@@ -191,15 +191,18 @@ end subroutine save
   !! tracking the particle or terminate it, as well as whether to
   !! record any output data as per selected reporting conditions.
   !<
-  subroutine check(this, particle, cell_defn)
+  subroutine check(this, particle, cell_defn, tmax)
     ! modules
     use TdisModule, only: endofsimulation, totim
     ! dummy
     class(MethodType), intent(inout) :: this
     type(ParticleType), pointer, intent(inout) :: particle
     type(CellDefnType), pointer, intent(inout) :: cell_defn
+    real(DP), intent(in) :: tmax
     ! local
     logical(LGP) :: dry_cell, dry_particle, no_exit_face, stop_zone, weak_sink
+    integer(I4B) :: i
+    real(DP) :: t
 
     dry_cell = this%fmi%ibdgwfsat0(cell_defn%icell) == 0
     dry_particle = particle%z > cell_defn%top
@@ -235,6 +238,8 @@ subroutine check(this, particle, cell_defn)
 
     if (dry_cell) then
       if (particle%idrymeth == 0) then
+        ! drop to cell bottom. handled by pass
+        ! to bottom method, nothing to do here
         no_exit_face = .false.
       else if (particle%idrymeth == 1) then
         ! stop
@@ -246,23 +251,38 @@ subroutine check(this, particle, cell_defn)
         ! stay
         no_exit_face = .false.
         particle%advancing = .false.
+
+        ! update tracking time to time
+        ! step end time and save record
         particle%ttrack = totim
+        call this%save(particle, reason=2)
+
+        ! record user tracking times
+        call this%tracktimes%advance()
+        if (this%tracktimes%any()) then
+          do i = this%tracktimes%selection(1), this%tracktimes%selection(2)
+            t = this%tracktimes%times(i)
+            if (t < particle%ttrack) cycle
+            if (t >= tmax) exit
+            particle%ttrack = t
+            call this%save(particle, reason=5)
+          end do
+        end if
+
         ! terminate if last period/step
         if (endofsimulation) then
           particle%istatus = 5
           call this%save(particle, reason=3)
           return
         end if
-        call this%save(particle, reason=2)
       end if
     else if (dry_particle .and. this%name /= "passtobottom") then
-      ! dry particle
       if (particle%idrymeth == 0) then
         ! drop to water table
         particle%z = cell_defn%top
         call this%save(particle, reason=1)
       else if (particle%idrymeth == 1) then
-        ! terminate
+        ! stop
         particle%advancing = .false.
         particle%istatus = 7
         call this%save(particle, reason=3)
@@ -271,6 +291,24 @@ subroutine check(this, particle, cell_defn)
         ! stay
         no_exit_face = .false.
         particle%advancing = .false.
+
+        ! update tracking time to time
+        ! step end time and save record
+        particle%ttrack = totim
+        call this%save(particle, reason=2)
+
+        ! record user tracking times
+        call this%tracktimes%advance()
+        if (this%tracktimes%any()) then
+          do i = this%tracktimes%selection(1), this%tracktimes%selection(2)
+            t = this%tracktimes%times(i)
+            if (t < particle%ttrack) cycle
+            if (t >= tmax) exit
+            particle%ttrack = t
+            call this%save(particle, reason=5)
+          end do
+        end if
+        
         return
       end if
     end if

src/Solution/ParticleTracker/Method.f90

@@ -47,7 +47,7 @@ subroutine apply_ptb(this, particle, tmax)
     real(DP), intent(in) :: tmax
 
     ! Check termination/reporting conditions
-    call this%check(particle, this%cell%defn)
+    call this%check(particle, this%cell%defn, tmax)
     if (.not. particle%advancing) return
 
     ! Pass to bottom face

src/Solution/ParticleTracker/MethodCellPassToBot.f90

@@ -130,7 +130,7 @@ subroutine apply_mcp(this, particle, tmax)
     select type (cell => this%cell)
     type is (CellRectType)
       ! Check termination/reporting conditions
-      call this%check(particle, cell%defn)
+      call this%check(particle, cell%defn, tmax)
       if (.not. particle%advancing) return
 
       ! Transform model coordinates to local cell coordinates

src/Solution/ParticleTracker/MethodCellPollock.f90

@@ -209,7 +209,7 @@ subroutine apply_mcpq(this, particle, tmax)
     select type (cell => this%cell)
     type is (CellRectQuadType)
       ! Check termination/reporting conditions
-      call this%check(particle, cell%defn)
+      call this%check(particle, cell%defn, tmax)
       if (.not. particle%advancing) return
 
       ! Transform model coordinates to local cell coordinates

src/Solution/ParticleTracker/MethodCellPollockQuad.f90

@@ -165,7 +165,7 @@ subroutine apply_mct(this, particle, tmax)
     select type (cell => this%cell)
     type is (CellPolyType)
       ! Check termination/reporting conditions
-      call this%check(particle, this%cell%defn)
+      call this%check(particle, this%cell%defn, tmax)
       if (.not. particle%advancing) return
 
       ! Number of vertices