adjoint_curl_curl.f90

Go to the documentation of this file.

 
!
module adjoint_curl_curl
  use num_types, only: rp
  use field, only: field_t
  use coefs, only: coef_t
  use gather_scatter, only: gs_t
  use gs_ops, only: gs_op_add
  use operators, only: opgrad, cdtp
  use field_math, only: field_sub2, field_cmult
  use neko_config, only: neko_bcknd_device
  use math, only: col2
  use device_math, only: device_col2
  use facet_normal, only: facet_normal_t
  use comm, only: neko_comm
  use mpi_f08, only: mpi_allreduce, mpi_in_place, mpi_logical, mpi_lor
  use utils, only: neko_error
  implicit none
  private
 
  public :: adjoint_curl_curl_load, adjoint_curl_curl_check
 
contains
 
  subroutine adjoint_curl_curl_check(prim_sym, adj_sym, cyclic)
    type(facet_normal_t), intent(in) :: prim_sym, adj_sym
    logical, intent(in) :: cyclic
    logical :: has_sym
 
    has_sym = prim_sym%marked_facet%size() .gt. 0 .or. &
         adj_sym%marked_facet%size() .gt. 0
    call mpi_allreduce(mpi_in_place, has_sym, 1, mpi_logical, mpi_lor, &
         neko_comm)
 
    if (has_sym) then
       call neko_error("Symmetry boundaries " // &
            "(case.fluid.boundary_conditions or " // &
            "case.adjoint_fluid.boundary_conditions) are not supported " // &
            "by the adjoint: the transpose of the curl-curl pressure " // &
            "term on symmetry surfaces is not implemented.")
    end if
 
    if (cyclic) then
       call neko_error("Cyclic boundaries (case.fluid.cyclic) are not " // &
            "supported by the adjoint: the transpose of the cyclic " // &
            "rotations in the curl-curl pressure term is not implemented.")
    end if
 
  end subroutine adjoint_curl_curl_check
 
  subroutine adjoint_curl_curl_load(load_x, load_y, load_z, p_adj, coef, gs, &
       mu, rho, work1, work2, work3, work4)
    type(field_t), intent(inout) :: load_x, load_y, load_z
    type(field_t), intent(in) :: p_adj
    type(coef_t), intent(in) :: coef
    type(gs_t), intent(inout) :: gs
    real(kind=rp), intent(in) :: mu, rho
    type(field_t), intent(inout) :: work1, work2, work3, work4
 
    ! g = Binv gs(G p), so that B g = M^T G p; held in load
    call opgrad(load_x%x, load_y%x, load_z%x, p_adj%x, coef)
    call adjoint_curl_curl_gs_binv(load_x, load_y, load_z, coef, gs)
 
    ! c = Binv gs(D_c^T B g), so that B c = M^T D_c^T M^T G p; held in
    ! work1..3
    call adjoint_curl_transpose(work1, work2, work3, load_x, load_y, load_z, &
         work4, coef)
    call adjoint_curl_curl_gs_binv(work1, work2, work3, coef, gs)
 
    ! L = (mu / rho) D_c^T B c
    call adjoint_curl_transpose(load_x, load_y, load_z, work1, work2, work3, &
         work4, coef)
    call field_cmult(load_x, mu / rho)
    call field_cmult(load_y, mu / rho)
    call field_cmult(load_z, mu / rho)
 
  end subroutine adjoint_curl_curl_load
 
  subroutine adjoint_curl_transpose(w1, w2, w3, y1, y2, y3, work, coef)
    type(field_t), intent(inout) :: w1, w2, w3
    type(field_t), intent(inout) :: y1, y2, y3
    type(field_t), intent(inout) :: work
    type(coef_t), intent(in) :: coef
 
    call cdtp(w1%x, y2%x, coef%drdz, coef%dsdz, coef%dtdz, coef)
    call cdtp(work%x, y3%x, coef%drdy, coef%dsdy, coef%dtdy, coef)
    call field_sub2(w1, work)
 
    call cdtp(w2%x, y3%x, coef%drdx, coef%dsdx, coef%dtdx, coef)
    call cdtp(work%x, y1%x, coef%drdz, coef%dsdz, coef%dtdz, coef)
    call field_sub2(w2, work)
 
    call cdtp(w3%x, y1%x, coef%drdy, coef%dsdy, coef%dtdy, coef)
    call cdtp(work%x, y2%x, coef%drdx, coef%dsdx, coef%dtdx, coef)
    call field_sub2(w3, work)
 
  end subroutine adjoint_curl_transpose
 
  subroutine adjoint_curl_curl_gs_binv(y1, y2, y3, coef, gs)
    type(field_t), intent(inout) :: y1, y2, y3
    type(coef_t), intent(in) :: coef
    type(gs_t), intent(inout) :: gs
    integer :: n
 
    n = y1%size()
 
    call gs%op(y1, gs_op_add)
    call gs%op(y2, gs_op_add)
    call gs%op(y3, gs_op_add)
 
    if (neko_bcknd_device .eq. 1) then
       call device_col2(y1%x_d, coef%Binv_d, n)
       call device_col2(y2%x_d, coef%Binv_d, n)
       call device_col2(y3%x_d, coef%Binv_d, n)
    else
       call col2(y1%x, coef%Binv, n)
       call col2(y2%x, coef%Binv, n)
       call col2(y3%x, coef%Binv, n)
    end if
 
  end subroutine adjoint_curl_curl_gs_binv
 
end module adjoint_curl_curl