minimum_dissipation_objective.f90

! Copyright (c) 2024-25, 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.
!
!
! I promise I'll write this document properly in the future...
!
! But the Borval Peterson (I think) paper had an objective function
! that had 2 terms, dissipation and this term they claimed represented
! out of plane stresses.
! I never really understood that extra term, I also don't think it
! applies to 3D cases, but everyone includes it anyway.
!
! It appears to me to be basically a heuristic penality that targets
! non-binary designs
!
! so let's call
!
! F = \int |\nabla u|^2  + K \int \chi \u^2
!
!      | dissipation |     |"lube term"|
!
! I say "lube term" because they said it came from lubrication theory...
! Anyway, we can change all this later (especially the names!)
 
! If the objective function \int |\nabla u|^2,
! the corresponding adjoint forcing is \int \nabla v \cdot \nabla u
!
! for the lube term, the adjoint forcing is \chi u
!
! This has always annoyed me...
! because now I see one objective and one constraint
!
module minimum_dissipation_objective
  use num_types, only: rp
  use field, only: field_t
  use field_math, only: field_col3, field_addcol3, field_cmult, field_add2s2, &
       field_copy
  use operators, only: grad
  use adjoint_fluid_pnpn, only: adjoint_fluid_pnpn_t
  use scratch_registry, only: neko_scratch_registry
  use adjoint_minimum_dissipation_source_term, only: &
       adjoint_minimum_dissipation_source_term_t
  use objective, only: objective_t
  use simulation_m, only: simulation_t
  use adjoint_fluid_scheme, only: adjoint_fluid_scheme_t
  use coefs, only: coef_t
  use field_registry, only: neko_field_registry
  use neko_config, only: neko_bcknd_device
  use math, only: glsc2, copy
  use device_math, only: device_copy, device_glsc2
  use design, only: design_t
  use brinkman_design, only: brinkman_design_t
  use point_zone, only: point_zone_t
  use mask_ops, only: mask_exterior_const
  use math_ext, only: glsc2_mask
  use utils, only: neko_error
  use json_module, only: json_file
  use json_utils, only: json_get_or_default
  implicit none
  private
 
  type, public, extends(objective_t) :: minimum_dissipation_objective_t
     private
 
     type(field_t), pointer :: u => null()
     type(field_t), pointer :: v => null()
     type(field_t), pointer :: w => null()
     type(coef_t), pointer :: c_xh => null()
 
     type(field_t), pointer :: adjoint_u => null()
     type(field_t), pointer :: adjoint_v => null()
     type(field_t), pointer :: adjoint_w => null()
 
   contains
     procedure, public, pass(this) :: init_json_sim => &
          minimum_dissipation_init_json_sim
     procedure, public, pass(this) :: init_from_attributes => &
          minimum_dissipation_init_attributes
     procedure, public, pass(this) :: free => minimum_dissipation_free
     procedure, public, pass(this) :: update_value => &
          minimum_dissipation_update_value
     procedure, public, pass(this) :: update_sensitivity => &
          minimum_dissipation_update_sensitivity
 
  end type minimum_dissipation_objective_t
 
contains
 
  subroutine minimum_dissipation_init_json_sim(this, json, design, simulation)
    class(minimum_dissipation_objective_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 :: name
    character(len=:), allocatable :: mask_name
    real(kind=rp) :: weight
 
    call json_get_or_default(json, "weight", weight, 1.0_rp)
    call json_get_or_default(json, "mask_name", mask_name, "")
    call json_get_or_default(json, "name", name, "Dissipation")
 
    call this%init_from_attributes(design, simulation, weight, name, mask_name)
  end subroutine minimum_dissipation_init_json_sim
 
  subroutine minimum_dissipation_init_attributes(this, design, simulation, &
       weight, name, mask_name)
    class(minimum_dissipation_objective_t), intent(inout) :: this
    class(design_t), intent(in) :: design
    type(simulation_t), target, intent(inout) :: simulation
    real(kind=rp), intent(in) :: weight
    character(len=*), intent(in) :: name
    character(len=*), intent(in) :: mask_name
 
    type(adjoint_minimum_dissipation_source_term_t) :: adjoint_forcing
 
    call this%init_base(name, design%size(), weight, mask_name)
 
    ! Save the simulation and design
    this%u => neko_field_registry%get_field('u')
    this%v => neko_field_registry%get_field('v')
    this%w => neko_field_registry%get_field('w')
    this%c_Xh => simulation%fluid%c_Xh
    this%adjoint_u => neko_field_registry%get_field('u_adj')
    this%adjoint_v => neko_field_registry%get_field('v_adj')
    this%adjoint_w => neko_field_registry%get_field('w_adj')
 
    ! you will need to init this!
    ! append a source term based on the minimum dissipation
    ! init the adjoint forcing term for the adjoint
    call adjoint_forcing%init_from_components( &
         simulation%adjoint_fluid%f_adj_x, &
         simulation%adjoint_fluid%f_adj_y, &
         simulation%adjoint_fluid%f_adj_z, &
         this%u, this%v, this%w, this%weight, &
         this%mask, this%has_mask, &
         this%c_Xh)
 
    ! append adjoint forcing term based on objective function
    select type (f => simulation%adjoint_fluid)
    type is (adjoint_fluid_pnpn_t)
       call f%source_term%add_source_term(adjoint_forcing)
    end select
 
  end subroutine minimum_dissipation_init_attributes
 
  subroutine minimum_dissipation_free(this)
    class(minimum_dissipation_objective_t), intent(inout) :: this
    call this%free_base()
 
    if (associated(this%u)) nullify(this%u)
    if (associated(this%v)) nullify(this%v)
    if (associated(this%w)) nullify(this%w)
    if (associated(this%c_Xh)) nullify(this%c_Xh)
 
    if (associated(this%adjoint_u)) nullify(this%adjoint_u)
    if (associated(this%adjoint_v)) nullify(this%adjoint_v)
    if (associated(this%adjoint_w)) nullify(this%adjoint_w)
 
  end subroutine minimum_dissipation_free
 
  subroutine minimum_dissipation_update_value(this, design)
    class(minimum_dissipation_objective_t), intent(inout) :: this
    class(design_t), intent(in) :: design
    type(field_t), pointer :: wo1, wo2, wo3, work
    type(field_t), pointer :: objective_field
    integer :: temp_indices(5)
    integer n
 
    call neko_scratch_registry%request_field(wo1, temp_indices(1))
    call neko_scratch_registry%request_field(wo2, temp_indices(2))
    call neko_scratch_registry%request_field(wo3, temp_indices(3))
    call neko_scratch_registry%request_field(objective_field, temp_indices(4))
    call neko_scratch_registry%request_field(work, temp_indices(5))
 
    ! update_value the objective function.
    call grad(wo1%x, wo2%x, wo3%x, this%u%x, this%c_Xh)
    call field_col3(objective_field, wo1, wo1)
    call field_addcol3(objective_field, wo2, wo2)
    call field_addcol3(objective_field, wo3, wo3)
 
    call grad(wo1%x, wo2%x, wo3%x, this%v%x, this%c_Xh)
    call field_addcol3(objective_field, wo1, wo1)
    call field_addcol3(objective_field, wo2, wo2)
    call field_addcol3(objective_field, wo3, wo3)
 
    call grad(wo1%x, wo2%x, wo3%x, this%w%x, this%c_Xh)
    call field_addcol3(objective_field, wo1, wo1)
    call field_addcol3(objective_field, wo2, wo2)
    call field_addcol3(objective_field, wo3, wo3)
 
    ! integrate the field
    n = wo1%size()
    if (this%has_mask) then
       if (neko_bcknd_device .eq. 1) then
          ! note, this could be done more elagantly by writing
          ! device_glsc2_mask
          call field_copy(work, objective_field)
          call mask_exterior_const(work, this%mask, 0.0_rp)
          this%value = device_glsc2(work%x_d, this%c_xh%B_d, n)
       else
          this%value = glsc2_mask(objective_field%x, this%C_Xh%b, &
               n, this%mask%mask%get(), this%mask%size)
       end if
    else
       if (neko_bcknd_device .eq. 1) then
          this%value = device_glsc2(objective_field%x_d, &
               this%C_Xh%b_d, n)
       else
          this%value = glsc2(objective_field%x, this%C_Xh%b, n)
       end if
    end if
 
    this%value = this%value * 0.5_rp
 
    call neko_scratch_registry%relinquish_field(temp_indices)
 
  end subroutine minimum_dissipation_update_value
 
  subroutine minimum_dissipation_update_sensitivity(this, design)
    class(minimum_dissipation_objective_t), intent(inout) :: this
    class(design_t), intent(in) :: design
    type(field_t), pointer :: work
    integer :: temp_indices(1)
 
    call neko_scratch_registry%request_field(work, temp_indices(1))
 
    ! here it should just be an inner product between the forward and adjoint
    call field_col3(work, this%u, this%adjoint_u)
    call field_addcol3(work, this%v, this%adjoint_v)
    call field_addcol3(work, this%w, this%adjoint_w)
    ! but negative
    call field_cmult(work, -1.0_rp)
 
    if (neko_bcknd_device .eq. 1) then
       call device_copy(this%sensitivity%x_d, work%x_d, this%sensitivity%size())
    else
       call copy(this%sensitivity%x, work%x, this%sensitivity%size())
    end if
 
    call neko_scratch_registry%relinquish_field(temp_indices)
 
  end subroutine minimum_dissipation_update_sensitivity
 
end module minimum_dissipation_objective