volume_constraint.f90

! Copyright (c) 2023, The Neko Authors
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
!
!   * Redistributions of source code must retain the above copyright
!     notice, this list of conditions and the following disclaimer.
!
!   * Redistributions in binary form must reproduce the above
!     copyright notice, this list of conditions and the following
!     disclaimer in the documentation and/or other materials provided
!     with the distribution.
!
!   * Neither the name of the authors nor the names of its
!     contributors may be used to endorse or promote products derived
!     from this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
! $V = \frac{1}{|\Omega_O|}\int_{\Omega_O} \tilde{\rho} d\Omega$
! Either
! $V < V_\text{max}$
! $V > V_\text{min}$
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 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 field_registry, only: neko_field_registry
  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
  use device_math, only: device_glsc2
  use vector_math, only: vector_cmult
  use math_ext, only: glsc2_mask
  use field_math, only: field_rone, field_copy
  use utils, only: neko_error
  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()
 
   contains
 
     procedure, public, pass(this) :: init_json_sim => &
          volume_constraint_init_json_sim
     procedure, public, pass(this) :: init_from_attributes => &
          volume_constraint_init_attributes
     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, 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, "mask_name", mask_name, "")
    call json_get_or_default(json, "name", name, "Volume constraint")
    call json_get_or_default(json, "is_max", is_max, .false.)
    call json_get(json, "limit", limit)
 
    call this%init_from_attributes(design, simulation, name, mask_name, &
         is_max, limit)
  end subroutine volume_constraint_init_json_sim
 
  subroutine volume_constraint_init_attributes(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
 
    real(kind=rp) :: volume
    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_field(work, temp_indices(1))
       call field_rone(work)
 
       if (neko_bcknd_device .eq. 1) then
          call mask_exterior_const(work, this%mask, 0.0_rp)
          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_field(temp_indices)
    else
       this%volume_domain = this%c_Xh%volume
    end if
 
    ! ------------------------------------------------------------------------ !
    ! Initialize the value of constraint
 
    ! Compute the volume of the design
    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 )
 
    ! ------------------------------------------------------------------------ !
    ! Initialize the sensitivity value
 
    this%sensitivity = 1.0_rp / this%volume_domain
 
    ! Invert the sign if it is a maximum constraint
    call vector_cmult(this%sensitivity, -1.0_rp)
 
    if (this%has_mask) then
       call mask_exterior_const(this%sensitivity, this%mask, 0.0_rp)
    end if
 
  end subroutine volume_constraint_init_attributes
 
  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
 
    ! Sensitivity is just a constant so it should not be updated
 
  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, design_indicator
    integer :: temp_indices(1)
 
    volume = 0.0_rp
    design_indicator => neko_field_registry%get_field("design_indicator")
 
    ! in the future we should be using the mapped design variable
    ! corresponding to this constraint!!!
    if (this%has_mask) then
 
       if (neko_bcknd_device .eq. 1) then
          call neko_scratch_registry%request_field(work, temp_indices(1))
          call field_copy(work, design_indicator)
          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_field(temp_indices)
       else
          volume = glsc2_mask(design_indicator%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(design_indicator%x_d, &
               this%c_xh%B_d, design%size())
       else
          volume = glsc2(design_indicator%x, &
               this%c_Xh%B, design%size())
       end if
 
    end if
 
  end function volume_brinkman_design
 
end module volume_constraint