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
  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_dudxyz, &
       opr_device_lambda2, opr_device_opgrad
  use space, only : space_t
  use coefs, only : coef_t
  use field, only : field_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
  use device, only : device_get_ptr
  use device_math, only : device_add2, device_cmult, device_copy, device_glsum, &
       device_cadd
  use scratch_registry, only : neko_scratch_registry
  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
  implicit none
  private
 
  public :: dudxyz, opgrad, ortho, cdtp, conv1, curl, cfl, cfl_compressible, &
       lambda2op, strain_rate, div, grad, set_convect_rst, runge_kutta
 
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 div(res, ux, uy, uz, 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) :: res
    real(kind=rp), dimension(coef%Xh%lx, coef%Xh%ly, coef%Xh%lz, &
         coef%msh%nelv), intent(in) :: ux, uy, uz
    type(field_t), pointer :: work
    integer :: ind
    type(c_ptr) :: res_d
 
    if (neko_bcknd_device .eq. 1) then
       res_d = device_get_ptr(res)
    end if
 
    call neko_scratch_registry%request_field(work, ind)
 
    ! 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
 
  subroutine grad(ux, uy, uz, u, coef)
    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
 
    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
 
  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, 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
       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
 
    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
       call opr_device_cdtp(dtx, x, dr, ds, dt, 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(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)
      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
 
  subroutine convect_scalar(du, u, c, 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)
    real(kind=rp), intent(inout) :: c(xh_gl%lxyz, coef_gl%msh%nelv, 3)
 
    if (neko_bcknd_sx .eq. 1) then
       call opr_sx_convect_scalar(du, u, c, 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, c, xh_gll, xh_gl, &
            coef_gll, coef_gl, gll_to_gl)
    else
       call opr_cpu_convect_scalar(du, u, c, 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, 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
       if (present(event)) then
          call opr_device_curl(w1, w2, w3, u1, u2, u3, &
               work1, work2, coef, event)
       else
          call opr_device_curl(w1, w2, w3, u1, u2, u3, work1, work2, coef)
       end if
    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)
    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
 
  !! 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(dt, max_wave_speed, 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) :: max_wave_speed
    real(kind=rp) :: cfl_compressible
    integer :: ierr, n
    type(field_t), pointer :: zero_vector
    integer :: ind
 
    n = xh%lx * xh%ly * xh%lz * nelv
 
    ! Request a scratch field for zero vector
    call neko_scratch_registry%request_field(zero_vector, ind)
 
    ! Initialize zero vector
    call field_rzero(zero_vector)
 
    ! Use incompressible CFL with max_wave_speed as u-component, zero v and w
    cfl_compressible = cfl(dt, max_wave_speed, zero_vector%x, zero_vector%x, xh, coef, nelv, gdim)
 
    ! Release the scratch field
    call neko_scratch_registry%relinquish_field(ind)
 
  end function cfl_compressible
 
  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) :: s13(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)
 
    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)
    end if
 
    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)
    end if
 
    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)
    end if
 
    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)
    end if
 
    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(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 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
 
  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
 
  subroutine set_convect_rst(cr, cs, ct, cx, cy, cz, Xh, coef)
    type(space_t), intent(inout) :: xh
    type(coef_t), intent(inout) :: coef
    real(kind=rp), dimension(Xh%lxyz, coef%msh%nelv), &
         intent(inout) :: cr, cs, ct
    real(kind=rp), dimension(Xh%lxyz, coef%msh%nelv), &
         intent(in) :: cx, cy, cz
 
    if (neko_bcknd_sx .eq. 1) then
       call opr_sx_set_convect_rst(cr, cs, ct, cx, cy, cz, xh, coef)
    else if (neko_bcknd_xsmm .eq. 1) then
       call opr_xsmm_set_convect_rst(cr, cs, ct, cx, cy, cz, xh, coef)
    else
       call opr_cpu_set_convect_rst(cr, cs, ct, cx, cy, cz, xh, coef)
    end if
 
  end subroutine set_convect_rst
 
  subroutine runge_kutta(phi, c_r1, c_r23, c_r4, 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
    real(kind=rp), dimension(n), intent(inout) :: phi
    real(kind=rp), dimension(3 * n_GL), intent(inout) :: c_r1, c_r23, c_r4
    real(kind=rp) :: c1, c2, c3
    ! Work Arrays
    real(kind=rp), dimension(n) :: u1, r1, r2, r3, r4
    real(kind=rp), dimension(n_GL) :: u1_gl
    integer :: i, e
 
    c1 = 1.
    c2 = -dtau/2.
    c3 = -dtau
 
    ! Stage 1:
    call invcol3 (u1, phi, coef%B, n)
    call gll_to_gl%map(u1_gl, u1, nel, xh_gl)
    call convect_scalar(r1, u1_gl, c_r1, xh_gll, xh_gl, coef, &
         coef_gl, gll_to_gl)
    call col2(r1, coef%B, n)
 
    ! Stage 2:
    call add3s2 (u1, phi, r1, c1, c2, n)
    call invcol2 (u1, coef%B, n)
    call gll_to_gl%map(u1_gl, u1, nel, xh_gl)
    call convect_scalar(r2, u1_gl, c_r23, xh_gll, xh_gl, coef, &
         coef_gl, gll_to_gl)
    call col2(r2, coef%B, n)
 
    ! Stage 3:
    call add3s2 (u1, phi, r2, c1, c2, n)
    call invcol2 (u1, coef%B, n)
    call gll_to_gl%map(u1_gl, u1, nel, xh_gl)
    call convect_scalar(r3, u1_gl, c_r23, xh_gll, xh_gl, coef, &
         coef_gl, gll_to_gl)
    call col2(r3, coef%B, n)
 
    ! Stage 4:
    call add3s2 (u1, phi, r3, c1, c3, n)
    call invcol2 (u1, coef%B, n)
    call gll_to_gl%map(u1_gl, u1, nel, xh_gl)
    call convect_scalar(r4, u1_gl, c_r4, xh_gll, xh_gl, coef, &
         coef_gl, gll_to_gl)
    call col2(r4, coef%B, n)
 
    c1 = -dtau/6.
    c2 = -dtau/3.
    do i = 1, n
       phi(i) = phi(i) + c1 * (r1(i) + r4(i)) + c2 * (r2(i) + r3(i))
    end do
 
  end subroutine runge_kutta
 
 
 
end module operators