problem.f90

Go to the documentation of this file.

 
 
module problem
  use num_types, only: rp
  use fld_file_output, only: fld_file_output_t
  use design, only: design_t
  use objective, only: objective_t, objective_wrapper_t, objective_factory
  use augmented_lagrangian_objective, only: augmented_lagrangian_objective_t
  use constraint, only: constraint_t, constraint_wrapper_t, constraint_factory
  use vector, only: vector_t
  use matrix, only: matrix_t
  use device, only: device_memcpy, host_to_device, device_to_host
  use neko_config, only: neko_bcknd_device
  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 device_math, only: device_copy
  use vector_math, only: vector_add2
 
  implicit none
  private
 
  type, public :: problem_t
     private
 
     integer :: n_design
     integer :: n_objectives
     integer :: n_constraints
 
     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
 
     ! ----------------------------------------------------------------------- !
     ! 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
 
     ! ----------------------------------------------------------------------- !
     ! 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
 
    this%n_design = design%size()
    this%n_objectives = 0
    this%n_constraints = 0
 
    ! minimum dissipation objective function
    call this%read_objectives(parameters, design, simulation)
 
    ! volume constraint
    call this%read_constraints(parameters, design, simulation)
 
  end subroutine problem_init
 
  subroutine problem_free(this)
    class(problem_t), intent(inout) :: this
    integer :: i
 
    ! 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
    class(objective_t), allocatable :: objective
    type(simulation_t), optional, intent(inout) :: simulation
 
    ! 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()
       call simulation%run_forward()
    end if
 
    call this%update_objectives(design)
    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)) call simulation%run_backward()
 
    call this%update_objective_sensitivities(design)
    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
 
  ! ========================================================================== !
  ! 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
 
  ! ========================================================================== !
  ! 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
 
    call all_objective_values%init(this%n_objectives)
    do i = 1, this%n_objectives
       all_objective_values%x(i) = this%objective_list(i)%objective%value
    end do
 
    if (neko_bcknd_device .eq. 1) then
       call device_memcpy(all_objective_values%x, all_objective_values%x_d, &
            this%n_objectives, host_to_device, sync = .true.)
    end if
 
  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
 
    call constraint_value%init(this%n_constraints)
    do i = 1, this%n_constraints
       constraint_value%x(i) = this%constraint_list(i)%constraint%value
    end do
 
    if (neko_bcknd_device .eq. 1) then
       call device_memcpy(constraint_value%x, constraint_value%x_d, &
            this%n_constraints, host_to_device, sync = .true.)
    end if
 
  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 sensitivity%init(this%n_design)
    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(in) :: this
    type(matrix_t), intent(inout) :: sensitivity
    type(vector_t) :: tmp
    integer :: i, j, n
 
    n = this%n_constraints * this%n_design
    call sensitivity%init(this%n_constraints, this%n_design)
    if (neko_bcknd_device .eq. 1) then
       call tmp%init(this%n_design)
    end if
 
    do i = 1, this%n_constraints
       if (neko_bcknd_device .eq. 1) then
          tmp = this%constraint_list(i)%constraint%sensitivity
          call device_memcpy(tmp%x, tmp%x_d, &
               this%n_design, device_to_host, sync = .true.)
          do j = 1, this%n_design
             sensitivity%x(i, j) = tmp%x(j)
          end do
       else
          do j = 1, this%n_design
             sensitivity%x(i, j) = &
                  this%constraint_list(i)%constraint%sensitivity%x(j)
          end do
       end if
    end do
 
    if (neko_bcknd_device .eq. 1) then
       call device_memcpy(sensitivity%x, sensitivity%x_d, n, &
            host_to_device, sync = .true.)
    end if
 
    ! Free the temporary vector
    if (neko_bcknd_device .eq. 1) then
       call tmp%free()
    end if
 
  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 thie iteration or KKT they can be
    ! appended to the begining 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