ParticleReleaseSchedule.f90

Go to the documentation of this file.

!> @brief Particle release scheduling.
module particlereleaseschedulemodule
 
  use arrayhandlersmodule, only: expandarray
  use constantsmodule, only: dzero, done, linelength
  use kindmodule, only: i4b, lgp, dp
  use mathutilmodule, only: is_close
  use timeselectmodule, only: timeselecttype
  use timestepselectmodule, only: timestepselecttype
 
  implicit none
  private
  public :: particlereleasescheduletype
  public :: create_release_schedule
 
  !> @brief Particle release scheduling utility.
  !!
  !! The release schedule composes a time selection object for any
  !! explicitly specified release times, and a time step selection
  !! object for release times specified in period/time step terms.
  !!
  !! Release time coincidence is computed within a given tolerance;
  !! times closer than the tolerance are merged into a single time.
  !!
  !! The release schedule must be refreshed each time step. This is
  !! achieved by calling `advance()`. After this, the `times` array
  !! is a debounced/consolidated schedule for the current time step.
  !<
  type :: particlereleasescheduletype
    real(dp), allocatable :: times(:) !< release times
    real(dp) :: tolerance !< release time coincidence tolerance
    type(timeselecttype), pointer :: time_select !< time selection
    type(timestepselecttype), pointer :: step_select !< time step selection
  contains
    procedure :: advance
    procedure :: any
    procedure :: count
    procedure :: destroy
    procedure :: log
    procedure :: schedule
  end type particlereleasescheduletype
 
contains
 
  !> @brief Create a new release schedule.
  function create_release_schedule(tolerance) result(schedule_)
    real(dp), intent(in) :: tolerance !< coincident release time tolerance
    type(particlereleasescheduletype), pointer :: schedule_ !< schedule pointer
 
    allocate (schedule_)
    allocate (schedule_%times(0))
    allocate (schedule_%time_select)
    allocate (schedule_%step_select)
    call schedule_%time_select%init()
    call schedule_%step_select%init()
    schedule_%tolerance = tolerance
  end function create_release_schedule
 
  !> @brief Deallocate the release schedule.
  subroutine destroy(this)
    class(particlereleasescheduletype), intent(inout) :: this !< this instance
 
    deallocate (this%times)
    call this%time_select%deallocate()
    call this%step_select%deallocate()
    deallocate (this%time_select)
    deallocate (this%step_select)
  end subroutine destroy
 
  !> @brief Write the release schedule to the given output unit.
  subroutine log(this, iout)
    class(particlereleasescheduletype), intent(inout) :: this !< this instance
    integer(I4B), intent(in) :: iout !< output unit
    character(len=*), parameter :: fmt = &
      &"(6x,A,': ',50(G0,' '))"
 
    if (this%any()) then
      write (iout, fmt) 'RELEASE SCHEDULE', this%times
    else
      write (iout, "(1x,a,1x,a)") 'NO RELEASES SCHEDULED'
    end if
  end subroutine log
 
  !> @brief Add a release time to the schedule.
  !!
  !! To schedule multiple release times at once, expand
  !! and populate the time selection object by hand. DO
  !! NOT attempt to manipulate the times array; this is
  !! a read-only property which the schedule maintains.
  !<
  subroutine schedule(this, trelease)
    class(particlereleasescheduletype), intent(inout) :: this
    real(DP), intent(in) :: trelease
    call expandarray(this%times)
    this%times(size(this%times)) = trelease
  end subroutine schedule
 
  !> @brief Refresh the schedule for the current time step.
  !!
  !! This involves several tasks: first, advance the time
  !! selection. Then, if period-block release setting lines
  !! are provided, reinitialize the time step selection for
  !! the given period. Finally, refresh the schedule array,
  !! deduplicating any times closer than the set tolerance.
  !!
  !! This routine is idempotent.
  !<
  subroutine advance(this, lines)
    use tdismodule, only: totimc, kstp, endofperiod
    class(particlereleasescheduletype), intent(inout) :: this
    character(len=LINELENGTH), intent(in), optional :: lines(:)
    integer(I4B) :: it, i
    real(DP) :: tprevious
    real(DP) :: trelease
 
    ! Advance the time selection.
    call this%time_select%advance()
 
    ! Reinitialize the time step selection if new
    ! period-block release settings are provided.
    if (present(lines)) then
      call this%step_select%init()
      do i = 1, size(lines)
        call this%step_select%read(lines(i))
      end do
    end if
 
    ! Reinitialize the release time schedule.
    if (allocated(this%times)) deallocate (this%times)
    allocate (this%times(0))
 
    tprevious = -done
    trelease = -done
 
    ! Add a release time configured by period-block
    ! settings, if one is scheduled this time step.
    if (this%step_select%is_selected(kstp, endofperiod=endofperiod)) then
      trelease = totimc
      call this%schedule(trelease)
      tprevious = trelease
    end if
 
    ! Schedule explicitly specified release times, up
    ! to the configured tolerance of coincidence
    if (this%time_select%any()) then
      do it = this%time_select%selection(1), this%time_select%selection(2)
        trelease = this%time_select%times(it)
        if (tprevious >= dzero .and. is_close( &
            tprevious, &
            trelease, &
            atol=this%tolerance)) cycle
 
        call this%schedule(trelease)
        tprevious = trelease
      end do
    end if
  end subroutine advance
 
  !> @brief Check if any releases are scheduled.
  !!
  !! Note: be sure to call advance() before calling this function,
  !! or the result may still be associated with a prior time step.
  !<
  logical function any(this) result(a)
    class(particlereleasescheduletype) :: this
    a = this%count() > 0
  end function any
 
  !> @brief Return the number of releases scheduled.
  !!
  !! Note: be sure to call advance() before calling this function,
  !! or the result may still be associated with a prior time step.
  !<
  integer function count(this) result(n)
    class(particlereleasescheduletype) :: this
    n = size(this%times)
  end function count
 
end module particlereleaseschedulemodule