simulation.f90

Go to the documentation of this file.

 
!
! Here, we simply march forward to steady state solutions
module simulation_m
  use case, only: case_t
  use user_access_singleton, only: neko_user_access
  use adjoint_case, only: adjoint_case_t
  use fluid_scheme_incompressible, only: fluid_scheme_incompressible_t
  use adjoint_fluid_scheme, only: adjoint_fluid_scheme_t
  use adjoint_fluid_pnpn, only: adjoint_fluid_pnpn_t
  use adjoint_scalars, only: adjoint_scalars_t
  use scalars, only: scalars_t
  use fluid_pnpn, only: fluid_pnpn_t
  use time_step_controller, only: time_step_controller_t
  use field_output, only: field_output_t
  use simcomp_executor, only: neko_simcomps
  use neko_ext, only: reset, reset_adjoint
  use json_file_module, only: json_file
  use json_utils, only: json_get, json_get_or_default
  use num_types, only: rp, sp, dp
  use mpi_f08, only: mpi_wtime
  use profiler, only: profiler_start, profiler_stop, &
       profiler_start_region, profiler_end_region
  use simulation_adjoint, only: simulation_adjoint_init, &
       simulation_adjoint_step, simulation_adjoint_finalize
  use simulation, only: simulation_init, simulation_step, simulation_finalize, &
       simulation_restart
  use state_recover, only: state_recover_t, state_recover_factory
  use runtime_stats, only: neko_rt_stats
  implicit none
  private
 
  type :: simulation_t
 
     type(case_t), public :: neko_case
     type(adjoint_case_t), public :: adjoint_case
     class(fluid_scheme_incompressible_t), public, pointer :: fluid => null()
     type(scalars_t), public, pointer :: scalars => null()
     class(adjoint_fluid_scheme_t), public, pointer :: adjoint_fluid => null()
     type(adjoint_scalars_t), public, pointer :: adjoint_scalars => null()
     type(field_output_t), public :: output_forward
     type(field_output_t), public :: output_adjoint
     character(len=:), allocatable :: forward_field_base_fname
     character(len=:), allocatable :: adjoint_field_base_fname
     integer :: current_design_iteration = 0
     logical :: unsteady = .false.
 
     logical :: have_scalar = .false.
     integer :: n_timesteps = 0
 
     ! ----------------------------------------------------------------------- !
     ! State recovery system
 
     class(state_recover_t), allocatable :: state_recover
 
   contains
     procedure, pass(this) :: init => simulation_initialize
     procedure, pass(this) :: free => simulation_free
     procedure, pass(this) :: run_forward => simulation_run_forward
     procedure, pass(this) :: run_backward => simulation_run_backward
     procedure, pass(this) :: reset => simulation_reset
     procedure, pass(this) :: set_output_counter => &
          simulation_set_output_counter
     procedure, pass(this) :: set_design_iteration => &
          simulation_set_design_iteration
     procedure, pass(this) :: write => simulation_write
     procedure, pass(this) :: write_forward => simulation_write_forward
     procedure, pass(this) :: write_adjoint => simulation_write_adjoint
 
  end type simulation_t
  public :: simulation_t
contains
 
  subroutine simulation_initialize(this, parameters)
    class(simulation_t), intent(inout), target :: this
    type(json_file), intent(inout) :: parameters
    type(json_file) :: state_recovery_params
    integer :: i, n_scalars
    character(len=:), allocatable :: output_directory, precision_s, file_format
    integer :: precision
    logical :: unsteady, subdivide
 
    ! initialize the primal Neko objects
    call this%neko_case%init(parameters)
    call neko_user_access%init(this%neko_case)
 
    call neko_rt_stats%init(parameters)
    call neko_simcomps%init(this%neko_case)
 
    ! initialize the adjoint
    call this%adjoint_case%init(this%neko_case)
 
    ! Capture Neko's own field output filenames before any design-iteration
    ! tag ever gets spliced into them (see `reset`/`neko_ext::reset`).
    this%forward_field_base_fname = &
         trim(this%neko_case%f_out%file_%get_base_fname())
    this%adjoint_field_base_fname = &
         trim(this%adjoint_case%f_out%file_%get_base_fname())
 
    ! Start the profiler
    call profiler_start
 
    select type (fluid => this%neko_case%fluid)
    type is (fluid_pnpn_t)
       this%fluid => fluid
    end select
 
    select type (adjoint_fluid => this%adjoint_case%fluid_adj)
    type is (adjoint_fluid_pnpn_t)
       this%adjoint_fluid => adjoint_fluid
    end select
 
    if (allocated(this%neko_case%scalars)) then
       this%scalars => this%neko_case%scalars
    end if
 
    if (allocated(this%adjoint_case%adjoint_scalars)) then
       this%adjoint_scalars => this%adjoint_case%adjoint_scalars
    end if
 
    !---------------------------------------------------------
    ! Initialize the output types
 
    ! Read settings for the output types, with some reasonable defaults
    call json_get_or_default(parameters, 'case.output_directory', &
         output_directory, '')
    call json_get_or_default(parameters, 'case.output_precision', precision_s, &
         'single')
    call json_get_or_default(parameters, 'case.fluid.output_format', &
         file_format, 'fld')
    call json_get_or_default(parameters, 'case.fluid.output_subdivide', &
         subdivide, .false.)
 
    if (trim(precision_s) .eq. 'double') then
       precision = dp
    else
       precision = sp
    end if
 
    ! Allocate the output type
    n_scalars = 0
    if (allocated(this%neko_case%scalars)) then
       n_scalars = size(this%neko_case%scalars%scalar_fields)
    end if
    call this%output_forward%init('forward_fields', 4 + n_scalars, &
         precision = precision, &
         path = trim(output_directory), &
         format = trim(file_format))
    call this%output_forward%file_%set_subdivide(subdivide)
 
    call this%output_forward%fields%assign(1, this%fluid%p)
    call this%output_forward%fields%assign(2, this%fluid%u)
    call this%output_forward%fields%assign(3, this%fluid%v)
    call this%output_forward%fields%assign(4, this%fluid%w)
 
    ! Assign all scalar fields
    if (allocated(this%neko_case%scalars)) then
       do i = 1, n_scalars
          call this%output_forward%fields%assign(4 + i, &
               this%scalars%scalar_fields(i)%scalar%s)
       end do
    end if
 
    ! Read settings for the output types, with some reasonable defaults
    call json_get_or_default(parameters, &
         'case.adjoint_fluid.output_format', file_format, 'fld')
    call json_get_or_default(parameters, &
         'case.adjoint_fluid.output_subdivide', subdivide, .false.)
 
    n_scalars = 0
    if (allocated(this%adjoint_case%adjoint_scalars)) then
       n_scalars = size(this%adjoint_case%adjoint_scalars%adjoint_scalar_fields)
    end if
    call this%output_adjoint%init('adjoint_fields', 4 + n_scalars, &
         precision = precision, &
         path = trim(output_directory), &
         format = trim(file_format))
    call this%output_adjoint%file_%set_subdivide(subdivide)
 
    call this%output_adjoint%fields%assign(1, this%adjoint_fluid%p_adj)
    call this%output_adjoint%fields%assign(2, this%adjoint_fluid%u_adj)
    call this%output_adjoint%fields%assign(3, this%adjoint_fluid%v_adj)
    call this%output_adjoint%fields%assign(4, this%adjoint_fluid%w_adj)
 
    ! Assign all scalar fields
    if (allocated(this%adjoint_case%adjoint_scalars)) then
       do i = 1, n_scalars
          call this%output_adjoint%fields%assign(4 + i, &
               this%adjoint_scalars%adjoint_scalar_fields(i)%s_adj)
       end do
    end if
 
    ! Check if the simulation is steady or unsteady
    call json_get_or_default(parameters, "unsteady", unsteady, .false.)
    this%unsteady = unsteady
 
    ! State recovery is only needed for unsteady runs.
    if (this%unsteady) then
       call json_get(parameters, 'state_recovery', state_recovery_params)
       call state_recover_factory(this%state_recover, this%neko_case, &
            state_recovery_params)
    end if
 
 
  end subroutine simulation_initialize
 
  subroutine simulation_free(this)
    class(simulation_t), intent(inout) :: this
 
    ! Stop the profiler
    call profiler_stop
 
    if (allocated(this%state_recover)) then
       call this%state_recover%free()
       deallocate(this%state_recover)
    end if
 
    ! Free the objects
    call this%neko_case%free()
    call this%adjoint_case%free()
    call this%output_forward%free()
    call this%output_adjoint%free()
 
    ! Nullify pointers
    nullify(this%fluid)
    nullify(this%scalars)
    nullify(this%adjoint_fluid)
    nullify(this%adjoint_scalars)
 
    ! Reset flags and counters
    this%unsteady = .false.
    this%have_scalar = .false.
    this%n_timesteps = 0
    this%current_design_iteration = 0
    if (allocated(this%forward_field_base_fname)) then
       deallocate(this%forward_field_base_fname)
    end if
    if (allocated(this%adjoint_field_base_fname)) then
       deallocate(this%adjoint_field_base_fname)
    end if
 
    ! Close global objects
    call neko_simcomps%free()
 
  end subroutine simulation_free
 
  subroutine simulation_run_forward(this)
    class(simulation_t), intent(inout) :: this
    type(time_step_controller_t) :: dt_controller
    real(kind=dp) :: loop_start
 
    call dt_controller%init(this%neko_case%params)
 
    call this%neko_case%time%reset()
    call simulation_init(this%neko_case, dt_controller)
 
    call profiler_start_region("Forward simulation")
    loop_start = mpi_wtime()
    this%n_timesteps = 0
    do while (.not. this%neko_case%time%is_done())
       this%n_timesteps = this%n_timesteps + 1
 
       call simulation_step(this%neko_case, dt_controller, loop_start)
 
       if (this%unsteady) then
          call this%state_recover%save()
       end if
    end do
    call profiler_end_region("Forward simulation")
 
    call simulation_finalize(this%neko_case)
 
  end subroutine simulation_run_forward
 
  subroutine simulation_run_backward(this)
    class(simulation_t), intent(inout) :: this
    type(time_step_controller_t) :: dt_controller
    real(kind=dp) :: loop_start
    real(kind=rp) :: cfl
    integer :: i
 
    call dt_controller%init(this%neko_case%params)
 
    call simulation_adjoint_init(this%adjoint_case, dt_controller)
 
    call profiler_start_region("Adjoint simulation")
    cfl = this%adjoint_case%fluid_adj%compute_cfl(this%adjoint_case%time%dt)
    loop_start = mpi_wtime()
    do i = this%n_timesteps, 1, -1
       if (this%unsteady) then
          call this%state_recover%restore(i)
       end if
 
       call simulation_adjoint_step(this%adjoint_case, dt_controller, cfl, &
            loop_start)
    end do
    call profiler_end_region("Adjoint simulation")
 
    call simulation_adjoint_finalize(this%adjoint_case)
 
  end subroutine simulation_run_backward
 
  subroutine simulation_reset(this)
    class(simulation_t), intent(inout) :: this
 
    call reset(this%neko_case, this%current_design_iteration, &
         this%forward_field_base_fname)
    call reset_adjoint(this%adjoint_case, this%neko_case, &
         this%current_design_iteration, this%adjoint_field_base_fname)
    if (this%unsteady) then
       call this%state_recover%reset()
    end if
 
  end subroutine simulation_reset
 
  subroutine simulation_set_design_iteration(this, iteration)
    class(simulation_t), intent(inout) :: this
    integer, intent(in) :: iteration
 
    this%current_design_iteration = iteration
 
  end subroutine simulation_set_design_iteration
 
  subroutine simulation_set_output_counter(this, idx)
    class(simulation_t), intent(inout) :: this
    integer, intent(in) :: idx
 
    call this%output_forward%set_counter(idx)
    call this%output_adjoint%set_counter(idx)
 
  end subroutine simulation_set_output_counter
 
  subroutine simulation_write(this, idx)
    class(simulation_t), intent(inout) :: this
    integer, intent(in) :: idx
 
    call this%output_forward%sample(real(idx, kind=rp))
    call this%output_adjoint%sample(real(idx, kind=rp))
 
  end subroutine simulation_write
 
  subroutine simulation_write_forward(this, idx)
    class(simulation_t), intent(inout) :: this
    integer, intent(in) :: idx
 
    call this%output_forward%sample(real(idx, kind=rp))
 
  end subroutine simulation_write_forward
 
  subroutine simulation_write_adjoint(this, idx)
    class(simulation_t), intent(inout) :: this
    integer, intent(in) :: idx
 
    call this%output_adjoint%sample(real(idx, kind=rp))
 
  end subroutine simulation_write_adjoint
 
end module simulation_m