base_functional.f90

Go to the documentation of this file.

 
 
module base_functional
  use design, only: design_t
  use json_module, only: json_file
  use json_utils, only: json_get
  use num_types, only: rp, dp, xp
  use point_zone, only: point_zone_t
  use simulation_m, only: simulation_t
  use time_state, only: time_state_t
  use vector, only: vector_t
  use utils, only: neko_error
  use scratch_registry, only: neko_scratch_registry
  use vector_math, only: vector_copy, vector_add3s2, vector_rzero
  implicit none
  private
 
  type, abstract, public :: base_functional_t
 
     character(len=:), allocatable :: name
 
     real(kind=rp) :: value = 0.0_dp
     type(vector_t) :: sensitivity
     logical :: has_mask = .false.
     class(point_zone_t), pointer :: mask => null()
 
     ! ----------------------------------------------------------------------- !
     ! Variables related to time averaging of the functional value and
     ! sensitivity
 
     real(kind=dp) :: start_time = 0.0_dp
     real(kind=dp) :: end_time = huge(0.0_dp)
     real(kind=xp) :: value_weight = 0.0_xp
     real(kind=xp) :: sensitivity_weight = 0.0_xp
 
   contains
 
     ! ----------------------------------------------------------------------- !
     ! Derived class interfaces
 
     generic :: init => init_json, init_json_sim
 
     procedure, pass(this) :: init_json => functional_init_json
     procedure, pass(this) :: init_json_sim => functional_init_json_sim
     procedure(functional_free), pass(this), deferred :: free
 
     procedure(functional_update_value), pass(this), deferred :: update_value
     procedure(functional_update_sensitivity), pass(this), deferred :: &
          update_sensitivity
 
     procedure, pass(this) :: get_name => functional_get_name
     procedure, pass(this) :: get_value => functional_get_value
     procedure, pass(this) :: get_sensitivity => functional_get_sensitivity
     procedure, pass(this) :: get_log_size => functional_get_log_size
     procedure, pass(this) :: get_log_headers => functional_get_log_headers
     procedure, pass(this) :: get_log_values => functional_get_log_values
     procedure, pass(this) :: reset_value => functional_reset_value
     procedure, pass(this) :: reset_sensitivity => functional_reset_sensitivity
     procedure, pass(this) :: accumulate_value => functional_accumulate_value
     procedure, pass(this) :: accumulate_sensitivity => &
          functional_accumulate_sensitivity
 
     procedure, private, pass(this) :: in_window => functional_in_window
  end type base_functional_t
 
  ! -------------------------------------------------------------------------- !
  ! Interface specifications for the derived types, these are the constructors
  ! for the different types of objective functions.
 
  abstract interface
 
     subroutine functional_free(this)
       import base_functional_t
       class(base_functional_t), intent(inout) :: this
     end subroutine functional_free
 
     subroutine functional_update_value(this, design)
       import base_functional_t, design_t
       class(base_functional_t), intent(inout) :: this
       class(design_t), intent(in) :: design
     end subroutine functional_update_value
 
     subroutine functional_update_sensitivity(this, design)
       import base_functional_t, design_t
       class(base_functional_t), intent(inout) :: this
       class(design_t), intent(in) :: design
     end subroutine functional_update_sensitivity
 
  end interface
 
contains
 
  subroutine functional_init_json(this, json, design)
    class(base_functional_t), intent(inout) :: this
    type(json_file), intent(inout) :: json
    class(design_t), intent(in) :: design
    character(len=:), allocatable :: type
 
    call json_get(json, 'type', type)
    call neko_error("Functional type: '" // type // &
         "' does not support initialization without simulation")
  end subroutine functional_init_json
 
  subroutine functional_init_json_sim(this, json, design, simulation)
    class(base_functional_t), intent(inout) :: this
    type(json_file), intent(inout) :: json
    class(design_t), intent(in) :: design
    type(simulation_t), target, intent(inout) :: simulation
    character(len=:), allocatable :: type
 
    call json_get(json, 'type', type)
    call neko_error("Functional type: '" // type // &
         "' does not support initialization with simulation")
  end subroutine functional_init_json_sim
 
  function functional_get_name(this) result(name)
    class(base_functional_t), intent(in) :: this
    character(len=:), allocatable :: name
 
    if (allocated(this%name)) then
       name = this%name
    else
       name = ""
    end if
  end function functional_get_name
 
  function functional_get_value(this) result(v)
    class(base_functional_t), intent(in) :: this
    real(kind=rp) :: v
 
    v = this%value
  end function functional_get_value
 
  subroutine functional_get_sensitivity(this, sensitivity)
    class(base_functional_t), intent(in) :: this
    type(vector_t), intent(inout) :: sensitivity
 
    sensitivity = this%sensitivity
  end subroutine functional_get_sensitivity
 
  function functional_get_log_size(this) result(n)
    class(base_functional_t), intent(in) :: this
    integer :: n
 
    n = 1
  end function functional_get_log_size
 
  subroutine functional_get_log_headers(this, headers)
    class(base_functional_t), intent(in) :: this
    character(len=*), intent(out) :: headers(:)
 
    if (size(headers) .eq. 0) return
    headers(1) = this%get_name()
  end subroutine functional_get_log_headers
 
  subroutine functional_get_log_values(this, values)
    class(base_functional_t), intent(in) :: this
    real(kind=rp), intent(out) :: values(:)
 
    if (size(values) .eq. 0) return
    values(1) = this%value
  end subroutine functional_get_log_values
 
  subroutine functional_reset_value(this)
    class(base_functional_t), intent(inout) :: this
 
    this%value = 0.0_dp
    this%value_weight = 0.0_dp
  end subroutine functional_reset_value
 
  subroutine functional_reset_sensitivity(this)
    class(base_functional_t), intent(inout) :: this
 
    call vector_rzero(this%sensitivity)
    this%sensitivity_weight = 0.0_dp
  end subroutine functional_reset_sensitivity
 
  subroutine functional_accumulate_value(this, design, time)
    class(base_functional_t), intent(inout) :: this
    class(design_t), intent(in) :: design
    type(time_state_t), intent(in) :: time
    real(kind=xp) :: value_new, value_old, weight_new, weight_old, dt
 
    if (.not. this%in_window(time)) return
 
    ! Store the old value and compute the new value.
    value_old = real(this%value, kind=xp)
    call this%update_value(design)
    value_new = real(this%value, kind=xp)
 
    ! Compute the weights for the old and new values.
    dt = real(time%dt, kind=xp)
    weight_new = dt / (this%value_weight + dt)
    weight_old = this%value_weight / (this%value_weight + dt)
 
    ! Rectangle rule; could potentially use higher order trapezoidal/Simpson
    ! etc, but this should suffice.
    this%value = real(value_old * weight_old + value_new * weight_new, kind=rp)
    this%value_weight = this%value_weight + dt
  end subroutine functional_accumulate_value
 
  subroutine functional_accumulate_sensitivity(this, design, time)
    class(base_functional_t), intent(inout) :: this
    class(design_t), intent(in) :: design
    type(time_state_t), intent(in) :: time
    real(kind=xp) :: weight_new, weight_old, dt
    type(vector_t), pointer :: sensitivity_new, sensitivity_old
    integer :: n, idx(2)
 
    if (.not. this%in_window(time)) return
 
    n = this%sensitivity%size()
    call neko_scratch_registry%request(sensitivity_new, idx(1), n, .false.)
    call neko_scratch_registry%request(sensitivity_old, idx(2), n, .false.)
 
    ! Store old sensitivity and compute new sensitivity.
    call vector_copy(sensitivity_old, this%sensitivity)
    call this%update_sensitivity(design)
    call vector_copy(sensitivity_new, this%sensitivity)
 
    ! Compute the weights for the old and new sensitivities.
    dt = real(time%dt, kind=xp)
    weight_new = dt / (this%sensitivity_weight + dt)
    weight_old = this%sensitivity_weight / (this%sensitivity_weight + dt)
 
    ! Rectangle rule; could potentially use higher order trapezoidal/Simpson
    ! etc, but this should suffice.
    call vector_add3s2(this%sensitivity, &
         sensitivity_new, sensitivity_old, weight_new, weight_old)
    this%sensitivity_weight = this%sensitivity_weight + dt
 
    nullify(sensitivity_new)
    nullify(sensitivity_old)
    call neko_scratch_registry%relinquish(idx)
  end subroutine functional_accumulate_sensitivity
 
  ! -------------------------------------------------------------------------- !
  ! Private helper methods
 
  pure function functional_in_window(this, time) result(inside)
    class(base_functional_t), intent(in) :: this
    type(time_state_t), intent(in) :: time
    logical :: inside
    real(kind=rp) :: tol
 
    tol = 1.0e-6_rp * time%dt
    inside = time%t .ge. this%start_time - tol .and. &
         time%t .le. this%end_time + tol
  end function functional_in_window
end module base_functional