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, i8
  use opr_cpu, only : opr_cpu_cfl, opr_cpu_curl, opr_cpu_opgrad, &
       opr_cpu_conv1, opr_cpu_convect_scalar, opr_cpu_cdtp, &
       opr_cpu_dudxyz, opr_cpu_lambda2, opr_cpu_set_convect_rst, &
       opr_cpu_rotate_cyc_r1, opr_cpu_rotate_cyc_r4
  use opr_sx, only : opr_sx_cfl, opr_sx_curl, opr_sx_opgrad, &
       opr_sx_conv1, opr_sx_convect_scalar, opr_sx_cdtp, &
       opr_sx_dudxyz, opr_sx_lambda2, opr_sx_set_convect_rst
  use opr_xsmm, only : opr_xsmm_cdtp, opr_xsmm_conv1, opr_xsmm_curl, &
       opr_xsmm_dudxyz, opr_xsmm_opgrad, &
       opr_xsmm_convect_scalar, opr_xsmm_set_convect_rst
  use opr_device, only : opr_device_cdtp, opr_device_cfl, opr_device_curl, &
       opr_device_conv1, opr_device_convect_scalar, opr_device_dudxyz, &
       opr_device_lambda2, opr_device_opgrad, opr_device_set_convect_rst, &
       opr_device_rotate_cyc
  use space, only : space_t
  use coefs, only : coef_t
  use field, only : field_t
  use field_list, only : field_list_t
  use field_math, only : field_rzero
  use interpolation, only : interpolator_t
  use math, only : glsum, cmult, add2, add3s2, cadd, copy, col2, invcol2, &
       invcol3, rzero, add5s4
  use device, only : device_get_ptr, device_map, device_free
  use device_math, only : device_add2, device_cmult, device_copy, device_cadd, &
       device_glsum, device_add3s2, device_invcol2, device_invcol3, &
       device_col2, device_add5s4
  use scratch_registry, only : neko_scratch_registry
  use vector, only : vector_t
  use comm, only : neko_comm, mpi_real_precision
  use mpi_f08, only : mpi_allreduce, mpi_in_place, mpi_max, mpi_sum
  use, intrinsic :: iso_c_binding, only : c_ptr
  use logger, only : neko_log
  implicit none
  private
 
  public :: dudxyz, opgrad, ortho, cdtp, conv1, curl, cfl, cfl_compressible, &
       lambda2op, strain_rate, div, grad, set_convect_rst, runge_kutta, &
       rotate_cyc
 
  interface dudxyz
     module procedure dudxyz_r4
     module procedure opr_device_dudxyz
     module procedure dudxyz_f
  end interface dudxyz
 
  interface div
     module procedure div_r4
     module procedure div_d
  end interface div
 
  interface grad
     module procedure grad_r4
     module procedure grad_d
  end interface grad
 
  interface cfl
     module procedure cfl_r4
     module procedure cfl_d
     module procedure cfl_f
  end interface cfl
 
  interface cfl_compressible
     module procedure cfl_compressible_r4
     module procedure cfl_compressible_d
     module procedure cfl_compressible_f
  end interface cfl_compressible
 
  interface rotate_cyc
     module procedure rotate_cyc_r1
     module procedure rotate_cyc_r4
     module procedure rotate_cyc_d
     module procedure rotate_cyc_f
  end interface rotate_cyc
 
  interface strain_rate
     module procedure strain_rate_r4
     module procedure strain_rate_d
     module procedure strain_rate_f
  end interface strain_rate
contains
 
  subroutine dudxyz_r4(du, u, dr, ds, dt, coef)
    real(kind=rp), contiguous, dimension(:,:,:,:), intent(inout) :: du
    real(kind=rp), contiguous, dimension(:,:,:,:), intent(in) :: u, dr, ds, dt
    type(coef_t), intent(in) :: coef
    type(c_ptr) :: du_d, u_d, dr_d, ds_d, dt_d
 
    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 neko_log%deprecated('Operator: dudxyz, implicit device', &
            '2.0.0', 'Please call opr_device_dudxyz instead.')
 
       du_d = device_get_ptr(du)
       u_d = device_get_ptr(u)
       dr_d = device_get_ptr(dr)
       ds_d = device_get_ptr(ds)
       dt_d = device_get_ptr(dt)
 
       call opr_device_dudxyz(du_d, u_d, dr_d, ds_d, dt_d, coef)
    else
       call opr_cpu_dudxyz(du, u, dr, ds, dt, coef)
    end if
 
  end subroutine dudxyz_r4
 
  subroutine dudxyz_f(du, u, dr, ds, dt, coef)
    type(field_t), intent(inout) :: du
    type(field_t), intent(in) :: u, dr, ds, dt
    type(coef_t), intent(in) :: coef
 
    if (neko_bcknd_sx .eq. 1) then
       call opr_sx_dudxyz(du%x, u%x, dr%x, ds%x, dt%x, coef)
    else if (neko_bcknd_xsmm .eq. 1) then
       call opr_xsmm_dudxyz(du%x, u%x, dr%x, ds%x, dt%x, coef)
    else if (neko_bcknd_device .eq. 1) then
       call opr_device_dudxyz(du%x_d, u%x_d, dr%x_d, ds%x_d, dt%x_d, coef)
    else
       call opr_cpu_dudxyz(du%x, u%x, dr%x, ds%x, dt%x, coef)
    end if
 
  end subroutine dudxyz_f
 
  subroutine div_r4(res, ux, uy, uz, coef)
    real(kind=rp), contiguous, dimension(:,:,:,:), intent(inout) :: res
    real(kind=rp), contiguous, dimension(:,:,:,:), intent(in) :: ux, uy, uz
    type(coef_t), intent(in), target :: coef
    type(field_t), pointer :: work
    integer :: ind
    type(c_ptr) :: res_d
 
    if (neko_bcknd_device .eq. 1) then
       call neko_log%deprecated('Operator: div, implicit device', &
            '2.0.0', 'Please call div_d instead.')
       res_d = device_get_ptr(res)
    end if
 
    call neko_scratch_registry%request_field(work, ind, .false.)
 
    ! Get dux / dx
    call dudxyz(res, ux, coef%drdx, coef%dsdx, coef%dtdx, coef)
 
    ! Get duy / dy
    call dudxyz(work%x, uy, coef%drdy, coef%dsdy, coef%dtdy, coef)
    if (neko_bcknd_device .eq. 1) then
       call device_add2(res_d, work%x_d, work%size())
    else
       call add2(res, work%x, work%size())
    end if
 
    ! Get dux / dz
    call dudxyz(work%x, uz, coef%drdz, coef%dsdz, coef%dtdz, coef)
    if (neko_bcknd_device .eq. 1) then
       call device_add2(res_d, work%x_d, work%size())
    else
       call add2(res, work%x, work%size())
    end if
 
    call neko_scratch_registry%relinquish_field(ind)
 
  end subroutine div_r4
 
  subroutine div_d(res_d, ux_d, uy_d, uz_d, coef)
    type(c_ptr), intent(inout) :: res_d
    type(c_ptr), intent(in) :: ux_d, uy_d, uz_d
    type(coef_t), intent(in), target :: coef
    type(field_t), pointer :: work
    integer :: ind
 
    call neko_scratch_registry%request_field(work, ind, .false.)
 
    ! Get dux / dx
    call dudxyz(res_d, ux_d, coef%drdx_d, coef%dsdx_d, coef%dtdx_d, coef)
 
    ! Get duy / dy
    call dudxyz(work%x_d, uy_d, coef%drdy_d, coef%dsdy_d, coef%dtdy_d, coef)
    call device_add2(res_d, work%x_d, work%size())
 
    ! Get duz / dz
    call dudxyz(work%x_d, uz_d, coef%drdz_d, coef%dsdz_d, coef%dtdz_d, coef)
    call device_add2(res_d, work%x_d, work%size())
 
    call neko_scratch_registry%relinquish_field(ind)
 
  end subroutine div_d
 
  subroutine grad_r4(ux, uy, uz, u, coef)
    type(coef_t), intent(in) :: coef
    real(kind=rp), contiguous, dimension(:,:,:,:), intent(inout) :: ux
    real(kind=rp), contiguous, dimension(:,:,:,:), intent(inout) :: uy
    real(kind=rp), contiguous, dimension(:,:,:,:), intent(inout) :: uz
    real(kind=rp), contiguous, dimension(:,:,:,:), intent(in) :: u
 
    call dudxyz(ux, u, coef%drdx, coef%dsdx, coef%dtdx, coef)
    call dudxyz(uy, u, coef%drdy, coef%dsdy, coef%dtdy, coef)
    call dudxyz(uz, u, coef%drdz, coef%dsdz, coef%dtdz, coef)
 
  end subroutine grad_r4
 
  subroutine grad_d(ux_d, uy_d, uz_d, u_d, coef)
    type(coef_t), intent(in) :: coef
    type(c_ptr), intent(inout) :: ux_d, uy_d, uz_d
    type(c_ptr), intent(in) :: u_d
 
    call dudxyz(ux_d, u_d, coef%drdx_d, coef%dsdx_d, coef%dtdx_d, coef)
    call dudxyz(uy_d, u_d, coef%drdy_d, coef%dsdy_d, coef%dtdy_d, coef)
    call dudxyz(uz_d, u_d, coef%drdz_d, coef%dsdz_d, coef%dtdz_d, coef)
 
  end subroutine grad_d
 
  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
    type(c_ptr) :: ux_d, uy_d, uz_d, u_d
 
    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
       ux_d = device_get_ptr(ux)
       uy_d = device_get_ptr(uy)
       uz_d = device_get_ptr(uz)
       u_d = device_get_ptr(u)
       call opr_device_opgrad(ux_d, uy_d, uz_d, u_d, coef)
    else
       call opr_cpu_opgrad(ux, uy, uz, u, coef, eblk_start, eblk_end)
    end if
 
  end subroutine opgrad
 
  subroutine ortho(x, glb_n_points, n)
    integer, intent(in) :: n
    integer(kind=i8), intent(in) :: glb_n_points
    real(kind=rp), dimension(n), intent(inout) :: x
    real(kind=rp) :: c
    type(c_ptr) :: x_d
 
    if (neko_bcknd_device .eq. 1) then
       call neko_log%deprecated('Operator: ortho, implicit device', &
            '2.0.0', 'Please call device_ortho instead.')
 
       x_d = device_get_ptr(x)
       c = device_glsum(x_d, n) / glb_n_points
       call device_cadd(x_d, -c, n)
    else
       c = glsum(x, n) / glb_n_points
       call cadd(x, -c, n)
    end if
 
  end subroutine ortho
 
  subroutine cdtp(dtx, x, dr, ds, dt, coef, es, ee)
    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
    integer, optional :: es, ee
    integer :: eblk_start, eblk_end
    type(c_ptr) :: dtx_d, x_d, dr_d, ds_d, dt_d
 
    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_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
       dtx_d = device_get_ptr(dtx)
       x_d = device_get_ptr(x)
       dr_d = device_get_ptr(dr)
       ds_d = device_get_ptr(ds)
       dt_d = device_get_ptr(dt)
       call opr_device_cdtp(dtx_d, x_d, dr_d, ds_d, dt_d, coef)
    else
       call opr_cpu_cdtp(dtx, x, dr, ds, dt, coef, eblk_start, eblk_end)
    end if
 
  end subroutine cdtp
 
  subroutine conv1(du, u, vx, vy, vz, Xh, coef, es, ee)
    type(space_t), intent(in) :: xh
    type(coef_t), intent(in) :: coef
    real(kind=rp), intent(inout) :: du(xh%lxyz, coef%msh%nelv)
    real(kind=rp), intent(in) :: u(xh%lx, xh%ly, xh%lz, coef%msh%nelv)
    real(kind=rp), intent(in) :: vx(xh%lx, xh%ly, xh%lz, coef%msh%nelv)
    real(kind=rp), intent(in) :: vy(xh%lx, xh%ly, xh%lz, coef%msh%nelv)
    real(kind=rp), intent(in) :: vz(xh%lx, xh%ly, xh%lz, coef%msh%nelv)
    integer, optional :: es, ee
    integer :: eblk_end, eblk_start
    type(c_ptr) :: du_d, u_d, vx_d, vy_d, vz_d
 
    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)
      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
         du_d = device_get_ptr(du)
         u_d = device_get_ptr(u)
         vx_d = device_get_ptr(vx)
         vy_d = device_get_ptr(vy)
         vz_d = device_get_ptr(vz)
         call opr_device_conv1(du_d, u_d, vx_d, vy_d, vz_d, 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
 
  subroutine convect_scalar(du, u, cr, cs, ct, Xh_GLL, Xh_GL, coef_GLL, &
       coef_GL, GLL_to_GL)
    type(space_t), intent(in) :: Xh_GL
    type(space_t), intent(in) :: Xh_GLL
    type(coef_t), intent(in) :: coef_GLL
    type(coef_t), intent(in) :: coef_GL
    type(interpolator_t), intent(inout) :: GLL_to_GL
    real(kind=rp), intent(inout) :: &
         du(xh_gll%lx, xh_gll%ly, xh_gll%lz, coef_gl%msh%nelv)
    real(kind=rp), intent(inout) :: &
         u(xh_gl%lx, xh_gl%lx, xh_gl%lx, coef_gl%msh%nelv)
    type(field_t), intent(inout) :: cr, cs, ct
    type(c_ptr) :: u_d
 
    if (neko_bcknd_sx .eq. 1) then
       call opr_sx_convect_scalar(du, u, cr%x, cs%x, ct%x, &
            xh_gll, xh_gl, coef_gll, coef_gl, gll_to_gl)
    else if (neko_bcknd_xsmm .eq. 1) then
       call opr_xsmm_convect_scalar(du, u, cr%x, cs%x, ct%x, &
            xh_gll, xh_gl, coef_gll, coef_gl, gll_to_gl)
    else if (neko_bcknd_device .eq. 1) then
       u_d = device_get_ptr(u)
       call opr_device_convect_scalar(du, u_d, cr%x_d, cs%x_d, ct%x_d, &
            xh_gll, xh_gl, coef_gll, coef_gl, gll_to_gl)
    else
       call opr_cpu_convect_scalar(du, u, cr%x, cs%x, ct%x, &
            xh_gll, xh_gl, coef_gll, coef_gl, gll_to_gl)
    end if
 
  end subroutine convect_scalar
 
  !! 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, event)
    type(field_t), intent(inout) :: w1
    type(field_t), intent(inout) :: w2
    type(field_t), intent(inout) :: w3
    type(field_t), intent(in) :: u1
    type(field_t), intent(in) :: u2
    type(field_t), intent(in) :: u3
    type(field_t), intent(inout) :: work1
    type(field_t), intent(inout) :: work2
    type(coef_t), intent(in) :: coef
    type(c_ptr), optional, intent(inout) :: event
 
    if (neko_bcknd_sx .eq. 1) then
       call opr_sx_curl(w1%x, w2%x, w3%x, u1%x, u2%x, u3%x, &
            work1%x, work2%x, coef)
    else if (neko_bcknd_xsmm .eq. 1) then
       call opr_xsmm_curl(w1%x, w2%x, w3%x, u1%x, u2%x, u3%x, &
            work1%x, work2%x, coef)
    else if (neko_bcknd_device .eq. 1) then
       call opr_device_curl(w1, w2, w3, u1, u2, u3, &
            work1, work2, coef, event)
    else
       call opr_cpu_curl(w1%x, w2%x, w3%x, u1%x, u2%x, u3%x, &
            work1%x, work2%x, 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_r4(dt, u, v, w, Xh, coef, nelv, gdim)
    real(kind=rp), intent(in) :: dt
    real(kind=rp), contiguous, dimension(:,:,:,:), intent(in) :: u, v, w
    type(space_t), intent(in) :: xh
    type(coef_t), intent(in) :: coef
    integer, intent(in) :: nelv, gdim
    real(kind=rp) :: cfl_r4
    integer :: ierr
    type(c_ptr) :: u_d, v_d, w_d
 
    if (neko_bcknd_sx .eq. 1) then
       cfl_r4 = opr_sx_cfl(dt, u, v, w, xh, coef, nelv)
    else if (neko_bcknd_device .eq. 1) then
       call neko_log%deprecated('Operator: cfl_r4, implicit device', &
            '2.0.0', 'Please call cfl_d instead.')
 
       u_d = device_get_ptr(u)
       v_d = device_get_ptr(v)
       w_d = device_get_ptr(w)
 
       cfl_r4 = opr_device_cfl(dt, u_d, v_d, w_d, xh, coef, nelv, gdim)
    else
       cfl_r4 = 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_r4, 1, &
            mpi_real_precision, mpi_max, neko_comm, ierr)
    end if
 
  end function cfl_r4
 
  function cfl_d(dt, u_d, v_d, w_d, Xh, coef, nelv, gdim)
    real(kind=rp), intent(in) :: dt
    type(c_ptr), intent(in) :: u_d, v_d, w_d
    type(space_t), intent(in) :: xh
    type(coef_t), intent(in) :: coef
    integer, intent(in) :: nelv, gdim
    real(kind=rp) :: cfl_d
    integer :: ierr
 
    cfl_d = opr_device_cfl(dt, u_d, v_d, w_d, xh, coef, nelv, gdim)
 
    if (.not. neko_device_mpi) then
       call mpi_allreduce(mpi_in_place, cfl_d, 1, &
            mpi_real_precision, mpi_max, neko_comm, ierr)
    end if
 
  end function cfl_d
 
  function cfl_f(dt, u, v, w, Xh, coef, nelv, gdim)
    real(kind=rp), intent(in) :: dt
    type(field_t), intent(in) :: u, v, w
    type(space_t), intent(in) :: xh
    type(coef_t), intent(in) :: coef
    integer, intent(in) :: nelv, gdim
    real(kind=rp) :: cfl_f
    integer :: ierr
 
    if (neko_bcknd_sx .eq. 1) then
       cfl_f = opr_sx_cfl(dt, u%x, v%x, w%x, xh, coef, nelv)
    else if (neko_bcknd_device .eq. 1) then
       cfl_f = opr_device_cfl(dt, u%x_d, v%x_d, w%x_d, xh, coef, nelv, gdim)
    else
       cfl_f = opr_cpu_cfl(dt, u%x, v%x, w%x, xh, coef, nelv, gdim)
    end if
 
    if (.not. neko_device_mpi) then
       call mpi_allreduce(mpi_in_place, cfl_f, 1, &
            mpi_real_precision, mpi_max, neko_comm, ierr)
    end if
 
  end function cfl_f
 
  !! Compute the CFL number for compressible flows
  !! @param dt The timestep.
  !! @param max_wave_speed The precomputed maximum wave speed field.
  !! @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_compressible_r4(dt, max_wave_speed, Xh, coef, nelv, gdim)
    real(kind=rp), intent(in) :: dt
    real(kind=rp), contiguous, dimension(:,:,:,:), intent(in) :: max_wave_speed
    type(space_t), intent(in) :: xh
    type(coef_t), intent(in) :: coef
    integer, intent(in) :: nelv, gdim
    real(kind=rp) :: cfl_compressible_r4
 
    cfl_compressible_r4 = cfl(dt, max_wave_speed, max_wave_speed, &
         max_wave_speed, xh, coef, nelv, gdim)
 
  end function cfl_compressible_r4
 
  !! Compute the CFL number for compressible flows
  !! @param dt The timestep.
  !! @param max_wave_speed The precomputed maximum wave speed field.
  !! @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_compressible_d(dt, max_wave_speed, Xh, coef, nelv, gdim)
    real(kind=rp), intent(in) :: dt
    type(c_ptr), intent(in) :: max_wave_speed
    type(space_t), intent(in) :: xh
    type(coef_t), intent(in) :: coef
    integer, intent(in) :: nelv, gdim
    real(kind=rp) :: cfl_compressible_d
 
    cfl_compressible_d = cfl(dt, max_wave_speed, max_wave_speed, &
         max_wave_speed, xh, coef, nelv, gdim)
 
  end function cfl_compressible_d
 
  !! Compute the CFL number for compressible flows
  !! @param dt The timestep.
  !! @param max_wave_speed The precomputed maximum wave speed field.
  !! @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_compressible_f(dt, max_wave_speed, Xh, coef, nelv, gdim)
    real(kind=rp), intent(in) :: dt
    type(field_t), intent(in) :: max_wave_speed
    type(space_t), intent(in) :: xh
    type(coef_t), intent(in) :: coef
    integer, intent(in) :: nelv, gdim
    real(kind=rp) :: cfl_compressible_f
 
    cfl_compressible_f = cfl(dt, max_wave_speed, max_wave_speed, &
         max_wave_speed, xh, coef, nelv, gdim)
 
  end function cfl_compressible_f
 
  subroutine strain_rate_r4(s11, s22, s33, s12, s13, s23, u, v, w, coef)
    real(kind=rp), contiguous, intent(inout) :: s11(:,:,:,:)
    real(kind=rp), contiguous, intent(inout) :: s22(:,:,:,:)
    real(kind=rp), contiguous, intent(inout) :: s33(:,:,:,:)
    real(kind=rp), contiguous, intent(inout) :: s12(:,:,:,:)
    real(kind=rp), contiguous, intent(inout) :: s13(:,:,:,:)
    real(kind=rp), contiguous, intent(inout) :: s23(:,:,:,:)
    real(kind=rp), contiguous, intent(in) :: u(:,:,:,:), v(:,:,:,:), w(:,:,:,:)
    type(coef_t), intent(in) :: coef
 
    type(c_ptr) :: s11_d, s22_d, s33_d, s12_d, s23_d, s13_d, u_d, v_d, w_d
 
    integer :: nelv, lxyz
 
    nelv = coef%msh%nelv
    lxyz = coef%Xh%lxyz
 
    if (neko_bcknd_device .eq. 1) then
       call neko_log%deprecated('Operator: strain_rate_r4, implicit device', &
            '2.0.0', 'Please call strain_rate_d instead.')
       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)
       u_d = device_get_ptr(u)
       v_d = device_get_ptr(v)
       w_d = device_get_ptr(w)
 
       call dudxyz(s12_d, u_d, coef%drdy_d, coef%dsdy_d, coef%dtdy_d, coef)
       call dudxyz(s11_d, v_d, coef%drdx_d, coef%dsdx_d, coef%dtdx_d, coef)
       call device_add2(s12_d, s11_d, nelv*lxyz)
 
       call dudxyz(s13_d, u_d, coef%drdz_d, coef%dsdz_d, coef%dtdz_d, coef)
       call dudxyz(s11_d, w_d, coef%drdx_d, coef%dsdx_d, coef%dtdx_d, coef)
       call device_add2(s13_d, s11_d, nelv*lxyz)
 
       call dudxyz(s23_d, v_d, coef%drdz_d, coef%dsdz_d, coef%dtdz_d, coef)
       call dudxyz(s11_d, w_d, coef%drdy_d, coef%dsdy_d, coef%dtdy_d, coef)
       call device_add2(s23_d, s11_d, nelv*lxyz)
 
       call dudxyz(s11_d, u_d, coef%drdx_d, coef%dsdx_d, coef%dtdx_d, coef)
       call dudxyz(s22_d, v_d, coef%drdy_d, coef%dsdy_d, coef%dtdy_d, coef)
       call dudxyz(s33_d, w_d, coef%drdz_d, coef%dsdz_d, coef%dtdz_d, coef)
       call device_cmult(s12_d, 0.5_rp, nelv*lxyz)
       call device_cmult(s13_d, 0.5_rp, nelv*lxyz)
       call device_cmult(s23_d, 0.5_rp, nelv*lxyz)
    else
       call dudxyz(s12, u, coef%drdy, coef%dsdy, coef%dtdy, coef)
       call dudxyz(s11, v, coef%drdx, coef%dsdx, coef%dtdx, coef)
       call add2(s12, s11, nelv*lxyz)
 
       call dudxyz(s13, u, coef%drdz, coef%dsdz, coef%dtdz, coef)
       call dudxyz(s11, w, coef%drdx, coef%dsdx, coef%dtdx, coef)
       call add2(s13, s11, nelv*lxyz)
 
       call dudxyz(s23, v, coef%drdz, coef%dsdz, coef%dtdz, coef)
       call dudxyz(s11, w, coef%drdy, coef%dsdy, coef%dtdy, coef)
       call add2(s23, s11, nelv*lxyz)
 
       call dudxyz(s11, u, coef%drdx, coef%dsdx, coef%dtdx, coef)
       call dudxyz(s22, v, coef%drdy, coef%dsdy, coef%dtdy, coef)
       call dudxyz(s33, w, coef%drdz, coef%dsdz, coef%dtdz, coef)
       call cmult(s12, 0.5_rp, nelv*lxyz)
       call cmult(s13, 0.5_rp, nelv*lxyz)
       call cmult(s23, 0.5_rp, nelv*lxyz)
    end if
 
  end subroutine strain_rate_r4
 
  subroutine strain_rate_d(s11_d, s22_d, s33_d, s12_d, s13_d, s23_d, &
       u_d, v_d, w_d, coef)
    type(c_ptr), intent(inout) :: s11_d, s22_d, s33_d, s12_d, s13_d, s23_d
    type(c_ptr), intent(in) :: u_d, v_d, w_d
    type(coef_t), intent(in) :: coef
 
    integer :: nelv, lxyz
 
    nelv = coef%msh%nelv
    lxyz = coef%Xh%lxyz
 
    call dudxyz(s12_d, u_d, coef%drdy_d, coef%dsdy_d, coef%dtdy_d, coef)
    call dudxyz(s11_d, v_d, coef%drdx_d, coef%dsdx_d, coef%dtdx_d, coef)
    call device_add2(s12_d, s11_d, nelv*lxyz)
 
    call dudxyz(s13_d, u_d, coef%drdz_d, coef%dsdz_d, coef%dtdz_d, coef)
    call dudxyz(s11_d, w_d, coef%drdx_d, coef%dsdx_d, coef%dtdx_d, coef)
    call device_add2(s13_d, s11_d, nelv*lxyz)
 
    call dudxyz(s23_d, v_d, coef%drdz_d, coef%dsdz_d, coef%dtdz_d, coef)
    call dudxyz(s11_d, w_d, coef%drdy_d, coef%dsdy_d, coef%dtdy_d, coef)
    call device_add2(s23_d, s11_d, nelv*lxyz)
 
    call dudxyz(s11_d, u_d, coef%drdx_d, coef%dsdx_d, coef%dtdx_d, coef)
    call dudxyz(s22_d, v_d, coef%drdy_d, coef%dsdy_d, coef%dtdy_d, coef)
    call dudxyz(s33_d, w_d, coef%drdz_d, coef%dsdz_d, coef%dtdz_d, coef)
    call device_cmult(s12_d, 0.5_rp, nelv*lxyz)
    call device_cmult(s13_d, 0.5_rp, nelv*lxyz)
    call device_cmult(s23_d, 0.5_rp, nelv*lxyz)
 
  end subroutine strain_rate_d
 
  subroutine strain_rate_f(s11, s22, s33, s12, s13, s23, u, v, w, coef)
    type(field_t), intent(inout) :: s11, s22, s33, s12, s13, s23
    type(field_t), intent(in) :: u, v, w
    type(coef_t), intent(in) :: coef
 
    integer :: n
    n = coef%Xh%lxyz * coef%msh%nelv
 
    if (neko_bcknd_device .eq. 1) then
       call dudxyz(s12%x_d, u%x_d, coef%drdy_d, coef%dsdy_d, coef%dtdy_d, coef)
       call dudxyz(s11%x_d, v%x_d, coef%drdx_d, coef%dsdx_d, coef%dtdx_d, coef)
       call device_add2(s12%x_d, s11%x_d, n)
 
       call dudxyz(s13%x_d, u%x_d, coef%drdz_d, coef%dsdz_d, coef%dtdz_d, coef)
       call dudxyz(s11%x_d, w%x_d, coef%drdx_d, coef%dsdx_d, coef%dtdx_d, coef)
       call device_add2(s13%x_d, s11%x_d, n)
 
       call dudxyz(s23%x_d, v%x_d, coef%drdz_d, coef%dsdz_d, coef%dtdz_d, coef)
       call dudxyz(s11%x_d, w%x_d, coef%drdy_d, coef%dsdy_d, coef%dtdy_d, coef)
       call device_add2(s23%x_d, s11%x_d, n)
 
       call dudxyz(s11%x_d, u%x_d, coef%drdx_d, coef%dsdx_d, coef%dtdx_d, coef)
       call dudxyz(s22%x_d, v%x_d, coef%drdy_d, coef%dsdy_d, coef%dtdy_d, coef)
       call dudxyz(s33%x_d, w%x_d, coef%drdz_d, coef%dsdz_d, coef%dtdz_d, coef)
       call device_cmult(s12%x_d, 0.5_rp, n)
       call device_cmult(s13%x_d, 0.5_rp, n)
       call device_cmult(s23%x_d, 0.5_rp, n)
    else
       call dudxyz(s12%x, u%x, coef%drdy, coef%dsdy, coef%dtdy, coef)
       call dudxyz(s11%x, v%x, coef%drdx, coef%dsdx, coef%dtdx, coef)
       call add2(s12%x, s11%x, n)
 
       call dudxyz(s13%x, u%x, coef%drdz, coef%dsdz, coef%dtdz, coef)
       call dudxyz(s11%x, w%x, coef%drdx, coef%dsdx, coef%dtdx, coef)
       call add2(s13%x, s11%x, n)
 
       call dudxyz(s23%x, v%x, coef%drdz, coef%dsdz, coef%dtdz, coef)
       call dudxyz(s11%x, w%x, coef%drdy, coef%dsdy, coef%dtdy, coef)
       call add2(s23%x, s11%x, n)
 
       call dudxyz(s11%x, u%x, coef%drdx, coef%dsdx, coef%dtdx, coef)
       call dudxyz(s22%x, v%x, coef%drdy, coef%dsdy, coef%dtdy, coef)
       call dudxyz(s33%x, w%x, coef%drdz, coef%dsdz, coef%dtdz, coef)
       call cmult(s12%x, 0.5_rp, n)
       call cmult(s13%x, 0.5_rp, n)
       call cmult(s23%x, 0.5_rp, n)
    end if
 
  end subroutine strain_rate_f
 
  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%x, u%x, v%x, w%x, coef)
    else if (neko_bcknd_device .eq. 1) then
       call opr_device_lambda2(lambda2%x_d, u%x_d, v%x_d, w%x_d, coef)
    else
       call opr_cpu_lambda2(lambda2%x, u%x, v%x, w%x, coef)
    end if
 
  end subroutine lambda2op
 
  subroutine set_convect_rst(cr, cs, ct, cx, cy, cz, Xh, coef)
    type(space_t), intent(inout) :: xh
    type(coef_t), intent(inout) :: coef
    type(field_t), intent(inout) :: cr, cs, ct
    real(kind=rp), dimension(Xh%lxyz, coef%msh%nelv), &
         intent(in) :: cx, cy, cz
    type(c_ptr) :: cx_d, cy_d, cz_d
 
    if (neko_bcknd_sx .eq. 1) then
       call opr_sx_set_convect_rst(cr%x, cs%x, ct%x, cx, cy, cz, xh, coef)
    else if (neko_bcknd_xsmm .eq. 1) then
       call opr_xsmm_set_convect_rst(cr%x, cs%x, ct%x, cx, cy, cz, xh, coef)
    else if (neko_bcknd_device .eq. 1) then
       cx_d = device_get_ptr(cx)
       cy_d = device_get_ptr(cy)
       cz_d = device_get_ptr(cz)
       call opr_device_set_convect_rst(cr%x_d, cs%x_d, ct%x_d, &
            cx_d, cy_d, cz_d, xh, coef)
    else
       call opr_cpu_set_convect_rst(cr%x, cs%x, ct%x, cx, cy, cz, xh, coef)
    end if
 
  end subroutine set_convect_rst
 
  subroutine runge_kutta(phi, conv_k1, conv_k23, conv_k4, Xh_GLL, Xh_GL, &
       coef, coef_GL, GLL_to_GL, tau, dtau, n, nel, n_GL)
    type(space_t), intent(in) :: xh_gll
    type(space_t), intent(inout) :: xh_gl
    type(coef_t), intent(in) :: coef
    type(coef_t), intent(inout) :: coef_gl
    type(interpolator_t) :: gll_to_gl
    real(kind=rp), intent(inout) :: tau, dtau
    integer, intent(in) :: n, nel, n_gl
    type(field_t), intent(inout) :: phi
    type(field_list_t) :: conv_k1, conv_k23, conv_k4
    real(kind=rp) :: c1, c2, c3
    type(field_t), pointer :: u1, k1, k2, k3, k4
    type(vector_t), pointer :: u1_gl
    integer :: ind(6), i, e
 
    call neko_scratch_registry%request_field(u1, ind(1), .false.)
    call neko_scratch_registry%request_field(k1, ind(2), .false.)
    call neko_scratch_registry%request_field(k2, ind(3), .false.)
    call neko_scratch_registry%request_field(k3, ind(4), .false.)
    call neko_scratch_registry%request_field(k4, ind(5), .false.)
    call neko_scratch_registry%request_vector(u1_gl, ind(6), n_gl, .false.)
 
    c1 = 1.0_rp
    c2 = -dtau / 2.0_rp
    c3 = -dtau
 
    if (neko_bcknd_device .eq. 1) then
 
       ! Stage 1:
       call device_invcol3(u1%x_d, phi%x_d, coef%B_d, n)
       call gll_to_gl%map(u1_gl%x, u1%x, nel, xh_gl)
       call convect_scalar(k1%x, u1_gl%x, conv_k1%items(1)%ptr, &
            conv_k1%items(2)%ptr, conv_k1%items(3)%ptr, &
            xh_gll, xh_gl, coef, coef_gl, gll_to_gl)
       call device_col2(k1%x_d, coef%B_d, n)
 
       ! Stage 2:
       call device_add3s2(u1%x_d, phi%x_d, k1%x_d, c1, c2, n)
       call device_invcol2(u1%x_d, coef%B_d, n)
       call gll_to_gl%map(u1_gl%x, u1%x, nel, xh_gl)
       call convect_scalar(k2%x, u1_gl%x, conv_k23%items(1)%ptr, &
            conv_k23%items(2)%ptr, conv_k23%items(3)%ptr, &
            xh_gll, xh_gl, coef, coef_gl, gll_to_gl)
       call device_col2(k2%x_d, coef%B_d, n)
 
       ! Stage 3:
       call device_add3s2(u1%x_d, phi%x_d, k2%x_d, c1, c2, n)
       call device_invcol2(u1%x_d, coef%B_d, n)
       call gll_to_gl%map(u1_gl%x, u1%x, nel, xh_gl)
       call convect_scalar(k3%x, u1_gl%x, conv_k23%items(1)%ptr, &
            conv_k23%items(2)%ptr, conv_k23%items(3)%ptr, &
            xh_gll, xh_gl, coef, coef_gl, gll_to_gl)
       call device_col2(k3%x_d, coef%B_d, n)
 
       ! Stage 4:
       call device_add3s2(u1%x_d, phi%x_d, k3%x_d, c1, c3, n)
       call device_invcol2(u1%x_d, coef%B_d, n)
       call gll_to_gl%map(u1_gl%x, u1%x, nel, xh_gl)
       call convect_scalar(k4%x, u1_gl%x, conv_k4%items(1)%ptr, &
            conv_k4%items(2)%ptr, conv_k4%items(3)%ptr, &
            xh_gll, xh_gl, coef, coef_gl, gll_to_gl)
       call device_col2(k4%x_d, coef%B_d, n)
 
       c1 = -dtau / 6.0_rp
       c2 = -dtau / 3.0_rp
 
       call device_add5s4(phi%x_d, k1%x_d, k2%x_d, k3%x_d, k4%x_d, &
            c1, c2, c2, c1, n)
 
    else
 
       ! Stage 1:
       call invcol3(u1%x, phi%x, coef%B, n)
       call gll_to_gl%map(u1_gl%x, u1%x, nel, xh_gl)
       call convect_scalar(k1%x, u1_gl%x, conv_k1%items(1)%ptr, &
            conv_k1%items(2)%ptr, conv_k1%items(3)%ptr, &
            xh_gll, xh_gl, coef, coef_gl, gll_to_gl)
       call col2(k1%x, coef%B, n)
 
       ! Stage 2:
       call add3s2(u1%x, phi%x, k1%x, c1, c2, n)
       call invcol2(u1%x, coef%B, n)
       call gll_to_gl%map(u1_gl%x, u1%x, nel, xh_gl)
       call convect_scalar(k2%x, u1_gl%x, conv_k23%items(1)%ptr, &
            conv_k23%items(2)%ptr, conv_k23%items(3)%ptr, &
            xh_gll, xh_gl, coef, coef_gl, gll_to_gl)
       call col2(k2%x, coef%B, n)
 
       ! Stage 3:
       call add3s2(u1%x, phi%x, k2%x, c1, c2, n)
       call invcol2(u1%x, coef%B, n)
       call gll_to_gl%map(u1_gl%x, u1%x, nel, xh_gl)
       call convect_scalar(k3%x, u1_gl%x, conv_k23%items(1)%ptr, &
            conv_k23%items(2)%ptr, conv_k23%items(3)%ptr, &
            xh_gll, xh_gl, coef, coef_gl, gll_to_gl)
       call col2(k3%x, coef%B, n)
 
       ! Stage 4:
       call add3s2(u1%x, phi%x, k3%x, c1, c3, n)
       call invcol2(u1%x, coef%B, n)
       call gll_to_gl%map(u1_gl%x, u1%x, nel, xh_gl)
       call convect_scalar(k4%x, u1_gl%x, conv_k4%items(1)%ptr, &
            conv_k4%items(2)%ptr, conv_k4%items(3)%ptr, &
            xh_gll, xh_gl, coef, coef_gl, gll_to_gl)
       call col2(k4%x, coef%B, n)
 
       c1 = -dtau / 6.0_rp
       c2 = -dtau / 3.0_rp
       call add5s4(phi%x, k1%x, k2%x, k3%x, k4%x, c1, c2, c2, c1, n)
    end if
 
    call neko_scratch_registry%relinquish(ind)
 
  end subroutine runge_kutta
 
  subroutine rotate_cyc_r1(vx, vy, vz, idir, coef)
    real(kind=rp), contiguous, dimension(:), intent(inout) :: vx, vy, vz
    integer, intent(in) :: idir
    type(coef_t), intent(in) :: coef
    type(c_ptr) :: vx_d, vy_d, vz_d
 
    if (coef%cyclic .and. coef%cyc_msk(0) .gt. 1) then
       if (neko_bcknd_device .eq. 1) then
          call neko_log%deprecated('Operator: rotate_cyc_r1, implicit device', &
               '2.0.0', 'Please call rotate_cyc_d instead.')
 
          vx_d = device_get_ptr(vx)
          vy_d = device_get_ptr(vy)
          vz_d = device_get_ptr(vz)
          call opr_device_rotate_cyc(vx_d, vy_d, vz_d, idir, coef)
       else
          call opr_cpu_rotate_cyc_r1(vx, vy, vz, idir, coef)
       end if
    end if
  end subroutine rotate_cyc_r1
 
  subroutine rotate_cyc_r4(vx, vy, vz, idir, coef)
    real(kind=rp), contiguous, dimension(:,:,:,:), intent(inout) :: vx, vy, vz
    integer, intent(in) :: idir
    type(coef_t), intent(in) :: coef
    type(c_ptr) :: vx_d, vy_d, vz_d
 
    if (coef%cyclic .and. coef%cyc_msk(0) .gt. 1) then
       if (neko_bcknd_device .eq. 1) then
          call neko_log%deprecated('Operator: rotate_cyc_r4, implicit device', &
               '2.0.0', 'Please call rotate_cyc_d instead.')
 
          vx_d = device_get_ptr(vx)
          vy_d = device_get_ptr(vy)
          vz_d = device_get_ptr(vz)
          call opr_device_rotate_cyc(vx_d, vy_d, vz_d, idir, coef)
       else
          call opr_cpu_rotate_cyc_r4(vx, vy, vz, idir, coef)
       end if
    end if
  end subroutine rotate_cyc_r4
 
  subroutine rotate_cyc_d(vx_d, vy_d, vz_d, idir, coef)
    type(c_ptr), intent(inout) :: vx_d, vy_d, vz_d
    integer, intent(in) :: idir
    type(coef_t), intent(in) :: coef
 
    if (coef%cyclic .and. coef%cyc_msk(0) .gt. 1) then
       call opr_device_rotate_cyc(vx_d, vy_d, vz_d, idir, coef)
    end if
  end subroutine rotate_cyc_d
 
  subroutine rotate_cyc_f(vx, vy, vz, idir, coef)
    type(field_t), intent(inout) :: vx, vy, vz
    integer, intent(in) :: idir
    type(coef_t), intent(in) :: coef
 
    if (coef%cyclic .and. coef%cyc_msk(0) .gt. 1) then
       if (neko_bcknd_device .eq. 1) then
          call opr_device_rotate_cyc(vx%x_d, vy%x_d, vz%x_d, idir, coef)
       else
          call opr_cpu_rotate_cyc_r4(vx%x, vy%x, vz%x, idir, coef)
       end if
    end if
  end subroutine rotate_cyc_f
 
end module operators