time_based_controller.f90

Go to the documentation of this file.

! Copyright (c) 2023, The Neko Authors
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
!
!   * Redistributions of source code must retain the above copyright
!     notice, this list of conditions and the following disclaimer.
!
!   * Redistributions in binary form must reproduce the above
!     copyright notice, this list of conditions and the following
!     disclaimer in the documentation and/or other materials provided
!     with the distribution.
!
!   * Neither the name of the authors nor the names of its
!     contributors may be used to endorse or promote products derived
!     from this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
module time_based_controller
  use num_types, only : dp, i8
  use utils, only : neko_error
  use time_state, only : time_state_t
  implicit none
  private
 
  real(kind=dp), public, parameter :: time_tol = 0.1_dp
 
  real(kind=dp), public, parameter :: span_tol = 1.0e-9_dp
 
  type, public :: time_based_controller_t
     real(kind=dp) :: frequency = 0.0_dp
     real(kind=dp) :: time_interval = 0.0_dp
     integer :: nsteps = 0
     real(kind=dp) :: start_time = 0.0_dp
     real(kind=dp) :: end_time = 0.0_dp
     integer :: nexecutions = 0
     logical :: never = .false.
     character(len=:), allocatable :: control_mode
     real(kind=dp) :: control_value = 0.0_dp
     real(kind=dp) :: anchor_time = 0.0_dp
     logical :: write_at_start = .true.
     integer(kind=i8) :: first_index = 0
     logical :: start_is_scheduled = .false.
     logical :: start_pending = .false.
     integer :: next_index = 0
     real(kind=dp) :: direction = 1.0_dp
     integer :: tstep_offset = 0
     integer :: last_tstep = -1
 
   contains
     procedure, pass(this) :: init => time_based_controller_init
     procedure, pass(this) :: free => time_based_controller_free
     procedure, pass(this) :: check => time_based_controller_check
     procedure, pass(this) :: register_execution => &
          time_based_controller_register_execution
     procedure, pass(this) :: set_counter => &
          time_based_controller_set_counter
     procedure, pass(this) :: next_time => time_based_controller_next_time
     procedure, pass(this) :: time_to_next => &
          time_based_controller_time_to_next
     procedure, pass(this) :: tolerance => time_based_controller_tolerance
  end type time_based_controller_t
 
contains
 
  subroutine time_based_controller_init(this, start_time, end_time, &
       control_mode, control_value, write_at_start, anchor_time, direction)
    class(time_based_controller_t), intent(inout) :: this
    real(kind=dp), intent(in) :: start_time
    real(kind=dp), intent(in) :: end_time
    character(len=*), intent(in) :: control_mode
    real(kind=dp), intent(in) :: control_value
    logical, intent(in), optional :: write_at_start
    real(kind=dp), intent(in), optional :: anchor_time
    real(kind=dp), intent(in), optional :: direction
    real(kind=dp) :: span, offset
 
    call this%free()
 
    this%start_time = start_time
    this%end_time = end_time
    this%control_mode = control_mode
    this%control_value = control_value
 
    if (present(write_at_start)) then
       this%write_at_start = write_at_start
    else
       this%write_at_start = .true.
    end if
 
    if (present(direction)) then
       this%direction = sign(1.0_dp, direction)
    else
       this%direction = sign(1.0_dp, end_time - start_time)
    end if
    ! Length of the interval covered, positive also for a run backwards.
    span = this%direction * (end_time - start_time)
 
    if (trim(control_mode) .eq. 'simulationtime') then
       if (control_value .le. 0.0_dp) then
          call neko_error("The output interval must be positive")
       end if
       this%time_interval = control_value
       this%frequency = 1.0_dp / this%time_interval
       this%nsteps = 0
    else if (trim(control_mode) .eq. 'nsamples') then
       if (control_value .le. 0.0_dp) then
          call neko_error("nsamples must be positive")
       end if
       if (span .le. 0.0_dp) then
          call neko_error("nsamples requires the output to start before the &
          &end of the simulation")
       end if
 
       this%frequency = control_value / span
       this%time_interval = 1.0_dp / this%frequency
       this%nsteps = 0
       ! The samples divide the interval from start_time to end_time, so
       ! they are counted from start_time.
       this%anchor_time = start_time
    else if (trim(control_mode) .eq. 'tsteps') then
       if (control_value .lt. 1.0_dp) then
          call neko_error("The output interval in time steps must be at &
          &least 1")
       end if
       this%nsteps = int(control_value)
       ! if the timestep will be variable, we cannot compute these.
       this%frequency = 0.0_dp
       this%time_interval = 0.0_dp
    else if (trim(control_mode) .eq. 'never') then
       this%never = .true.
    else
       call neko_error("The control parameter must be simulationtime, nsamples&
       & tsteps, or never, but received "//trim(control_mode))
    end if
 
    if (present(anchor_time)) this%anchor_time = anchor_time
 
    ! Find the index of the first scheduled time, and whether start_time is
    ! itself one of the scheduled times.
    this%first_index = 0
    this%start_is_scheduled = .true.
    this%start_pending = .false.
 
    if (this%time_interval .gt. 0.0_dp) then
       offset = this%direction * (start_time - this%anchor_time) / &
            this%time_interval
       this%first_index = ceiling(offset - span_tol * max(1.0_dp, &
            abs(offset)), kind = i8)
       this%start_is_scheduled = abs(real(this%first_index, dp) - offset) &
            .le. span_tol * max(1.0_dp, abs(offset))
       if (this%start_is_scheduled .and. .not. this%write_at_start) then
          this%first_index = this%first_index + 1_i8
       end if
       this%start_pending = this%write_at_start .and. &
            .not. this%start_is_scheduled
    end if
 
    if (this%nsteps .gt. 0 .and. .not. this%write_at_start) then
       ! The first scheduled step of a step based schedule is the first step
       ! of the run. Skip it, as start_time is skipped above.
       this%next_index = 1
    else
       this%next_index = 0
    end if
 
  end subroutine time_based_controller_init
 
  subroutine time_based_controller_free(this)
    class(time_based_controller_t), intent(inout) :: this
 
    if (allocated(this%control_mode)) then
       deallocate(this%control_mode)
    end if
 
    this%frequency = 0.0_dp
    this%time_interval = 0.0_dp
    this%nsteps = 0
    this%start_time = 0.0_dp
    this%end_time = 0.0_dp
    this%nexecutions = 0
    this%never = .false.
    this%control_value = 0.0_dp
    this%write_at_start = .true.
    this%anchor_time = 0.0_dp
    this%first_index = 0
    this%start_is_scheduled = .false.
    this%start_pending = .false.
    this%next_index = 0
    this%direction = 1.0_dp
    this%tstep_offset = 0
    this%last_tstep = -1
  end subroutine time_based_controller_free
 
  pure function time_based_controller_check(this, time, force) &
       result(check)
    class(time_based_controller_t), intent(in) :: this
    type(time_state_t), intent(in) :: time
    logical, intent(in), optional :: force
    logical :: check
    logical :: ifforce
    real(kind=dp) :: progress, tol, t_next, t_start, t_end
    integer :: nstep
 
    if (present(force)) then
       ifforce = force
    else
       ifforce = .false.
    end if
 
    check = .false.
 
    ! At most one execution per time step.
    if (this%last_tstep .eq. time%tstep) return
 
    ! Nothing is scheduled, but an execution can still be forced.
    if (this%never .and. .not. ifforce) return
 
    progress = this%direction * (time%t - this%anchor_time)
    t_start = this%direction * (this%start_time - this%anchor_time)
    tol = this%tolerance(time%dt)
 
    ! Nothing is executed before start_time.
    if (progress .lt. t_start - tol) return
 
    if (ifforce) then
       check = .true.
    else if (this%nsteps .gt. 0) then
       nstep = time%tstep - this%tstep_offset
       check = nstep .ge. this%next_index * this%nsteps
    else if (this%start_pending) then
       ! The execution at start_time, which is not one of the scheduled
       ! times.
       check = .true.
    else
       t_end = this%direction * (this%end_time - this%anchor_time)
       t_next = real(this%first_index + this%next_index, dp) * &
            this%time_interval
       ! Nothing is scheduled after end_time. The condition is on the
       ! scheduled time and not on the current time, so a step that
       ! overshoots end_time still performs the execution scheduled for it.
       if (t_next .gt. t_end + span_tol * max(abs(t_end), &
            this%time_interval)) return
       check = progress .ge. t_next - tol
    end if
 
  end function time_based_controller_check
 
  pure function next_index_after(this, time) result(index)
    class(time_based_controller_t), intent(in) :: this
    type(time_state_t), intent(in) :: time
    integer :: index
    real(kind=dp) :: progress, tol
 
    if (this%nsteps .gt. 0) then
       index = (time%tstep - this%tstep_offset) / this%nsteps + 1
    else if (this%time_interval .gt. 0.0_dp) then
       progress = this%direction * (time%t - this%anchor_time)
       tol = this%tolerance(time%dt)
       index = int(floor((progress + tol) / this%time_interval, kind = i8) &
            - this%first_index) + 1
    else
       ! Nothing is scheduled.
       index = this%next_index
    end if
 
    ! The result is not limited from below by the current `next_index` on
    ! purpose: a scheduled execution always gives a larger index, and a
    ! forced execution before the next scheduled time must leave
    ! `next_index` unchanged.
    index = max(index, 0)
 
  end function next_index_after
 
  subroutine time_based_controller_register_execution(this, time)
    class(time_based_controller_t), intent(inout) :: this
    type(time_state_t), intent(in) :: time
 
    this%nexecutions = this%nexecutions + 1
    this%start_pending = .false.
    this%next_index = next_index_after(this, time)
    this%last_tstep = time%tstep
 
  end subroutine time_based_controller_register_execution
 
  subroutine time_based_controller_set_counter(this, time)
    class(time_based_controller_t), intent(inout) :: this
    type(time_state_t), intent(in) :: time
    real(kind=dp) :: progress, tol, dt, t_start, t_first
    integer :: n_passed
 
    if (this%never) return
 
    if (this%nsteps .gt. 0) then
       ! `tstep` is not stored in the checkpoint, so a step based schedule
       ! cannot be continued. Start a new one from the restart step, without
       ! executing at the restart step itself.
       this%tstep_offset = time%tstep
       this%next_index = 1
       this%start_pending = .false.
       this%last_tstep = time%tstep
       return
    end if
 
    ! The size of the step that produced the checkpoint, which is the step
    ! `check` used at that time, and not the step of the new run.
    dt = time%dt
    if (abs(time%dtlag(1)) .gt. 0.0_dp) dt = time%dtlag(1)
 
    progress = this%direction * (time%t - this%anchor_time)
    t_start = this%direction * (this%start_time - this%anchor_time)
    tol = this%tolerance(dt)
 
    if (progress .lt. t_start - tol) then
       ! The run has not reached start_time.
       this%next_index = 0
       this%nexecutions = 0
    else
       n_passed = int(floor((progress + tol) / this%time_interval, &
            kind = i8) - this%first_index) + 1
       this%next_index = max(n_passed, 0)
       this%nexecutions = this%next_index
       if (this%start_pending) then
          ! The execution at start_time was performed by the previous run.
          ! It is counted separately unless the first scheduled time was
          ! within tolerance of start_time, in which case the same execution
          ! was registered for both (see next_index_after).
          t_first = real(this%first_index, dp) * this%time_interval
          if (t_first .gt. t_start + tol) then
             this%nexecutions = this%nexecutions + 1
          end if
       end if
       this%start_pending = .false.
    end if
 
    ! Anything scheduled for the restart time was executed by the run that
    ! wrote the checkpoint.
    this%last_tstep = time%tstep
 
  end subroutine time_based_controller_set_counter
 
  pure function time_based_controller_tolerance(this, dt) result(tol)
    class(time_based_controller_t), intent(in) :: this
    real(kind=dp), intent(in) :: dt
    real(kind=dp) :: tol
 
    tol = time_tol * abs(dt)
    if (this%time_interval .gt. 0.0_dp) then
       tol = min(tol, time_tol * this%time_interval)
    end if
 
  end function time_based_controller_tolerance
 
  pure function time_based_controller_next_time(this) result(t)
    class(time_based_controller_t), intent(in) :: this
    real(kind=dp) :: t
 
    if (this%never .or. this%nsteps .gt. 0) then
       t = this%end_time
    else if (this%start_pending) then
       t = this%start_time
    else
       t = this%anchor_time + this%direction * &
            real(this%first_index + this%next_index, dp) * this%time_interval
    end if
 
  end function time_based_controller_next_time
 
  pure function time_based_controller_time_to_next(this, time) &
       result(t_to_next)
    class(time_based_controller_t), intent(in) :: this
    type(time_state_t), intent(in) :: time
    real(kind=dp) :: t_to_next
    real(kind=dp) :: progress, tol, t_next, t_end
    integer(kind=i8) :: k
 
    t_to_next = huge(0.0_dp)
 
    ! A step based or disabled control has no scheduled times: it executes
    ! after a number of steps, or only when forced, wherever the steps end.
    if (this%never .or. this%nsteps .gt. 0) return
 
    ! All times are measured as progress from the anchor in the direction of
    ! the run, so that they increase also for a run backwards in time. The
    ! tolerance only absorbs round-off in the accumulated time.
    progress = this%direction * (time%t - this%anchor_time)
    tol = span_tol * max(abs(progress), this%time_interval)
 
    if (this%start_pending) then
       ! The execution at start_time is still to come. It is not one of the
       ! scheduled times, e.g. for a component with a start_time of its own
       ! between two of them, so it is landed on separately.
       t_next = this%direction * (this%start_time - this%anchor_time)
    else
       ! The next scheduled time that has not been executed yet.
       t_next = real(this%first_index + this%next_index, dp) * &
            this%time_interval
    end if
 
    if (t_next .le. progress + tol) then
       ! The execution is due already but has not been performed, which
       ! happens at the start of the run for the simulation components: they
       ! are not executed before the time loop. It is performed at the end of
       ! the coming step, wherever that step ends, so the step is fitted to
       ! the first scheduled time after the current time instead.
       k = max(this%first_index, &
            floor((progress + tol) / this%time_interval, kind = i8) + 1_i8)
       t_next = real(k, dp) * this%time_interval
    end if
 
    ! A scheduled time after end_time is never executed. The end itself is
    ! landed on by the time step controller.
    t_end = this%direction * (this%end_time - this%anchor_time)
    if (t_next .gt. t_end + span_tol * max(abs(t_end), this%time_interval)) &
         return
 
    t_to_next = t_next - progress
 
  end function time_based_controller_time_to_next
 
end module time_based_controller