TimeSelect.f90

Go to the documentation of this file.

!> @brief Specify times for some event to occur.
module timeselectmodule
 
  use kindmodule, only: dp, i4b, lgp
  use constantsmodule, only: dzero, done
  use arrayhandlersmodule, only: expandarray
  use errorutilmodule, only: pstop
  use sortmodule, only: qsort
 
  implicit none
  public :: timeselecttype
 
  !> @brief Represents a series of instants at which some event should occur.
  !!
  !! Maintains an array of configured times which can be sliced to match e.g.
  !! the current period & time step. Slicing can be performed manually, with
  !! the select() routine, or automatically, with the advance() routine, for
  !! a convenient view onto the applicable subset of the complete time array.
  !!
  !! Array storage can be expanded manually. Note: array expansion must take
  !! place before selection; when expand() is called the selection is wiped.
  !! Alternatively, the extend() routine will automatically expand the array
  !! and sort it.
  !!
  !! Most use cases likely assume a strictly increasing time selection; this
  !! can be checked with increasing(). Note that the sort() routine does not
  !! check for duplicates, and should usually be followed by an increasing()
  !! check before the time selection is used.
  !<
  type :: timeselecttype
    real(dp), allocatable :: times(:)
    integer(I4B) :: selection(2)
  contains
    procedure :: deallocate
    procedure :: expand
    procedure :: init
    procedure :: increasing
    procedure :: log
    procedure :: select
    procedure :: advance
    procedure :: any
    procedure :: count
    procedure :: sort
    procedure :: extend
  end type timeselecttype
 
contains
 
  !> @brief Deallocate the time selection object.
  subroutine deallocate (this)
    class(timeselecttype) :: this
    deallocate (this%times)
  end subroutine deallocate
 
  !> @brief Expand capacity by the given amount. Resets the current slice.
  subroutine expand(this, increment)
    class(timeselecttype) :: this
    integer(I4B), optional, intent(in) :: increment
 
    call expandarray(this%times, increment=increment)
    this%selection = (/1, size(this%times)/)
  end subroutine expand
 
  !> @brief Initialize or clear the time selection object.
  subroutine init(this)
    class(timeselecttype) :: this
 
    if (allocated(this%times)) deallocate (this%times)
    allocate (this%times(0))
    this%selection = (/0, 0/)
  end subroutine
 
  !> @brief Determine if times strictly increase.
  !!
  !! Returns true if the times array strictly increases,
  !! as well as if the times array is empty, or not yet
  !! allocated. Note that this function operates on the
  !! entire times array, not the current selection. Note
  !! also that this function conducts exact comparisons;
  !! deduplication with tolerance must be done manually.
  !<
  function increasing(this) result(inc)
    class(timeselecttype) :: this
    logical(LGP) :: inc
    integer(I4B) :: i
    real(dp) :: l, t
 
    inc = .true.
    if (.not. allocated(this%times)) return
    do i = 1, size(this%times)
      t = this%times(i)
      if (i /= 1) then
        if (l >= t) then
          inc = .false.
          return
        end if
      end if
      l = t
    end do
  end function increasing
 
  !> @brief Show the current time selection, if any.
  subroutine log(this, iout, verb)
    ! dummy
    class(timeselecttype) :: this !< this instance
    integer(I4B), intent(in) :: iout !< output unit
    character(len=*), intent(in) :: verb !< selection name
    ! formats
    character(len=*), parameter :: fmt = &
      &"(6x,'THE FOLLOWING TIMES WILL BE ',A,': ',50(G0,' '))"
 
    if (this%any()) then
      write (iout, fmt) verb, this%times(this%selection(1):this%selection(2))
    else
      write (iout, "(a,1x,a)") 'NO TIMES WILL BE', verb
    end if
  end subroutine log
 
  !> @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 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
  !! multiple consuming components (e.g., subdomain particle tracking
  !! solutions) to share the object efficiently, provided all proceed
  !! 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 select(this, t0, t1, changed)
    ! dummy
    class(timeselecttype) :: this
    real(DP), intent(in) :: t0, t1
    logical(LGP), intent(inout), optional :: changed
    ! local
    integer(I4B) :: i, i0, i1
    integer(I4B) :: l, u, lp, up
    real(DP) :: t
 
    ! by default, need to iterate over all times
    i0 = 1
    i1 = size(this%times)
 
    ! if no times fall within the slice, set to [-1, -1]
    l = -1
    u = -1
 
    ! previous bounding indices
    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
    ! time simultaneous with or later than the slice's start.
    ! The upper doesn't need to change if it indexes the last
    ! time before or simultaneous with the slice's end.
    if (lp > 0 .and. up > 0) then
      if (lp > 1) then
        if (this%times(lp - 1) < t0 .and. &
            this%times(lp) >= t0) then
          l = lp
          i0 = l
        end if
      end if
      if (up > 1 .and. up < i1) then
        if (this%times(up + 1) > t1 .and. &
            this%times(up) <= t1) then
          u = up
          i1 = u
        end if
      end if
      if (l == lp .and. u == up) then
        this%selection = (/l, u/)
        if (present(changed)) changed = .false.
        return
      end if
    end if
 
    ! recompute bounding indices if needed
    do i = i0, i1
      t = this%times(i)
      if (l < 0 .and. t >= t0 .and. t <= t1) l = i
      if (l > 0 .and. t <= t1) u = i
    end do
    this%selection = (/l, u/)
    if (present(changed)) changed = l /= lp .or. u /= up
 
  end subroutine
 
  !> @brief Update the selection to the current time step.
  subroutine 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 advance
 
  !> @brief Check if any times are currently selected.
  !!
  !! Indicates whether any times are selected for the
  !! current time step.
  !!
  !! Note that this routine does NOT indicate whether
  !! the times array has nonzero size; use the size
  !! intrinsic for that.
  !<
  function any(this) result(a)
    class(timeselecttype) :: this
    logical(LGP) :: a
 
    a = all(this%selection > 0)
  end function any
 
  !> @brief Return the number of times currently selected.
  !!
  !! Returns the number of times selected for the current
  !! time step.
  !!
  !! Note that this routine does NOT return the total size
  !! of the times array; use the size intrinsic for that.
  !<
  function count(this) result(n)
    class(timeselecttype) :: this
    integer(I4B) :: n
 
    if (this%any()) then
      n = this%selection(2) - this%selection(1)
    else
      n = 0
    end if
  end function count
 
  !> @brief Sort the time selection in increasing order.
  !!
  !! Note that this routine does NOT remove duplicate times.
  !! Call increasing() to check for duplicates in the array.
  !<
  subroutine sort(this)
    class(timeselecttype) :: this
    integer(I4B), allocatable :: indx(:)
 
    allocate (indx(size(this%times)))
    call qsort(indx, this%times)
    deallocate (indx)
  end subroutine sort
 
  !> @brief Extend the time selection with the given array.
  !!
  !! This routine sorts the selection after appending the
  !! elements of the given array, but users should likely
  !! still call increasing() to check for duplicate times.
  !<
  subroutine extend(this, a)
    class(timeselecttype) :: this
    real(DP) :: a(:)
 
    this%times = [this%times, a]
    call this%sort()
  end subroutine extend
 
end module timeselectmodule