problem.f90

Go to the documentation of this file.

 
 
module problem
  use num_types, only: rp, dp
  use design, only: design_t
  use objective, only: objective_t, objective_wrapper_t, objective_factory
  use constraint, only: constraint_t, constraint_wrapper_t, constraint_factory
  use augmented_lagrangian_objective, only: augmented_lagrangian_objective_t
  use vector, only: vector_t
  use matrix, only: matrix_t
  use device, only: host_to_device, device_to_host
  use json_module, only: json_file
  use json_utils, only: json_extract_item, json_get, json_get_or_default
  use simulation_m, only: simulation_t
  use logger, only: neko_log
  use math, only: copy
  use vector_math, only: vector_add2, vector_cfill
  use time_step_controller, only: time_step_controller_t
  use simulation_adjoint, only: simulation_adjoint_init, &
       simulation_adjoint_step, simulation_adjoint_finalize
  use simulation, only: simulation_init, simulation_step, simulation_finalize
  use mpi_f08, only: mpi_wtime
  use profiler, only: profiler_start_region, profiler_end_region
 
  implicit none
  private
 
  type, public :: problem_t
     private
 
     integer :: n_design = 0
     integer :: n_objectives = 0
     integer :: n_constraints = 0
 
     class(objective_wrapper_t), allocatable, dimension(:) :: objective_list
     class(constraint_wrapper_t), allocatable, dimension(:) :: constraint_list
 
   contains
 
     ! ----------------------------------------------------------------------- !
     ! Interfaces
 
     procedure, pass(this), public :: init => problem_init
     procedure, pass(this), public :: free => problem_free
 
     procedure, pass(this), public :: compute => problem_compute
 
     procedure, pass(this), public :: compute_sensitivity => &
          problem_compute_sensitivity
     procedure, pass(this), public :: run_forward_unsteady => &
          problem_run_forward_unsteady
     procedure, pass(this), public :: run_backward_unsteady => &
          problem_run_backward_unsteady
     ! ----------------------------------------------------------------------- !
     ! Base class methods
 
     procedure, pass(this), public :: read_objectives => problem_read_objectives
     procedure, pass(this), public :: read_constraints => &
          problem_read_constraints
 
     ! ----------------------------------------------------------------------- !
     ! Actual methods
 
     procedure, pass(this), public :: write => problem_write
 
     procedure, pass(this), public :: add_objective => problem_add_objective
     procedure, pass(this), public :: add_constraint => problem_add_constraint
 
     ! ----------------------------------------------------------------------- !
     ! Internal Updater methods
 
     procedure, pass(this) :: update_objectives => &
          problem_update_objectives
     procedure, pass(this) :: update_constraints => &
          problem_update_constraints
     procedure, pass(this) :: update_objective_sensitivities => &
          problem_update_objective_sensitivities
     procedure, pass(this) :: update_constraint_sensitivities => &
          problem_update_constraint_sensitivities
 
     procedure, pass(this) :: reset_objectives => &
          problem_reset_objectives
     procedure, pass(this) :: reset_constraints => &
          problem_reset_constraints
     procedure, pass(this) :: reset_objective_sensitivities => &
          problem_reset_objective_sensitivities
     procedure, pass(this) :: reset_constraint_sensitivities => &
          problem_reset_constraint_sensitivities
 
     procedure, pass(this) :: accumulate_objectives => &
          problem_accumulate_objectives
     procedure, pass(this) :: accumulate_constraints => &
          problem_accumulate_constraints
     procedure, pass(this) :: accumulate_objective_sensitivities => &
          problem_accumulate_objective_sensitivities
     procedure, pass(this) :: accumulate_constraint_sensitivities => &
          problem_accumulate_constraint_sensitivities
 
     ! ----------------------------------------------------------------------- !
     ! Public Getters
 
     procedure, pass(this), public :: get_objective_value => &
          problem_get_objective_value
     procedure, pass(this), public :: get_all_objective_values => &
          problem_get_all_objective_values
     procedure, pass(this), public :: get_constraint_values => &
          problem_get_constraint_values
     procedure, pass(this), public :: get_objective_sensitivities => &
          problem_get_objective_sensitivities
     procedure, pass(this), public :: get_constraint_sensitivities => &
          problem_get_constraint_sensitivities
 
     procedure, pass(this) :: get_n_objectives => problem_get_num_objectives
     procedure, pass(this) :: get_n_constraints => problem_get_num_constraints
 
     procedure, pass(this) :: get_log_header => problem_get_log_header
 
  end type problem_t
 
contains
 
  ! ========================================================================== !
  ! Base class methods
 
  subroutine problem_init(this, parameters, design, simulation)
    class(problem_t), intent(inout) :: this
    type(json_file), intent(inout) :: parameters
    class(design_t), intent(in) :: design
    type(simulation_t), optional, intent(inout) :: simulation
 
    call this%free()
 
    this%n_design = design%size()
 
    ! Read the objectives and constraints
    call this%read_objectives(parameters, design, simulation)
    call this%read_constraints(parameters, design, simulation)
 
  end subroutine problem_init
 
  subroutine problem_free(this)
    class(problem_t), intent(inout) :: this
    integer :: i
 
    this%n_design = 0
    this%n_objectives = 0
    this%n_constraints = 0
 
    ! Free the objective list
    if (allocated(this%objective_list)) then
       do i = 1, size(this%objective_list)
          call this%objective_list(i)%free()
       end do
       deallocate(this%objective_list)
    end if
 
    ! Free the constraint list
    if (allocated(this%constraint_list)) then
       do i = 1, size(this%constraint_list)
          call this%constraint_list(i)%free()
       end do
       deallocate(this%constraint_list)
    end if
  end subroutine problem_free
 
  subroutine problem_write(this, idx)
    class(problem_t), intent(inout) :: this
    integer, intent(in) :: idx
 
  end subroutine problem_write
 
  ! ========================================================================== !
  ! Handling constraints and objectives
 
  subroutine problem_read_objectives(this, parameters, design, simulation)
    class(problem_t), intent(inout) :: this
    type(json_file), intent(inout) :: parameters
    class(design_t), intent(in) :: design
    type(simulation_t), optional, intent(inout) :: simulation
    class(objective_t), allocatable :: objective
 
    ! A single objective term as its own json_file.
    character(len=:), allocatable :: path, type
    type(json_file) :: objective_json
    integer :: n_objectives, i
    logical :: dealias
 
    call neko_log%section("Reading objectives")
 
    ! Get the number of objectives.
    path = "optimization.objectives"
    if (parameters%valid_path(path)) then
       call parameters%info(path, n_children = n_objectives)
 
       ! Grab a single parameters entry and create a constraint from it.
       do i = 1, n_objectives
          call json_extract_item(parameters, path, i, objective_json)
          call json_get(objective_json, "type", type)
          call neko_log%message(type)
 
          call objective_factory(objective, objective_json, design, simulation)
          call this%add_objective(objective)
       end do
    end if
 
    if (present(simulation)) then
       if (allocated(objective)) deallocate(objective)
       allocate(augmented_lagrangian_objective_t::objective)
       select type(alo => objective)
       class is (augmented_lagrangian_objective_t)
          call json_get_or_default(parameters, &
               "adjoint_fluid.dealias_sensitivity", dealias, .true.)
          call alo%init_from_attributes(design, simulation, weight = 1.0_rp, &
               name = "Augmented Lagrangian", mask_name = "", &
               dealias = dealias)
       end select
       call this%add_objective(objective)
    end if
 
    call neko_log%end_section()
 
  end subroutine problem_read_objectives
 
  subroutine problem_read_constraints(this, parameters, design, simulation)
    class(problem_t), intent(inout) :: this
    type(json_file), intent(inout) :: parameters
    class(design_t), intent(in) :: design
    class(constraint_t), allocatable :: constraint
    type(simulation_t), optional, intent(inout) :: simulation
 
    ! A single constraint term as its own json_file.
    character(len=:), allocatable :: path, type
    type(json_file) :: constraint_json
    integer :: n_constraints, i
 
    call neko_log%section("Reading constraints")
 
    ! Get the number of constraints.
    path = "optimization.constraints"
 
    if (parameters%valid_path(path)) then
       call parameters%info(path, n_children = n_constraints)
 
       ! Grab a single parameters entry and create a constraint from it.
       do i = 1, n_constraints
          call json_extract_item(parameters, path, i, constraint_json)
          call json_get(constraint_json, "type", type)
          call neko_log%message(type)
 
          call constraint_factory(constraint, constraint_json, design, simulation)
          call this%add_constraint(constraint)
       end do
    end if
 
    call neko_log%end_section()
 
  end subroutine problem_read_constraints
 
  subroutine problem_add_objective(this, objective)
    class(problem_t), intent(inout) :: this
    class(objective_t), allocatable, intent(inout) :: objective
    class(objective_wrapper_t), allocatable, dimension(:) :: temp_list
    integer :: i, n
 
    n = 0
    if (allocated(this%objective_list)) then
       n = size(this%objective_list)
       call move_alloc(this%objective_list, temp_list)
       allocate(this%objective_list(n + 1))
       if (allocated(temp_list)) then
          do i = 1, n
             call move_alloc(temp_list(i)%objective, &
                  this%objective_list(i)%objective)
          end do
       end if
    else
       allocate(this%objective_list(1))
    end if
 
    call move_alloc(objective, this%objective_list(n + 1)%objective)
    this%n_objectives = n + 1
  end subroutine problem_add_objective
 
  subroutine problem_add_constraint(this, constraint)
    class(problem_t), intent(inout) :: this
    class(constraint_t), allocatable, intent(inout) :: constraint
    class(constraint_wrapper_t), allocatable, dimension(:) :: temp_list
    integer :: i, n
 
    n = 0
    if (allocated(this%constraint_list)) then
       n = size(this%constraint_list)
       call move_alloc(this%constraint_list, temp_list)
       allocate(this%constraint_list(n + 1))
       if (allocated(temp_list)) then
          do i = 1, n
             call move_alloc(temp_list(i)%constraint, &
                  this%constraint_list(i)%constraint)
          end do
       end if
    else
       allocate(this%constraint_list(1))
    end if
 
    call move_alloc(constraint, this%constraint_list(n + 1)%constraint)
    this%n_constraints = n + 1
  end subroutine problem_add_constraint
 
  ! ========================================================================== !
  ! Problem part computation
 
  subroutine problem_compute(this, design, simulation)
    class(problem_t), intent(inout) :: this
    class(design_t), intent(inout) :: design
    class(simulation_t), optional, intent(inout) :: simulation
 
    if (present(simulation)) then
       call simulation%reset()
       if (simulation%unsteady) then
          ! Objective value accumulated
          call this%run_forward_unsteady(simulation, design)
       else
          call simulation%run_forward()
          ! Compute objective value on steady field
          call this%update_objectives(design)
       end if
    else
       call this%update_objectives(design)
    end if
 
    call this%update_constraints(design)
 
  end subroutine problem_compute
 
  subroutine problem_compute_sensitivity(this, design, simulation)
    class(problem_t), intent(inout) :: this
    class(design_t), intent(inout) :: design
    class(simulation_t), optional, intent(inout) :: simulation
 
    type(vector_t) :: objective_sensitivity
 
    if (present(simulation)) then
       if (simulation%unsteady) then
          ! Objective sensitivity accumulated
          call this%run_backward_unsteady(simulation, design)
       else
          call simulation%run_backward()
          ! Compute sensitivity on steady field
          call this%update_objective_sensitivities(design)
       end if
    else
       call this%update_objective_sensitivities(design)
    end if
 
    call this%update_constraint_sensitivities(design)
 
    call objective_sensitivity%init(this%n_design)
    call this%get_objective_sensitivities(objective_sensitivity)
 
    call design%map_backward(objective_sensitivity)
 
    call objective_sensitivity%free()
  end subroutine problem_compute_sensitivity
 
  subroutine problem_run_forward_unsteady(this, simulation, design)
    class(problem_t), intent(inout) :: this
    class(simulation_t), intent(inout) :: simulation
    class(design_t), intent(inout) :: design
    type(time_step_controller_t) :: dt_controller
    real(kind=dp) :: loop_start
 
    call dt_controller%init(simulation%neko_case%params)
 
    call simulation%reset()
    call simulation_init(simulation%neko_case, dt_controller)
 
    ! Reset the objective value to zero
    call this%reset_objectives()
 
    call profiler_start_region("Forward simulation")
    loop_start = mpi_wtime()
    simulation%n_timesteps = 0
    do while (simulation%neko_case%time%t .lt. simulation%neko_case%time%end_time)
       simulation%n_timesteps = simulation%n_timesteps + 1
       ! step forward
       call simulation_step(simulation%neko_case, dt_controller, loop_start)
       ! accumulate objective value
       call this%accumulate_objectives(design, simulation%adjoint_case%time%dt)
       ! save a checkpoint
       call simulation%checkpoint%save(simulation%neko_case)
    end do
    call profiler_end_region("Forward simulation")
 
    call simulation_finalize(simulation%neko_case)
 
  end subroutine problem_run_forward_unsteady
 
  subroutine problem_run_backward_unsteady(this, simulation, design)
    class(problem_t), intent(inout) :: this
    class(simulation_t), intent(inout) :: simulation
    class(design_t), intent(inout) :: design
    type(time_step_controller_t) :: dt_controller
    real(kind=dp) :: loop_start
    real(kind=rp) :: cfl
    real(kind=rp) :: total_time
    integer :: i
 
    call dt_controller%init(simulation%neko_case%params)
 
    call simulation_adjoint_init(simulation%adjoint_case, dt_controller)
 
    ! Reset the sensitivity value to zero
    call this%reset_objective_sensitivities()
 
    cfl = simulation%adjoint_case%fluid_adj%compute_cfl(simulation%adjoint_case%time%dt)
    loop_start = mpi_wtime()
 
    total_time = simulation%n_timesteps * simulation%adjoint_case%time%dt
 
    call profiler_start_region("Adjoint simulation")
 
    ! TODO. IC's need to be handled rather carefully, we should take a
    ! checkpoint at i = 0. However, 99% of the time we use an initial condition
    ! for the fluid of u=0, so this doesn't matter. Then we use u_adj = 0 on
    ! the other end.
    !
    ! we have:
    !  - n    time steps to compute
    !  - n+1  fields to consider
    !  - n-1  non-zero contributions to u * u_adj
    !
    !              non-zero
    !             |--------|
    !  primal  o--x--x--x--x--x
    !          x--x--x--x--x--o  adjoint
    !          ^              ^
    !         u=0          u_adj=0
 
    do i = simulation%n_timesteps, 1, -1
       ! restore primal field
       call simulation%checkpoint%restore(simulation%neko_case, i)
       ! accumulate objective sensitivity
       call this%accumulate_objective_sensitivities(design, &
            simulation%adjoint_case%time%dt)
       ! step the adjoint backwards
       call simulation_adjoint_step(simulation%adjoint_case, dt_controller, &
            cfl, loop_start, total_time)
    end do
 
    call profiler_end_region("Forward simulation")
 
    call simulation_adjoint_finalize(simulation%adjoint_case)
 
  end subroutine problem_run_backward_unsteady
 
  ! ========================================================================== !
  ! Update the objectives and constraints
 
  subroutine problem_update_objectives(this, design)
    class(problem_t), intent(inout) :: this
    class(design_t), intent(in) :: design
    integer :: i
 
    do i = 1, this%n_objectives
       call this%objective_list(i)%objective%update_value(design)
    end do
  end subroutine problem_update_objectives
 
  subroutine problem_update_constraints(this, design)
    class(problem_t), intent(inout) :: this
    class(design_t), intent(in) :: design
    integer :: i
 
    do i = 1, this%n_constraints
       call this%constraint_list(i)%constraint%update_value(design)
    end do
  end subroutine problem_update_constraints
 
  subroutine problem_update_objective_sensitivities(this, design)
    class(problem_t), intent(inout) :: this
    class(design_t), intent(in) :: design
    integer :: i
 
    do i = 1, this%n_objectives
       call this%objective_list(i)%objective%update_sensitivity(design)
    end do
  end subroutine problem_update_objective_sensitivities
 
  subroutine problem_update_constraint_sensitivities(this, design)
    class(problem_t), intent(inout) :: this
    class(design_t), intent(in) :: design
    integer :: i
 
    do i = 1, this%n_constraints
       call this%constraint_list(i)%constraint%update_sensitivity(design)
    end do
  end subroutine problem_update_constraint_sensitivities
 
  ! ========================================================================== !
  ! Reset the objectives and constraints
 
  subroutine problem_reset_objectives(this)
    class(problem_t), intent(inout) :: this
    integer :: i
 
    do i = 1, this%n_objectives
       call this%objective_list(i)%objective%reset_value()
    end do
  end subroutine problem_reset_objectives
 
  subroutine problem_reset_constraints(this)
    class(problem_t), intent(inout) :: this
    integer :: i
 
    do i = 1, this%n_constraints
       call this%constraint_list(i)%constraint%reset_value()
    end do
  end subroutine problem_reset_constraints
 
  subroutine problem_reset_objective_sensitivities(this)
    class(problem_t), intent(inout) :: this
    integer :: i
 
    do i = 1, this%n_objectives
       call this%objective_list(i)%objective%reset_sensitivity()
    end do
  end subroutine problem_reset_objective_sensitivities
 
  subroutine problem_reset_constraint_sensitivities(this)
    class(problem_t), intent(inout) :: this
    integer :: i
 
    do i = 1, this%n_constraints
       call this%constraint_list(i)%constraint%reset_sensitivity()
    end do
  end subroutine problem_reset_constraint_sensitivities
 
  ! ========================================================================== !
  ! Accumulate the objectives and constraints
 
  subroutine problem_accumulate_objectives(this, design, dt)
    class(problem_t), intent(inout) :: this
    class(design_t), intent(in) :: design
    real(kind=rp), intent(in) :: dt
    integer :: i
 
    do i = 1, this%n_objectives
       call this%objective_list(i)%objective%accumulate_value(design, dt)
    end do
  end subroutine problem_accumulate_objectives
 
  subroutine problem_accumulate_constraints(this, design, dt)
    class(problem_t), intent(inout) :: this
    class(design_t), intent(in) :: design
    real(kind=rp), intent(in) :: dt
    integer :: i
 
    do i = 1, this%n_constraints
       call this%constraint_list(i)%constraint%accumulate_value(design, dt)
    end do
  end subroutine problem_accumulate_constraints
 
  subroutine problem_accumulate_objective_sensitivities(this, design, dt)
    class(problem_t), intent(inout) :: this
    class(design_t), intent(in) :: design
    real(kind=rp), intent(in) :: dt
    integer :: i
 
    do i = 1, this%n_objectives
       call this%objective_list(i)%objective%accumulate_sensitivity(design, dt)
    end do
  end subroutine problem_accumulate_objective_sensitivities
 
  subroutine problem_accumulate_constraint_sensitivities(this, design, dt)
    class(problem_t), intent(inout) :: this
    class(design_t), intent(in) :: design
    real(kind=rp), intent(in) :: dt
    integer :: i
 
    do i = 1, this%n_constraints
       call this%constraint_list(i)%constraint%accumulate_sensitivity(design, dt)
    end do
  end subroutine problem_accumulate_constraint_sensitivities
 
  ! ========================================================================== !
  ! Problem part getters
 
  subroutine problem_get_objective_value(this, objective_value)
    class(problem_t), intent(in) :: this
    real(kind=rp), intent(out) :: objective_value
    integer :: i
 
    objective_value = 0.0_rp
    do i = 1, this%n_objectives
       objective_value = objective_value + &
            this%objective_list(i)%objective%get_weight() * &
            this%objective_list(i)%objective%get_value()
    end do
 
  end subroutine problem_get_objective_value
 
  subroutine problem_get_all_objective_values(this, all_objective_values)
    class(problem_t), intent(in) :: this
    type(vector_t), intent(inout) :: all_objective_values
    integer :: i
 
    do i = 1, this%n_objectives
       all_objective_values%x(i) = this%objective_list(i)%objective%value
    end do
 
    call all_objective_values%copy_from(host_to_device, sync = .true.)
 
  end subroutine problem_get_all_objective_values
 
  subroutine problem_get_constraint_values(this, constraint_value)
    class(problem_t), intent(in) :: this
    type(vector_t), intent(inout) :: constraint_value
    integer :: i
 
    do i = 1, this%n_constraints
       constraint_value%x(i) = this%constraint_list(i)%constraint%value
    end do
 
    call constraint_value%copy_from(host_to_device, sync = .true.)
 
  end subroutine problem_get_constraint_values
 
  subroutine problem_get_objective_sensitivities(this, sensitivity)
    class(problem_t), intent(in) :: this
    type(vector_t), intent(inout) :: sensitivity
    integer :: i
 
    call vector_cfill(sensitivity, 0.0_rp)
    do i = 1, this%n_objectives
       call vector_add2(sensitivity, &
            this%objective_list(i)%objective%sensitivity)
    end do
 
  end subroutine problem_get_objective_sensitivities
 
  subroutine problem_get_constraint_sensitivities(this, sensitivity)
    class(problem_t), intent(inout) :: this
    type(matrix_t), target, intent(inout) :: sensitivity
    real(kind=rp), pointer :: row(:)
    integer :: i
 
    ! Copy all constraint sensitivities to host, sync on last one
    do i = 1, this%n_constraints
       call this%constraint_list(i)%constraint%sensitivity%copy_from( &
            device_to_host, sync = i .eq. this%n_constraints)
    end do
 
    do i = 1, this%n_constraints
       row(1:this%n_design) => sensitivity%x(i, :)
 
       call copy(row, this%constraint_list(i)%constraint%sensitivity%x, &
            this%n_design)
    end do
 
    call sensitivity%copy_from(host_to_device, sync = .true.)
 
  end subroutine problem_get_constraint_sensitivities
 
  ! ========================================================================== !
  ! Simple getters
 
  pure function problem_get_num_objectives(this) result(n)
    class(problem_t), intent(in) :: this
    integer :: n
 
    n = this%n_objectives
  end function problem_get_num_objectives
 
  pure function problem_get_num_constraints(this) result(n)
    class(problem_t), intent(in) :: this
    integer :: n
 
    n = this%n_constraints
  end function problem_get_num_constraints
 
  function problem_get_log_header(this) result(buff)
    class(problem_t), intent(in) :: this
    character(len=1024) :: buff
    character(len=50) :: mini_buff
    integer :: i
 
    ! When it comes to multi-objective optimization
    ! (handled in the way that we do) we want to know the value of each
    ! objective individually, not just the combined effect.
    !
    ! my vision is:
    !
    !      | Total F | F_1 | F_2 | ... | F_n | C_1 | C_2 | ... | C_m |
    !
    ! And then if we also want things like this iteration or KKT they can be
    ! appended to the beginning or end of this by the optimizer.
    !
    ! iter | Total F | F_1 | F_2 | ... | F_n | C_1 | C_2 | ... | C_m | KKT
    buff = "Total objective function"
    do i = 1, this%get_n_objectives()
       mini_buff = ""
       write(mini_buff, '(", ", A)') this%objective_list(i)%objective%name
       buff = trim(buff) // trim(mini_buff)
    end do
 
    do i = 1, this%get_n_constraints()
       mini_buff = ""
       write(mini_buff, '(", ", A)') &
            this%constraint_list(i)%constraint%name
       buff = trim(buff) // trim(mini_buff)
    end do
 
  end function problem_get_log_header
end module problem