simulation.f90

Go to the documentation of this file.

! Copyright (c) 2020-2021, 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 simulation
  use mpi_f08, only : mpi_wtime
  use case, only : case_t
  use checkpoint, only : chkp_t
  use num_types, only : rp, dp
  use time_scheme_controller, only : time_scheme_controller_t
  use file, only : file_t
  use logger, only : log_size, neko_log
  use jobctrl, only : jobctrl_time_limit
  use profiler, only : profiler_start, profiler_stop, &
       profiler_start_region, profiler_end_region
  use simcomp_executor, only : neko_simcomps
  use json_utils, only : json_get, json_get_or_default
  use time_state, only : time_state_t
  use time_step_controller, only : time_step_controller_t
  implicit none
  private
 
  interface simulation_restart
     module procedure case_restart_from_parameters, &
          case_restart_from_checkpoint
  end interface simulation_restart
 
  public :: simulation_init, simulation_step, simulation_finalize, &
       simulation_restart
 
contains
 
  subroutine simulation_init(C, dt_controller)
    type(case_t), intent(inout) :: c
    type(time_step_controller_t), intent(inout) :: dt_controller
    character(len=LOG_SIZE) :: log_buf
    logical :: found
    character(len=:), allocatable :: restart_file
 
    ! Restart the case if needed
    call c%params%get('case.restart_file', restart_file, found)
    if (found .and. len_trim(restart_file) .gt. 0) then
       ! Restart the case
       call simulation_restart(c)
    end if
 
    ! Write the initial logging message
    call neko_log%section('Starting simulation')
    write(log_buf, '(A, E15.7,A,E15.7,A)') &
         'T  : [', c%time%t, ', ', c%time%end_time, ']'
    call neko_log%message(log_buf)
    if (.not. dt_controller%is_variable_dt) then
       write(log_buf, '(A, E15.7)') 'dt :  ', c%time%dt
    else
       write(log_buf, '(A, E15.7)') 'CFL :  ', dt_controller%cfl_trg
    end if
    call neko_log%message(log_buf)
    if (dt_controller%exact_output_time) then
       call neko_log%message('dt is fitted to the output times')
    end if
 
    ! Execute outputs and user-init before time loop
    call neko_log%section('Preprocessing')
    call c%user%initialize(c%time)
 
    ! With a variable time step, set the first step before the initial
    ! output.
    if (dt_controller%is_variable_dt) then
       c%time%dt = dt_controller%first_dt(c%time, &
            c%fluid%compute_cfl(c%time%dt))
    end if
    call c%output_controller%execute(c%time)
 
    call neko_log%end_section()
    call neko_log%newline()
 
  end subroutine simulation_init
 
  subroutine simulation_finalize(C)
    type(case_t), intent(inout) :: c
    logical :: output_at_end
 
    ! Run a final output if specified in the json
    call json_get_or_default(c%params, 'case.output_at_end', &
         output_at_end, .true.)
    call c%output_controller%execute(c%time, output_at_end)
 
    if (.not. (output_at_end) .and. .not. c%time%is_done()) then
       call simulation_joblimit_chkp(c, c%time%t)
    end if
 
    ! Finalize the user modules
    call c%user%finalize(c%time)
 
    call neko_log%end_section('Normal end.')
 
  end subroutine simulation_finalize
 
  subroutine simulation_step(C, dt_controller, tstep_loop_start_time)
    type(case_t), intent(inout) :: c
    type(time_step_controller_t), intent(inout) :: dt_controller
    real(kind=dp), optional, intent(in) :: tstep_loop_start_time
    real(kind=dp) :: start_time, end_time, tstep_start_time
    real(kind=dp) :: cfl, dt_previous
    character(len=LOG_SIZE) :: log_buf
 
    ! Setup the time step, and start time
    call profiler_start_region('Time-Step')
    start_time = mpi_wtime()
    tstep_start_time = start_time
 
    ! Compute the next time step size
    dt_previous = c%time%dt
    cfl = c%fluid%compute_cfl(c%time%dt)
    call dt_controller%set_dt(c%time, cfl)
 
    ! Shorten the step to land exactly on the next sampling or output time
    if (dt_controller%exact_output_time) then
       call dt_controller%land(c%time, &
            min(c%output_controller%time_to_next(c%time), &
            neko_simcomps%time_to_next(c%time)))
    end if
 
    ! The CFL number of the step about to be taken, for the log
    if (c%time%dt .ne. dt_previous) cfl = c%fluid%compute_cfl(c%time%dt)
 
    ! Advance time step from t to t+dt and print the status
    call simulation_settime(c%time, c%fluid%ext_bdf)
    call c%time%status()
    call neko_log%begin()
 
    call neko_log%section('Preprocessing')
 
    write(log_buf, "(A,F8.4,2x,A,1P,E14.7)") 'CFL:', cfl, 'dt:', c%time%dt
    call neko_log%message(log_buf)
 
    ! Run the preprocessing
    call c%user%preprocess(c%time)
    call neko_simcomps%preprocess(c%time)
    call neko_log%end_section()
 
    ! Fluid step
    call neko_log%section('Fluid')
    start_time = mpi_wtime()
    call c%fluid%step(c%time, dt_controller)
    end_time = mpi_wtime()
    write(log_buf, '(A,3X,E15.7)') 'Fluid step time (s):', end_time - start_time
    call neko_log%end_section(log_buf)
 
    ! Scalar step
    if (allocated(c%scalars)) then
       start_time = mpi_wtime()
       call neko_log%section('Scalar')
       call c%scalars%step(c%time, c%fluid%ext_bdf, dt_controller)
       end_time = mpi_wtime()
       write(log_buf, '(A,6X,E15.7)') 'Scalar step time:', end_time - start_time
       call neko_log%end_section(log_buf)
    end if
 
    ! Postprocessing
    call neko_log%section('Postprocessing')
 
    ! Execute all simulation components
    call neko_simcomps%compute(c%time)
 
    ! Run the user compute routine
    call c%user%compute(c%time)
 
    ! Run any IO needed.
    call c%output_controller%execute(c%time)
 
    call neko_log%end_section()
 
    ! End the step and print summary
    end_time = mpi_wtime()
    call neko_log%section('Step summary')
    write(log_buf, '(A20,I8,A6,E15.7)') &
         'Total time for step ', c%time%tstep, ' (s): ', &
         end_time-tstep_start_time
    call neko_log%message(log_buf)
 
    if (present(tstep_loop_start_time)) then
       write(log_buf, '(A34,E15.7)') &
            'Total elapsed time (s):           ', &
            end_time - tstep_loop_start_time
       call neko_log%message(log_buf)
    end if
 
    call neko_log%end_section()
    call neko_log%end()
    call profiler_end_region
 
  end subroutine simulation_step
 
  subroutine simulation_settime(time, ext_bdf)
    type(time_state_t), intent(inout) :: time
    type(time_scheme_controller_t), intent(inout), allocatable :: ext_bdf
    integer :: i
 
    ! Always shift the time lag arrays (needed by entropy viscosity etc.)
    do i = 10, 2, -1
       time%tlag(i) = time%tlag(i-1)
       time%dtlag(i) = time%dtlag(i-1)
    end do
 
    time%dtlag(1) = time%dt
    time%tlag(1) = time%t
 
    if (allocated(ext_bdf)) then
       if (ext_bdf%ndiff .eq. 0) then
          time%dtlag(2) = time%dt
          time%tlag(2) = time%t
       end if
 
       call ext_bdf%set_coeffs(real(time%dtlag, kind=rp))
    end if
 
    time%tstep = time%tstep + 1
    time%t = time%t + time%dt
 
  end subroutine simulation_settime
 
  subroutine case_restart_from_parameters(C)
    type(case_t), intent(inout) :: C
    type(file_t) :: chkpf, previous_meshf
    character(len=LOG_SIZE) :: log_buf
    character(len=:), allocatable :: restart_file
    character(len=:), allocatable :: restart_mesh_file
    real(kind=rp) :: tol
    logical :: found, check_cont
    integer :: i
 
    call json_get(c%params, 'case.restart_file', restart_file)
    call json_get_or_default(c%params, 'case.restart_mesh_file', &
         restart_mesh_file, "")
 
    if (restart_mesh_file .ne. "") then
       call previous_meshf%init(trim(restart_mesh_file))
       call previous_meshf%read(c%chkp%previous_mesh)
 
       call json_get_or_default(c%params, 'case.mesh2mesh_tolerance', &
            c%chkp%mesh2mesh_tol, 1e-6_dp)
    end if
 
    call neko_log%section('Restarting from checkpoint')
    write(log_buf, '(A,A)') 'File :   ', trim(restart_file)
    call neko_log%message(log_buf)
 
    call chkpf%init(trim(restart_file))
    call chkpf%read(c%chkp)
 
    call case_restart_from_checkpoint(c, c%chkp)
 
    ! Free the previous mesh
    call c%chkp%previous_mesh%free()
 
    write(log_buf, '(A,E15.7)') 'Time : ', c%time%t
    call neko_log%message(log_buf)
    call neko_log%end_section()
 
  end subroutine case_restart_from_parameters
 
  subroutine case_restart_from_checkpoint(C, chkp)
    type(case_t), intent(inout) :: C
    type(chkp_t), intent(inout) :: chkp
    character(len=LOG_SIZE) :: log_buf
    integer :: i
 
    ! Restart the time state and BDF coefficients
    call chkp%set_time_state(c%time)
    do i = 1, size(c%time%dtlag)
       call c%fluid%ext_bdf%set_coeffs(real(c%time%dtlag, kind=rp))
    end do
 
    ! Restart the fluid and scalars
    call c%fluid%restart(chkp)
    if (allocated(c%scalars)) call c%scalars%restart(chkp)
 
    ! Restart the output controller
    call c%output_controller%set_counter(c%time)
 
    ! Restart the simulation components
    call neko_simcomps%restart(c%time)
 
  end subroutine case_restart_from_checkpoint
 
  subroutine simulation_joblimit_chkp(C, t)
    type(case_t), intent(inout) :: C
    real(kind=dp), intent(inout) :: t
    type(file_t) :: chkpf
    character(len=:), allocatable :: chkp_format
    character(len=LOG_SIZE) :: log_buf
    character(len=10) :: format_str
    logical :: found
 
    call json_get_or_default(c%params, 'case.checkpoint_format', chkp_format, &
         'default')
    call c%chkp%sync_host()
    if (chkp_format .eq. 'hdf5') then
       format_str = '.h5'
    else
       format_str = '.chkp'
    end if
    call chkpf%init(c%output_directory // 'joblimit' // trim(format_str))
    write(log_buf, '(A)') '! saving checkpoint >>>'
    call neko_log%message(log_buf)
    call neko_log%flush()
    call chkpf%write(c%chkp, t)
    call neko_log%flush()
 
  end subroutine simulation_joblimit_chkp
 
end module simulation