drag_torque.f90

Go to the documentation of this file.

! Copyright (c) 2008-2020, UCHICAGO ARGONNE, LLC.
! Copyright (c) 2025, The Neko Authors
!
! The UChicago Argonne, LLC as Operator of Argonne National
! Laboratory holds copyright in the Software. The copyright holder
! reserves all rights except those expressly granted to licensees,
! and U.S. Government license rights.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
!
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the disclaimer below.
!
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions and the disclaimer (as noted below)
! in the documentation and/or other materials provided with the
! distribution.
!
! 3. Neither the name of ANL 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
! UCHICAGO ARGONNE, LLC, THE U.S. DEPARTMENT OF
! ENERGY 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.
!
! Additional BSD Notice
! ---------------------
! 1. This notice is required to be provided under our contract with
! the U.S. Department of Energy (DOE). This work was produced at
! Argonne National Laboratory under Contract
! No. DE-AC02-06CH11357 with the DOE.
!
! 2. Neither the United States Government nor UCHICAGO ARGONNE,
! LLC nor any of their employees, makes any warranty,
! express or implied, or assumes any liability or responsibility for the
! accuracy, completeness, or usefulness of any information, apparatus,
! product, or process disclosed, or represents that its use would not
! infringe privately-owned rights.
!
! 3. Also, reference herein to any specific commercial products, process,
! or services by trade name, trademark, manufacturer or otherwise does
! not necessarily constitute or imply its endorsement, recommendation,
! or favoring by the United States Government or UCHICAGO ARGONNE LLC.
! The views and opinions of authors expressed
! herein do not necessarily state or reflect those of the United States
! Government or UCHICAGO ARGONNE, LLC, and shall
! not be used for advertising or product endorsement purposes.
!
 
module drag_torque
  use field, only : field_t
  use coefs, only : coef_t
  use vector, only : vector_t
  use facet_zone, only : facet_zone_t
  use math, only : rzero, col3, vdot3, col2
  use space, only : space_t
  use num_types, only : rp
  use utils, only : nonlinear_index, neko_error
  use iso_c_binding, only : c_ptr
  use neko_config, only : neko_bcknd_device
  use device, only : host_to_device
  use device_math, only : device_cmult, device_col2, &
       device_col3, device_vdot3, device_rzero
  use comm, only : neko_comm, mpi_real_precision
  use mpi_f08, only : mpi_allreduce, mpi_in_place, mpi_sum
  implicit none
  private
  public :: drag_torque_zone, drag_torque_facet, drag_torque_pt, &
       setup_normals, calc_force_array, device_calc_force_array
 
contains
  subroutine drag_torque_zone(dgtq, tstep, zone, center, s11, s22, s33, s12, &
       s13, s23, p, coef, visc)
    integer, intent(in) :: tstep
    type(facet_zone_t) :: zone
    type(coef_t), intent(inout) :: coef
    real(kind=rp), intent(inout) :: s11(coef%Xh%lx, coef%Xh%ly, coef%Xh%lz, &
         coef%msh%nelv)
    real(kind=rp), intent(inout) :: s22(coef%Xh%lx, coef%Xh%ly, coef%Xh%lz, &
         coef%msh%nelv)
    real(kind=rp), intent(inout) :: s33(coef%Xh%lx, coef%Xh%ly, coef%Xh%lz, &
         coef%msh%nelv)
    real(kind=rp), intent(inout) :: s12(coef%Xh%lx, coef%Xh%ly, coef%Xh%lz, &
         coef%msh%nelv)
    real(kind=rp), intent(inout) :: s13(coef%Xh%lx, coef%Xh%ly, coef%Xh%lz, &
         coef%msh%nelv)
    real(kind=rp), intent(inout) :: s23(coef%Xh%lx, coef%Xh%ly, coef%Xh%lz, &
         coef%msh%nelv)
    type(field_t), intent(inout) :: p
    real(kind=rp), intent(in) :: visc, center(3)
    real(kind=rp) :: dgtq(3,4)
    real(kind=rp) :: dragpx = 0.0_rp ! pressure
    real(kind=rp) :: dragpy = 0.0_rp
    real(kind=rp) :: dragpz = 0.0_rp
    real(kind=rp) :: dragvx = 0.0_rp ! viscous
    real(kind=rp) :: dragvy = 0.0_rp
    real(kind=rp) :: dragvz = 0.0_rp
    real(kind=rp) :: torqpx = 0.0_rp ! pressure
    real(kind=rp) :: torqpy = 0.0_rp
    real(kind=rp) :: torqpz = 0.0_rp
    real(kind=rp) :: torqvx = 0.0_rp ! viscous
    real(kind=rp) :: torqvy = 0.0_rp
    real(kind=rp) :: torqvz = 0.0_rp
    real(kind=rp) :: dragx, dragy, dragz
    integer :: ie, ifc, mem, ierr
    dragx = 0.0
    dragy = 0.0
    dragz = 0.0
 
!
!     Fill up viscous array w/ default
!
    dragpx = 0.0
    dragpy = 0.0
    dragpz = 0.0
    dragvx = 0.0
    dragvy = 0.0
    dragvz = 0.0
    do mem = 1, zone%size
       ie = zone%facet_el(mem)%x(2)
       ifc = zone%facet_el(mem)%x(1)
       call drag_torque_facet(dgtq, coef%dof%x, coef%dof%y, coef%dof%z, &
            center, s11, s22, s33, s12, s13, s23, p%x, visc, ifc, ie, coef, &
            coef%Xh)
 
       dragpx = dragpx + dgtq(1,1) ! pressure
       dragpy = dragpy + dgtq(2,1)
       dragpz = dragpz + dgtq(3,1)
 
       dragvx = dragvx + dgtq(1,2) ! viscous
       dragvy = dragvy + dgtq(2,2)
       dragvz = dragvz + dgtq(3,2)
 
       torqpx = torqpx + dgtq(1,3) ! pressure
       torqpy = torqpy + dgtq(2,3)
       torqpz = torqpz + dgtq(3,3)
 
       torqvx = torqvx + dgtq(1,4) ! viscous
       torqvy = torqvy + dgtq(2,4)
       torqvz = torqvz + dgtq(3,4)
    end do
!
!     Sum contributions from all processors
!
    call mpi_allreduce(mpi_in_place, dragpx, 1, &
         mpi_real_precision, mpi_sum, neko_comm, ierr)
    call mpi_allreduce(mpi_in_place, dragpy, 1, &
         mpi_real_precision, mpi_sum, neko_comm, ierr)
    call mpi_allreduce(mpi_in_place, dragpz, 1, &
         mpi_real_precision, mpi_sum, neko_comm, ierr)
    call mpi_allreduce(mpi_in_place, dragvx, 1, &
         mpi_real_precision, mpi_sum, neko_comm, ierr)
    call mpi_allreduce(mpi_in_place, dragvy, 1, &
         mpi_real_precision, mpi_sum, neko_comm, ierr)
    call mpi_allreduce(mpi_in_place, dragvz, 1, &
         mpi_real_precision, mpi_sum, neko_comm, ierr)
    !Torque
    call mpi_allreduce(mpi_in_place, torqpx, 1, &
         mpi_real_precision, mpi_sum, neko_comm, ierr)
    call mpi_allreduce(mpi_in_place, torqpy, 1, &
         mpi_real_precision, mpi_sum, neko_comm, ierr)
    call mpi_allreduce(mpi_in_place, torqpz, 1, &
         mpi_real_precision, mpi_sum, neko_comm, ierr)
    call mpi_allreduce(mpi_in_place, torqvx, 1, &
         mpi_real_precision, mpi_sum, neko_comm, ierr)
    call mpi_allreduce(mpi_in_place, torqvy, 1, &
         mpi_real_precision, mpi_sum, neko_comm, ierr)
    call mpi_allreduce(mpi_in_place, torqvz, 1, &
         mpi_real_precision, mpi_sum, neko_comm, ierr)
 
    dgtq(1,1) = dragpx ! pressure
    dgtq(2,1) = dragpy
    dgtq(3,1) = dragpz
 
    dgtq(1,2) = dragvx ! viscous
    dgtq(2,2) = dragvy
    dgtq(3,2) = dragvz
 
    dgtq(1,3) = torqpx ! pressure
    dgtq(2,3) = torqpy
    dgtq(3,3) = torqpz
 
    dgtq(1,4) = torqvx ! viscous
    dgtq(2,4) = torqvy
    dgtq(3,4) = torqvz
 
  end subroutine drag_torque_zone
 
  subroutine drag_torque_facet(dgtq, xm0, ym0, zm0, center, s11, s22, s33, &
       s12, s13, s23, pm1, visc, f, e, coef, Xh)
    type(coef_t), intent(in) :: coef
    type(space_t), intent(in) :: xh
    real(kind=rp), intent(out) :: dgtq(3,4)
    real(kind=rp), intent(in) :: center(3)
    real(kind=rp), intent(in) :: xm0(xh%lx, xh%ly, xh%lz, coef%msh%nelv)
    real(kind=rp), intent(in) :: ym0(xh%lx, xh%ly, xh%lz, coef%msh%nelv)
    real(kind=rp), intent(in) :: zm0(xh%lx, xh%ly, xh%lz, coef%msh%nelv)
    real(kind=rp), intent(in) :: s11(xh%lx, xh%ly, xh%lz, coef%msh%nelv)
    real(kind=rp), intent(in) :: s22(xh%lx, xh%ly, xh%lz, coef%msh%nelv)
    real(kind=rp), intent(in) :: s33(xh%lx, xh%ly, xh%lz, coef%msh%nelv)
    real(kind=rp), intent(in) :: s12(xh%lx, xh%ly, xh%lz, coef%msh%nelv)
    real(kind=rp), intent(in) :: s13(xh%lx, xh%ly, xh%lz, coef%msh%nelv)
    real(kind=rp), intent(in) :: s23(xh%lx, xh%ly, xh%lz, coef%msh%nelv)
    real(kind=rp), intent(in) :: pm1(xh%lx, xh%ly, xh%lz, coef%msh%nelv)
    real(kind=rp), intent(in) :: visc
    integer, intent(in) :: f, e
    integer :: pf, i, j1, j2
    real(kind=rp) :: n1, n2, n3, a, v, dgtq_i(3,4)
    integer :: skpdat(6,6), nx, ny, nz
    integer :: js1
    integer :: jf1
    integer :: jskip1
    integer :: js2
    integer :: jf2
    integer :: jskip2
    real(kind=rp) :: s11_, s12_, s22_, s13_, s23_, s33_
 
 
    nx = xh%lx
    ny = xh%ly
    nz = xh%lz
    skpdat(1,1) = 1
    skpdat(2,1) = nx*(ny-1) + 1
    skpdat(3,1) = nx
    skpdat(4,1) = 1
    skpdat(5,1) = ny*(nz-1) + 1
    skpdat(6,1) = ny
 
    skpdat(1,2) = 1 + (nx-1)
    skpdat(2,2) = nx*(ny-1)+1 + (nx-1)
    skpdat(3,2) = nx
    skpdat(4,2) = 1
    skpdat(5,2) = ny*(nz-1)+1
    skpdat(6,2) = ny
 
    skpdat(1,3) = 1
    skpdat(2,3) = nx
    skpdat(3,3) = 1
    skpdat(4,3) = 1
    skpdat(5,3) = ny*(nz-1)+1
    skpdat(6,3) = ny
 
    skpdat(1,4) = 1 + nx*(ny-1)
    skpdat(2,4) = nx + nx*(ny-1)
    skpdat(3,4) = 1
    skpdat(4,4) = 1
    skpdat(5,4) = ny*(nz-1)+1
    skpdat(6,4) = ny
 
    skpdat(1,5) = 1
    skpdat(2,5) = nx
    skpdat(3,5) = 1
    skpdat(4,5) = 1
    skpdat(5,5) = ny
    skpdat(6,5) = 1
 
    skpdat(1,6) = 1 + nx*ny*(nz-1)
    skpdat(2,6) = nx + nx*ny*(nz-1)
    skpdat(3,6) = 1
    skpdat(4,6) = 1
    skpdat(5,6) = ny
    skpdat(6,6) = 1
    pf = f
    js1 = skpdat(1, pf)
    jf1 = skpdat(2, pf)
    jskip1 = skpdat(3, pf)
    js2 = skpdat(4, pf)
    jf2 = skpdat(5, pf)
    jskip2 = skpdat(6, pf)
    call rzero(dgtq,12)
    i = 0
    a = 0
    do j2 = js2, jf2, jskip2
       do j1 = js1, jf1, jskip1
          i = i + 1
          n1 = coef%nx(i,1,f,e) * coef%area(i,1,f,e)
          n2 = coef%ny(i,1,f,e) * coef%area(i,1,f,e)
          n3 = coef%nz(i,1,f,e) * coef%area(i,1,f,e)
          a = a + coef%area(i,1,f,e)
          v = visc
          s11_ = s11(j1, j2, 1, e)
          s12_ = s12(j1, j2, 1, e)
          s22_ = s22(j1, j2, 1, e)
          s13_ = s13(j1, j2, 1, e)
          s23_ = s23(j1, j2, 1, e)
          s33_ = s33(j1, j2, 1, e)
          call drag_torque_pt(dgtq_i, xm0(j1, j2, 1, e), ym0(j1, j2, 1, e), &
               zm0(j1, j2, 1, e), center, s11_, s22_, s33_, s12_, s13_, s23_, &
               pm1(j1, j2, 1, e), n1, n2, n3, v)
          dgtq = dgtq + dgtq_i
       end do
    end do
  end subroutine drag_torque_facet
 
  subroutine drag_torque_pt(dgtq, x, y, z, center, s11, s22, s33, s12, s13, &
       s23, p, n1, n2, n3, v)
    real(kind=rp), intent(inout) :: dgtq(3, 4)
    real(kind=rp), intent(in) :: x
    real(kind=rp), intent(in) :: y
    real(kind=rp), intent(in) :: z
    real(kind=rp), intent(in) :: p
    real(kind=rp), intent(in) :: v
    real(kind=rp), intent(in) :: n1, n2, n3, center(3)
    real(kind=rp), intent(in) :: s11, s12, s22, s13, s23, s33
    real(kind=rp) :: s21, s31, s32, r1, r2, r3
 
    call rzero(dgtq, 12)
    s21 = s12
    s32 = s23
    s31 = s13
    !pressure drag
    dgtq(1,1) = p*n1
    dgtq(2,1) = p*n2
    dgtq(3,1) = p*n3
    ! viscous drag
    dgtq(1,2) = -2*v*(s11*n1 + s12*n2 + s13*n3)
    dgtq(2,2) = -2*v*(s21*n1 + s22*n2 + s23*n3)
    dgtq(3,2) = -2*v*(s31*n1 + s32*n2 + s33*n3)
    r1 = x - center(1)
    r2 = y - center(2)
    r3 = z - center(3)
    !pressure torque
    dgtq(1,3) = r2*dgtq(3,1) - r3*dgtq(2,1)
    dgtq(2,3) = r3*dgtq(1,1) - r1*dgtq(3,1)
    dgtq(3,3) = r1*dgtq(2,1) - r2*dgtq(1,1)
    !viscous torque
    dgtq(1,4) = r2*dgtq(3,2) - r3*dgtq(2,2)
    dgtq(2,4) = r3*dgtq(1,2) - r1*dgtq(3,2)
    dgtq(3,4) = r1*dgtq(2,2) - r2*dgtq(1,2)
  end subroutine drag_torque_pt
 
  subroutine calc_force_array(force1, force2, force3, force4, force5, force6,&
       s11, s22, s33, s12, s13, s23,&
       p, n1, n2, n3, mu, n_pts)
    integer :: n_pts
    real(kind=rp), intent(inout), dimension(n_pts) :: force1, force2, force3
    real(kind=rp), intent(inout), dimension(n_pts) :: force4, force5, force6
    real(kind=rp), intent(in) :: p(n_pts)
    real(kind=rp), intent(in) :: mu(n_pts)
    real(kind=rp), intent(in) :: n1(n_pts), n2(n_pts), n3(n_pts)
    real(kind=rp), intent(in), dimension(n_pts) :: s11, s12, s22, s13, s23, s33
    real(kind=rp) :: v2(n_pts)
    call rzero(force4, n_pts)
    call rzero(force5, n_pts)
    call rzero(force6, n_pts)
    !pressure force
    call col3(force1, p, n1, n_pts)
    call col3(force2, p, n2, n_pts)
    call col3(force3, p, n3, n_pts)
    ! viscous force
    v2 = -2.0_rp*mu
    call vdot3(force4, s11, s12, s13, n1, n2, n3, n_pts)
    call vdot3(force5, s12, s22, s23, n1, n2, n3, n_pts)
    call vdot3(force6, s13, s23, s33, n1, n2, n3, n_pts)
    call col2(force4, v2, n_pts)
    call col2(force5, v2, n_pts)
    call col2(force6, v2, n_pts)
  end subroutine calc_force_array
 
  subroutine device_calc_force_array(force1, force2, force3,&
       force4, force5, force6,&
       s11, s22, s33, s12, s13, s23,&
       p, n1, n2, n3, mu, n_pts)
    integer :: n_pts
    type(c_ptr) :: force1, force2, force3
    type(c_ptr) :: force4, force5, force6
    type(c_ptr) :: p, n1, n2, n3
    type(c_ptr) :: s11, s12, s22, s13, s23, s33
    type(c_ptr) :: mu
    if (neko_bcknd_device .eq. 1) then
       call device_rzero(force4, n_pts)
       call device_rzero(force5, n_pts)
       call device_rzero(force6, n_pts)
       !pressure force
       call device_col3(force1, p, n1, n_pts)
       call device_col3(force2, p, n2, n_pts)
       call device_col3(force3, p, n3, n_pts)
 
       call device_vdot3(force4, s11, s12, s13, n1, n2, n3, n_pts)
       call device_vdot3(force5, s12, s22, s23, n1, n2, n3, n_pts)
       call device_vdot3(force6, s13, s23, s33, n1, n2, n3, n_pts)
       call device_cmult(force4, -2.0_rp, n_pts)
       call device_cmult(force5, -2.0_rp, n_pts)
       call device_cmult(force6, -2.0_rp, n_pts)
       call device_col2(force4, mu, n_pts)
       call device_col2(force5, mu, n_pts)
       call device_col2(force6, mu, n_pts)
    else
       call neko_error('error in drag_torque, no device bcklnd configured')
    end if
  end subroutine device_calc_force_array
 
 
  subroutine setup_normals(coef, mask, facets, n1, n2, n3, n_pts)
    type(coef_t) :: coef
    integer :: n_pts
    type(vector_t), intent(inout) :: n1, n2, n3
    integer :: mask(0:n_pts), facets(0:n_pts), fid, idx(4)
    real(kind=rp) :: normal(3), area
    integer :: i
 
    do i = 1, n_pts
       fid = facets(i)
       if ( (fid .lt. 1) .or. (fid .gt. 6)) then
          call neko_error('setup_normals: invalid facet id. ' // &
               'Finalize the bc with only_facets = .true.')
       end if
       idx = nonlinear_index(mask(i), coef%Xh%lx, coef%Xh%lx, coef%Xh%lx)
       normal = coef%get_normal(idx(1), idx(2), idx(3), idx(4), fid)
       area = coef%get_area(idx(1), idx(2), idx(3), idx(4), fid)
       n1%x(i) = normal(1)*area
       n2%x(i) = normal(2)*area
       n3%x(i) = normal(3)*area
    end do
 
    if (neko_bcknd_device .eq. 1) then
       call n1%copy_from(host_to_device, .false.)
       call n2%copy_from(host_to_device, .false.)
       call n3%copy_from(host_to_device, .true.)
    end if
  end subroutine setup_normals
 
end module drag_torque