volume_constraint.f90

Go to the documentation of this file.

 
!
module volume_constraint
  use constraint, only: constraint_t
 
  use design, only: design_t
  use brinkman_design, only: brinkman_design_t
  use simulation_m, only: simulation_t
  use mapping_handler, only: mapping_handler_t
  use num_types, only: rp
  use coefs, only: coef_t
  use json_module, only: json_file
  use json_utils, only: json_get, json_get_or_default
  use field, only: field_t
  use scratch_registry, only: neko_scratch_registry
  use neko_config, only: neko_bcknd_device
  use mask_ops, only: mask_exterior_const
  use math, only: glsc2, copy, cmult
  use device_math, only: device_glsc2, device_copy, device_cmult
  use vector_math, only: vector_cmult
  use math_ext, only: glsc2_mask
  use field_math, only: field_rone, field_copy, field_cmult, field_cfill
  use utils, only: neko_error
  use vector, only: vector_t
  use neko_ext, only: field_to_vector
  implicit none
  private
 
  type, public, extends(constraint_t) :: volume_constraint_t
     private
 
     logical :: is_max
     real(kind=rp) :: limit
     real(kind=rp) :: volume_domain
     class(coef_t), pointer :: c_xh => null()
     type(mapping_handler_t) :: mapping
     logical :: if_mapping = .false.
 
   contains
 
     procedure, public, pass(this) :: init_json_sim => &
          volume_constraint_init_json_sim
     procedure, public, pass(this) :: init_from_components => &
          volume_constraint_init_from_components
     procedure, public, pass(this) :: free => volume_constraint_free
     procedure, public, pass(this) :: update_value => &
          volume_constraint_update_value
     procedure, public, pass(this) :: update_sensitivity => &
          volume_constraint_update_sensitivity
     procedure, public, pass(this) :: get_log_size => &
          volume_constraint_get_log_size
     procedure, public, pass(this) :: get_log_headers => &
          volume_constraint_get_log_headers
     procedure, public, pass(this) :: get_log_values => &
          volume_constraint_get_log_values
 
     procedure, private, pass(this) :: compute_volume
 
  end type volume_constraint_t
 
contains
 
  subroutine volume_constraint_init_json_sim(this, json, design, simulation)
    class(volume_constraint_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 :: mask_name
    character(len=:), allocatable :: name
    logical :: is_max
    real(kind=rp) :: limit
 
    call json_get_or_default(json, "name", name, "Volume constraint")
    call json_get_or_default(json, "mask_name", mask_name, "")
    call json_get_or_default(json, "is_max", is_max, .false.)
    call json_get(json, "limit", limit)
 
    call this%init_from_components(design, simulation, name, mask_name, &
         is_max, limit)
 
    ! check if we have a mapping
    if ('mapping' .in. json) then
       ! Initialize the mapper
       call this%mapping%init_base(this%c_Xh)
       call this%mapping%add(json, 'mapping')
       this%if_mapping = .true.
       ! recompute value after mapping
       call this%update_value(design)
    end if
 
  end subroutine volume_constraint_init_json_sim
 
  subroutine volume_constraint_init_from_components(this, design, simulation, &
       name, mask_name, is_max, limit)
    class(volume_constraint_t), intent(inout) :: this
    class(design_t), intent(in) :: design
    type(simulation_t), target, intent(inout) :: simulation
    character(len=*), intent(in) :: mask_name
    character(len=*), intent(in) :: name
    logical, intent(in) :: is_max
    real(kind=rp), intent(in) :: limit
 
    type(field_t), pointer :: work
    integer :: temp_indices(1)
 
    ! Initialize the base class
    call this%init_base(name, design%size(), mask_name)
 
    ! Store the attributes
    this%is_max = is_max
    this%limit = limit
    this%c_Xh => simulation%neko_case%fluid%c_Xh
 
    ! Now we can extract the mask/has_mask from the design
    if (this%has_mask) then
 
       ! calculate the volume of the optimization domain
       call neko_scratch_registry%request(work, temp_indices(1), .false.)
       call field_rone(work)
       call mask_exterior_const(work, this%mask, 0.0_rp)
 
       if (neko_bcknd_device .eq. 1) then
          this%volume_domain = device_glsc2(work%x_d, this%c_xh%B_d, &
               work%size())
       else
          this%volume_domain = glsc2_mask(work%x, this%c_Xh%B, &
               design%size(), this%mask%mask%get(), this%mask%size)
       end if
 
       call neko_scratch_registry%relinquish(temp_indices)
    else
       this%volume_domain = this%c_Xh%volume
    end if
 
    ! ------------------------------------------------------------------------ !
    ! Initialize the value of constraint
 
    call this%update_value(design)
 
    ! ------------------------------------------------------------------------ !
    ! Initialize the sensitivity value
 
    if (neko_bcknd_device .eq. 1) then
       call device_copy(this%sensitivity%x_d, this%c_xh%B_d, design%size())
       call device_cmult(this%sensitivity%x_d, -1.0_rp / this%volume_domain, &
            design%size())
 
       if (this%is_max) then
          call device_cmult(this%sensitivity%x_d, -1.0_rp, design%size())
       end if
    else
       call copy(this%sensitivity%x, this%c_Xh%B, design%size())
       call cmult(this%sensitivity%x, -1.0_rp / this%volume_domain, &
            design%size())
 
       if (this%is_max) then
          call vector_cmult(this%sensitivity, -1.0_rp)
       end if
    end if
 
    if (this%has_mask) then
       call mask_exterior_const(this%sensitivity, this%mask, 0.0_rp)
    end if
 
  end subroutine volume_constraint_init_from_components
 
  subroutine volume_constraint_free(this)
    class(volume_constraint_t), intent(inout) :: this
 
    call this%free_base()
  end subroutine volume_constraint_free
 
  subroutine volume_constraint_update_value(this, design)
    class(volume_constraint_t), intent(inout) :: this
    class(design_t), intent(in) :: design
    real(kind=rp) :: volume
 
    volume = this%compute_volume(design)
 
    ! Compute the distance to the target volume
    this%value = this%limit - volume / this%volume_domain
 
    ! Invert the sign if it is a maximum constraint
    if (this%is_max) this%value = -this%value
 
  end subroutine volume_constraint_update_value
 
  subroutine volume_constraint_update_sensitivity(this, design)
    class(volume_constraint_t), intent(inout) :: this
    class(design_t), intent(in) :: design
    type(field_t), pointer :: unmapped, mapped
    integer :: temp_indices(2)
 
    if (this%if_mapping) then
       ! Recompute and map backward
       call neko_scratch_registry%request(unmapped, temp_indices(1), &
            .false.)
       call neko_scratch_registry%request(mapped, temp_indices(2), &
            .false.)
       ! The mapping will handle the mass matrix
       call field_cfill(unmapped, -1.0_rp / this%volume_domain)
       if (this%is_max) then
          call field_cmult(unmapped, -1.0_rp)
       end if
       ! mask
       if (this%has_mask) then
          call mask_exterior_const(unmapped, this%mask, 0.0_rp)
       end if
       ! map backwards
       call this%mapping%apply_backward(mapped, unmapped)
       call field_to_vector(this%sensitivity, mapped)
 
       call neko_scratch_registry%relinquish(temp_indices)
 
    else
       ! Sensitivity is just a constant so it should not be updated
    end if
 
  end subroutine volume_constraint_update_sensitivity
 
  ! ========================================================================== !
  ! The actual volume computations for different types of designs
 
  function compute_volume(this, design) result(volume)
    class(volume_constraint_t), intent(inout) :: this
    class(design_t), intent(in) :: design
    real(kind=rp) :: volume
 
    volume = 0.0_rp
    select type (design)
    type is (brinkman_design_t)
       volume = volume_brinkman_design(this, design)
 
    class default
       call neko_error('Volume constraint only works with brinkman_design')
    end select
 
  end function compute_volume
 
  function volume_brinkman_design(this, design) result(volume)
    class(volume_constraint_t), intent(inout) :: this
    type(brinkman_design_t), intent(in) :: design
    real(kind=rp) :: volume
    type(field_t), pointer :: work
    type(vector_t), pointer :: values, unmapped_values
    integer :: temp_indices, ind_value, ind_um_value
 
    call neko_scratch_registry%request(values, ind_value, design%size(), &
         .false.)
 
    volume = 0.0_rp
    if (this%if_mapping) then
       call neko_scratch_registry%request(unmapped_values, ind_um_value, &
            design%size(), .false.)
       call design%get_values(unmapped_values)
       call this%mapping%apply_forward(values, unmapped_values)
       call neko_scratch_registry%relinquish(ind_um_value)
    else
       call design%get_values(values)
    end if
 
    if (this%has_mask) then
 
       if (neko_bcknd_device .eq. 1) then
          call neko_scratch_registry%request(work, temp_indices, .false.)
          call device_copy(work%x_d, values%x_d, design%size())
          call mask_exterior_const(work, this%mask, 0.0_rp)
 
          volume = device_glsc2(work%x_d, this%c_xh%B_d, design%size())
 
          call neko_scratch_registry%relinquish(temp_indices)
       else
          volume = glsc2_mask(values%x, this%c_Xh%B, design%size(), &
               this%mask%mask%get(), this%mask%size)
       end if
 
    else
 
       if (neko_bcknd_device .eq. 1) then
          volume = device_glsc2(values%x_d, this%c_xh%B_d, design%size())
       else
          volume = glsc2(values%x, this%c_Xh%B, design%size())
       end if
 
    end if
 
    call neko_scratch_registry%relinquish(ind_value)
 
  end function volume_brinkman_design
 
  function volume_constraint_get_log_size(this) result(n)
    class(volume_constraint_t), intent(in) :: this
    integer :: n
 
    n = 2
  end function volume_constraint_get_log_size
 
  subroutine volume_constraint_get_log_headers(this, headers)
    class(volume_constraint_t), intent(in) :: this
    character(len=*), intent(out) :: headers(:)
    character(len=64) :: prefix
 
    if (size(headers) .lt. 1) return
    prefix = trim(this%name)
    headers(1) = prefix
    if (size(headers) .lt. 2) return
    headers(2) = trim(prefix) // '.volume'
 
  end subroutine volume_constraint_get_log_headers
 
  subroutine volume_constraint_get_log_values(this, values)
    class(volume_constraint_t), intent(in) :: this
    real(kind=rp), intent(out) :: values(:)
    real(kind=rp) :: volume_ratio
 
    if (size(values) .lt. 1) return
    if (this%is_max) then
       volume_ratio = this%limit + this%value
    else
       volume_ratio = this%limit - this%value
    end if
 
    values(1) = this%value
    if (size(values) .lt. 2) return
    values(2) = volume_ratio
 
  end subroutine volume_constraint_get_log_values
 
end module volume_constraint