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.
!
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 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, 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, 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_field(work, temp_indices(1))
       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_field(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_field(unmapped, temp_indices(1))
       call neko_scratch_registry%request_field(mapped, temp_indices(2))
       ! 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_field(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) :: values, unmapped_values
    integer :: temp_indices(1)
 
    volume = 0.0_rp
    if (this%if_mapping) then
       call design%get_values(unmapped_values)
       call values%init(unmapped_values%size())
       call this%mapping%apply_forward(values, unmapped_values)
       call unmapped_values%free()
    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_field(work, temp_indices(1))
          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_field(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 values%free()
 
  end function volume_brinkman_design
 
end module volume_constraint