simulation_adjoint.f90

! 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_adjoint
  use mpi_f08, only: mpi_wtime
  use case, only: case_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 field, only: field_t
  use profiler, only: profiler_start, profiler_stop, &
       profiler_start_region, profiler_end_region
  use json_utils, only: json_get_or_default
  use time_state, only : time_state_t
  use time_step_controller, only: time_step_controller_t
  use adjoint_case, only: adjoint_case_t
  implicit none
  private
 
  public :: simulation_adjoint_init, simulation_adjoint_step, &
       simulation_adjoint_finalize, simulation_adjoint_restart
 
contains
 
 
  subroutine simulation_adjoint_init(C, dt_controller)
    type(adjoint_case_t), intent(inout) :: c
    type(time_step_controller_t), intent(inout) :: dt_controller
    character(len=LOG_SIZE) :: log_buf
 
    ! Write the initial logging message
    call neko_log%section('Adjoint 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)
 
    ! Call stats, samplers and user-init before time loop
    call neko_log%section('Postprocessing')
    call c%output_controller%execute(c%time)
 
    call c%case%user%initialize(c%time)
    call neko_log%end_section()
    call neko_log%newline()
 
  end subroutine simulation_adjoint_init
 
  subroutine simulation_adjoint_finalize(C)
    type(adjoint_case_t), intent(inout) :: c
    logical :: output_at_end
 
    ! Run a final output if specified in the json
    call json_get_or_default(c%case%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. c%time%t .lt. c%time%end_time) then
       call simulation_adjoint_joblimit_chkp(c, c%time%t)
    end if
 
    ! Finalize the user modules
    call c%case%user%finalize(c%time)
 
    call neko_log%end_section('Normal end.')
 
  end subroutine simulation_adjoint_finalize
 
  subroutine simulation_adjoint_step(C, dt_controller, cfl, tstep_loop_start_time)
    type(adjoint_case_t), intent(inout) :: c
    real(kind=rp), intent(inout) :: cfl
    type(time_step_controller_t), intent(inout) :: dt_controller
    real(kind=dp), intent(in) :: tstep_loop_start_time
    real(kind=dp) :: start_time, end_time, tstep_start_time
    character(len=LOG_SIZE) :: log_buf
 
    ! Setup the time step, and start time
    call profiler_start_region('Time-Step Adjoint')
    c%time%tstep = c%time%tstep + 1
    start_time = mpi_wtime()
    tstep_start_time = start_time
 
    ! Compute the next time step
    ! NOTE. we should be wary here since CFL is based on the convective velocity
    ! not the adjoint velocity
    cfl = c%fluid_adj%compute_cfl(c%time%dt)
    call dt_controller%set_dt(c%time, cfl)
    if (dt_controller%is_variable_dt) cfl = c%fluid_adj%compute_cfl(c%time%dt)
 
    ! Calculate the cfl after the possibly varied dt
    ! cfl = C%fluid_adj%compute_cfl(C%time%dt)
 
    ! Advance time step from t to t+dt and print the status
    call simulation_settime(c%time, c%fluid_adj%ext_bdf)
    call c%time%status()
    call neko_log%begin()
 
    write(log_buf, '(A,E15.7,1x,A,E15.7)') 'CFL:', cfl, 'dt:', c%time%dt
    call neko_log%message(log_buf)
 
    ! Scalar step
    ! (Note that for the adjoint we should the adjoint_scalars first)
    if (allocated(c%adjoint_scalars)) then
       start_time = mpi_wtime()
       call neko_log%section('Adjoint scalar')
       call c%adjoint_scalars%step(c%time, &
            c%case%fluid%ext_bdf, dt_controller)
       end_time = mpi_wtime()
       write(log_buf, '(A,E15.7)') &
            'Scalar step time:      ', end_time-start_time
       call neko_log%end_section(log_buf)
    end if
 
    ! Fluid step
    call neko_log%section('Adjoint fluid')
    call c%fluid_adj%step(c%time, dt_controller)
    end_time = mpi_wtime()
    write(log_buf, '(A,E15.7)') &
         'Fluid step time (s):   ', end_time-start_time
    call neko_log%end_section(log_buf)
 
    ! Postprocessing
    call neko_log%section('Postprocessing')
 
    ! 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, '(A,I8,A,E15.7)') &
         'Total time for step ', c%time%tstep, ' (s): ', end_time-tstep_start_time
    call neko_log%message(log_buf)
    write(log_buf, '(A,E15.7)') &
         'Total elapsed time (s):           ', end_time-tstep_loop_start_time
    call neko_log%message(log_buf)
 
    call neko_log%end_section()
    call neko_log%end()
    call profiler_end_region
 
  end subroutine simulation_adjoint_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
 
    if (allocated(ext_bdf)) then
       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 (ext_bdf%ndiff .eq. 0) then
          time%dtlag(2) = time%dt
          time%tlag(2) = time%t
       end if
 
       call ext_bdf%set_coeffs(time%dtlag)
    end if
 
    time%t = time%t + time%dt
 
  end subroutine simulation_settime
 
  subroutine simulation_adjoint_restart(C)
    type(adjoint_case_t), intent(inout) :: c
    integer :: i
    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
 
    call c%case%params%get('case.restart_file', restart_file, found)
    call c%case%params%get('case.restart_mesh_file', restart_mesh_file, found)
 
    if (found) then
       call previous_meshf%init(trim(restart_mesh_file))
       call previous_meshf%read(c%fluid_adj%chkp%previous_mesh)
    end if
 
    call c%case%params%get('case.mesh2mesh_tolerance', tol, found)
 
    if (found) c%case%fluid%chkp%mesh2mesh_tol = tol
 
    call chkpf%init(trim(restart_file))
    call chkpf%read(c%fluid_adj%chkp)
    c%time%dtlag = c%fluid_adj%chkp%dtlag
    c%time%tlag = c%fluid_adj%chkp%tlag
 
    ! Free the previous mesh, dont need it anymore
    do i = 1, size(c%time%dtlag)
       call c%case%fluid%ext_bdf%set_coeffs(c%time%dtlag)
    end do
 
    call c%fluid_adj%restart(c%case%chkp)
    call c%case%fluid%chkp%previous_mesh%free()
    if (allocated(c%adjoint_scalars)) then
       call c%adjoint_scalars%restart(c%case%chkp)
    end if
 
    c%time%t = real(c%case%fluid%chkp%restart_time(), kind=rp)
    call neko_log%section('Restarting from checkpoint')
    write(log_buf, '(A,A)') 'File :   ', trim(restart_file)
    call neko_log%message(log_buf)
    write(log_buf, '(A,E15.7)') 'Time : ', c%time%t
    call neko_log%message(log_buf)
    call neko_log%end_section()
 
    call c%output_controller%set_counter(c%time)
  end subroutine simulation_adjoint_restart
 
  subroutine simulation_adjoint_joblimit_chkp(C, t)
    type(adjoint_case_t), intent(inout) :: c
    real(kind=rp), 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 c%case%params%get('case.checkpoint_format', chkp_format, found)
    call c%case%fluid%chkp%sync_host()
    format_str = '.chkp'
    if (found) then
       if (chkp_format .eq. 'hdf5') then
          format_str = '.h5'
       end if
    end if
    call chkpf%init(c%case%output_directory // 'joblimit' // trim(format_str))
    call chkpf%write(c%case%fluid%chkp, t)
    write(log_buf, '(A)') '! saving checkpoint >>>'
    call neko_log%message(log_buf)
 
  end subroutine simulation_adjoint_joblimit_chkp
 
end module simulation_adjoint