lube_term_objective.f90

Go to the documentation of this file.

 
!
!
! 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 lube_term_objective
  use objective, only: objective_t
  use design, only: design_t
  use brinkman_design, only: brinkman_design_t
  use simulation_m, only: simulation_t
  use adjoint_lube_source_term, only: adjoint_lube_source_term_t
  use adjoint_fluid_pnpn, only: adjoint_fluid_pnpn_t
  use num_types, only: rp
  use field, only: field_t
  use scratch_registry, only: neko_scratch_registry, scratch_registry_t
  use neko_config, only: neko_bcknd_device
  use mask_ops, only: mask_exterior_const, compute_masked_volume
  use utils, only: neko_error
  use json_module, only: json_file
  use json_utils, only: json_get_or_default
  use registry, only: neko_registry
  use interpolation, only: interpolator_t
  use space, only: space_t, gl
  use coefs, only: coef_t
  use math, only: glsc2, copy, col2, invcol2
  use device_math, only: device_copy, device_glsc2, device_col2, device_invcol2
  use math_ext, only: glsc2_mask
  use field_math, only: field_col3, field_addcol3, field_cmult, field_col2
  implicit none
  private
 
  type, public, extends(objective_t) :: lube_term_objective_t
     private
 
     type(field_t), pointer :: u => null()
     type(field_t), pointer :: v => null()
     type(field_t), pointer :: w => null()
     type(field_t), pointer :: brinkman_amplitude => null()
     logical :: dealias_sensitivity
     logical :: dealias_forcing
     real(kind=rp) :: volume
 
     ! ---- everything GLL ----
     type(coef_t), pointer :: c_xh_gll
     type(space_t), pointer :: xh_gll
 
     ! ---- everything for GL ----
     type(space_t), pointer :: xh_gl
     type(coef_t), pointer :: c_xh_gl
     type(interpolator_t), pointer :: gll_to_gl
     type(scratch_registry_t), pointer :: scratch_gl
     integer :: gdim
 
 
   contains
 
     procedure, public, pass(this) :: init_json_sim => lube_term_init_json_sim
     procedure, public, pass(this) :: init_from_attributes => &
          lube_term_init_attributes
     procedure, public, pass(this) :: free => lube_term_free
     procedure, public, pass(this) :: update_value => &
          lube_term_update_value
     procedure, public, pass(this) :: update_sensitivity => &
          lube_term_update_sensitivity
 
  end type lube_term_objective_t
 
contains
 
  subroutine lube_term_init_json_sim(this, json, design, simulation)
    class(lube_term_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 :: mask_name
    character(len=:), allocatable :: name
    real(kind=rp) :: weight
    logical :: dealias_sensitivity, dealias_forcing
 
    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, "Out of plane stresses")
    call json_get_or_default(json, "dealias_sensitivity", &
         dealias_sensitivity, .true.)
    call json_get_or_default(json, "dealias_forcing", &
         dealias_forcing, .true.)
 
    call this%init_from_attributes(design, simulation, weight, name, &
         mask_name, dealias_sensitivity, dealias_forcing)
  end subroutine lube_term_init_json_sim
 
  subroutine lube_term_init_attributes(this, design, simulation, weight, &
       name, mask_name, dealias_sensitivity, dealias_forcing)
    class(lube_term_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) :: mask_name
    character(len=*), intent(in) :: name
    logical, intent(in) :: dealias_sensitivity
    logical, intent(in) :: dealias_forcing
    type(adjoint_lube_source_term_t) :: lube_term
 
    ! Call the base initializer
    call this%init_base(name, design%size(), weight, mask_name)
 
    this%dealias_forcing = dealias_forcing
    this%dealias_sensitivity = dealias_sensitivity
 
    ! Grab the brinkman amplitude for the lube term
    select type (design)
    type is (brinkman_design_t)
       this%brinkman_amplitude => &
            neko_registry%get_field("brinkman_amplitude")
 
 
    class default
       call neko_error('Minimum dissipation only works with brinkman_design')
    end select
 
    this%u => neko_registry%get_field('u')
    this%v => neko_registry%get_field('v')
    this%w => neko_registry%get_field('w')
 
    ! GLL
    this%c_Xh_GLL => simulation%neko_case%fluid%c_Xh
    this%Xh_GLL => simulation%neko_case%fluid%c_Xh%Xh
    this%gdim = this%c_Xh_GLL%msh%gdim
 
    ! GL
    this%c_Xh_GL => simulation%adjoint_case%fluid_adj%c_Xh_GL
    this%Xh_GL => this%c_Xh_GL%Xh
 
    ! GLL to GL
    this%GLL_to_GL => simulation%adjoint_case%fluid_adj%GLL_to_GL
 
    this%scratch_GL => simulation%adjoint_case%fluid_adj%scratch_GL
 
    ! compute the volume of the objective domain
    if (this%has_mask) then
       this%volume = compute_masked_volume(this%mask, this%c_Xh_GLL)
    else
       this%volume = this%c_Xh_GLL%volume
    end if
 
    ! if we have the lube term we need to initialize and append that too
 
    associate(f_adj_x => simulation%adjoint_fluid%f_adj_x, &
         f_adj_y => simulation%adjoint_fluid%f_adj_y, &
         f_adj_z => simulation%adjoint_fluid%f_adj_z)
 
      call lube_term%init_from_components(f_adj_x, f_adj_y, f_adj_z, design, &
           this%weight, this%u, this%v, this%w, this%mask, &
           this%has_mask, this%c_Xh_GLL, this%c_Xh_GL, this%GLL_to_GL, &
           this%dealias_forcing, this%volume, &
           this%scratch_GL, this%gdim)
 
    end associate
 
    ! 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(lube_term)
    end select
 
  end subroutine lube_term_init_attributes
 
  subroutine lube_term_free(this)
    class(lube_term_objective_t), intent(inout) :: this
    call this%free_base()
 
    this%u => null()
    this%v => null()
    this%w => null()
    this%c_Xh_GLL => null()
    this%brinkman_amplitude => null()
    nullify(this%c_Xh_GL)
    nullify(this%Xh_GL)
    nullify(this%Xh_GLL)
    nullify(this%GLL_to_GL)
    nullify(this%scratch_GL)
 
  end subroutine lube_term_free
 
  subroutine lube_term_update_value(this, design)
    class(lube_term_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), .false.)
 
    call field_col3(work, this%u, this%u)
    call field_addcol3(work, this%v, this%v)
    if (this%gdim .eq. 3) then
       call field_addcol3(work, this%w, this%w)
    end if
    call field_col2(work, this%brinkman_amplitude)
 
    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 mask_exterior_const(work, this%mask, 0.0_rp)
          this%value = device_glsc2(work%x_d, this%c_Xh_GLL%B_d, design%size())
       else
          this%value = glsc2_mask(work%x, this%c_Xh_GLL%B, design%size(), &
               this%mask%mask%get(), this%mask%size)
       end if
    else
       if (neko_bcknd_device .eq. 1) then
          this%value = device_glsc2(work%x_d, this%c_Xh_GLL%B_d, design%size())
       else
          this%value = glsc2(work%x, this%c_Xh_GLL%B, design%size())
       end if
    end if
    this%value = 0.5_rp * this%value / this%volume
 
    call neko_scratch_registry%relinquish_field(temp_indices)
 
  end subroutine lube_term_update_value
 
  subroutine lube_term_update_sensitivity(this, design)
    class(lube_term_objective_t), intent(inout) :: this
    class(design_t), intent(in) :: design
    type(field_t), pointer :: work
    integer :: temp_indices(1)
    integer :: n_GL, nel
    type(field_t), pointer :: accumulate, fld_GL
    integer :: temp_indices_GL(2)
 
    ! if we have the lube term we also get an extra term in the sensitivity
 
    call neko_scratch_registry%request_field(work, temp_indices(1), .false.)
 
    if(this%dealias_sensitivity) then
       nel = this%c_Xh_GLL%msh%nelv
       n_gl = nel * this%Xh_GL%lxyz
       call this%scratch_GL%request_field(accumulate, temp_indices_gl(1), &
            .false.)
       call this%scratch_GL%request_field(fld_gl, temp_indices_gl(2), .false.)
 
       call this%GLL_to_GL%map(fld_gl%x, this%u%x, nel, this%Xh_GL)
       call field_col3(accumulate, fld_gl, fld_gl)
       call this%GLL_to_GL%map(fld_gl%x, this%v%x, nel, this%Xh_GL)
       call field_addcol3(accumulate, fld_gl, fld_gl)
       if (this%gdim .eq. 3) then
          call this%GLL_to_GL%map(fld_gl%x, this%w%x, nel, this%Xh_GL)
          call field_addcol3(accumulate, fld_gl, fld_gl)
       end if
       ! scale
       call field_cmult(accumulate, this%weight * 0.5_rp / this%volume)
 
       ! Evaluate term on GL and preempt the GLL premultiplication
       if (neko_bcknd_device .eq. 1) then
          call device_col2(accumulate%x_d, this%c_Xh_GL%B_d, n_gl)
          call this%GLL_to_GL%map(work%x, accumulate%x, nel, this%Xh_GLL)
          call device_invcol2(work%x_d, this%c_Xh_GLL%B_d, work%size())
       else
          call col2(accumulate%x, this%c_Xh_GL%B, n_gl)
          call this%GLL_to_GL%map(work%x, accumulate%x, nel, this%Xh_GLL)
          call invcol2(work%x, this%c_Xh_GLL%B, work%size())
       end if
 
       call this%scratch_GL%relinquish_field(temp_indices_gl)
 
    else
       call field_col3(work, this%u, this%u)
       call field_addcol3(work, this%v, this%v)
       if (this%gdim .eq. 3) then
          call field_addcol3(work, this%w, this%w)
       end if
       ! scale
       call field_cmult(work, this%weight * 0.5_rp / this%volume)
    end if
 
    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 lube_term_update_sensitivity
 
end module lube_term_objective