scale_meshfield_spectral_transform.F90

Go to the documentation of this file.

!> module FElib / Data / Utility
!!
!! @par Description
!!           A module for providing utility routines associated with spectral transformation for DG fields
!! @author Yuta Kawai, Team SCALE
!!
!<
#include "scaleFElib.h"
module scale_meshfield_spectral_transform
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_const, only: &
    pi => const_pi
  use scale_prc, only: &
    prc_abort
  
  use scale_polynomial, only: &
    polynomial_genlagrangepoly
  use scale_element_base, only: &
    elementbase1d, elementbase2d, elementbase3d
  use scale_element_line, only: &
      lineelement
  use scale_localmesh_1d, only: localmesh1d
  use scale_localmesh_2d, only: localmesh2d
  use scale_localmesh_3d, only: localmesh3d
  use scale_mesh_base1d, only: meshbase1d
  use scale_mesh_base2d, only: meshbase2d
 
  use scale_localmeshfield_base, only: &
    localmeshfield1d,      &
    localmeshfieldbaselist
  use scale_meshfield_base, only: &
    meshfield1d, meshfield2d, &
    meshfield1dlist, meshfield2dlist
 
  use scale_fast_fourier_transform, only: &
    fastfouriertransform1d
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  !
 
  type, public :: meshfield_spetraltransformbase
    integer :: eval_type_id
    integer :: ndim
    integer :: var_num
 
    real(rp), allocatable :: fftintrpmat(:,:)
    real(rp), allocatable :: fft_xi(:)
    integer :: nsampleptperelem
 
    real(rp), allocatable :: intintrpmat(:,:)
    real(rp), allocatable :: intweight(:)
    real(rp), allocatable :: intxi(:)
    integer :: nintglpt
  end type meshfield_spetraltransformbase
 
  !> A derived type for spectral transform of MeshField1D
  type, public, extends(meshfield_spetraltransformbase) :: meshfield_spetraltransform1d
    integer :: kall
    integer :: ks, ke
 
    real(rp) :: xmin_gl, xmax_gl
    real(rp) :: delx
 
    real(rp), allocatable :: k(:)
    real(rp), allocatable :: spectral_coef(:,:,:)
 
    type(fastfouriertransform1d) :: fft
  contains
    procedure :: init => meshfield_spetraltransform1d_init
    procedure :: final => meshfield_spetraltransform1d_final
    procedure :: transform => meshfield_spetraltransform1d_transform
  end type meshfield_spetraltransform1d
 
  !> A derived type for spectral transform of MeshField2D
  type, public, extends(meshfield_spetraltransformbase) :: meshfield_spetraltransform2d
    integer :: kall
    integer :: ks, ke
    integer :: lall
    integer :: ls, le
 
    real(rp) :: xmin_gl, xmax_gl
    real(rp) :: delx
    real(rp) :: ymin_gl, ymax_gl
    real(rp) :: dely
 
    integer :: nprcx    
    integer :: nprcy
    integer :: negx, negy
 
    integer :: nsampleptperelem1d
 
    real(rp), allocatable :: k(:)
    real(rp), allocatable :: l(:)
    real(rp), allocatable :: spectral_coef(:,:,:,:)
 
    type(fastfouriertransform1d) :: fft_x
    type(fastfouriertransform1d) :: fft_y
  contains
    procedure :: init => meshfield_spetraltransform2d_init
    procedure :: final => meshfield_spetraltransform2d_final
    procedure :: transform => meshfield_spetraltransform2d_transform
  end type meshfield_spetraltransform2d
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  integer, public, parameter :: st_evaltype_sample_uniform_pts = 1
  integer, public, parameter :: st_evaltype_l2projection_1     = 2
  integer, public, parameter :: st_evaltype_l2projection_2     = 3
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
 
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  !-----------------------------------------------------------------------------
  
contains
!- Base type
 
!OCL SERIAL
  subroutine meshfield_spetraltransformbase_init( this, eval_type, ndim, var_num, NintGLpt, NsamplePtPerElem )
    implicit none
    class(meshfield_spetraltransformbase), intent(inout) :: this
    integer, intent(in) :: eval_type
    integer, intent(in) :: ndim
    integer, intent(in) :: var_num
    integer, intent(in), optional :: nintglpt
    integer, intent(in), optional :: nsampleptperelem
    !--------------------------------------------
    
    this%eval_type_id = eval_type
    call check_eval_type_id( eval_type )
 
    this%ndim = ndim
    this%var_num = var_num
 
    if (     this%eval_type_id == st_evaltype_l2projection_1 &
        .or. this%eval_type_id == st_evaltype_l2projection_2 ) then
      if ( .not. present(nintglpt) ) then
        log_info("MeshField_SpetralTransformBase_Init",*) "The order of Gaussian quadrature should be given. Check!"
        call prc_abort
      end if
      this%NintGLPt = nintglpt
    else
      this%NintGLPt = -1
    end if
 
    if ( this%eval_type_id == st_evaltype_sample_uniform_pts ) then
      if ( .not. present(nsampleptperelem) ) then
        log_info("MeshField_SpetralTransformBase_Init",*) "NsamplePtPerElem shouled be given. Check!"
      else
        this%NsamplePtPerElem = nsampleptperelem
      end if
    else
      this%NsamplePtPerElem = -1
    end if
    log_info("MeshField_SpetralTransformBase_Init",*) "NintGLPt=", this%NintGLpt, "NsamplePtPerElem=", this%NsamplePtPerElem
 
    return
  end subroutine meshfield_spetraltransformbase_init
 
  subroutine check_eval_type_id( eval_type ) 
    implicit none
    integer, intent(in) :: eval_type
    !-------------------------
    select case(eval_type)
!    case(ST_EVALTYPE_L2PROJECTION_1,ST_EVALTYPE_L2PROJECTION_2,ST_EVALTYPE_SAMPLE_UNIFORM_PTS)
    case( st_evaltype_l2projection_1, st_evaltype_sample_uniform_pts )
    case default
      log_info("MeshField_SpetralTransformBase_Init",*) "Unsupported evaluation type is specified. Check!"
      call prc_abort
    end select
    return
  end subroutine check_eval_type_id
 
!OCL SERIAL
  subroutine meshfield_spetraltransformbase_final( this )
    implicit none
    class(meshfield_spetraltransformbase), intent(inout) :: this
    !--------------------------------------------
 
    if ( this%NintGLpt > 0 ) &
      deallocate( this%IntIntrpMat, this%IntXi, this%IntWeight )
    
    if ( this%NsamplePtPerElem > 0 ) &
      deallocate( this%FFTIntrpMat, this%FFT_xi )
 
    return
  end subroutine meshfield_spetraltransformbase_final  
 
!- 1D
 
!OCL SERIAL
  subroutine meshfield_spetraltransform1d_init( this, eval_type, ks, ke, mesh1D, var_num, &
    GLQuadOrd, NsamplePtPerElem )
    implicit none
    class(meshfield_spetraltransform1d), intent(inout) :: this
    integer, intent(in) :: eval_type
    integer, intent(in) :: ks, ke
    class(meshbase1d), intent(in), target :: mesh1d
    integer, intent(in) :: var_num
    integer, optional :: glquadord
    integer, optional :: nsampleptperelem
 
    integer :: i
    integer :: lx
    real(rp) :: dx
 
    class(localmesh1d), pointer :: lmesh
    class(elementbase1d), pointer :: elem
 
    type(lineelement) :: elem_dummy
    !--------------------------------------------
 
    lmesh => mesh1d%lcmesh_list(1)
    elem => lmesh%refElem1D
 
    call meshfield_spetraltransformbase_init( this, eval_type, 1, var_num, &
      glquadord, nsampleptperelem )
 
    this%ks = ks; this%ke = ke
    this%kall = ke - ks + 1
 
    this%xmin_gl = mesh1d%xmin_gl
    this%xmax_gl = mesh1d%xmax_gl
    this%delx = 1.0_rp / real(mesh1d%NeG,kind=rp)
 
    allocate( this%k(ks:ke) )
    lx = this%xmax_gl - this%xmin_gl
    do i=ks, ke
      this%k(i) = i * 2.0_rp * pi / lx
    end do
    
    allocate( this%spectral_coef(ks:ke,2,var_num) )
 
    call elem_dummy%Init(elem%PolyOrder, .false.)
 
    if ( this%NintGLpt > 0 ) then
      allocate( this%IntIntrpMat(this%NintGLpt,elem%Np) )
      allocate( this%IntXi(this%NintGLpt), this%IntWeight(this%NintGLpt) )
 
      this%IntIntrpMat(:,:) = elem_dummy%GenIntGaussLegendreIntrpMat( glquadord, this%IntWeight, this%IntXi )
    end if
 
    if ( this%NsamplePtPerElem > 0 ) then
      allocate( this%FFTIntrpMat(this%NsamplePtPerElem,elem%Np) )
      allocate( this%FFT_Xi(this%NsamplePtPerElem) )
 
      dx = 2.0_rp / real(this%NsamplePtPerElem,kind=rp)
      do i=1, this%NsamplePtPerElem
        this%FFT_Xi(i) = -1.0_rp + real(i - 0.5_rp,kind=rp) * dx
      end do
      this%FFTIntrpMat(:,:) = polynomial_genlagrangepoly( elem_dummy%PolyOrder, elem_dummy%x1, this%FFT_xi )
 
      call this%fft%Init( nsampleptperelem * mesh1d%NeG )
    end if
 
    call elem_dummy%Final()
 
    return
  end subroutine meshfield_spetraltransform1d_init
 
!OCL SERIAL
  subroutine meshfield_spetraltransform1d_final( this )
    implicit none
    class(meshfield_spetraltransform1d), intent(inout) :: this
    !--------------------------------------------
 
    call meshfield_spetraltransformbase_final( this )
 
    if ( this%NsamplePtPerElem > 0 ) then
      call this%fft%Final()
    end if
 
    deallocate( this%k, this%spectral_coef )
    
    return
  end subroutine meshfield_spetraltransform1d_final
 
!OCL SERIAL
  subroutine meshfield_spetraltransform1d_transform( this, q_list, mesh_num )
    implicit none
    class(meshfield_spetraltransform1d), intent(inout) :: this
    integer, intent(in) :: mesh_num
    type(meshfield1dlist), target :: q_list(this%var_num,mesh_num)
 
    class(meshbase1d), pointer :: mesh1d
    class(localmesh1d), pointer :: lmesh
    !--------------------------------------------
 
    mesh1d => q_list(1,1)%ptr%mesh
    lmesh => mesh1d%lcmesh_list(1)
    
    select case( this%eval_type_id )
    case( st_evaltype_sample_uniform_pts ) 
      call spectral_transform1d_dft( this%spectral_coef, & ! (out)
        q_list, this%var_num, mesh_num,                                        & ! (in)
        this%ks, this%ke, lmesh%refElem1D%Np, lmesh%Ne, this%NsamplePtPerElem, & ! (in)
        this%fft, this%FFTIntrpMat, this%FFT_xi,                               & ! (in)
        this%xmin_gl, this%delx, this%xmax_gl - this%xmin_gl                   ) ! (in)
    case( st_evaltype_l2projection_1)
      call spectral_transform1d_l2projection( this%spectral_coef,    & ! (out)
        q_list, this%var_num, mesh_num,                              & ! (in)
        this%ks, this%ke, lmesh%refElem1D%Np, lmesh%Ne,              & ! (in)
        this%IntIntrpMat, this%IntXi, this%IntWeight, this%NintGLPt, & ! (in)
        this%xmin_gl, this%delx, this%xmax_gl - this%xmin_gl         ) ! (in)
    end select
    
    return
  end subroutine meshfield_spetraltransform1d_transform
 
!- 2D
 
!OCL SERIAL
  subroutine meshfield_spetraltransform2d_init( this, eval_type, ks, ke, ls, le, mesh2D, var_num, &
    GLQuadOrd, NsamplePtPerElem1D )
    use scale_element_line, only: &
      lineelement
    use scale_mesh_rectdom2d, only: &
      meshrectdom2d
    implicit none
    class(meshfield_spetraltransform2d), intent(inout) :: this
    integer, intent(in) :: eval_type
    integer, intent(in) :: ks, ke
    integer, intent(in) :: ls, le
    class(meshbase2d), intent(in), target :: mesh2d
    integer, intent(in) :: var_num
    integer, optional :: glquadord
    integer, optional :: nsampleptperelem1d
 
    integer :: i, j
    integer :: lx, ly
    real(rp) :: dx
 
    class(localmesh2d), pointer :: lmesh
    class(elementbase2d), pointer :: elem
 
    type(lineelement) :: elem_dummy
    class(meshbase2d), pointer :: ptr_mesh2d
    !--------------------------------------------
 
    lmesh => mesh2d%lcmesh_list(1)
    elem => lmesh%refElem2D
 
    call meshfield_spetraltransformbase_init( this, eval_type, 1, var_num, glquadord, nsampleptperelem1d )
 
    this%ks = ks; this%ke = ke
    this%kall = ke - ks + 1
 
    this%ls = ls; this%le = le
    this%lall = le - ls + 1
 
    select type(ptr_mesh2d => mesh2d)
    class is (meshrectdom2d)
      this%xmin_gl = ptr_mesh2d%xmin_gl
      this%xmax_gl = ptr_mesh2d%xmax_gl
      this%delx = ( this%xmax_gl - this%xmin_gl ) / real(ptr_mesh2d%NeGX, kind=rp)
 
      this%ymin_gl = ptr_mesh2d%ymin_gl
      this%ymax_gl = ptr_mesh2d%ymax_gl
      this%dely = ( this%ymax_gl - this%ymin_gl ) / real(ptr_mesh2d%NeGY, kind=rp)
 
      this%NeGX = ptr_mesh2d%NeGX
      this%NeGY = ptr_mesh2d%NeGY
      this%NprcX = ptr_mesh2d%NprcX
      this%NprcY = ptr_mesh2d%NprcY
    class default
      log_info('MeshField_SpetralTransform2D_Init',*) 'Unexpected mesh type is given. Check!'
      call prc_abort
    end select
 
    allocate( this%k(ks:ke) )
    lx = this%xmax_gl - this%xmin_gl
    do i=ks, ke
      this%k(i) = i * 2.0_rp * pi / lx
    end do
 
    allocate( this%l(ls:le) )
    ly = this%ymax_gl - this%ymin_gl
    do j=ls, le
      this%l(j) = j * 2.0_rp * pi / ly
    end do
    
    allocate( this%spectral_coef(ks:ke,ls:le,2,var_num) )
 
    call elem_dummy%Init(elem%PolyOrder, .false.)
 
    if ( this%NintGLpt > 0 ) then
      allocate( this%IntIntrpMat(this%NintGLpt,elem%Nfp) )
      allocate( this%IntXi(this%NintGLpt), this%IntWeight(this%NintGLpt) )
 
      this%IntIntrpMat(:,:) = elem_dummy%GenIntGaussLegendreIntrpMat( glquadord, this%IntWeight, this%IntXi )
    end if
 
    if ( this%NsamplePtPerElem > 0 ) then
      this%NsamplePtPerElem1D = nsampleptperelem1d
      allocate( this%FFTIntrpMat(this%NsamplePtPerElem1D,elem%Nfp) )
      allocate( this%FFT_Xi(this%NsamplePtPerElem1D) )
 
      dx = 2.0_rp / real(this%NsamplePtPerElem1D,kind=rp)
      do i=1, this%NsamplePtPerElem1D
        this%FFT_Xi(i) = -1.0_rp + real(i - 0.5_rp,kind=rp) * dx
      end do
      this%FFTIntrpMat(:,:) = polynomial_genlagrangepoly( elem_dummy%PolyOrder, elem_dummy%x1, this%FFT_xi )
 
      call this%fft_x%Init( this%NsamplePtPerElem1D * this%NeGX )
      call this%fft_y%Init( this%NsamplePtPerElem1D * this%NeGY )
    end if
 
    call elem_dummy%Final()
    return
  end subroutine meshfield_spetraltransform2d_init
 
!OCL SERIAL
  subroutine meshfield_spetraltransform2d_final( this )
    implicit none
    class(meshfield_spetraltransform2d), intent(inout) :: this
    !--------------------------------------------
 
    call meshfield_spetraltransformbase_final( this )
 
    if ( this%NsamplePtPerElem > 0 ) then
      call this%fft_x%Final()
      call this%fft_y%Final()
    end if
    
    deallocate( this%k, this%l, this%spectral_coef )
    return
  end subroutine meshfield_spetraltransform2d_final
 
!OCL SERIAL
  subroutine meshfield_spetraltransform2d_transform( this, q_list, mesh_num_x, mesh_num_y )
    implicit none
    class(meshfield_spetraltransform2d), intent(inout) :: this
    integer, intent(in) :: mesh_num_x, mesh_num_y
    type(meshfield2dlist), target :: q_list(this%var_num,mesh_num_x,mesh_num_y)
 
    class(meshbase2d), pointer :: mesh2d
    class(localmesh2d), pointer :: lmesh
    !--------------------------------------------
 
    mesh2d => q_list(1,1,1)%ptr%mesh
    lmesh => mesh2d%lcmesh_list(1)
    
    select case( this%eval_type_id )
    case( st_evaltype_sample_uniform_pts )
      call spectral_transform2d_dft( this%spectral_coef, &
        q_list, this%var_num, mesh_num_x, mesh_num_y,                                   & ! (in)
        this%ks, this%ke,  this%ls, this%le, lmesh%refElem2D%Nfp, lmesh%NeX, lmesh%NeY, & ! (in)
        this%NeGY, this%NprcY, this%NsamplePtPerElem1D,                                 & ! (in)
        this%fft_x, this%fft_y, this%FFTIntrpMat )
    case( st_evaltype_l2projection_1)
      call spectral_transform2d_l2projection( this%spectral_coef,                       & ! (out)
        q_list, this%var_num, mesh_num_x*mesh_num_y,                                    & ! (in)
        this%ks, this%ke,  this%ls, this%le, lmesh%refElem2D%Nfp, lmesh%NeX, lmesh%NeY, & ! (in)
        this%IntIntrpMat, this%IntXi, this%IntWeight, this%NintGLPt,                    & ! (in)
        this%xmin_gl, this%delx, this%xmax_gl - this%xmin_gl,                           & ! (in)
        this%ymin_gl, this%dely, this%ymax_gl - this%ymin_gl                            ) ! (in)
    end select
    
    return
  end subroutine meshfield_spetraltransform2d_transform
 
!--- Private subroutines (1D)
 
!OCL SERIAL
  subroutine spectral_transform1d_dft( spectral_coef, &
    q_list, var_num, mesh_num, ks, ke, Np, Ne, NsamplePerElem, &
    fft, FFTIntrpMat, FFT_xi, &
    xmin_gl, delx, Lx )
 
    use mpi
    use scale_prc, only: &
      prc_local_comm_world    
    implicit none
    integer, intent(in) :: ks, ke
    integer, intent(in) :: np
    integer, intent(in) :: ne
    integer, intent(in) :: nsampleperelem
    integer, intent(in) :: var_num
    integer, intent(in) :: mesh_num
    real(rp), intent(out) :: spectral_coef(ks:ke,2,var_num)
    type(meshfield1dlist), target :: q_list(var_num,mesh_num)
    class(fastfouriertransform1d), intent(in) :: fft
    real(rp), intent(in) :: fftintrpmat(nsampleperelem,np)
    real(rp), intent(in) :: fft_xi(nsampleperelem)
    real(rp), intent(in) :: xmin_gl
    real(rp), intent(in) :: delx
    real(rp), intent(in) :: lx
 
    real(rp) :: g_q(nsampleperelem,ne,mesh_num,var_num)
    real(rp) :: q_tmp(np)
    complex(RP) :: s_q(nsampleperelem*ne*mesh_num,var_num)
 
    integer :: m
    integer :: v
    integer :: kel
 
    integer :: nall
    integer :: kk
 
    real(rp) :: x_(nsampleperelem)
    !-----------------------------------------------
 
    !$omp parallel do collapse(3) private(v,m,kel,q_tmp, x_)
    do v=1, var_num
    do m=1, mesh_num
    do kel=1, ne
      q_tmp(:) = q_list(v,m)%ptr%local(1)%val(:,kel)
      g_q(:,kel,m,v) = matmul( fftintrpmat, q_tmp )
    end do
    end do
    end do
 
    !$omp parallel do
    do v=1, var_num
      call fft%Forward_real( g_q(:,:,:,v), s_q(:,v) )
    end do
 
    nall = nsampleperelem * ne * mesh_num
    !$omp parallel do private(kk)
    do v=1, var_num
      do kk=1, nall/2+1
        spectral_coef(kk-1,1,v) = real(s_q(kk,v))
        spectral_coef(kk-1,2,v) = aimag(s_q(kk,v))
      end do
      do kk=nall/2+1, nall
        spectral_coef(kk-1-nall,1,v) = real(s_q(kk,v))
        spectral_coef(kk-1-nall,2,v) = aimag(s_q(kk,v))
      end do
    end do
    return
  end subroutine spectral_transform1d_dft
  
!OCL SERIAL
  subroutine spectral_transform1d_l2projection( spectral_coef, &
    q_list, var_num, mesh_num, ks, ke, Np, Ne, IntIntrpMat, IntXi, IntWeight, NintGLPt, &
    xmin_gl, delx, Lx )
 
    use mpi
    use scale_prc, only: &
      prc_local_comm_world 
    implicit none
    integer, intent(in) :: ks, ke
    integer, intent(in) :: np
    integer, intent(in) :: ne
    integer, intent(in) :: var_num
    integer, intent(in) :: mesh_num
    class(meshfield1dlist), target :: q_list(var_num,mesh_num)
    real(rp), intent(out) :: spectral_coef(ks:ke,2,var_num)
    integer, intent(in) :: nintglpt
    real(rp), intent(in) :: intintrpmat(nintglpt,np)
    real(rp), intent(in) :: intxi(nintglpt)
    real(rp), intent(in) :: intweight(nintglpt)
    real(rp), intent(in) :: xmin_gl
    real(rp), intent(in) :: delx
    real(rp), intent(in) :: lx
 
    real(rp), allocatable  :: s_coef_lc(:,:,:)
    real(rp), allocatable  :: s_coef(:,:,:)
    real(rp), allocatable :: q_tmp(:,:,:)
 
    integer :: vec_size
    integer :: km
 
    class(meshbase1d), pointer :: mesh1d
    class(localmesh1d), pointer :: lmesh
    integer :: ldom
    integer :: meshid
    type(localmeshfieldbaselist) :: lcfields(var_num)
 
    integer :: kel, v
    integer :: k, kk
    integer :: ierr
    !-----------------------------------------------
 
    vec_size = var_num
    km = ke
 
    allocate( s_coef_lc(vec_size,2,0:km) )
    allocate( s_coef(vec_size,2,0:km) )
    
    !-
    mesh1d => q_list(1,1)%ptr%mesh
    lmesh => mesh1d%lcmesh_list(1)
 
    !-
 
    allocate( q_tmp(np,vec_size,ne) )
 
    s_coef_lc(:,:,:) = 0.0_rp
    
    do meshid=1, mesh_num
      mesh1d => q_list(1,meshid)%ptr%mesh
 
      do ldom=1, mesh1d%LOCAL_MESH_NUM
        do v=1, var_num
          lcfields(v)%ptr => q_list(v,meshid)%ptr%local(ldom)
        end do
        !$omp parallel do collapse(2)
        do kel=lmesh%NeS, lmesh%NeE
        do v=1, vec_size
          q_tmp(:,v,kel) = lcfields(v)%ptr%val(:,kel)
        end do
        end do      
        call spectral_transform1d_l2projection_lc( s_coef_lc, &
          q_tmp, 0, km, np, ne, vec_size, intintrpmat, intxi, intweight, nintglpt, &
          (lmesh%xmin - xmin_gl)/lx-0.5_rp, delx/lx  )
      end do
    end do
 
    ! global sum
    call mpi_allreduce( s_coef_lc, s_coef, vec_size * (km+1) * 2, &
      mpi_double_precision, mpi_sum, prc_local_comm_world, ierr   )
 
    !-
    !$omp parallel do private(kk)
    do v=1, var_num
    do k=ks, ke
      kk = abs(k)
      spectral_coef(k,1,v) = s_coef(v,1,kk)
      spectral_coef(k,2,v) = s_coef(v,2,kk)
    end do
    end do
 
    return
  end subroutine spectral_transform1d_l2projection
!OCL SERIAL
  subroutine spectral_transform1d_l2projection_lc( spectral_coef, &
    q, ks, ke, Np, Ne, vec_size, IntIntrpMat, IntXi, IntWeight, NintGLPt, &
    xmin_lc, delx  )
    implicit none
    integer, intent(in) :: ks, ke
    integer, intent(in) :: np
    integer, intent(in) :: ne
    integer, intent(in) :: vec_size
    real(rp), intent(in) :: q(np,vec_size,ne)
    real(rp), intent(inout) :: spectral_coef(vec_size,2,ks:ke)
    integer, intent(in) :: nintglpt
    real(rp), intent(in) :: intintrpmat(nintglpt,np)
    real(rp), intent(in) :: intxi(nintglpt)
    real(rp), intent(in) :: intweight(nintglpt)
    real(rp), intent(in) :: xmin_lc
    real(rp), intent(in) :: delx
 
    integer :: k
    integer :: kel
    integer :: v
 
    real(rp) :: q_intrp(nintglpt,vec_size,ne)
    real(rp) :: int_x(nintglpt,ne)
 
    real(rp) :: phi(nintglpt)
    real(rp) :: cos_kx(nintglpt)
    real(rp) :: sin_kx(nintglpt)
 
    real(rp) :: int_w(nintglpt)
    !-------------------------------
 
    int_w(:) = 0.5_rp * intweight(:) * delx
 
    !$omp parallel private(kel,k,v,cos_kx,sin_kx,phi)
    !$omp do
    do kel=1, ne
      int_x(:,kel) = 2.0_rp * pi * ( xmin_lc + delx * ( ( kel - 1.0_rp ) + 0.5_rp * ( 1.0_rp + intxi(:) ) ) )
      q_intrp(:,:,kel) = matmul(intintrpmat, q(:,:,kel))
    end do
    !$omp do
    do k=ks, ke
      do kel=1, ne
        phi(:) = mod( k * int_x(:,kel), 2.0_rp * pi )
        cos_kx(:) = cos( phi )
        sin_kx(:) = sin( phi )
 
        do v=1, vec_size
          spectral_coef(v,1,k) = spectral_coef(v,1,k) + sum(int_w(:) * cos_kx(:) * q_intrp(:,v,kel))
          spectral_coef(v,2,k) = spectral_coef(v,2,k) - sum(int_w(:) * sin_kx(:) * q_intrp(:,v,kel))
        end do
      end do
    end do
    !$omp end parallel
    return
  end subroutine spectral_transform1d_l2projection_lc
 
!--- Private subroutines (2D)
 
!OCL SERIAL
  subroutine spectral_transform2d_dft( spectral_coef, &
    q_list, var_num, mesh_num_x, mesh_num_y, ks, ke, ls, le, Np1D, NeX, NeY, NeGY, NprcY, NsamplePerElem1D, &
    fft_x, fft_y, FFTIntrpMat )
 
    use mpi
    use scale_prc, only: &
      prc_local_comm_world    
    implicit none
    integer, intent(in) :: ks, ke
    integer, intent(in) :: ls, le
    integer, intent(in) :: np1d
    integer, intent(in) :: nex
    integer, intent(in) :: ney, negy, nprcy
    integer, intent(in) :: nsampleperelem1d
    integer, intent(in) :: var_num
    integer, intent(in) :: mesh_num_x
    integer, intent(in) :: mesh_num_y
    real(rp), intent(out) :: spectral_coef(ks:ke,ls:le,2,var_num)
    class(meshfield2dlist), target :: q_list(var_num,mesh_num_x,mesh_num_y)
    class(fastfouriertransform1d), intent(in) :: fft_x
    class(fastfouriertransform1d), intent(in) :: fft_y
    real(rp), intent(in) :: fftintrpmat(nsampleperelem1d,np1d)
 
    real(rp) :: fftintrpmat_tr(np1d,nsampleperelem1d)
 
    real(rp) :: g_q(nsampleperelem1d*nex*mesh_num_x,nsampleperelem1d*ney*mesh_num_y*var_num)
    complex(RP) :: s_qx(nsampleperelem1d*nex*mesh_num_x,size(g_q,2))
    complex(RP) :: g_qy(nsampleperelem1d*ney*mesh_num_y*nprcy,var_num,nsampleperelem1d*nex*mesh_num_x/nprcy)
    complex(RP) :: s_q_lc(nsampleperelem1d*negy,var_num,nsampleperelem1d*nex*mesh_num_x/nprcy)
    complex(RP) :: s_q(nsampleperelem1d*negy,var_num,nsampleperelem1d*nex*mesh_num_x)
 
    integer :: v
    integer :: xdim
    integer :: prcy
    integer :: j
 
    integer :: nall_x, nall_y
    integer :: kk, ll, kk_os
 
    integer :: sendcount, recvcount
    complex(RP) :: sendbuf(nsampleperelem1d*ney*mesh_num_y*var_num,nsampleperelem1d*nex*mesh_num_x)
    complex(RP) :: recvbuf(nsampleperelem1d*ney*mesh_num_y,var_num,nsampleperelem1d*nex*mesh_num_x/nprcy,nprcy)
    integer :: nx, local_nx
    integer :: ny, local_ny
    integer :: ierr
    !-----------------------------------------------
 
    !-
    nall_x = nsampleperelem1d * nex * mesh_num_x
    nall_y = nsampleperelem1d * negy
 
    !-
    fftintrpmat_tr(:,:) = transpose(fftintrpmat)
    call spectral_transform2d_dft_sampling( g_q, &
      q_list, var_num, mesh_num_x, mesh_num_y, np1d, nex, ney, nsampleperelem1d, &
      fftintrpmat_tr )
   
    !-
    !$omp parallel do
    do v=1, size(g_q,2)
      call fft_x%Forward( g_q(:,v), s_qx(:,v) )
    end do
 
    !-
    local_nx = nall_x
    nx = local_nx
    
    local_ny = nsampleperelem1d * ney * mesh_num_y * var_num
    ny = nsampleperelem1d * negy * var_num
    
    sendbuf(:,:) = transpose(s_qx)
    sendcount = local_ny * (nx / nprcy)
    recvcount = sendcount
    call mpi_alltoall( sendbuf, sendcount, mpi_double_complex, &
                       recvbuf, recvcount, mpi_double_complex, prc_local_comm_world, ierr )
 
    local_ny = nsampleperelem1d * ney * mesh_num_y                     
    do v=1, var_num
    do prcy=1, nprcy
    do xdim=1, nsampleperelem1d*nex*mesh_num_x/nprcy
    do j=1, local_ny
      g_qy(j+(prcy-1)*local_ny,v,xdim) = recvbuf(j,v,xdim,prcy)
    end do
    end do
    end do
    end do
 
    !-
    !$omp parallel do collapse(2)
    do v=1, var_num
    do xdim=1, nx/nprcy
      call fft_y%Forward( g_qy(:,v,xdim), s_q_lc(:,v,xdim) )
    end do
    end do
 
    sendcount = size(s_q_lc)
    recvcount = size(s_q_lc)
    call mpi_allgather( s_q_lc, sendcount, mpi_double_complex, &
                           s_q, recvcount, mpi_double_complex, prc_local_comm_world, ierr )
 
    !-
    !$omp parallel do private(v,kk,ll)
    do v=1, var_num
      do ll=1, nall_y/2+1
        do kk=1, nall_x/2+1
          spectral_coef(kk-1,ll-1,1,v) = real(s_q(ll,v,kk))
          spectral_coef(kk-1,ll-1,2,v) = aimag(s_q(ll,v,kk))
        end do
        do kk=nall_x/2+1, nall_x
          spectral_coef(kk-1-nall_x,ll-1,1,v) = real(s_q(ll,v,kk))
          spectral_coef(kk-1-nall_x,ll-1,2,v) = aimag(s_q(ll,v,kk))
        end do
      end do
      do ll=nall_y/2+1, nall_y
        do kk=1, nall_x/2+1
          spectral_coef(kk-1,ll-1-nall_y,1,v) = real(s_q(ll,v,kk))
          spectral_coef(kk-1,ll-1-nall_y,2,v) = aimag(s_q(ll,v,kk))
        end do
        do kk=nall_x/2+1, nall_x
          spectral_coef(kk-1-nall_x,ll-1-nall_y,1,v) = real(s_q(ll,v,kk))
          spectral_coef(kk-1-nall_x,ll-1-nall_y,2,v) = aimag(s_q(ll,v,kk))
        end do
      end do
    end do
    return
  end subroutine spectral_transform2d_dft
 
!OCL SERIAL
  subroutine spectral_transform2d_dft_sampling( g_q, &
    q_list, var_num, mesh_num_x, mesh_num_y, Np1D, NeX, NeY, NsamplePerElem1D, &
    FFTIntrpMat_tr )
 
    use mpi
    use scale_prc, only: &
      prc_local_comm_world    
    implicit none
    integer, intent(in) :: np1d
    integer, intent(in) :: nex, ney
    integer, intent(in) :: nsampleperelem1d
    integer, intent(in) :: var_num
    integer, intent(in) :: mesh_num_x
    integer, intent(in) :: mesh_num_y
    real(rp), intent(out) :: g_q(nsampleperelem1d,nex,mesh_num_x,nsampleperelem1d,ney,mesh_num_y,var_num)
    class(meshfield2dlist), intent(in) :: q_list(var_num,mesh_num_x,mesh_num_y)
    real(rp), intent(in) :: fftintrpmat_tr(np1d,nsampleperelem1d)
 
    real(rp) :: q_tmp(np1d,np1d)
    real(rp) :: q_intrp1(nsampleperelem1d,np1d)
    real(rp) :: q_intrp2(nsampleperelem1d,nsampleperelem1d)
 
    integer :: mx, my, v
    integer :: kel_x, kel_y, kel
    integer :: p, p1, pp1, p2, pp2
    !-----------------------------------------------
 
    !$omp parallel do private(mx,my,v, kel_x,kel_y,kel, p,p1,pp1,p2,pp2, q_tmp,q_intrp1,q_intrp2) collapse(5)
    do v=1, var_num
    do my=1, mesh_num_y
    do mx=1, mesh_num_x
    do kel_y=1, ney
    do kel_x=1, nex
      kel = kel_x + (kel_y-1)*nex
      do p2=1, np1d
      do p1=1, np1d
        p = p1 + (p2-1)*np1d
        q_tmp(p1,p2) = q_list(v,mx,my)%ptr%local(1)%val(p,kel)
      end do
      end do
 
      q_intrp1(:,:) = 0.0_rp
      do p2=1, np1d
      do p1=1, nsampleperelem1d
      do pp1=1, np1d
        q_intrp1(p1,p2) = q_intrp1(p1,p2) + fftintrpmat_tr(pp1,p1) * q_tmp(pp1,p2)
      end do
      end do
      end do      
      q_intrp2(:,:) = 0.0_rp
      do p2=1, nsampleperelem1d
      do pp2=1, np1d
      do p1=1, nsampleperelem1d
        q_intrp2(p1,p2) = q_intrp2(p1,p2) + fftintrpmat_tr(pp2,p2) * q_intrp1(p1,pp2)
      end do
      end do
      end do      
 
      do p2=1, nsampleperelem1d
        g_q(:,kel_x,mx,p2,kel_y,my,v) = q_intrp2(:,p2)
      end do
    end do
    end do
    end do
    end do
    end do
    return
  end subroutine spectral_transform2d_dft_sampling
 
!OCL SERIAL
  subroutine spectral_transform2d_l2projection( spectral_coef, &
    q_list, var_num, mesh_num, ks, ke, ls, le, Np1D, NeX, NeY, &
    IntIntrpMat1D, IntXi1D, IntWeight1D, NintGLPt1D,           &
    xmin_gl, delx, Lx, ymin_gl, dely, Ly )
 
    use mpi
    use scale_prc, only: &
      prc_local_comm_world
    implicit none
    integer, intent(in) :: ks, ke
    integer, intent(in) :: ls, le
    integer, intent(in) :: np1d
    integer, intent(in) :: nex, ney
    integer, intent(in) :: var_num
    integer, intent(in) :: mesh_num
    type(meshfield2dlist), target :: q_list(var_num,mesh_num)
    real(rp), intent(out) :: spectral_coef(ks:ke,ls:le,2,var_num)
    integer, intent(in) :: nintglpt1d
    real(rp), intent(in) :: intintrpmat1d(nintglpt1d,np1d)
    real(rp), intent(in) :: intxi1d(nintglpt1d)
    real(rp), intent(in) :: intweight1d(nintglpt1d)
    real(rp), intent(in) :: xmin_gl
    real(rp), intent(in) :: delx
    real(rp), intent(in) :: lx
    real(rp), intent(in) :: ymin_gl
    real(rp), intent(in) :: dely
    real(rp), intent(in) :: ly
 
    real(rp), allocatable  :: s_coef_lc(:,:,:,:)
    real(rp), allocatable  :: s_coef(:,:,:,:)
    real(rp), allocatable :: q_tmp(:,:,:)
 
    real(rp) :: intintrpmat1d_tr(np1d,nintglpt1d)
 
    integer :: vec_size
    integer :: km
 
    class(meshbase2d), pointer :: mesh2d
    class(localmesh2d), pointer :: lmesh
    integer :: ldom
    integer :: meshid
 
    integer :: kel, v
    integer :: k, kk, l
    integer :: ierr
 
    type(localmeshfieldbaselist) :: lcfields(var_num)
    !-----------------------------------------------
 
    vec_size = var_num
    km = ke
 
    allocate( s_coef_lc(vec_size,2,0:km,ls:le) )
    allocate( s_coef(vec_size,2,0:km,ls:le) )
    
    !-
    mesh2d => q_list(1,1)%ptr%mesh
    lmesh => mesh2d%lcmesh_list(1)
 
    intintrpmat1d_tr(:,:) = transpose(intintrpmat1d)
    !-
 
    allocate( q_tmp(np1d**2,vec_size,nex*ney) )
 
    s_coef_lc(:,:,:,:) = 0.0_rp
    
    do meshid=1, mesh_num
      mesh2d => q_list(1,meshid)%ptr%mesh
      do ldom=1, mesh2d%LOCAL_MESH_NUM
        lmesh => mesh2d%lcmesh_list(ldom)
 
        do v=1, var_num
          lcfields(v)%ptr => q_list(v,meshid)%ptr%local(ldom)
        end do
        !$omp parallel do collapse(2)
        do kel=lmesh%NeS, lmesh%NeE
        do v=1, vec_size
          q_tmp(:,v,kel) = lcfields(v)%ptr%val(:,kel)
        end do
        end do      
        call spectral_transform2d_l2projection_lc( s_coef_lc, &
          q_tmp, 0, km, ls, le, np1d, nex, ney, vec_size,     &
          intintrpmat1d_tr, intxi1d, intweight1d, nintglpt1d, &
          (lmesh%xmin - xmin_gl)/lx-0.5_rp, delx/lx,          &
          (lmesh%ymin - ymin_gl)/ly-0.5_rp, dely/ly           )
      end do
    end do
 
    ! global sum
    call mpi_allreduce( s_coef_lc, s_coef, vec_size * (km+1)*(le-ls+1) * 2, &
      mpi_double_precision, mpi_sum, prc_local_comm_world, ierr   )
 
    !-
    !$omp parallel do collapse(2) private(kk)
    do v=1, var_num
    do l=ls, le
      do k=ks, ke
        kk = abs(k)
        if ( k >= 0 ) then
          spectral_coef(k,l,1,v) = s_coef(v,1,kk,l)
          spectral_coef(k,l,2,v) = s_coef(v,2,kk,l)
        else
          spectral_coef(k,l,1,v) = s_coef(v,1,kk,-l)
          spectral_coef(k,l,2,v) = - s_coef(v,2,kk,-l)
        end if
      end do
    end do
    end do
 
    return
  end subroutine spectral_transform2d_l2projection
!OCL SERIAL
  subroutine spectral_transform2d_l2projection_lc( spectral_coef, &
    q, ks, ke, ls, le, Np1D, NeX, NeY, vec_size, IntIntrpMat1D_tr, IntXi1D, IntWeight1D, NintGLpt1D, &
    xmin_lc, delx, ymin_lc, dely  )
    implicit none
    integer, intent(in) :: ks, ke
    integer, intent(in) :: ls, le
    integer, intent(in) :: np1d
    integer, intent(in) :: nex
    integer, intent(in) :: ney
    integer, intent(in) :: vec_size
    real(rp), intent(in) :: q(np1d,np1d,vec_size,nex,ney)
    real(rp), intent(inout) :: spectral_coef(vec_size,2,ks:ke,ls:le)
    integer, intent(in) :: nintglpt1d
    real(rp), intent(in) :: intintrpmat1d_tr(np1d,nintglpt1d)
    real(rp), intent(in) :: intxi1d(nintglpt1d)
    real(rp), intent(in) :: intweight1d(nintglpt1d)
    real(rp), intent(in) :: xmin_lc
    real(rp), intent(in) :: delx
    real(rp), intent(in) :: ymin_lc
    real(rp), intent(in) :: dely
 
    integer :: k, l
    integer :: kel_x, kel_y
    integer :: p, px, py
    integer :: v
 
    real(rp) :: q_intrp_tmp(nintglpt1d,np1d,vec_size)
    real(rp) :: q_intrp_tmp2(nintglpt1d**2,vec_size)
    real(rp) :: q_intrp(nintglpt1d**2,vec_size,nex,ney)
    real(rp) :: int_x(nintglpt1d,nex)
    real(rp) :: int_y(nintglpt1d,ney)
 
    real(rp) :: phi(nintglpt1d**2)
    real(rp) :: cos_phi(nintglpt1d**2)
    real(rp) :: sin_phi(nintglpt1d**2)
 
    real(rp) :: int_w(nintglpt1d**2)
    !-------------------------------
 
    !$omp parallel private(kel_x,kel_y,k,l,v,p,px,py,cos_phi,sin_phi,phi,q_intrp_tmp,q_intrp_tmp2)
    !$omp do
    do py=1, nintglpt1d
    do px=1, nintglpt1d
      p = px + (py-1)*nintglpt1d
      int_w(p) = 0.25_rp * delx * dely * intweight1d(px) * intweight1d(py)
    end do
    end do
    !$omp do
    do kel_x=1, nex
      int_x(:,kel_x) = 2.0_rp * pi * ( xmin_lc + delx * ( ( kel_x - 1.0_rp ) + 0.5_rp * ( 1.0_rp + intxi1d(:) ) ) )
    end do
    !$omp do
    do kel_y=1, ney
      int_y(:,kel_y) = 2.0_rp * pi * ( ymin_lc + dely * ( ( kel_y - 1.0_rp ) + 0.5_rp * ( 1.0_rp + intxi1d(:) ) ) )
    end do
    !$omp do collapse(2)
    do kel_y=1, ney
    do kel_x=1, nex
      q_intrp_tmp(:,:,:) = 0.0_rp
      do v=1, vec_size
      do py=1, np1d
      do px=1, nintglpt1d
        q_intrp_tmp(px,py,v) = q_intrp_tmp(px,py,v) &
          + sum(intintrpmat1d_tr(:,px) * q(:,py,v,kel_x,kel_y))
      end do
      end do
      end do
 
      q_intrp_tmp2(:,:) = 0.0_rp
      do py=1, nintglpt1d
      do v=1, vec_size
      do px=1, nintglpt1d
        p = px + (py-1)*nintglpt1d
        q_intrp_tmp2(p,v) = q_intrp_tmp2(p,v) &
          + sum(intintrpmat1d_tr(:,py) * q_intrp_tmp(px,:,v))
      end do
      end do
      end do
      q_intrp(:,:,kel_x,kel_y) = q_intrp_tmp2(:,:)
    end do
    end do
    !$omp do collapse(2)
    do l=ls, le
    do k=ks, ke
      do kel_y=1, ney
      do kel_x=1, nex
        do py=1, nintglpt1d
        do px=1, nintglpt1d
          p = px + (py-1)*nintglpt1d
          phi(p) = mod( k * int_x(px,kel_x) + l * int_y(py,kel_y), 2.0_rp * pi )
        end do
        end do
        cos_phi(:) = cos( phi(:) )
        sin_phi(:) = sin( phi(:) )
 
        do v=1, vec_size
          spectral_coef(v,1,k,l) = spectral_coef(v,1,k,l) + sum(int_w(:) * cos_phi(:) * q_intrp(:,v,kel_x,kel_y))
          spectral_coef(v,2,k,l) = spectral_coef(v,2,k,l) - sum(int_w(:) * sin_phi(:) * q_intrp(:,v,kel_x,kel_y))
        end do
      end do
      end do
    end do
    end do
    !$omp end parallel
    return
  end subroutine spectral_transform2d_l2projection_lc
 
end module scale_meshfield_spectral_transform