operators.f90

Go to the documentation of this file.

! Copyright (c) 2020-2024, 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 operators
  use neko_config, only : neko_bcknd_sx, neko_bcknd_device, neko_bcknd_xsmm,&
                          neko_device_mpi
  use num_types, only : rp
  use opr_cpu, only : opr_cpu_cfl, opr_cpu_curl, opr_cpu_opgrad, opr_cpu_conv1,&
                      opr_cpu_cdtp, opr_cpu_dudxyz, opr_cpu_lambda2
  use opr_sx, only : opr_sx_cfl, opr_sx_curl, opr_sx_dudxyz, opr_sx_opgrad, &
                     opr_sx_cdtp, opr_sx_conv1, opr_sx_lambda2
  use opr_xsmm, only : opr_xsmm_cdtp, opr_xsmm_conv1, opr_xsmm_curl, &
                       opr_xsmm_dudxyz, opr_xsmm_opgrad
  use opr_device, only : opr_device_cdtp, opr_device_cfl, opr_device_curl, &
                         opr_device_conv1, opr_device_dudxyz, &
                         opr_device_lambda2, opr_device_opgrad
  use space, only : space_t
  use coefs, only : coef_t
  use field, only : field_t
  use math, only : glsum, cmult, add2, cadd
  use device, only : c_ptr, device_get_ptr
  use device_math, only : device_add2, device_cmult
  use comm
  implicit none
  private
 
  public :: dudxyz, opgrad, ortho, cdtp, conv1, curl, cfl,&
            lambda2op, strain_rate
  
contains
 
  subroutine dudxyz (du, u, dr, ds, dt, coef)
    type(coef_t), intent(in), target :: coef
    real(kind=rp), dimension(coef%Xh%lx,coef%Xh%ly,coef%Xh%lz,coef%msh%nelv), &
         intent(inout) ::  du
    real(kind=rp), dimension(coef%Xh%lx,coef%Xh%ly,coef%Xh%lz,coef%msh%nelv), &
         intent(in) ::  u, dr, ds, dt
 
    if (neko_bcknd_sx .eq. 1) then
       call opr_sx_dudxyz(du, u, dr, ds, dt, coef)
    else if (neko_bcknd_xsmm .eq. 1) then
       call opr_xsmm_dudxyz(du, u, dr, ds, dt, coef)
    else if (neko_bcknd_device .eq. 1) then
       call opr_device_dudxyz(du, u, dr, ds, dt, coef)
    else
       call opr_cpu_dudxyz(du, u, dr, ds, dt, coef)
    end if
 
  end subroutine dudxyz
 
  subroutine opgrad(ux, uy, uz, u, coef, es, ee)
    type(coef_t), intent(in) :: coef
    real(kind=rp), dimension(coef%Xh%lxyz,coef%msh%nelv), intent(inout) :: ux
    real(kind=rp), dimension(coef%Xh%lxyz,coef%msh%nelv), intent(inout) :: uy
    real(kind=rp), dimension(coef%Xh%lxyz,coef%msh%nelv), intent(inout) :: uz
    real(kind=rp), dimension(coef%Xh%lxyz,coef%msh%nelv), intent(in) :: u
    integer, optional :: es, ee
    integer :: eblk_start, eblk_end
 
    if (present(es)) then
       eblk_start = es
    else
       eblk_start = 1
    end if
 
    if (present(ee)) then
       eblk_end = ee
    else
       eblk_end = coef%msh%nelv
    end if
 
    if (neko_bcknd_sx .eq. 1) then
       call opr_sx_opgrad(ux, uy, uz, u, coef)
    else if (neko_bcknd_xsmm .eq. 1) then
       call opr_xsmm_opgrad(ux, uy, uz, u, coef)
    else if (neko_bcknd_device .eq. 1) then
       call opr_device_opgrad(ux, uy, uz, u, coef)
    else
       call opr_cpu_opgrad(ux, uy, uz, u, coef, eblk_start, eblk_end)
    end if
 
  end subroutine opgrad
 
  subroutine ortho(x, n, glb_n)
    integer, intent(in) :: n
    integer, intent(in) :: glb_n
    real(kind=rp), dimension(n), intent(inout) :: x
    real(kind=rp) :: rlam
 
    rlam = glsum(x, n)/glb_n
    call cadd(x, -rlam, n)
 
  end subroutine ortho
 
  subroutine cdtp (dtx, x, dr, ds, dt, coef)
    type(coef_t), intent(in) :: coef
    real(kind=rp), dimension(coef%Xh%lxyz,coef%msh%nelv), intent(inout) :: dtx
    real(kind=rp), dimension(coef%Xh%lxyz,coef%msh%nelv), intent(inout) :: x
    real(kind=rp), dimension(coef%Xh%lxyz,coef%msh%nelv), intent(in) :: dr
    real(kind=rp), dimension(coef%Xh%lxyz,coef%msh%nelv), intent(in) :: ds
    real(kind=rp), dimension(coef%Xh%lxyz,coef%msh%nelv), intent(in) :: dt
 
    if (neko_bcknd_sx .eq. 1) then
       call opr_sx_cdtp(dtx, x, dr, ds, dt, coef)
    else if (neko_bcknd_xsmm .eq. 1) then
       call opr_xsmm_cdtp(dtx, x, dr, ds, dt, coef)
    else if (neko_bcknd_device .eq. 1) then
       call opr_device_cdtp(dtx, x, dr, ds, dt, coef)
    else
       call opr_cpu_cdtp(dtx, x, dr, ds, dt, coef)
    end if
 
  end subroutine cdtp
 
  subroutine conv1(du, u, vx, vy, vz, Xh, coef, es, ee)
    type(space_t), intent(inout) :: xh
    type(coef_t), intent(inout) :: coef
    real(kind=rp), intent(inout) :: du(xh%lxyz,coef%msh%nelv)
    real(kind=rp), intent(inout) :: u(xh%lx,xh%ly,xh%lz,coef%msh%nelv)
    real(kind=rp), intent(inout) :: vx(xh%lx,xh%ly,xh%lz,coef%msh%nelv)
    real(kind=rp), intent(inout) :: vy(xh%lx,xh%ly,xh%lz,coef%msh%nelv)
    real(kind=rp), intent(inout) :: vz(xh%lx,xh%ly,xh%lz,coef%msh%nelv)
    integer, optional :: es, ee
    integer :: eblk_end, eblk_start
 
    associate(nelv => coef%msh%nelv, gdim => coef%msh%gdim)
      if (present(es)) then
         eblk_start = es
      else
         eblk_start = 1
      end if
 
      if (present(ee)) then
         eblk_end = ee
      else
         eblk_end = coef%msh%nelv
      end if
 
      if (neko_bcknd_sx .eq. 1) then
         call opr_sx_conv1(du, u, vx, vy, vz, xh, coef, nelv, gdim)
      else if (neko_bcknd_xsmm .eq. 1) then
         call opr_xsmm_conv1(du, u, vx, vy, vz, xh, coef, nelv, gdim)
      else if (neko_bcknd_device .eq. 1) then
         call opr_device_conv1(du, u, vx, vy, vz, xh, coef, nelv, gdim)
      else
         call opr_cpu_conv1(du, u, vx, vy, vz, xh, coef, eblk_start, eblk_end)
      end if
    end associate
 
  end subroutine conv1
 
  !! Compute the curl fo a vector field.
  !! @param w1 Will store the x component of the curl.
  !! @param w2 Will store the y component of the curl.
  !! @param w3 Will store the z component of the curl.
  !! @param u1 The x component of the vector field.
  !! @param u2 The y component of the vector field.
  !! @param u3 The z component of the vector field.
  !! @param work1 A temporary array for computations.
  !! @param work2 A temporary array for computations.
  !! @param coef The SEM coefficients.
  subroutine curl(w1, w2, w3, u1, u2, u3, work1, work2, coef)
    type(field_t), intent(inout) :: w1
    type(field_t), intent(inout) :: w2
    type(field_t), intent(inout) :: w3
    type(field_t), intent(inout) :: u1
    type(field_t), intent(inout) :: u2
    type(field_t), intent(inout) :: u3
    type(field_t), intent(inout) :: work1
    type(field_t), intent(inout) :: work2
    type(coef_t), intent(in)  :: coef
 
    if (neko_bcknd_sx .eq. 1) then
       call opr_sx_curl(w1, w2, w3, u1, u2, u3, work1, work2, coef)
    else if (neko_bcknd_xsmm .eq. 1) then
       call opr_xsmm_curl(w1, w2, w3, u1, u2, u3, work1, work2, coef)
    else if (neko_bcknd_device .eq. 1) then
       call opr_device_curl(w1, w2, w3, u1, u2, u3, work1, work2, coef)
    else
       call opr_cpu_curl(w1, w2, w3, u1, u2, u3, work1, work2, coef)
    end if
 
  end subroutine curl
 
  !! Compute the CFL number
  !! @param dt The timestep.
  !! @param u The x component of velocity.
  !! @param v The y component of velocity.
  !! @param w The z component of velocity.
  !! @param Xh The SEM function space.
  !! @param coef The SEM coefficients.
  !! @param nelv The total number of elements.
  !! @param gdim Number of geometric dimensions.
  function cfl(dt, u, v, w, Xh, coef, nelv, gdim)
    type(space_t), intent(in) :: xh
    type(coef_t), intent(in) :: coef
    integer, intent(in) :: nelv, gdim
    real(kind=rp), intent(in) :: dt
    real(kind=rp), dimension(Xh%lx,Xh%ly,Xh%lz,nelv), intent(in) ::  u, v, w
    real(kind=rp) :: cfl
    integer :: ierr
 
    if (neko_bcknd_sx .eq. 1) then
       cfl = opr_sx_cfl(dt, u, v, w, xh, coef, nelv, gdim)
    else if (neko_bcknd_device .eq. 1) then
       cfl = opr_device_cfl(dt, u, v, w, xh, coef, nelv, gdim)
    else
       cfl = opr_cpu_cfl(dt, u, v, w, xh, coef, nelv, gdim)
    end if
 
    if (.not. neko_device_mpi) then
       call mpi_allreduce(mpi_in_place, cfl, 1, &
            mpi_real_precision, mpi_max, neko_comm, ierr)
    end if
 
  end function cfl
 
  subroutine strain_rate(s11, s22, s33, s12, s13, s23, u, v, w, coef)
    type(field_t), intent(in) :: u, v, w
    type(coef_t), intent(in) :: coef
    real(kind=rp), intent(inout) :: s11(u%Xh%lx, u%Xh%ly, u%Xh%lz, u%msh%nelv)
    real(kind=rp), intent(inout) :: s22(u%Xh%lx, u%Xh%ly, u%Xh%lz, u%msh%nelv)
    real(kind=rp), intent(inout) :: s33(u%Xh%lx, u%Xh%ly, u%Xh%lz, u%msh%nelv)
    real(kind=rp), intent(inout) :: s12(u%Xh%lx, u%Xh%ly, u%Xh%lz, u%msh%nelv)
    real(kind=rp), intent(inout) :: s23(u%Xh%lx, u%Xh%ly, u%Xh%lz, u%msh%nelv)
    real(kind=rp), intent(inout) :: s13(u%Xh%lx, u%Xh%ly, u%Xh%lz, u%msh%nelv)
 
    type(c_ptr) :: s11_d, s22_d, s33_d, s12_d, s23_d, s13_d
 
    integer :: nelv, lxyz
 
    if (neko_bcknd_device .eq. 1) then
       s11_d = device_get_ptr(s11)
       s22_d = device_get_ptr(s22)
       s33_d = device_get_ptr(s33)
       s12_d = device_get_ptr(s12)
       s23_d = device_get_ptr(s23)
       s13_d = device_get_ptr(s13)
    endif
 
    nelv = u%msh%nelv
    lxyz = u%Xh%lxyz
 
    ! we use s11 as a work array here
    call dudxyz (s12, u%x, coef%drdy, coef%dsdy, coef%dtdy, coef)
    call dudxyz (s11, v%x, coef%drdx, coef%dsdx, coef%dtdx, coef)
    if (neko_bcknd_device .eq. 1) then
       call device_add2(s12_d, s11_d, nelv*lxyz)
    else
       call add2(s12, s11, nelv*lxyz)
    endif
 
    call dudxyz (s13, u%x, coef%drdz, coef%dsdz, coef%dtdz, coef)
    call dudxyz (s11, w%x, coef%drdx, coef%dsdx, coef%dtdx, coef)
    if (neko_bcknd_device .eq. 1) then
       call device_add2(s13_d, s11_d, nelv*lxyz)
    else
       call add2(s13, s11, nelv*lxyz)
    endif
 
    call dudxyz (s23, v%x, coef%drdz, coef%dsdz, coef%dtdz, coef)
    call dudxyz (s11, w%x, coef%drdy, coef%dsdy, coef%dtdy, coef)
    if (neko_bcknd_device .eq. 1) then
       call device_add2(s23_d, s11_d, nelv*lxyz)
    else
       call add2(s23, s11, nelv*lxyz)
    endif
 
    call dudxyz (s11, u%x, coef%drdx, coef%dsdx, coef%dtdx, coef)
    call dudxyz (s22, v%x, coef%drdy, coef%dsdy, coef%dtdy, coef)
    call dudxyz (s33, w%x, coef%drdz, coef%dsdz, coef%dtdz, coef)
 
    if (neko_bcknd_device .eq. 1) then
       call device_cmult(s11_d, 2.0_rp, nelv*lxyz)
       call device_cmult(s22_d, 2.0_rp, nelv*lxyz)
       call device_cmult(s33_d, 2.0_rp, nelv*lxyz)
    else
       call cmult(s11, 2.0_rp, nelv*lxyz)
       call cmult(s22, 2.0_rp, nelv*lxyz)
       call cmult(s33, 2.0_rp, nelv*lxyz)
    endif
 
  end subroutine strain_rate
 
  subroutine lambda2op(lambda2, u, v, w, coef)
    type(coef_t), intent(in) :: coef
    type(field_t), intent(inout) :: lambda2
    type(field_t), intent(in) :: u, v, w
 
    if (neko_bcknd_sx .eq. 1) then
       call opr_sx_lambda2(lambda2, u, v, w, coef)
    else if (neko_bcknd_device .eq. 1) then
       call opr_device_lambda2(lambda2, u, v, w, coef)
    else
       call opr_cpu_lambda2(lambda2, u, v, w, coef)
    end if
 
  end subroutine lambda2op
 
end module operators