add try_advance() routine to timeselect, simplify subcel methods

autotest/TestTimeSelect.f90

@@ -25,15 +25,15 @@ subroutine test_is_increasing(error)
 
     ! increasing
     ts%times = (/0.0_DP, 1.0_DP, 2.0_DP/)
-    call check(error, ts%is_increasing())
+    call check(error, ts%increasing())
 
     ! not decreasing
     ts%times = (/0.0_DP, 0.0_DP, 2.0_DP/)
-    call check(error,.not. ts%is_increasing())
+    call check(error,.not. ts%increasing())
 
     ! decreasing
     ts%times = (/2.0_DP, 1.0_DP, 0.0_DP/)
-    call check(error,.not. ts%is_increasing())
+    call check(error,.not. ts%increasing())
   end subroutine
 
   subroutine test_slice(error)
@@ -49,63 +49,63 @@ subroutine test_slice(error)
       "expected size 3, got"//to_string(size(ts%times)))
 
     ! empty slice
-    call ts%set_slice(1.1_DP, 1.9_DP)
+    call ts%select(1.1_DP, 1.9_DP)
     call check( &
       error, &
-      ts%slice(1) == -1 .and. ts%slice(2) == -1, &
+      ts%selection(1) == -1 .and. ts%selection(2) == -1, &
       "empty slice failed, got ["// &
-      to_string(ts%slice(1))//","//to_string(ts%slice(2))//"]")
+      to_string(ts%selection(1))//","//to_string(ts%selection(2))//"]")
 
     ! single-item slice
-    call ts%set_slice(0.5_DP, 1.5_DP)
+    call ts%select(0.5_DP, 1.5_DP)
     call check( &
       error, &
-      ts%slice(1) == 2 .and. ts%slice(2) == 2, &
+      ts%selection(1) == 2 .and. ts%selection(2) == 2, &
       "1-item slice failed, got ["// &
-      to_string(ts%slice(1))//","//to_string(ts%slice(2))//"]")
+      to_string(ts%selection(1))//","//to_string(ts%selection(2))//"]")
 
     ! multi-item slice
     changed = .false.
-    call ts%set_slice(0.5_DP, 2.5_DP, changed=changed)
+    call ts%select(0.5_DP, 2.5_DP, changed=changed)
     call check(error, changed)
     call check( &
       error, &
-      ts%slice(1) == 2 .and. ts%slice(2) == 3, &
+      ts%selection(1) == 2 .and. ts%selection(2) == 3, &
       "2-item slice failed, got ["// &
-      to_string(ts%slice(1))//","//to_string(ts%slice(2))//"]")
+      to_string(ts%selection(1))//","//to_string(ts%selection(2))//"]")
 
     ! no-change
-    call ts%set_slice(0.1_DP, 2.5_DP, changed=changed)
+    call ts%select(0.1_DP, 2.5_DP, changed=changed)
     call check(error,.not. changed)
     call check( &
       error, &
-      ts%slice(1) == 2 .and. ts%slice(2) == 3, &
+      ts%selection(1) == 2 .and. ts%selection(2) == 3, &
       "2-item slice failed, got ["// &
-      to_string(ts%slice(1))//","//to_string(ts%slice(2))//"]")
+      to_string(ts%selection(1))//","//to_string(ts%selection(2))//"]")
 
     ! lower bound equal to a time value
-    call ts%set_slice(0.0_DP, 2.5_DP)
+    call ts%select(0.0_DP, 2.5_DP)
     call check( &
       error, &
-      ts%slice(1) == 1 .and. ts%slice(2) == 3, &
+      ts%selection(1) == 1 .and. ts%selection(2) == 3, &
       "lb eq slice failed, got [" &
-      //to_string(ts%slice(1))//","//to_string(ts%slice(2))//"]")
+      //to_string(ts%selection(1))//","//to_string(ts%selection(2))//"]")
 
     ! upper bound equal to a time value
-    call ts%set_slice(-0.5_DP, 2.0_DP)
+    call ts%select(-0.5_DP, 2.0_DP)
     call check( &
       error, &
-      ts%slice(1) == 1 .and. ts%slice(2) == 3, &
+      ts%selection(1) == 1 .and. ts%selection(2) == 3, &
       "ub eq slice failed, got [" &
-      //to_string(ts%slice(1))//","//to_string(ts%slice(2))//"]")
+      //to_string(ts%selection(1))//","//to_string(ts%selection(2))//"]")
 
     ! both bounds equal to a time value
-    call ts%set_slice(0.0_DP, 2.0_DP)
+    call ts%select(0.0_DP, 2.0_DP)
     call check( &
       error, &
-      ts%slice(1) == 1 .and. ts%slice(2) == 3, &
+      ts%selection(1) == 1 .and. ts%selection(2) == 3, &
       "lb ub eq slice failed, got [" &
-      //to_string(ts%slice(1))//","//to_string(ts%slice(2))//"]")
+      //to_string(ts%selection(1))//","//to_string(ts%selection(2))//"]")
 
   end subroutine test_slice
 end module TestTimeSelect

src/Model/ModelUtilities/TimeSelect.f90

@@ -10,20 +10,21 @@ module TimeSelectModule
 
   !> @brief Represents a series of instants at which some event should occur.
   !!
-  !! Supports slicing e.g. to filter times within a given period & time step.
+  !! Supports selection e.g. to filter times in a selected period & time step.
   !! Array storage can be expanded as needed. Note: array expansion must take
-  !! place before slicing; whenever expand() is invoked, the slice is cleared.
-  !! The series is assumed to strictly increase, is_increasing() checks this.
+  !! place before selection; when expand() is invoked the selection is cleared.
+  !! The time series is assumed to strictly increase, increasing() checks this.
   !<
   type :: TimeSelectType
     real(DP), allocatable :: times(:)
-    integer(I4B) :: slice(2)
+    integer(I4B) :: selection(2)
   contains
     procedure :: destroy
     procedure :: expand
     procedure :: init
-    procedure :: is_increasing
-    procedure :: set_slice
+    procedure :: increasing
+    procedure :: select
+    procedure :: try_advance
   end type TimeSelectType
 
 contains
@@ -39,19 +40,19 @@ subroutine expand(this, increment)
     class(TimeSelectType) :: this
     integer(I4B), optional, intent(in) :: increment
     call ExpandArray(this%times, increment=increment)
-    this%slice = (/1, size(this%times)/)
+    this%selection = (/1, size(this%times)/)
   end subroutine expand
 
   !> @brief Initialize or clear the time selection object.
   subroutine init(this)
     class(TimeSelectType) :: this
     if (.not. allocated(this%times)) allocate (this%times(0))
-    this%slice = (/0, 0/)
+    this%selection = (/0, 0/)
   end subroutine
 
   !> @brief Determine if times strictly increase.
   !! Returns true if empty or not yet allocated.
-  function is_increasing(this) result(inc)
+  function increasing(this) result(inc)
     class(TimeSelectType) :: this
     logical(LGP) :: inc
     integer(I4B) :: i
@@ -69,17 +70,17 @@ function is_increasing(this) result(inc)
       end if
       l = t
     end do
-  end function is_increasing
+  end function increasing
 
-  !> @brief Slice the time selection between t0 and t1 (inclusive).
+  !> @brief Select times between t0 and t1 (inclusive).
   !!
   !! Finds and stores the index of the first time at the same instant
   !! as or following the start time, and of the last time at the same
   !! instant as or preceding the end time. Allows filtering the times
   !! for e.g. a particular stress period and time step. Array indices
   !! are assumed to start at 1. If no times are found to fall within
-  !! the slice (i.e. the slice falls entirely between two consecutive
-  !! times or beyond the min/max range), indices are set to [-1, -1].
+  !! the selection (i.e. it falls entirely between two consecutive
+  !! times or beyond the time range), indices are set to [-1, -1].
   !!
   !! The given start and end times are first checked against currently
   !! stored indices to avoid recalculating them if possible, allowing
@@ -88,7 +89,7 @@ end function is_increasing
   !! through stress periods and time steps in lockstep, i.e. they all
   !! solve any given period/step before any will proceed to the next.
   !<
-  subroutine set_slice(this, t0, t1, changed)
+  subroutine select(this, t0, t1, changed)
     ! -- dummy
     class(TimeSelectType) :: this
     real(DP), intent(in) :: t0, t1
@@ -107,8 +108,8 @@ subroutine set_slice(this, t0, t1, changed)
     u = -1
 
     ! -- previous bounding indices
-    lp = this%slice(1)
-    up = this%slice(2)
+    lp = this%selection(1)
+    up = this%selection(2)
 
     ! -- Check if we can reuse either the lower or upper bound.
     !    The lower doesn't need to change if it indexes the 1st
@@ -131,7 +132,7 @@ subroutine set_slice(this, t0, t1, changed)
         end if
       end if
       if (l == lp .and. u == up) then
-        this%slice = (/l, u/)
+        this%selection = (/l, u/)
         if (present(changed)) changed = .false.
         return
       end if
@@ -143,9 +144,24 @@ subroutine set_slice(this, t0, t1, changed)
       if (l < 0 .and. t >= t0 .and. t <= t1) l = i
       if (l > 0 .and. t <= t1) u = i
     end do
-    this%slice = (/l, u/)
+    this%selection = (/l, u/)
     if (present(changed)) changed = l /= lp .or. u /= up
 
   end subroutine
 
+  !> @brief Update the selection to match the current time step.
+  subroutine try_advance(this)
+    ! -- modules
+    use TdisModule, only: kper, kstp, nper, nstp, totimc, delt
+    ! -- dummy
+    class(TimeSelectType) :: this
+    ! -- local
+    real(DP) :: l, u
+    l = minval(this%times)
+    u = maxval(this%times)
+    if (.not. (kper == 1 .and. kstp == 1)) l = totimc
+    if (.not. (kper == nper .and. kstp == nstp(kper))) u = totimc + delt
+    call this%select(l, u)
+  end subroutine try_advance
+
 end module TimeSelectModule

src/Model/ParticleTracking/prt1oc1.f90

@@ -304,7 +304,7 @@ subroutine prt_oc_read_options(this)
             call this%tracktimes%expand()
             this%tracktimes%times(size(this%tracktimes%times)) = dval
           end do ttloop
-          if (.not. this%tracktimes%is_increasing()) then
+          if (.not. this%tracktimes%increasing()) then
             errmsg = "TRACK TIMES MUST STRICTLY INCREASE"
             call store_error(errmsg)
             call this%parser%StoreErrorUnit(terminate=.true.)
@@ -328,7 +328,7 @@ subroutine prt_oc_read_options(this)
             read (line, '(f30.0)') dval
             this%tracktimes%times(i) = dval
           end do
-          if (.not. this%tracktimes%is_increasing()) then
+          if (.not. this%tracktimes%increasing()) then
             errmsg = "TRACK TIMES MUST STRICTLY INCREASE"
             call store_error(errmsg)
             call this%parser%StoreErrorUnit(terminate=.true.)

src/Model/ParticleTracking/prt1prp1.f90

@@ -725,7 +725,7 @@ subroutine prp_options(this, option, found)
         call this%releasetimes%expand()
         this%releasetimes%times(size(this%releasetimes%times)) = dval
       end do rtloop
-      if (.not. this%releasetimes%is_increasing()) then
+      if (.not. this%releasetimes%increasing()) then
         errmsg = "RELEASE TIMES MUST STRICTLY INCREASE"
         call store_error(errmsg)
         call this%parser%StoreErrorUnit(terminate=.true.)
@@ -749,7 +749,7 @@ subroutine prp_options(this, option, found)
         read (line, '(f30.0)') dval
         this%releasetimes%times(i) = dval
       end do
-      if (.not. this%releasetimes%is_increasing()) then
+      if (.not. this%releasetimes%increasing()) then
         errmsg = "RELEASE TIMES MUST STRICTLY INCREASE"
         call store_error(errmsg)
         call this%parser%StoreErrorUnit(terminate=.true.)

src/Solution/ParticleTracker/MethodSubcellPollock.f90

@@ -83,7 +83,7 @@ end subroutine apply_msp
   subroutine track_subcell(this, subcell, particle, tmax)
     ! modules
     use ParticleModule, only: get_particle_id
-    use TdisModule, only: kper, kstp, nper, nstp, totimc, delt
+    use TdisModule, only: kper, kstp
     ! dummy
     class(MethodSubcellPollockType), intent(inout) :: this
     class(SubcellRectType), intent(in) :: subcell
@@ -92,11 +92,12 @@ subroutine track_subcell(this, subcell, particle, tmax)
     ! local
     doubleprecision :: vx, dvxdx, vy, dvydy, vz, dvzdz
     doubleprecision :: dtexitx, dtexity, dtexitz, dtexit, texit, dt, t, t0
-    doubleprecision :: x, y, z, l, u
+    doubleprecision :: x, y, z
     integer :: statusVX, statusVY, statusVZ, i
     doubleprecision :: initialX, initialY, initialZ
     integer :: exitFace
     integer :: reason
+    integer(I4B) :: tslice(2) !< user-time slice for the current time step
 
     reason = -1
 
@@ -157,19 +158,14 @@ subroutine track_subcell(this, subcell, particle, tmax)
     texit = particle%ttrack + dtexit
     t0 = particle%ttrack
 
-    if (all(this%tracktimes%slice > 0)) then
-      ! -- Select user tracking times to solve. If this is the first time step
-      !    of the simulation, include all times before it begins; if it is the
-      !    last time step, include all times after it ends. Otherwise take the
-      !    times within the current period and time step only.
-      l = minval(this%tracktimes%times)
-      u = maxval(this%tracktimes%times)
-      if (.not. (kper == 1 .and. kstp == 1)) l = totimc
-      if (.not. (kper == nper .and. kstp == nstp(kper))) u = totimc + delt
-      call this%tracktimes%set_slice(l, u)
-
-      ! -- Solve for tracking times in this slice
-      do i = this%tracktimes%slice(1), this%tracktimes%slice(2)
+    ! -- Select user tracking times to solve. If this is the first time step
+    !    of the simulation, include all times before it begins; if it is the
+    !    last time step include all times after it ends, otherwise the times
+    !    within the current period and time step only.
+    call this%tracktimes%try_advance()
+    tslice = this%tracktimes%selection
+    if (all(tslice > 0)) then
+      do i = tslice(1), tslice(2)
         t = this%tracktimes%times(i)
         if (t >= particle%ttrack .and. t < texit .and. t < tmax) then
           dt = t - t0

src/Solution/ParticleTracker/MethodSubcellTernary.f90

@@ -66,7 +66,7 @@ end subroutine apply_mst
   !> @brief Track a particle across a triangular subcell using the ternary method
   subroutine track_subcell(this, subcell, particle, tmax)
     ! modules
-    use TdisModule, only: kper, kstp, nper, nstp, totimc, delt
+    use TdisModule, only: kper, kstp
     ! dummy
     class(MethodSubcellTernaryType), intent(inout) :: this
     class(SubcellTriType), intent(in) :: subcell
@@ -82,12 +82,13 @@ subroutine track_subcell(this, subcell, particle, tmax)
     double precision :: rot(2, 2), res(2), loc(2)
     double precision :: alp, bet, alp0, bet0, alp1, bet1, alp2, bet2, alpi, beti
     double precision :: vzbot, vztop, vzi, vziodz, az, dtexitz
-    double precision :: dt, t, t0, dtexitxy, texit, x, y, z, l, u
+    double precision :: dt, t, t0, dtexitxy, texit, x, y, z
     integer :: izstatus, itopbotexit
     integer :: ntmax, nsave, isolv, itrifaceenter, itrifaceexit
     double precision :: tol, step, dtexit, alpexit, betexit
     integer :: ntdebug ! kluge
     integer :: reason, i
+    integer(I4B) :: tslice(2)
 
     lbary = .true. ! kluge
     ntmax = 10000
@@ -196,19 +197,14 @@ subroutine track_subcell(this, subcell, particle, tmax)
     texit = particle%ttrack + dtexit
     t0 = particle%ttrack
 
-    if (all(this%tracktimes%slice > 0)) then
-      ! -- Select user tracking times to solve. If this is the first time step
-      !    of the simulation, include all times before it begins; if it is the
-      !    last time step, include all times after it ends. Otherwise take the
-      !    times within the current period and time step only.
-      l = minval(this%tracktimes%times)
-      u = maxval(this%tracktimes%times)
-      if (.not. (kper == 1 .and. kstp == 1)) l = totimc
-      if (.not. (kper == nper .and. kstp == nstp(kper))) u = totimc + delt
-      call this%tracktimes%set_slice(l, u)
-
-      ! -- Solve for tracking times in this slice
-      do i = this%tracktimes%slice(1), this%tracktimes%slice(2)
+    ! -- Select user tracking times to solve. If this is the first time step
+    !    of the simulation, include all times before it begins; if it is the
+    !    last time step, include all times after it ends. Otherwise take the
+    !    times within the current period and time step only.
+    call this%tracktimes%try_advance()
+    tslice = this%tracktimes%selection
+    if (all(tslice > 0)) then
+      do i = tslice(1), tslice(2)
         t = this%tracktimes%times(i)
         if (t >= particle%ttrack .and. t < texit .and. t < tmax) then
           dt = t - t0