simulation_adjoint.f90

Go to the documentation of this file.

 
!
module simulation_adjoint
  use mpi_f08, only: mpi_wtime
  use neko_config, only: neko_bcknd_device
  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 profiler, only: 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
  use scratch_registry, only: neko_scratch_registry
  use device_math, only: device_glsc3
  use math, only: glsc3
  use vector, only: vector_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 simulation_adjoint_norm_output(c, 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, final_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=rp), optional, intent(in) :: final_time
    real(kind=rp) :: t_bkp
    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)
    ! for cosmetic reasons we want the simulation to run backwards
    if (present(final_time)) then
       t_bkp = c%time%t
       c%time%t = final_time - t_bkp
    end if
    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')
 
    ! Correct the time so the output fields are the same as the primal
    if (present(final_time)) then
       c%time%t = t_bkp
    end if
 
    ! Run any IO needed.
    call c%output_controller%execute(c%time)
    call simulation_adjoint_norm_output(c, 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)
    if (c%norm_output_enabled) then
       call c%norm_output_ctrl%set_counter(c%time)
    end if
  end subroutine simulation_adjoint_restart
 
  subroutine simulation_adjoint_norm_output(C, time_output)
    type(adjoint_case_t), intent(inout) :: c
    type(time_state_t), intent(in) :: time_output
    type(vector_t), pointer :: data_line
    real(kind=rp) :: norm_l2
    integer :: n, idx
 
    if (.not. c%norm_output_enabled) return
    if (.not. c%norm_output_ctrl%check(time_output)) return
 
    n = c%fluid_adj%c_Xh%dof%size()
    if (neko_bcknd_device .eq. 1) then
       norm_l2 = device_glsc3(c%fluid_adj%u_adj%x_d, &
            c%fluid_adj%u_adj%x_d, c%fluid_adj%c_Xh%B_d, n) + &
            device_glsc3(c%fluid_adj%v_adj%x_d, &
            c%fluid_adj%v_adj%x_d, c%fluid_adj%c_Xh%B_d, n) + &
            device_glsc3(c%fluid_adj%w_adj%x_d, &
            c%fluid_adj%w_adj%x_d, c%fluid_adj%c_Xh%B_d, n)
    else
       norm_l2 = glsc3(c%fluid_adj%u_adj%x, c%fluid_adj%u_adj%x, &
            c%fluid_adj%c_Xh%B, n) + &
            glsc3(c%fluid_adj%v_adj%x, c%fluid_adj%v_adj%x, &
            c%fluid_adj%c_Xh%B, n) + &
            glsc3(c%fluid_adj%w_adj%x, c%fluid_adj%w_adj%x, &
            c%fluid_adj%c_Xh%B, n)
    end if
 
    norm_l2 = sqrt(norm_l2) / c%fluid_adj%c_Xh%volume
 
    call neko_scratch_registry%request(data_line, idx, 1, .false.)
    data_line%x = [norm_l2]
    call c%norm_output_file%write(data_line, time_output%t)
    call neko_scratch_registry%relinquish(idx)
    call c%norm_output_ctrl%register_execution()
  end subroutine simulation_adjoint_norm_output
 
  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