fix vertical termination conditions

autotest/test_prt_dry.py

@@ -369,10 +369,10 @@ def check_output(idx, test, snapshot):
     hds_path = gwf_ws / hds_file
     trackcsv_file = prt_name + ".csv"
     trackcsv_path = prt_ws / trackcsv_file
-    mf6pathlines = pd.read_csv(trackcsv_path)
-    startpts = mf6pathlines[mf6pathlines.ireason == 0]
+    pls = pd.read_csv(trackcsv_path)
+    strtpts = pls[pls.ireason == 0]
 
-    assert snapshot == mf6pathlines.drop(["name", "icell"], axis=1).round(3).to_records(
+    assert snapshot == pls.drop(["name", "icell"], axis=1).round(3).to_records(
         index=False
     )
 
@@ -438,7 +438,7 @@ def plot_pathlines_and_timeseries(
             plt.scatter(pathlines["x"], pathlines["y"])
             mm.plot_pathline(pathlines, layer="all", colors=["blue"], lw=0.75)
             ax.set_title(plottitle, fontsize=12)
-            ax.scatter(startpts.x, startpts.y)
+            ax.scatter(strtpts.x, strtpts.y)
             from shapely.geometry import Polygon
 
             cellids = [105, 85, 125, 104, 106]
@@ -468,7 +468,7 @@ def plot_pathlines_and_timeseries(
         ibd[ilay, irow, icol] = 2
         ibd = np.ma.masked_equal(ibd, 0)
 
-        plot_pathlines_and_timeseries(ax, gwf.modelgrid, ibd, mf6pathlines, None, name)
+        plot_pathlines_and_timeseries(ax, gwf.modelgrid, ibd, pls, None, name)
 
     plot_3d = False
     if plot_3d:
@@ -482,7 +482,7 @@ def plot_pathlines_and_timeseries(
         vert_exag = 1
         vtk = Vtk(model=gwf, binary=False, vertical_exageration=vert_exag, smooth=False)
         vtk.add_model(gwf)
-        vtk.add_pathline_points(mf6pathlines.to_records(index=False))
+        vtk.add_pathline_points(pls.to_records(index=False))
         gwf_mesh, prt_mesh = vtk.to_pyvista()
         axes = pv.Axes(show_actor=False, actor_scale=2.0, line_width=5)
         p = pv.Plotter(window_size=[700, 700])

autotest/test_prt_disv1.py → autotest/test_prt_newton_disv1.py

@@ -6,13 +6,7 @@
 refined cells, instead of the new ternary method
 which applies more generally to polygonal cells.
 
-The simulation includes a single stress period
-with multiple time steps. This serves to test
-whether PRT properly solves trajectories over
-"internal" time steps, i.e. within the step's
-slice of simulation time, as well as extending
-tracking to termination or particle stop times
-during the simulation's final time step.
+The simulation uses the Newton formulation.
 
 Several cases are provided:
     - default: No user-specified tracking times, MP7 in pathline mode.

src/Model/ParticleTracking/prt-prp.f90

@@ -443,7 +443,7 @@ subroutine initialize_particle(this, particle, ip, trelease)
     real(DP), intent(in) :: trelease !< release time
     ! local
     integer(I4B) :: irow, icol, ilay, icpl
-    integer(I4B) :: ic, icu
+    integer(I4B) :: ic, icu, ic_old
     real(DP) :: x, y, z
     real(DP) :: top, bot, hds
 
@@ -461,7 +461,7 @@ subroutine initialize_particle(this, particle, ip, trelease)
     particle%irpt = ip
     particle%istopweaksink = this%istopweaksink
     particle%istopzone = this%istopzone
-    particle%idry = this%idry
+    particle%idrymeth = this%idry
     particle%icu = icu
 
     select type (dis => this%dis)
@@ -478,9 +478,14 @@ subroutine initialize_particle(this, particle, ip, trelease)
     ! to the top-most active cell beneath it if drape is
     ! enabled, or else terminate permanently unreleased.
     if (this%ibound(ic) == 0) then
-      if (this%idrape > 0) &
+      ic_old = ic
+      if (this%idrape > 0) then
         call this%dis%highest_active(ic, this%ibound)
-      if (this%ibound(ic) == 0) particle%istatus = 8
+        if (ic == ic_old .or. this%ibound(ic) == 0) &
+          particle%istatus = 8
+      else
+        particle%istatus = 8
+      end if
     end if
 
     ! Load coordinates and transform if needed
@@ -928,7 +933,7 @@ subroutine prp_read_packagedata(this)
 
         if (n < 1 .or. n > this%nreleasepoints) then
           write (errmsg, '(a,1x,i0,a)') &
-            'Release point number must be greater than 0 and less than ', &
+            'Release point number must be greater than 0 and less than', &
             'or equal to', this%nreleasepoints, '.'
           call store_error(errmsg)
           cycle

src/Model/ParticleTracking/prt.f90

@@ -932,15 +932,15 @@ subroutine prt_solve(this)
 
           ! -- If particle is permanently unreleased, record its initial/terminal state
           if (particle%istatus == 8) &
-            call this%method%save(particle, reason=3) ! reason=3: termination
+            call this%method%save(particle, reason=3)
 
           ! If particle is inactive or not yet to be released, cycle
           if (particle%istatus > 1) cycle
 
           ! If particle released this time step, record its initial state
           particle%istatus = 1
           if (particle%trelease >= totimc) &
-            call this%method%save(particle, reason=0) ! reason=0: release
+            call this%method%save(particle, reason=0)
 
           ! Maximum time is the end of the time step or the particle
           ! stop time, whichever comes first, unless it's the final

src/Solution/ParticleTracker/CellDefn.f90

@@ -17,6 +17,7 @@ module CellDefnModule
     integer(I4B), public :: npolyverts !< number of vertices for cell polygon
     real(DP), public :: porosity !< cell porosity
     real(DP), public :: retfactor !< cell retardation factor
+    integer(I4B), public :: ilay !< layer number
     integer(I4B), public :: izone !< cell zone number
     integer(I4B), public :: iweaksink !< weak sink indicator
     integer(I4B), public :: inoexitface !< no exit face indicator

src/Solution/ParticleTracker/Method.f90

@@ -7,6 +7,8 @@ module MethodModule
   use SubcellModule, only: SubcellType
   use ParticleModule
   use BaseDisModule, only: DisBaseType
+  use DisModule, only: DisType
+  use DisvModule, only: DisvType
   use PrtFmiModule, only: PrtFmiType
   use CellModule, only: CellType
   use CellDefnModule, only: CellDefnType
@@ -29,7 +31,7 @@ module MethodModule
   !! depending on cell geometry (implementing the strategy pattern).
   !<
   type, abstract :: MethodType
-    character(len=40), pointer, public :: type !< method name
+    character(len=40), pointer, public :: name !< method name
     logical(LGP), public :: delegates !< whether the method delegates
     type(PrtFmiType), pointer, public :: fmi => null() !< ptr to fmi
     class(CellType), pointer, public :: cell => null() !< ptr to the current cell
@@ -52,7 +54,7 @@ module MethodModule
     procedure :: save
     procedure :: track
     procedure :: try_pass
-    procedure :: prepare
+    procedure :: check
   end type MethodType
 
   abstract interface
@@ -96,7 +98,8 @@ subroutine init(this, fmi, cell, subcell, trackctl, tracktimes, &
     if (present(retfactor)) this%retfactor => retfactor
   end subroutine init
 
-  !> @brief Track particle through subdomains
+  !> @brief Track the particle over domains of the given
+  ! level until the particle terminates or tmax elapses.
   recursive subroutine track(this, particle, level, tmax)
     ! dummy
     class(MethodType), intent(inout) :: this
@@ -108,6 +111,7 @@ recursive subroutine track(this, particle, level, tmax)
     integer(I4B) :: nextlevel
     class(methodType), pointer :: submethod
 
+    ! Advance the particle over subdomains
     advancing = .true.
     nextlevel = level + 1
     do while (advancing)
@@ -132,8 +136,9 @@ subroutine try_pass(this, particle, nextlevel, advancing)
       ! otherwise pass the particle to the next subdomain.
       ! if that leaves it on a boundary, stop advancing.
       call this%pass(particle)
-      if (particle%iboundary(nextlevel - 1) .ne. 0) &
+      if (particle%iboundary(nextlevel - 1) .ne. 0) then
         advancing = .false.
+      end if
     end if
   end subroutine try_pass
 
@@ -180,26 +185,24 @@ subroutine save(this, particle, reason)
       end if
     end if
 
-    ! Save the particle's state to any registered tracking output files
-    call this%trackctl%save(particle, kper=per, &
-                            kstp=stp, reason=reason)
+    call this%trackctl%save(particle, kper=per, kstp=stp, reason=reason)
   end subroutine save
 
-  !> @brief Prepare to apply the tracking method to the particle.
+  !> @brief Check reporting/terminating conditions before tracking.
   !!
   !! Check a number of conditions determining whether to continue
-  !! tracking the particle or terminate it. This includes a check
-  !! for any reporting conditions as well.
+  !! tracking the particle or terminate it, as well as whether to
+  !! record any output data as per selected reporting conditions.
   !<
-  subroutine prepare(this, particle, cell_defn)
+  subroutine check(this, particle, cell_defn)
+    ! 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
     ! local
-    logical(LGP) :: cell_dry
-    logical(LGP) :: locn_dry
-    integer(I4B) :: ic
+    logical(LGP) :: no_exit_face
 
     ! stop zone
     particle%izone = cell_defn%izone
@@ -210,40 +213,57 @@ subroutine prepare(this, particle, cell_defn)
       return
     end if
 
-    ! dry
-    cell_dry = this%fmi%ibdgwfsat0(cell_defn%icell) == 0
-    locn_dry = particle%z > cell_defn%top
-    if (cell_dry .or. locn_dry) then
-      if (particle%idry == 0) then
-        ! drop
-        if (cell_dry) then
-          ! if no active cell underneath position, terminate
-          ic = particle%idomain(2)
-          call this%fmi%dis%highest_active(ic, this%fmi%ibound)
-          if (this%fmi%ibound(ic) == 0) then
-            particle%advancing = .false.
-            particle%istatus = 7
-            call this%save(particle, reason=3)
-            return
-          end if
-        else if (locn_dry) then
-          particle%z = cell_defn%top
-          call this%save(particle, reason=1)
-        end if
-      else if (particle%idry == 1) then
+    ! no exit face
+    no_exit_face = cell_defn%inoexitface > 0
+
+    ! dry cell
+    if (this%fmi%ibdgwfsat0(cell_defn%icell) == 0) then
+      if (particle%idrymeth == 0) then
+        no_exit_face = .false.
+      else if (particle%idrymeth == 1) then
         ! stop
         particle%advancing = .false.
         particle%istatus = 7
         call this%save(particle, reason=3)
         return
-      else if (particle%idry == 2) then
+      else if (particle%idrymeth == 2) then
         ! stay
+        no_exit_face = .false.
         particle%advancing = .false.
+        particle%ttrack = totim
+        ! 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 (this%name /= "passtobottom" .and. particle%z > cell_defn%top) then
+      ! dry particle
+      if (particle%idrymeth == 0) then
+        ! drop to water table
+        no_exit_face = .false.
+        particle%z = cell_defn%top
+        call this%save(particle, reason=1)
+      else if (particle%idrymeth == 1) then
+        ! terminate
+        particle%advancing = .false.
+        particle%istatus = 7
+        call this%save(particle, reason=3)
+        return
+      else if (particle%idrymeth == 2) then
+        ! stay
+        no_exit_face = .false.
+        particle%advancing = .false.
+        return
+      end if
+    end if
 
-      ! cell with no exit face (mutually exclusive with
-      ! dry because a dry cell will have no exit face)
-    else if (cell_defn%inoexitface > 0) then
+    ! cell with no exit face. handle after
+    ! dry conditions because dry cells will
+    ! also have no exit face.
+    if (no_exit_face) then
       particle%advancing = .false.
       particle%istatus = 5
       call this%save(particle, reason=3)
@@ -261,6 +281,6 @@ subroutine prepare(this, particle, cell_defn)
         return
       end if
     end if
-  end subroutine prepare
+  end subroutine check
 
 end module MethodModule

src/Solution/ParticleTracker/MethodCellPassToBot.f90

@@ -9,6 +9,8 @@ module MethodCellPassToBotModule
   use CellModule, only: CellType
   use SubcellModule, only: SubcellType
   use TrackControlModule, only: TrackControlType
+  use DisModule, only: DisType
+  use DisvModule, only: DisvType
   implicit none
 
   private
@@ -28,15 +30,15 @@ module MethodCellPassToBotModule
   subroutine create_method_cell_ptb(method)
     type(MethodCellPassToBotType), pointer :: method
     allocate (method)
-    allocate (method%type)
-    method%type = "passtobottom"
+    allocate (method%name)
+    method%name = "passtobottom"
     method%delegates = .false.
   end subroutine create_method_cell_ptb
 
   !> @brief Deallocate the pass-to-bottom tracking method
   subroutine deallocate (this)
     class(MethodCellPassToBotType), intent(inout) :: this
-    deallocate (this%type)
+    deallocate (this%name)
   end subroutine deallocate
 
   !> @brief Pass particle vertically and instantaneously to the cell bottom
@@ -45,11 +47,32 @@ subroutine apply_ptb(this, particle, tmax)
     class(MethodCellPassToBotType), intent(inout) :: this
     type(ParticleType), pointer, intent(inout) :: particle
     real(DP), intent(in) :: tmax
+    ! local
+    integer(I4B) :: nlay, reason
 
-    call this%prepare(particle, this%cell%defn)
+    ! Check termination/reporting conditions
+    call this%check(particle, this%cell%defn)
     if (.not. particle%advancing) return
+
+    ! Pass to bottom face
     particle%z = this%cell%defn%bot
     particle%iboundary(2) = this%cell%defn%npolyverts + 2
+
+    ! Terminate if bottom layer
+    select type (dis => this%fmi%dis)
+    type is (DisType)
+      nlay = dis%nlay
+    type is (DisvType)
+      nlay = dis%nlay
+    end select
+    if (this%cell%defn%ilay == nlay) then
+      particle%istatus = 5
+      reason = 3
+    else
+      reason = 1
+    end if
+
+    ! Save datum
     call this%save(particle, reason=1)
   end subroutine apply_ptb