time_step_controller.f90

Go to the documentation of this file.

! Copyright (c) 2022, 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_step_controller
  use num_types, only : dp, i8
  use logger, only : neko_log, log_size
  use utils, only : neko_error
  use json_module, only : json_file
  use json_utils, only : json_get_or_default, json_get_or_lookup_or_default
  use time_state, only : time_state_t
  use math, only : abscmp
  use comm, only : neko_global_comm, is_mpmd
  use mpi_f08, only : mpi_min, mpi_max, mpi_in_place, mpi_allreduce, &
       mpi_double_precision, mpi_integer
  implicit none
  private
 
  real(kind=dp), public, parameter :: landing_tol = 1.0e-6_dp
 
  real(kind=dp), parameter :: max_landing_steps = 1.0e15_dp
 
  type, public :: time_step_controller_t
     logical :: is_variable_dt = .false.
     real(kind=dp) :: cfl_trg = 0.0_dp
     real(kind=dp) :: cfl_avg = 0.0_dp
     real(kind=dp) :: init_dt = huge(0.0_dp)
     real(kind=dp) :: max_dt = 0.0_dp
     real(kind=dp) :: min_dt = 0.0_dp
     integer :: max_update_frequency = 0
     integer :: min_update_frequency = 0
     integer :: dt_last_change = -1
     real(kind=dp) :: alpha = 0.0_dp 
     real(kind=dp) :: max_dt_increase_factor = 0.0_dp
     real(kind=dp) :: min_dt_decrease_factor = 0.0_dp
     real(kind=dp) :: dev_tol = 0.0_dp
     logical :: exact_output_time = .false.
     real(kind=dp) :: dt_previous = 0.0_dp
   contains
     procedure, pass(this) :: init => time_step_controller_init
     procedure, pass(this) :: set_dt => time_step_controller_set_dt
     procedure, pass(this) :: first_dt => time_step_controller_first_dt
     procedure, pass(this) :: land => time_step_controller_land
 
  end type time_step_controller_t
 
contains
 
  subroutine time_step_controller_init(this, params)
    class(time_step_controller_t), intent(inout) :: this
    type(json_file), intent(inout) :: params
    integer :: flags_any(2), flags_all(2), ierr
 
    this%dt_last_change = -1
    call json_get_or_default(params, 'variable_timestep', &
         this%is_variable_dt, .false.)
    if (this%is_variable_dt) then
       call json_get_or_lookup_or_default(params, 'target_cfl', &
            this%cfl_trg, 0.4_dp)
       call json_get_or_lookup_or_default(params, 'timestep', &
            this%init_dt, huge(0.0_dp))
       call json_get_or_lookup_or_default(params, 'max_timestep', &
            this%max_dt, huge(0.0_dp))
       call json_get_or_lookup_or_default(params, 'min_timestep', &
            this%min_dt, 0.0_dp)
       call json_get_or_lookup_or_default(params, 'max_update_frequency',&
            this%max_update_frequency, 0)
       call json_get_or_lookup_or_default(params, 'min_update_frequency',&
            this%min_update_frequency, huge(0))
       call json_get_or_lookup_or_default(params, 'cfl_running_avg_coeff', &
            this%alpha, 0.5_dp)
       call json_get_or_lookup_or_default(params, 'max_dt_increase_factor', &
            this%max_dt_increase_factor, 1.2_dp)
       call json_get_or_lookup_or_default(params, 'min_dt_decrease_factor', &
            this%min_dt_decrease_factor, 0.5_dp)
       call json_get_or_lookup_or_default(params, 'cfl_deviation_tolerance', &
            this%dev_tol, 0.2_dp)
    end if
 
    call json_get_or_default(params, 'exact_output_time', &
         this%exact_output_time, .false.)
    if (this%exact_output_time .and. .not. this%is_variable_dt) then
       call neko_error('exact_output_time requires variable_timestep')
    end if
 
    ! A variable time step takes a collective at every step in an MPMD run,
    ! and the landing on the scheduled times another one, so the coupled
    ! simulations have to agree on both, or their collectives would not
    ! match up.
    if (is_mpmd()) then
       flags_any = merge(1, 0, [this%is_variable_dt, this%exact_output_time])
       flags_all = flags_any
       call mpi_allreduce(mpi_in_place, flags_any, 2, mpi_integer, &
            mpi_max, neko_global_comm, ierr)
       call mpi_allreduce(mpi_in_place, flags_all, 2, mpi_integer, &
            mpi_min, neko_global_comm, ierr)
       if (any(flags_any .ne. flags_all)) then
          call neko_error('variable_timestep and exact_output_time must &
          &each be set in all the coupled cases of an MPMD run, or in &
          &none')
       end if
    end if
 
  end subroutine time_step_controller_init
 
  pure function time_step_controller_first_dt(this, time, cfl) result(dt)
    class(time_step_controller_t), intent(in) :: this
    type(time_state_t), intent(in) :: time
    real(kind=dp), intent(in) :: cfl
    real(kind=dp) :: dt
 
    dt = this%init_dt
    if (cfl .gt. 0.0_dp) dt = min(this%cfl_trg / cfl * time%dt, dt)
    dt = max(min(dt, this%max_dt), this%min_dt)
 
  end function time_step_controller_first_dt
 
  subroutine time_step_controller_set_dt(this, time, cfl)
    class(time_step_controller_t), intent(inout) :: this
    type(time_state_t), intent(inout) :: time
    real(kind=dp), intent(in) :: cfl
    real(kind=dp) :: dt_old, scaling_factor, global_min_dt
    character(len=LOG_SIZE) :: log_buf
    integer :: ierr
 
    this%dt_previous = time%dt
 
    ! Check if variable dt is requested
    if (.not. this%is_variable_dt) return
 
    ! Reset the average cfl if it is the first time step since the last change
    if (this%dt_last_change .eq. 0) then
       this%cfl_avg = cfl
    end if
 
    if (this%dt_last_change .eq. -1) then
 
       ! Set the first dt for desired cfl, or use the provided initial dt if it
       ! is smaller. Then clamp between max and min dt if provided.
       time%dt = this%first_dt(time, cfl)
       this%dt_last_change = 0
       this%cfl_avg = cfl
       ! There is no step before the first one
       this%dt_previous = time%dt
 
    else
       ! Calculate the average of cfl over the desired interval
       this%cfl_avg = this%alpha * cfl + (1 - this%alpha) * this%cfl_avg
 
       if (abs(this%cfl_avg - this%cfl_trg) .ge. this%dev_tol * this%cfl_trg &
            .and. this%dt_last_change .ge. this%max_update_frequency &
            .or. this%dt_last_change .ge. this%min_update_frequency) then
 
          ! Test for an exact zero cfl first, so that it is never divided by.
          if (abscmp(cfl, 0.0_dp)) then
             ! no limit on dt, so increase it as much as allowed
             scaling_factor = this%max_dt_increase_factor
          else if (this%cfl_trg/cfl .ge. 1) then
             ! increase of time step
             scaling_factor = min(this%max_dt_increase_factor, this%cfl_trg/cfl)
          else
             ! reduction of time step
             scaling_factor = max(this%min_dt_decrease_factor, this%cfl_trg/cfl)
          end if
 
          dt_old = time%dt
          time%dt = scaling_factor * dt_old
          time%dt = max(min(time%dt, this%max_dt), this%min_dt)
 
          write(log_buf, '(A,E15.7,1x,A,E15.7)') &
               'Average CFL:', this%cfl_avg, &
               'Target  CFL:', this%cfl_trg
          call neko_log%message(log_buf)
 
          write(log_buf, '(A,E15.7,1x,A,E15.7)') 'Old dt:', dt_old, &
               'New dt:', time%dt
          call neko_log%message(log_buf)
 
          this%dt_last_change = 0
 
       else
          this%dt_last_change = this%dt_last_change + 1
       end if
    end if
 
    ! If running in mpmd, the new dt is the minimum across simulations
    if (is_mpmd()) then
       global_min_dt = time%dt
       call mpi_allreduce(mpi_in_place, global_min_dt, 1, &
            mpi_double_precision, mpi_min, neko_global_comm, ierr)
 
       ! If my dt is larger that the global min, mark a change
       if (time%dt .gt. global_min_dt) then
          time%dt = global_min_dt
          this%dt_last_change = 0
       end if
    end if
 
 
  end subroutine time_step_controller_set_dt
 
  subroutine time_step_controller_land(this, time, time_to_next)
    class(time_step_controller_t), intent(inout) :: this
    type(time_state_t), intent(inout) :: time
    real(kind=dp), intent(in) :: time_to_next
    real(kind=dp) :: direction, dt_prev, lo, remaining, remaining_end, x
    integer(kind=i8) :: n
    integer :: ierr
 
    if (.not. (this%exact_output_time .and. this%is_variable_dt)) return
 
    direction = sign(1.0_dp, time%dt)
    dt_prev = abs(this%dt_previous)
    lo = this%min_dt_decrease_factor * dt_prev
 
    ! The end, or the next scheduled time unless it is too close
    remaining_end = direction * (time%end_time - time%t)
    if (remaining_end .le. time_to_next + lo) then
       remaining = remaining_end
    else if (time_to_next .ge. lo) then
       remaining = time_to_next
    else
       remaining = 0.0_dp
    end if
 
    x = remaining / abs(time%dt)
    if (remaining .gt. 0.0_dp .and. x .lt. max_landing_steps) then
       ! Candidate count: x proposed steps fit into the remaining time,
       ! round down, but take at least one step
       n = max(1_i8, floor(x, kind = i8))
       ! Round up instead if x is closer to n + 1 in ratio, that is, if
       ! (n + 1) / x < x / n. The boundary is the geometric mean
       ! sqrt(n (n + 1)), not n + 0.5, so that a step twice too long and one
       ! half as long count as equally far off
       if (x * x .gt. real(n, dp) * real(n + 1_i8, dp)) n = n + 1_i8
       ! Lower bound on the step, min_dt_decrease_factor times the previous
       ! one: at most remaining / lo steps, and at least one. The tolerance
       ! keeps round-off from dropping a count that fits exactly
       n = min(n, max(1_i8, floor(remaining / lo + landing_tol, kind = i8)))
       ! Upper bound on the step, max_dt_increase_factor times the previous
       ! one or max_timestep, whichever is smaller: at least that many
       ! steps. Applied last, so it wins if the two bounds conflict
       n = max(n, ceiling(remaining / min(this%max_dt_increase_factor * &
            dt_prev, this%max_dt) - landing_tol, kind = i8))
       ! The fitted step: n equal steps end exactly on the target
       time%dt = direction * remaining / real(n, dp)
    end if
 
    ! In an MPMD run the simulations advance in lockstep
    if (is_mpmd()) then
       call mpi_allreduce(mpi_in_place, time%dt, 1, mpi_double_precision, &
            mpi_min, neko_global_comm, ierr)
    end if
 
    ! The projection spaces watch for a change of the step taken
    if (abs(abs(time%dt) - dt_prev) .gt. landing_tol * dt_prev) then
       this%dt_last_change = 0
    end if
 
  end subroutine time_step_controller_land
 
end module time_step_controller