fdm.f90

Go to the documentation of this file.

! Copyright (c) 2008-2020, UCHICAGO ARGONNE, LLC.
!
! 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 fdm
  use neko_config, only : neko_bcknd_device, neko_bcknd_sx, neko_bcknd_xsmm
  use num_types, only : rp, sp, dp, qp
  use mesh, only : mesh_t
  use space, only : space_t
  use dofmap, only : dofmap_t
  use gather_scatter, only : gs_t, gs_op_add
  use fdm_cpu, only : fdm_do_fast_cpu
  use fdm_device, only : fdm_do_fast_device
  use fdm_sx, only : fdm_do_fast_sx
  use fdm_xsmm, only : fdm_do_fast_xsmm
  use utils, only : neko_error, neko_warning
  use comm, only : pe_rank
  use math, only : vlmax
  use device, only : glb_cmd_queue, device_to_host, host_to_device, &
       device_memcpy, device_map, device_free
  use fast3d, only : semhat
  use tensor, only : trsp
  use math, only : rzero, row_zero
  use, intrinsic :: iso_c_binding
  implicit none
  private
 
  type, public :: fdm_t
     real(kind=rp), allocatable :: s(:,:,:,:)
     real(kind=rp), allocatable :: d(:,:)
     type(c_ptr) :: s_d = c_null_ptr
     type(c_ptr) :: d_d = c_null_ptr
     real(kind=rp), allocatable :: len_lr(:), len_ls(:), len_lt(:)
     real(kind=rp), allocatable :: len_mr(:), len_ms(:), len_mt(:)
     real(kind=rp), allocatable :: len_rr(:), len_rs(:), len_rt(:)
     real(kind=rp), allocatable :: swplen(:,:,:,:)
     type(c_ptr) :: swplen_d = c_null_ptr
     type(space_t), pointer :: xh => null()
     type(dofmap_t), pointer :: dof => null()
     type(gs_t), pointer :: gs_h => null()
     type(mesh_t), pointer :: msh => null()
   contains
     procedure, pass(this) :: init => fdm_init
     procedure, pass(this) :: free => fdm_free
     procedure, pass(this) :: compute => fdm_compute
  end type fdm_t
 
  interface sygv
     module procedure sp_sygv, dp_sygv, qp_sygv
  end interface sygv
 
contains
 
  subroutine fdm_init(this, Xh, dof, gs_h)
    class(fdm_t), intent(inout) :: this
    type(space_t), target, intent(inout) :: Xh
    type(dofmap_t), target, intent(in) :: dof
    type(gs_t), target, intent(inout) :: gs_h
    !We only really use ah, bh
    real(kind=rp), dimension((Xh%lx)**2) :: ah, bh, ch, dh, zh
    real(kind=rp), dimension((Xh%lx)**2) :: dph, jph, bgl, zglhat, dgl, jgl, wh
    integer :: nl, n, nelv
 
    n = xh%lx - 1 !Polynomnial degree
    nl = xh%lx + 2 !Schwarz!
    nelv = dof%msh%nelv
    call fdm_free(this)
    allocate(this%s(nl*nl, 2, dof%msh%gdim, dof%msh%nelv))
    allocate(this%d(nl**3, dof%msh%nelv))
    allocate(this%swplen(xh%lx, xh%lx, xh%lx, dof%msh%nelv))
    allocate(this%len_lr(nelv), this%len_ls(nelv), this%len_lt(nelv))
    allocate(this%len_mr(nelv), this%len_ms(nelv), this%len_mt(nelv))
    allocate(this%len_rr(nelv), this%len_rs(nelv), this%len_rt(nelv))
 
    ! Zeroing here enables easier debugging since then
    ! MPI messages in GS are deterministic
    call rzero(this%swplen, xh%lxyz * dof%msh%nelv)
 
    if (neko_bcknd_device .eq. 1) then
       call device_map(this%s, this%s_d, nl*nl*2*dof%msh%gdim*dof%msh%nelv)
       call device_map(this%d, this%d_d, nl**dof%msh%gdim*dof%msh%nelv)
       call device_map(this%swplen, this%swplen_d, xh%lxyz*dof%msh%nelv)
    end if
 
    call semhat(ah, bh, ch, dh, zh, dph, jph, bgl, zglhat, dgl, jgl, n, wh)
    this%Xh => xh
    this%dof => dof
    this%gs_h => gs_h
    this%msh => dof%msh
 
    call swap_lengths(this, dof%x, dof%y, dof%z, dof%msh%nelv, dof%msh%gdim)
 
    call fdm_setup_fast(this, ah, bh, nl, n)
 
    if (neko_bcknd_device .eq. 1) then
       call device_memcpy(this%s, this%s_d, &
            nl*nl*2*dof%msh%gdim*dof%msh%nelv, host_to_device, sync = .false.)
       call device_memcpy(this%d, this%d_d, &
            nl**dof%msh%gdim*dof%msh%nelv, host_to_device, sync = .false.)
       call device_memcpy(this%swplen, this%swplen_d, &
            xh%lxyz*dof%msh%nelv, host_to_device, sync = .false.)
    end if
  end subroutine fdm_init
 
  subroutine swap_lengths(this, x, y, z, nelv, gdim)
    type(fdm_t), intent(inout) :: this
    integer, intent(in) :: gdim, nelv
    real(kind=rp), dimension(this%Xh%lxyz, nelv), intent(in) :: x, y, z
    integer :: j, k, e, n2, nz0, nzn, nx, lx1, n
 
    associate(l => this%swplen, xh => this%Xh, &
         llr => this%len_lr, lls => this%len_ls, llt => this%len_lt, &
         lmr => this%len_mr, lms => this%len_ms, lmt => this%len_mt, &
         lrr => this%len_rr, lrs => this%len_rs, lrt => this%len_rt)
      lx1 = this%Xh%lx
      n2 = lx1 - 1
      nz0 = 1
      nzn = 1
      nx = lx1 - 2
      if (gdim .eq. 3) then
         nz0 = 0
         nzn = n2
      end if
      call plane_space(lmr, lms, lmt, 0, n2, xh%wx, x, y, z, &
           nx, n2, nz0, nzn, nelv, gdim)
      n = n2 + 1
      if (gdim .eq. 3) then
         do e = 1, nelv
            do j = 2, n2
               do k = 2, n2
                  l(1, k, j, e) = lmr(e)
                  l(n, k, j, e) = lmr(e)
                  l(k, 1, j, e) = lms(e)
                  l(k, n, j, e) = lms(e)
                  l(k, j, 1, e) = lmt(e)
                  l(k, j, n, e) = lmt(e)
               end do
            end do
         end do
         if (neko_bcknd_device .eq. 1) then
            call device_memcpy(l, this%swplen_d, this%dof%size(), &
                 host_to_device, sync = .false.)
            call this%gs_h%op(l, this%dof%size(), gs_op_add)
            call device_memcpy(l, this%swplen_d, this%dof%size(), &
                 device_to_host, sync = .true.)
         else
            call this%gs_h%op(l, this%dof%size(), gs_op_add)
         end if
 
         do e = 1, nelv
            llr(e) = l(1, 2, 2, e) - lmr(e)
            lrr(e) = l(n, 2, 2, e) - lmr(e)
            lls(e) = l(2, 1, 2, e) - lms(e)
            lrs(e) = l(2, n, 2, e) - lms(e)
            llt(e) = l(2, 2, 1, e) - lmt(e)
            lrt(e) = l(2, 2, n, e) - lmt(e)
         end do
      else
         do e = 1, nelv
            do j = 2, n2
               l(1, j, 1, e) = lmr(e)
               l(n, j, 1, e) = lmr(e)
               l(j, 1, 1, e) = lms(e)
               l(j, n, 1, e) = lms(e)
            end do
         end do
 
         if (neko_bcknd_device .eq. 1) then
            call device_memcpy(l, this%swplen_d, this%dof%size(), &
                 host_to_device, sync = .false.)
            call this%gs_h%op(l, this%dof%size(), gs_op_add)
            call device_memcpy(l, this%swplen_d, this%dof%size(), &
                 device_to_host, sync = .true.)
         else
            call this%gs_h%op(l, this%dof%size(), gs_op_add)
         end if
 
         do e = 1, nelv
            llr(e) = l(1, 2, 1, e) - lmr(e)
            lrr(e) = l(n, 2, 1, e) - lmr(e)
            lls(e) = l(2, 1, 1, e) - lms(e)
            lrs(e) = l(2, n, 1, e) - lms(e)
         end do
      end if
    end associate
  end subroutine swap_lengths
 
  subroutine plane_space(lr, ls, lt, i1, i2, w, x, y, z, &
       nx, nxn, nz0, nzn, nelv, gdim)
    integer, intent(in) :: nxn, nzn, i1, i2, nelv, gdim, nx, nz0
    real(kind=rp), intent(inout) :: lr(nelv), ls(nelv), lt(nelv)
    real(kind=rp), intent(inout) :: w(nx)
    real(kind=rp), intent(in) :: x(0:nxn, 0:nxn, nz0:nzn, nelv)
    real(kind=rp), intent(in) :: y(0:nxn, 0:nxn, nz0:nzn, nelv)
    real(kind=rp), intent(in) :: z(0:nxn, 0:nxn, nz0:nzn, nelv)
    real(kind=rp) :: lr2, ls2, lt2, weight, wsum
    integer :: ny, nz, j1, k1, j2, k2, i, j, k, ie
    ny = nx
    nz = nx
    j1 = i1
    k1 = i1
    j2 = i2
    k2 = i2
    !   Now, for each element, compute lr,ls,lt between specified planes
    do ie = 1, nelv
       if (gdim .eq. 3) then
          lr2 = 0d0
          wsum = 0d0
          do k = 1, nz
             do j = 1, ny
                weight = w(j)*w(k)
                lr2 = lr2 + weight /&
                     ( (x(i2, j, k, ie) - x(i1, j, k, ie))**2&
                     + (y(i2, j, k, ie) - y(i1, j, k, ie))**2&
                     + (z(i2, j, k, ie) - z(i1, j, k, ie))**2 )
                wsum = wsum + weight
             end do
          end do
          lr2 = lr2/wsum
          lr(ie) = 1d0/sqrt(lr2)
          ls2 = 0d0
          wsum = 0d0
          do k = 1, nz
             do i = 1, nx
                weight = w(i)*w(k)
                ls2 = ls2 + weight / &
                     ( (x(i, j2, k, ie) - x(i, j1, k, ie))**2 &
                     + (y(i, j2, k, ie) - y(i, j1, k, ie))**2 &
                     + (z(i, j2, k, ie) - z(i, j1, k, ie))**2 )
                wsum = wsum + weight
             end do
          end do
          ls2 = ls2/wsum
          ls(ie) = 1d0/sqrt(ls2)
          lt2 = 0d0
          wsum = 0d0
          do j = 1, ny
             do i = 1, nx
                weight = w(i)*w(j)
                lt2 = lt2 + weight / &
                     ( (x(i, j, k2, ie) - x(i, j, k1, ie))**2 &
                     + (y(i, j, k2, ie) - y(i, j, k1, ie))**2 &
                     + (z(i, j, k2, ie) - z(i, j, k1, ie))**2 )
                wsum = wsum + weight
             end do
          end do
          lt2 = lt2/wsum
          lt(ie) = 1d0/sqrt(lt2)
       else ! 2D
          lr2 = 0d0
          wsum = 0d0
          do j = 1, ny
             weight = w(j)
             lr2 = lr2 + weight / &
                  ( (x(i2, j, 1, ie) - x(i1, j, 1, ie))**2 &
                  + (y(i2, j, 1, ie) - y(i1, j, 1, ie))**2 )
             wsum = wsum + weight
          end do
          lr2 = lr2/wsum
          lr(ie) = 1d0/sqrt(lr2)
          ls2 = 0d0
          wsum = 0d0
          do i = 1, nx
             weight = w(i)
             ls2 = ls2 + weight / &
                  ( (x(i, j2, 1, ie) - x(i, j1, 1, ie))**2 &
                  + (y(i, j2, 1, ie) - y(i, j1, 1, ie))**2 )
             wsum = wsum + weight
          end do
          ls2 = ls2/wsum
          ls(ie) = 1d0/sqrt(ls2)
       end if
    end do
    ie = 1014
  end subroutine plane_space
 
  subroutine fdm_setup_fast(this, ah, bh, nl, n)
    integer, intent(in) :: nl, n
    type(fdm_t), intent(inout) :: this
    real(kind=rp), intent(inout) :: ah(n+1, n+1), bh(n+1)
    real(kind=rp), dimension(2*this%Xh%lx + 4) :: lr, ls, lt
    integer :: i, j, k
    integer :: ie, il, nr, ns, nt
    integer :: lbr, rbr, lbs, rbs, lbt, rbt
    real(kind=rp) :: eps, diag
 
    associate(s => this%s, d => this%d, &
         llr => this%len_lr, lls => this%len_ls, llt => this%len_lt, &
         lmr => this%len_mr, lms => this%len_ms, lmt => this%len_mt, &
         lrr => this%len_rr, lrs => this%len_rs, lrt => this%len_rt)
      do ie = 1, this%dof%msh%nelv
         lbr = this%dof%msh%facet_type(1, ie)
         rbr = this%dof%msh%facet_type(2, ie)
         lbs = this%dof%msh%facet_type(3, ie)
         rbs = this%dof%msh%facet_type(4, ie)
         lbt = this%dof%msh%facet_type(5, ie)
         rbt = this%dof%msh%facet_type(6, ie)
 
         nr = nl
         ns = nl
         nt = nl
         call fdm_setup_fast1d(s(1, 1, 1, ie), lr, nr, lbr, rbr, &
              llr(ie), lmr(ie), lrr(ie), ah, bh, n)
         call fdm_setup_fast1d(s(1, 1, 2, ie), ls, ns, lbs, rbs, &
              lls(ie), lms(ie), lrs(ie), ah, bh, n)
         if (this%dof%msh%gdim .eq. 3) then
            call fdm_setup_fast1d(s(1, 1, 3, ie), lt, nt, lbt, rbt, &
                 llt(ie), lmt(ie), lrt(ie), ah, bh, n)
         end if
 
         il = 1
         if (.not. this%dof%msh%gdim .eq. 3) then
            eps = 1d-5 * (vlmax(lr(2), nr - 2) + vlmax(ls(2), ns - 2))
            do j = 1, ns
               do i = 1, nr
                  diag = lr(i) + ls(j)
                  if (diag .gt. eps) then
                     d(il, ie) = 1.0_rp / diag
                  else
                     d(il, ie) = 0.0_rp
                  end if
                  il = il + 1
               end do
            end do
         else
            eps = 1d-5 * (vlmax(lr(2), nr - 2) + &
                 vlmax(ls(2), ns - 2) + vlmax(lt(2), nt - 2))
            do k = 1, nt
               do j = 1, ns
                  do i = 1, nr
                     diag = lr(i) + ls(j) + lt(k)
                     if (diag .gt. eps) then
                        d(il, ie) = 1.0_rp / diag
                     else
                        d(il, ie) = 0.0_rp
                     end if
                     il = il + 1
                  end do
               end do
            end do
         end if
      end do
    end associate
 
  end subroutine fdm_setup_fast
 
  subroutine fdm_setup_fast1d(s, lam, nl, lbc, rbc, ll, lm, lr, ah, bh, n)
    integer, intent(in) :: nl, lbc, rbc, n
    real(kind=rp), intent(inout) :: s(nl, nl, 2), lam(nl), ll, lm, lr
    real(kind=rp), intent(inout) :: ah(0:n, 0:n), bh(0:n)
    integer :: lx1, lxm
    real(kind=rp) :: b(2*(n+3)**2)
 
    lx1 = n + 1
    lxm = lx1 + 2
 
    call fdm_setup_fast1d_a(s, lbc, rbc, ll, lm, lr, ah, n)
    call fdm_setup_fast1d_b(b, lbc, rbc, ll, lm, lr, bh, n)
    call generalev(s, b, lam, nl, lx1)
    if (lbc .gt. 0) call row_zero(s, nl, nl, 1)
    if (lbc .eq. 1) call row_zero(s, nl, nl, 2)
    if (rbc .gt. 0) call row_zero(s, nl, nl, nl)
    if (rbc .eq. 1) call row_zero(s, nl, nl, nl-1)
 
    call trsp(s(1, 1, 2), nl, s, nl)
 
  end subroutine fdm_setup_fast1d
 
  subroutine generalev(a, b, lam, n, lx)
    integer, intent(in) :: n, lx
    real(kind=rp), intent(inout) :: a(n, n), b(n, n), lam(n)
    integer :: lbw, lw
    real(kind=rp) :: bw(4*(lx+2)**3)
 
    lbw = 4*(lx+2)**3
    lw = n*n
    call sygv(a, b, lam, n, lx, bw, lbw)
 
  end subroutine generalev
 
  subroutine sp_sygv(a, b, lam, n, lx, bw, lbw)
    integer, intent(in) :: n, lx, lbw
    real(kind=sp), intent(inout) :: a(n, n), b(n, n), lam(n)
    real(kind=sp) :: bw(4*(lx+2)**3)
    integer :: info = 0
    call ssygv(1, 'V', 'U', n, a, n, b, n, lam, bw, lbw, info)
  end subroutine sp_sygv
 
  subroutine dp_sygv(a, b, lam, n, lx, bw, lbw)
    integer, intent(in) :: n, lx, lbw
    real(kind=dp), intent(inout) :: a(n, n), b(n, n), lam(n)
    real(kind=dp) :: bw(4*(lx+2)**3)
    integer :: info = 0
    call dsygv(1, 'V', 'U', n, a, n, b, n, lam, bw, lbw, info)
  end subroutine dp_sygv
 
  subroutine qp_sygv(a, b, lam, n, lx, bw, lbw)
    integer, intent(in) :: n, lx, lbw
    real(kind=qp), intent(inout) :: a(n, n), b(n, n), lam(n)
    real(kind=dp) :: a2(n, n), b2(n, n), lam2(n)
    real(kind=qp) :: bw(4*(lx+2)**3)
    real(kind=dp) :: bw2(4*(lx+2)**3)
    integer :: info = 0
 
    a2 = real(a, dp)
    b2 = real(b, dp)
    lam2 = real(lam, dp)
    call dsygv(1, 'V', 'U', n, a2, n, b2, n, lam2, bw2, lbw, info)
    a = real(a2, qp)
    b = real(b2, qp)
    lam = real(lam2, qp)
    if (pe_rank .eq. 0) then
       call neko_warning('Real precision choice not supported for fdm, ' // &
            'treating it as double')
    end if
 
  end subroutine qp_sygv
 
  subroutine fdm_setup_fast1d_a(a, lbc, rbc, ll, lm, lr, ah, n)
    integer, intent(in) ::lbc, rbc, n
    real(kind=rp), intent(inout) :: a(0:n+2, 0:n+2), ll, lm, lr
    real(kind=rp), intent(inout) :: ah(0:n, 0:n)
    real(kind=rp) :: fac
    integer :: i, j, i0, i1
 
    i0 = 0
    if (lbc .eq. 1) i0 = 1
    i1 = n
    if (rbc .eq. 1) i1 = n - 1
 
    call rzero(a, (n+3) * (n+3))
 
    fac = 2.0_rp / lm
    a(1, 1) = 1.0_rp
    a(n+1, n+1) = 1.0-rp
 
    do j = i0, i1
       do i = i0, i1
          a(i+1, j+1) = fac * ah(i, j)
       end do
    end do
 
    if (lbc .eq. 0) then
       fac = 2.0_rp / ll
       a(0, 0) = fac * ah(n-1, n-1)
       a(1, 0) = fac * ah(n, n-1)
       a(0, 1) = fac * ah(n-1, n)
       a(1, 1) = a(1, 1) + fac * ah(n, n)
    else
       a(0, 0) = 1.0_rp
    end if
 
    if (rbc .eq. 0) then
       fac = 2.0_rp / lr
       a(n+1, n+1) = a(n+1, n+1) + fac * ah(0, 0)
       a(n+2, n+1) = fac * ah(1, 0)
       a(n+1, n+2) = fac * ah(0, 1)
       a(n+2, n+2) = fac * ah(1, 1)
    else
       a(n+2, n+2) = 1.0_rp
    end if
 
  end subroutine fdm_setup_fast1d_a
 
  subroutine fdm_setup_fast1d_b(b, lbc, rbc, ll, lm, lr, bh, n)
    integer, intent(in) :: lbc, rbc, n
    real(kind=rp), intent(inout) :: b(0:n+2, 0:n+2), ll, lm, lr
    real(kind=rp), intent(inout) :: bh(0:n)
    real(kind=rp) :: fac
    integer :: i, i0, i1
 
    i0 = 0
    if (lbc .eq. 1) i0 = 1
    i1 = n
    if (rbc .eq. 1) i1 = n - 1
 
    call rzero(b, (n + 3) * (n + 3))
 
    fac = 0.5_rp * lm
    b(1, 1) = 1.0_rp
    b(n+1, n+1) = 1.0_rp
 
    do i = i0, i1
       b(i+1, i+1) = fac * bh(i)
    end do
 
    if (lbc .eq. 0) then
       fac = 0.5_rp * ll
       b(0, 0) = fac * bh(n-1)
       b(1, 1) = b(1, 1) + fac * bh(n)
    else
       b(0, 0) = 1.0_rp
    end if
 
    if (rbc .eq. 0) then
       fac = 0.5_rp * lr
       b(n+1, n+1) = b(n+1, n+1) + fac * bh(0)
       b(n+2, n+2) = fac * bh(1)
    else
       b(n+2, n+2) = 1.0_rp
    end if
 
  end subroutine fdm_setup_fast1d_b
 
  subroutine fdm_free(this)
    class(fdm_t), intent(inout) :: this
 
    if (allocated(this%s)) then
       deallocate(this%s)
    end if
 
    if (allocated(this%d)) then
       deallocate(this%d)
    end if
 
    if (allocated(this%len_lr)) then
       deallocate(this%len_lr)
    end if
 
    if (allocated(this%len_ls)) then
       deallocate(this%len_ls)
    end if
 
    if (allocated(this%len_lt)) then
       deallocate(this%len_lt)
    end if
 
    if (allocated(this%len_mr)) then
       deallocate(this%len_mr)
    end if
 
    if (allocated(this%len_ms)) then
       deallocate(this%len_ms)
    end if
 
    if (allocated(this%len_mt)) then
       deallocate(this%len_mt)
    end if
 
    if (allocated(this%len_rr)) then
       deallocate(this%len_rr)
    end if
 
    if (allocated(this%len_rs)) then
       deallocate(this%len_rs)
    end if
 
    if (allocated(this%len_rt)) then
       deallocate(this%len_rt)
    end if
 
    if (allocated(this%swplen)) then
       deallocate(this%swplen)
    end if
 
    nullify(this%Xh)
    nullify(this%dof)
    nullify(this%gs_h)
    nullify(this%msh)
 
    if (c_associated(this%s_d)) then
       call device_free(this%s_d)
    end if
    if (c_associated(this%d_d)) then
       call device_free(this%d_d)
    end if
    if (c_associated(this%swplen_d)) then
       call device_free(this%swplen_d)
    end if
 
  end subroutine fdm_free
 
  subroutine fdm_compute(this, e, r, stream)
    class(fdm_t), intent(inout) :: this
    real(kind=rp), dimension((this%Xh%lx + 2)**3, this%msh%nelv), &
         intent(inout) :: e, r
    type(c_ptr), optional :: stream
    type(c_ptr) :: strm
 
    if (present(stream)) then
       strm = stream
    else
       strm = glb_cmd_queue
    end if
 
    if (neko_bcknd_sx .eq. 1) then
       call fdm_do_fast_sx(e, r, this%s, this%d, &
            this%Xh%lx+2, this%msh%gdim, this%msh%nelv)
    else if (neko_bcknd_xsmm .eq. 1) then
       call fdm_do_fast_xsmm(e, r, this%s, this%d, &
            this%Xh%lx+2, this%msh%gdim, this%msh%nelv)
    else if (neko_bcknd_device .eq. 1) then
       call fdm_do_fast_device(e, r, this%s, this%d, &
            this%Xh%lx+2, this%msh%gdim, this%msh%nelv, strm)
    else
       call fdm_do_fast_cpu(e, r, this%s, this%d, &
            this%Xh%lx+2, this%msh%gdim, this%msh%nelv)
    end if
 
  end subroutine fdm_compute
 
 
end module fdm