coef.f90

Go to the documentation of this file.

! Copyright (c) 2020-2026, 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 coefs
  use gather_scatter, only : gs_t
  use gs_ops, only : gs_op_add
  use neko_config, only : neko_bcknd_device, neko_bcknd_opencl
  use num_types, only : rp
  use dofmap, only : dofmap_t
  use space, only : space_t
  use math, only : rone, invcol1, addcol3, subcol3, copy, &
       chsign, rzero, invers2, glsum, neko_eps
  use mesh, only : mesh_t
  use device_math, only : device_rone, device_invcol1, &
       device_glsum, device_copy
  use device_coef, only : device_coef_generate_geo, &
       device_coef_generate_dxydrst, device_coef_generate_mass, &
       device_coef_generate_area_and_normal
  use mxm_wrapper, only : mxm
  use device, only : device_map, device_memcpy, device_to_host, &
       host_to_device, device_unmap
  use utils, only : index_is_on_facet, linear_index, &
       neko_error
  use comm, only : neko_comm
  use neko_config, only : neko_bcknd_device
  use mpi_f08, only : mpi_allreduce, mpi_integer, mpi_sum
  use, intrinsic :: iso_c_binding
  implicit none
  private
 
  type, public :: coef_t
     real(kind=rp), allocatable :: g11(:,:,:,:)
     real(kind=rp), allocatable :: g22(:,:,:,:)
     real(kind=rp), allocatable :: g33(:,:,:,:)
     real(kind=rp), allocatable :: g12(:,:,:,:)
     real(kind=rp), allocatable :: g13(:,:,:,:)
     real(kind=rp), allocatable :: g23(:,:,:,:)
 
     real(kind=rp), allocatable :: g11_compressed(:,:,:,:)
     real(kind=rp), allocatable :: g22_compressed(:,:,:,:)
     real(kind=rp), allocatable :: g33_compressed(:,:,:,:)
     real(kind=rp), allocatable :: g12_compressed(:,:,:,:)
     real(kind=rp), allocatable :: g13_compressed(:,:,:,:)
     real(kind=rp), allocatable :: g23_compressed(:,:,:,:)
     integer, allocatable :: compression_inds(:)
 
     real(kind=rp), allocatable :: mult(:,:,:,:) 
     real(kind=rp), allocatable :: dxdr(:,:,:,:), dydr(:,:,:,:), dzdr(:,:,:,:)
     real(kind=rp), allocatable :: dxds(:,:,:,:), dyds(:,:,:,:), dzds(:,:,:,:)
     real(kind=rp), allocatable :: dxdt(:,:,:,:), dydt(:,:,:,:), dzdt(:,:,:,:)
     real(kind=rp), allocatable :: drdx(:,:,:,:), drdy(:,:,:,:), drdz(:,:,:,:)
     real(kind=rp), allocatable :: dsdx(:,:,:,:), dsdy(:,:,:,:), dsdz(:,:,:,:)
     real(kind=rp), allocatable :: dtdx(:,:,:,:), dtdy(:,:,:,:), dtdz(:,:,:,:)
 
     real(kind=rp), allocatable :: h1(:,:,:,:) 
     real(kind=rp), allocatable :: h2(:,:,:,:) 
     logical :: ifh2
 
     real(kind=rp), allocatable :: jac(:,:,:,:) 
     real(kind=rp), allocatable :: jacinv(:,:,:,:) 
     real(kind=rp), allocatable :: b(:,:,:,:) 
     real(kind=rp), allocatable :: binv(:,:,:,:) 
     real(kind=rp), pointer :: blag(:,:,:,:) => null() 
     real(kind=rp), pointer :: blaglag(:,:,:,:) => null() 
     real(kind=rp), allocatable :: area(:,:,:,:) 
     real(kind=rp), allocatable :: nx(:,:,:,:) 
     real(kind=rp), allocatable :: ny(:,:,:,:) 
     real(kind=rp), allocatable :: nz(:,:,:,:) 
     logical :: cyclic = .false.
     integer, allocatable :: cyc_msk(:)
     real(kind=rp), allocatable :: r11(:) 
     real(kind=rp), allocatable :: r12(:) 
 
     !! True if geometric metrics have been initialized
     logical, private :: coef_metrics_initialized = .false.
 
 
     real(kind=rp) :: volume
 
     type(space_t), pointer :: xh => null()
     type(mesh_t), pointer :: msh => null()
     type(dofmap_t), pointer :: dof => null()
     type(gs_t), pointer :: gs_h=> null()
 
     !
     ! Device pointers (if present)
     !
 
     type(c_ptr) :: g11_d = c_null_ptr
     type(c_ptr) :: g22_d = c_null_ptr
     type(c_ptr) :: g33_d = c_null_ptr
     type(c_ptr) :: g12_d = c_null_ptr
     type(c_ptr) :: g13_d = c_null_ptr
     type(c_ptr) :: g23_d = c_null_ptr
     type(c_ptr) :: dxdr_d = c_null_ptr
     type(c_ptr) :: dydr_d = c_null_ptr
     type(c_ptr) :: dzdr_d = c_null_ptr
     type(c_ptr) :: dxds_d = c_null_ptr
     type(c_ptr) :: dyds_d = c_null_ptr
     type(c_ptr) :: dzds_d = c_null_ptr
     type(c_ptr) :: dxdt_d = c_null_ptr
     type(c_ptr) :: dydt_d = c_null_ptr
     type(c_ptr) :: dzdt_d = c_null_ptr
     type(c_ptr) :: drdx_d = c_null_ptr
     type(c_ptr) :: drdy_d = c_null_ptr
     type(c_ptr) :: drdz_d = c_null_ptr
     type(c_ptr) :: dsdx_d = c_null_ptr
     type(c_ptr) :: dsdy_d = c_null_ptr
     type(c_ptr) :: dsdz_d = c_null_ptr
     type(c_ptr) :: dtdx_d = c_null_ptr
     type(c_ptr) :: dtdy_d = c_null_ptr
     type(c_ptr) :: dtdz_d = c_null_ptr
     type(c_ptr) :: mult_d = c_null_ptr
     type(c_ptr) :: h1_d = c_null_ptr
     type(c_ptr) :: h2_d = c_null_ptr
     type(c_ptr) :: jac_d = c_null_ptr
     type(c_ptr) :: jacinv_d = c_null_ptr
     type(c_ptr) :: b_d = c_null_ptr
     type(c_ptr) :: blag_d = c_null_ptr
     type(c_ptr) :: blaglag_d = c_null_ptr
     type(c_ptr) :: binv_d = c_null_ptr
     type(c_ptr) :: area_d = c_null_ptr
     type(c_ptr) :: nx_d = c_null_ptr
     type(c_ptr) :: ny_d = c_null_ptr
     type(c_ptr) :: nz_d = c_null_ptr
     type(c_ptr) :: cyc_msk_d = c_null_ptr
     type(c_ptr) :: r11_d = c_null_ptr
     type(c_ptr) :: r12_d = c_null_ptr
 
 
 
   contains
     procedure, private, pass(this) :: init_empty => coef_init_empty
     procedure, private, pass(this) :: init_all => coef_init_all
     procedure, pass(this) :: free => coef_free
     procedure, pass(this) :: get_normal => coef_get_normal
     procedure, pass(this) :: get_area => coef_get_area
     procedure, pass(this) :: generate_cyclic_bc => coef_generate_cyclic_bc
     procedure, pass(this) :: recompute_metrics => coef_recompute_metrics
     procedure, pass(this) :: enable_b_history => coef_enable_lagged_mass
     procedure, pass(this) :: update_b_history => coef_update_lagged_mass
     generic :: init => init_empty, init_all
  end type coef_t
 
contains
 
  subroutine coef_init_empty(this, Xh, msh)
    class(coef_t), intent(inout) :: this
    type(space_t), intent(inout), target :: Xh
    type(mesh_t), intent(inout), target :: msh
    integer :: n
    call this%free()
    this%msh => msh
    this%Xh => xh
 
    allocate(this%drdx(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
    allocate(this%dsdx(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
    allocate(this%dtdx(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
 
    allocate(this%drdy(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
    allocate(this%dsdy(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
    allocate(this%dtdy(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
 
    allocate(this%drdz(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
    allocate(this%dsdz(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
    allocate(this%dtdz(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
 
 
    !
    ! Setup device memory (if present)
    !
 
    n = this%Xh%lx * this%Xh%ly * this%Xh%lz * this%msh%nelv
    if (neko_bcknd_device .eq. 1) then
 
       call device_map(this%drdx, this%drdx_d, n)
       call device_map(this%drdy, this%drdy_d, n)
       call device_map(this%drdz, this%drdz_d, n)
 
       call device_map(this%dsdx, this%dsdx_d, n)
       call device_map(this%dsdy, this%dsdy_d, n)
       call device_map(this%dsdz, this%dsdz_d, n)
 
       call device_map(this%dtdx, this%dtdx_d, n)
       call device_map(this%dtdy, this%dtdy_d, n)
       call device_map(this%dtdz, this%dtdz_d, n)
 
    end if
 
  end subroutine coef_init_empty
 
  subroutine coef_init_all(this, gs_h)
    class(coef_t), intent(inout), target :: this
    type(gs_t), intent(inout), target :: gs_h
    integer :: n, m, ncyc
    call this%free()
 
    this%msh => gs_h%dofmap%msh
    this%Xh => gs_h%dofmap%Xh
    this%dof => gs_h%dofmap
    this%gs_h => gs_h
 
    !
    ! Allocate arrays for geometric data
    !
    allocate(this%G11(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
    allocate(this%G22(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
    allocate(this%G33(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
    allocate(this%G12(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
    allocate(this%G13(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
    allocate(this%G23(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
 
    allocate(this%dxdr(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
    allocate(this%dxds(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
    allocate(this%dxdt(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
 
    allocate(this%dydr(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
    allocate(this%dyds(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
    allocate(this%dydt(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
 
    allocate(this%dzdr(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
    allocate(this%dzds(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
    allocate(this%dzdt(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
 
    allocate(this%drdx(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
    allocate(this%dsdx(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
    allocate(this%dtdx(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
 
    allocate(this%drdy(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
    allocate(this%dsdy(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
    allocate(this%dtdy(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
 
    allocate(this%drdz(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
    allocate(this%dsdz(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
    allocate(this%dtdz(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
 
    allocate(this%jac(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
    allocate(this%jacinv(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
 
    allocate(this%area(this%Xh%lx, this%Xh%ly, 6, this%msh%nelv))
    allocate(this%nx(this%Xh%lx, this%Xh%ly, 6, this%msh%nelv))
    allocate(this%ny(this%Xh%lx, this%Xh%ly, 6, this%msh%nelv))
    allocate(this%nz(this%Xh%lx, this%Xh%ly, 6, this%msh%nelv))
 
    allocate(this%B(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
    allocate(this%Binv(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
 
    ! We do this so in a static simulation we don't allocate extra memory
    this%Blag => this%B
    this%Blaglag => this%B
 
    allocate(this%h1(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
    allocate(this%h2(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
 
    allocate(this%mult(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
 
 
    !
    ! Setup device memory (if present)
    !
 
    n = this%Xh%lx * this%Xh%ly * this%Xh%lz * this%msh%nelv
    if (neko_bcknd_device .eq. 1) then
       call device_map(this%G11, this%G11_d, n)
       call device_map(this%G22, this%G22_d, n)
       call device_map(this%G33, this%G33_d, n)
       call device_map(this%G12, this%G12_d, n)
       call device_map(this%G13, this%G13_d, n)
       call device_map(this%G23, this%G23_d, n)
 
       call device_map(this%dxdr, this%dxdr_d, n)
       call device_map(this%dydr, this%dydr_d, n)
       call device_map(this%dzdr, this%dzdr_d, n)
 
       call device_map(this%dxds, this%dxds_d, n)
       call device_map(this%dyds, this%dyds_d, n)
       call device_map(this%dzds, this%dzds_d, n)
 
       call device_map(this%dxdt, this%dxdt_d, n)
       call device_map(this%dydt, this%dydt_d, n)
       call device_map(this%dzdt, this%dzdt_d, n)
 
       call device_map(this%drdx, this%drdx_d, n)
       call device_map(this%drdy, this%drdy_d, n)
       call device_map(this%drdz, this%drdz_d, n)
 
       call device_map(this%dsdx, this%dsdx_d, n)
       call device_map(this%dsdy, this%dsdy_d, n)
       call device_map(this%dsdz, this%dsdz_d, n)
 
       call device_map(this%dtdx, this%dtdx_d, n)
       call device_map(this%dtdy, this%dtdy_d, n)
       call device_map(this%dtdz, this%dtdz_d, n)
 
       call device_map(this%mult, this%mult_d, n)
       call device_map(this%h1, this%h1_d, n)
       call device_map(this%h2, this%h2_d, n)
 
       call device_map(this%jac, this%jac_d, n)
       call device_map(this%jacinv, this%jacinv_d, n)
       call device_map(this%B, this%B_d, n)
       call device_map(this%Binv, this%Binv_d, n)
 
       this%Blag_d = this%B_d
       this%Blaglag_d = this%B_d
 
       m = this%Xh%lx * this%Xh%ly * 6 * this%msh%nelv
 
       call device_map(this%area, this%area_d, m)
       call device_map(this%nx, this%nx_d, m)
       call device_map(this%ny, this%ny_d, m)
       call device_map(this%nz, this%nz_d, m)
 
    end if
 
    call coef_generate_dxyzdrst(this)
 
    call coef_generate_geo(this)
 
    ! call coef_generate_geo_compressed(this)
 
    call coef_generate_area_and_normal(this)
 
    call coef_generate_mass(this)
 
    this%coef_metrics_initialized = .true.
 
 
    ! This is a placeholder, just for now
    ! We can probably find a prettier solution
    if (neko_bcknd_device .eq. 1) then
       call device_rone(this%h1_d, n)
       call device_rone(this%h2_d, n)
       call device_memcpy(this%h1, this%h1_d, n, &
            device_to_host, sync = .false.)
       call device_memcpy(this%h2, this%h2_d, n, &
            device_to_host, sync = .false.)
    else
       call rone(this%h1,n)
       call rone(this%h2,n)
    end if
 
    this%ifh2 = .false.
 
    !
    ! Set up multiplicity
    !
    if (neko_bcknd_device .eq. 1) then
       call device_rone(this%mult_d, n)
    else
       call rone(this%mult, n)
    end if
 
    call gs_h%op(this%mult, n, gs_op_add)
 
    if (neko_bcknd_device .eq. 1) then
       call device_invcol1(this%mult_d, n)
       call device_memcpy(this%mult, this%mult_d, n, &
            device_to_host, sync = .true.)
    else
       call invcol1(this%mult, n)
    end if
 
    ncyc = this%msh%periodic%size * this%Xh%lx * this%Xh%lx
    allocate(this%cyc_msk(0:ncyc))
    this%cyc_msk(0) = ncyc + 1
    if (ncyc .gt. 0) then
       allocate(this%R11(ncyc))
       allocate(this%R12(ncyc))
 
       call rone(this%R11, ncyc)
       call rzero(this%R12, ncyc)
 
       if (neko_bcknd_device .eq. 1) then
          call device_map(this%cyc_msk, this%cyc_msk_d, ncyc+1)
          call device_map(this%R11, this%R11_d, ncyc)
          call device_map(this%R12, this%R12_d, ncyc)
 
          call device_memcpy(this%cyc_msk, this%cyc_msk_d, ncyc+1, &
               host_to_device, sync = .false.)
          call device_memcpy(this%R11, this%R11_d, ncyc, &
               host_to_device, sync = .false.)
          call device_memcpy(this%R12, this%R12_d, ncyc, &
               host_to_device, sync = .false.)
       end if
 
    end if
  end subroutine coef_init_all
 
  subroutine coef_free(this)
    class(coef_t), intent(inout), target :: this
 
    if (allocated(this%G11)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%G11, this%G11_d)
       deallocate(this%G11)
    end if
 
    if (allocated(this%G22)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%G22, this%G22_d)
       deallocate(this%G22)
    end if
 
    if (allocated(this%G33)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%G33, this%G33_d)
       deallocate(this%G33)
    end if
 
    if (allocated(this%G12)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%G12, this%G12_d)
       deallocate(this%G12)
    end if
 
    if (allocated(this%G13)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%G13, this%G13_d)
       deallocate(this%G13)
    end if
 
    if (allocated(this%G23)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%G23, this%G23_d)
       deallocate(this%G23)
    end if
 
    if (allocated(this%G11_compressed)) then
       deallocate(this%G11_compressed)
    end if
 
    if (allocated(this%compression_inds)) then
       deallocate(this%compression_inds)
    end if
 
    if (allocated(this%G22_compressed)) then
       deallocate(this%G22_compressed)
    end if
 
    if (allocated(this%G33_compressed)) then
       deallocate(this%G33_compressed)
    end if
 
    if (allocated(this%G12_compressed)) then
       deallocate(this%G12_compressed)
    end if
 
    if (allocated(this%G13_compressed)) then
       deallocate(this%G13_compressed)
    end if
 
    if (allocated(this%G23_compressed)) then
       deallocate(this%G23_compressed)
    end if
 
    if (allocated(this%mult)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%mult, this%mult_d)
       deallocate(this%mult)
    end if
 
    if (associated(this%Blag) .and. &
         .not. associated(this%Blag, this%B)) then
       if (c_associated(this%Blag_d) .and. &
            .not. c_associated(this%Blag_d, this%B_d)) then
          call device_unmap(this%Blag, this%Blag_d)
       end if
       deallocate(this%Blag)
    end if
    nullify(this%Blag)
 
    if (associated(this%Blaglag) .and. &
         .not. associated(this%Blaglag, this%B)) then
       if (c_associated(this%Blaglag_d) .and. &
            .not. c_associated(this%Blaglag_d, this%B_d)) then
          call device_unmap(this%Blaglag, this%Blaglag_d)
       end if
       deallocate(this%Blaglag)
    end if
    nullify(this%Blaglag)
 
    if (c_associated(this%Blag_d) .and. &
         .not. c_associated(this%Blag_d, this%B_d)) then
       this%Blag_d = c_null_ptr
    end if
    this%Blag_d = c_null_ptr
 
    if (c_associated(this%Blaglag_d) .and. &
         .not. c_associated(this%Blaglag_d, this%B_d)) then
       this%Blaglag_d = c_null_ptr
    end if
    this%Blaglag_d = c_null_ptr
 
    if (allocated(this%B)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%B, this%B_d)
       deallocate(this%B)
    end if
 
    if (allocated(this%Binv)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%Binv, this%Binv_d)
       deallocate(this%Binv)
    end if
 
    if (allocated(this%dxdr)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%dxdr, this%dxdr_d)
       deallocate(this%dxdr)
    end if
 
    if (allocated(this%dxds)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%dxds, this%dxds_d)
       deallocate(this%dxds)
    end if
 
    if (allocated(this%dxdt)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%dxdt, this%dxdt_d)
       deallocate(this%dxdt)
    end if
 
    if (allocated(this%dydr)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%dydr, this%dydr_d)
       deallocate(this%dydr)
    end if
 
    if (allocated(this%dyds)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%dyds, this%dyds_d)
       deallocate(this%dyds)
    end if
 
    if (allocated(this%dydt)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%dydt, this%dydt_d)
       deallocate(this%dydt)
    end if
 
    if (allocated(this%dzdr)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%dzdr, this%dzdr_d)
       deallocate(this%dzdr)
    end if
 
    if (allocated(this%dzds)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%dzds, this%dzds_d)
       deallocate(this%dzds)
    end if
 
    if (allocated(this%dzdt)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%dzdt, this%dzdt_d)
       deallocate(this%dzdt)
    end if
 
    if (allocated(this%drdx)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%drdx, this%drdx_d)
       deallocate(this%drdx)
    end if
 
    if (allocated(this%dsdx)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%dsdx, this%dsdx_d)
       deallocate(this%dsdx)
    end if
 
    if (allocated(this%dtdx)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%dtdx, this%dtdx_d)
       deallocate(this%dtdx)
    end if
 
    if (allocated(this%drdy)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%drdy, this%drdy_d)
       deallocate(this%drdy)
    end if
 
    if (allocated(this%dsdy)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%dsdy, this%dsdy_d)
       deallocate(this%dsdy)
    end if
 
    if (allocated(this%dtdy)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%dtdy, this%dtdy_d)
       deallocate(this%dtdy)
    end if
 
    if (allocated(this%drdz)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%drdz, this%drdz_d)
       deallocate(this%drdz)
    end if
 
    if (allocated(this%dsdz)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%dsdz, this%dsdz_d)
       deallocate(this%dsdz)
    end if
 
    if (allocated(this%dtdz)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%dtdz, this%dtdz_d)
       deallocate(this%dtdz)
    end if
 
    if (allocated(this%jac)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%jac, this%jac_d)
       deallocate(this%jac)
    end if
 
    if (allocated(this%jacinv)) then
       if (neko_bcknd_device .eq. 1) then
          call device_unmap(this%jacinv, this%jacinv_d)
       end if
       deallocate(this%jacinv)
    end if
 
    if (allocated(this%h1)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%h1, this%h1_d)
       deallocate(this%h1)
    end if
 
    if (allocated(this%h2)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%h2, this%h2_d)
       deallocate(this%h2)
    end if
 
    if (allocated(this%area)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%area, this%area_d)
       deallocate(this%area)
    end if
 
    if (allocated(this%nx)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%nx, this%nx_d)
       deallocate(this%nx)
    end if
 
    if (allocated(this%ny)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%ny, this%ny_d)
       deallocate(this%ny)
    end if
 
    if (allocated(this%nz)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%nz, this%nz_d)
       deallocate(this%nz)
    end if
 
    if (allocated(this%cyc_msk)) then
       if (neko_bcknd_device .eq. 1) then
          call device_unmap(this%cyc_msk, this%cyc_msk_d)
       end if
       deallocate(this%cyc_msk)
    end if
 
    if (allocated(this%R11)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%R11, this%R11_d)
       deallocate(this%R11)
    end if
 
    if (allocated(this%R12)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%R12, this%R12_d)
       deallocate(this%R12)
    end if
 
 
    nullify(this%msh)
    nullify(this%Xh)
    nullify(this%dof)
    nullify(this%gs_h)
 
  end subroutine coef_free
 
  subroutine coef_generate_dxyzdrst(c)
    type(coef_t), intent(inout) :: c
    integer :: e, i, lxy, lyz, ntot
 
    lxy = c%Xh%lx*c%Xh%ly
    lyz = c%Xh%ly*c%Xh%lz
    ntot = c%dof%size()
 
    associate(drdx => c%drdx, drdy => c%drdy, drdz => c%drdz, &
         dsdx => c%dsdx, dsdy => c%dsdy, dsdz => c%dsdz, &
         dtdx => c%dtdx, dtdy => c%dtdy, dtdz => c%dtdz, &
         dxdr => c%dxdr, dydr => c%dydr, dzdr => c%dzdr, &
         dxds => c%dxds, dyds => c%dyds, dzds => c%dzds, &
         dxdt => c%dxdt, dydt => c%dydt, dzdt => c%dzdt, &
         dx => c%Xh%dx, dy => c%Xh%dy, dz => c%Xh%dz, &
         x => c%dof%x, y => c%dof%y, z => c%dof%z, &
         lx => c%Xh%lx, ly => c%Xh%ly, lz => c%Xh%lz, &
         dyt => c%Xh%dyt, dzt => c%Xh%dzt, &
         jacinv => c%jacinv, jac => c%jac)
 
      if (neko_bcknd_device .eq. 1) then
 
         call device_coef_generate_dxydrst(c%drdx_d, c%drdy_d, c%drdz_d, &
              c%dsdx_d, c%dsdy_d, c%dsdz_d, c%dtdx_d, c%dtdy_d, c%dtdz_d, &
              c%dxdr_d, c%dydr_d, c%dzdr_d, c%dxds_d, c%dyds_d, c%dzds_d, &
              c%dxdt_d, c%dydt_d, c%dzdt_d, c%Xh%dx_d, c%Xh%dy_d, c%Xh%dz_d, &
              c%dof%x_d, c%dof%y_d, c%dof%z_d, c%jacinv_d, c%jac_d, &
              c%Xh%lx, c%msh%nelv)
 
         ! copy to host only at initialization.
         if (.not. c%coef_metrics_initialized) then
            call device_memcpy(dxdr, c%dxdr_d, ntot, device_to_host, &
                 sync = .false.)
            call device_memcpy(dydr, c%dydr_d, ntot, device_to_host, &
                 sync = .false.)
            call device_memcpy(dzdr, c%dzdr_d, ntot, device_to_host, &
                 sync = .false.)
            call device_memcpy(dxds, c%dxds_d, ntot, device_to_host, &
                 sync = .false.)
            call device_memcpy(dyds, c%dyds_d, ntot, device_to_host, &
                 sync = .false.)
            call device_memcpy(dzds, c%dzds_d, ntot, device_to_host, &
                 sync = .false.)
            call device_memcpy(dxdt, c%dxdt_d, ntot, device_to_host, &
                 sync = .false.)
            call device_memcpy(dydt, c%dydt_d, ntot, device_to_host, &
                 sync = .false.)
            call device_memcpy(dzdt, c%dzdt_d, ntot, device_to_host, &
                 sync = .false.)
            call device_memcpy(drdx, c%drdx_d, ntot, device_to_host, &
                 sync = .false.)
            call device_memcpy(drdy, c%drdy_d, ntot, device_to_host, &
                 sync = .false.)
            call device_memcpy(drdz, c%drdz_d, ntot, device_to_host, &
                 sync = .false.)
            call device_memcpy(dsdx, c%dsdx_d, ntot, device_to_host, &
                 sync = .false.)
            call device_memcpy(dsdy, c%dsdy_d, ntot, device_to_host, &
                 sync = .false.)
            call device_memcpy(dsdz, c%dsdz_d, ntot, device_to_host, &
                 sync = .false.)
            call device_memcpy(dtdx, c%dtdx_d, ntot, device_to_host, &
                 sync = .false.)
            call device_memcpy(dtdy, c%dtdy_d, ntot, device_to_host, &
                 sync = .false.)
            call device_memcpy(dtdz, c%dtdz_d, ntot, device_to_host, &
                 sync = .false.)
            call device_memcpy(jac, c%jac_d, ntot, device_to_host, &
                 sync = .false.)
            call device_memcpy(jacinv, c%jacinv_d, ntot, device_to_host, &
                 sync = .true.)
         end if
 
      else
         !$omp parallel do private(i)
         do e = 1, c%msh%nelv
            call mxm(dx, lx, x(1,1,1,e), lx, dxdr(1,1,1,e), lyz)
            call mxm(dx, lx, y(1,1,1,e), lx, dydr(1,1,1,e), lyz)
            call mxm(dx, lx, z(1,1,1,e), lx, dzdr(1,1,1,e), lyz)
 
            do i = 1, lz
               call mxm(x(1,1,i,e), lx, dyt, ly, dxds(1,1,i,e), ly)
               call mxm(y(1,1,i,e), lx, dyt, ly, dyds(1,1,i,e), ly)
               call mxm(z(1,1,i,e), lx, dyt, ly, dzds(1,1,i,e), ly)
            end do
 
            ! We actually take 2d into account, wow, need to do that for the rest.
            if (c%msh%gdim .eq. 3) then
               call mxm(x(1,1,1,e), lxy, dzt, lz, dxdt(1,1,1,e), lz)
               call mxm(y(1,1,1,e), lxy, dzt, lz, dydt(1,1,1,e), lz)
               call mxm(z(1,1,1,e), lxy, dzt, lz, dzdt(1,1,1,e), lz)
            else
               call rzero(dxdt(1,1,1,e), lxy)
               call rzero(dydt(1,1,1,e), lxy)
               call rone(dzdt(1,1,1,e), lxy)
            end if
         end do
         !$omp end parallel do
 
         if (c%msh%gdim .eq. 2) then
            call rzero (jac, ntot)
            call addcol3 (jac, dxdr, dyds, ntot)
            call subcol3 (jac, dxds, dydr, ntot)
            call copy (drdx, dyds, ntot)
            call copy (drdy, dxds, ntot)
            call chsign (drdy, ntot)
            call copy (dsdx, dydr, ntot)
            call chsign (dsdx, ntot)
            call copy (dsdy, dxdr, ntot)
            call rzero (drdz, ntot)
            call rzero (dsdz, ntot)
            call rone (dtdz, ntot)
         else
            !$omp parallel private(i)
            !$omp do simd
            do i = 1, ntot
               c%jac(i, 1, 1, 1) = 0.0_rp
            end do
            !$omp end do simd
            !$omp do simd
            do i = 1, ntot
               c%jac(i, 1, 1, 1) = c%jac(i, 1, 1, 1) + ( c%dxdr(i, 1, 1, 1) &
                    * c%dyds(i, 1, 1, 1) * c%dzdt(i, 1, 1, 1) )
 
               c%jac(i, 1, 1, 1) = c%jac(i, 1, 1, 1) + ( c%dxdt(i, 1, 1, 1) &
                    * c%dydr(i, 1, 1, 1) * c%dzds(i, 1, 1, 1) )
 
               c%jac(i, 1, 1, 1) = c%jac(i, 1, 1, 1) + ( c%dxds(i, 1, 1, 1) &
                    * c%dydt(i, 1, 1, 1) * c%dzdr(i, 1, 1, 1) )
            end do
            !$omp end do simd
            !$omp do simd
            do i = 1, ntot
               c%jac(i, 1, 1, 1) = c%jac(i, 1, 1, 1) - ( c%dxdr(i, 1, 1, 1) &
                    * c%dydt(i, 1, 1, 1) * c%dzds(i, 1, 1, 1) )
 
               c%jac(i, 1, 1, 1) = c%jac(i, 1, 1, 1) - ( c%dxds(i, 1, 1, 1) &
                    * c%dydr(i, 1, 1, 1) * c%dzdt(i, 1, 1, 1) )
 
               c%jac(i, 1, 1, 1) = c%jac(i, 1, 1, 1) - ( c%dxdt(i, 1, 1, 1) &
                    * c%dyds(i, 1, 1, 1) * c%dzdr(i, 1, 1, 1) )
            end do
            !$omp end do simd
            !$omp do simd
            do i = 1, ntot
               c%drdx(i, 1, 1, 1) = c%dyds(i, 1, 1, 1) * c%dzdt(i, 1, 1, 1) &
                    - c%dydt(i, 1, 1, 1) * c%dzds(i, 1, 1, 1)
 
               c%drdy(i, 1, 1, 1) = c%dxdt(i, 1, 1, 1) * c%dzds(i, 1, 1, 1) &
                    - c%dxds(i, 1, 1, 1) * c%dzdt(i, 1, 1, 1)
 
               c%drdz(i, 1, 1, 1) = c%dxds(i, 1, 1, 1) * c%dydt(i, 1, 1, 1) &
                    - c%dxdt(i, 1, 1, 1) * c%dyds(i, 1, 1, 1)
            end do
            !$omp end do simd
            !$omp do simd
            do i = 1, ntot
               c%dsdx(i, 1, 1, 1) = c%dydt(i, 1, 1, 1) * c%dzdr(i, 1, 1, 1) &
                    - c%dydr(i, 1, 1, 1) * c%dzdt(i, 1, 1, 1)
 
               c%dsdy(i, 1, 1, 1) = c%dxdr(i, 1, 1, 1) * c%dzdt(i, 1, 1, 1) &
                    - c%dxdt(i, 1, 1, 1) * c%dzdr(i, 1, 1, 1)
 
               c%dsdz(i, 1, 1, 1) = c%dxdt(i, 1, 1, 1) * c%dydr(i, 1, 1, 1) &
                    - c%dxdr(i, 1, 1, 1) * c%dydt(i, 1, 1, 1)
            end do
            !$omp end do simd
            !$omp do simd
            do i = 1, ntot
               c%dtdx(i, 1, 1, 1) = c%dydr(i, 1, 1, 1) * c%dzds(i, 1, 1, 1) &
                    - c%dyds(i, 1, 1, 1) * c%dzdr(i, 1, 1, 1)
 
               c%dtdy(i, 1, 1, 1) = c%dxds(i, 1, 1, 1) * c%dzdr(i, 1, 1, 1) &
                    - c%dxdr(i, 1, 1, 1) * c%dzds(i, 1, 1, 1)
 
               c%dtdz(i, 1, 1, 1) = c%dxdr(i, 1, 1, 1) * c%dyds(i, 1, 1, 1) &
                    - c%dxds(i, 1, 1, 1) * c%dydr(i, 1, 1, 1)
            end do
            !$omp end do simd
            !$omp end parallel
         end if
         call invers2(jacinv, jac, ntot)
      end if
    end associate
 
  end subroutine coef_generate_dxyzdrst
 
  subroutine coef_generate_geo(c)
    type(coef_t), intent(inout) :: c
    integer :: e, i, lxyz, ntot
 
    lxyz = c%Xh%lx * c%Xh%ly * c%Xh%lz
    ntot = c%dof%size()
 
    if (neko_bcknd_device .eq. 1) then
 
       call device_coef_generate_geo(c%G11_d, c%G12_d, c%G13_d, &
            c%G22_d, c%G23_d, c%G33_d, &
            c%drdx_d, c%drdy_d, c%drdz_d, &
            c%dsdx_d, c%dsdy_d, c%dsdz_d, &
            c%dtdx_d, c%dtdy_d, c%dtdz_d, &
            c%jacinv_d, c%Xh%w3_d, c%msh%nelv, &
            c%Xh%lx, c%msh%gdim)
 
       ! copy to host only at initialization.
       if (.not. c%coef_metrics_initialized) then
          call device_memcpy(c%G11, c%G11_d, ntot, device_to_host, &
               sync = .false.)
          call device_memcpy(c%G22, c%G22_d, ntot, device_to_host, &
               sync = .false.)
          call device_memcpy(c%G33, c%G33_d, ntot, device_to_host, &
               sync = .false.)
          call device_memcpy(c%G12, c%G12_d, ntot, device_to_host, &
               sync = .false.)
          call device_memcpy(c%G13, c%G13_d, ntot, device_to_host, &
               sync = .false.)
          call device_memcpy(c%G23, c%G23_d, ntot, device_to_host, &
               sync = .true.)
       end if
 
    else
       if (c%msh%gdim .eq. 2) then
 
          do i = 1, ntot
             c%G11(i, 1, 1, 1) = c%drdx(i, 1, 1, 1) * c%drdx(i, 1, 1, 1) &
                  + c%drdy(i, 1, 1, 1) * c%drdy(i, 1, 1, 1)
 
             c%G22(i, 1, 1, 1) = c%dsdx(i, 1, 1, 1) * c%dsdx(i, 1, 1, 1) &
                  + c%dsdy(i, 1, 1, 1) * c%dsdy(i, 1, 1, 1)
 
             c%G12(i, 1, 1, 1) = c%drdx(i, 1, 1, 1) * c%dsdx(i, 1, 1, 1) &
                  + c%drdy(i, 1, 1, 1) * c%dsdy(i, 1, 1, 1)
          end do
 
          do i = 1, ntot
             c%G11(i, 1, 1, 1) = c%G11(i, 1, 1, 1) * c%jacinv(i, 1, 1, 1)
             c%G22(i, 1, 1, 1) = c%G22(i, 1, 1, 1) * c%jacinv(i, 1, 1, 1)
             c%G12(i, 1, 1, 1) = c%G12(i, 1, 1, 1) * c%jacinv(i, 1, 1, 1)
             c%G33(i, 1, 1, 1) = 0.0_rp
             c%G13(i, 1, 1, 1) = 0.0_rp
             c%G23(i, 1, 1, 1) = 0.0_rp
          end do
 
          do concurrent(e = 1:c%msh%nelv)
             do concurrent(i = 1:lxyz)
                c%G11(i,1,1,e) = c%G11(i,1,1,e) * c%Xh%w3(i,1,1)
                c%G22(i,1,1,e) = c%G22(i,1,1,e) * c%Xh%w3(i,1,1)
                c%G12(i,1,1,e) = c%G12(i,1,1,e) * c%Xh%w3(i,1,1)
             end do
          end do
 
       else
          !$omp parallel private(i)
          !$omp do
          do i = 1, ntot
             c%G11(i, 1, 1, 1) = c%drdx(i, 1, 1, 1) * c%drdx(i, 1, 1, 1) &
                  + c%drdy(i, 1, 1, 1) * c%drdy(i, 1, 1, 1) &
                  + c%drdz(i, 1, 1, 1) * c%drdz(i, 1, 1, 1)
 
             c%G22(i, 1, 1, 1) = c%dsdx(i, 1, 1, 1) * c%dsdx(i, 1, 1, 1) &
                  + c%dsdy(i, 1, 1, 1) * c%dsdy(i, 1, 1, 1) &
                  + c%dsdz(i, 1, 1, 1) * c%dsdz(i, 1, 1, 1)
 
             c%G33(i, 1, 1, 1) = c%dtdx(i, 1, 1, 1) * c%dtdx(i, 1, 1, 1) &
                  + c%dtdy(i, 1, 1, 1) * c%dtdy(i, 1, 1, 1) &
                  + c%dtdz(i, 1, 1, 1) * c%dtdz(i, 1, 1, 1)
          end do
          !$omp end do
          !$omp do
          do i = 1, ntot
             c%G11(i, 1, 1, 1) = c%G11(i, 1, 1, 1) * c%jacinv(i, 1, 1, 1)
             c%G22(i, 1, 1, 1) = c%G22(i, 1, 1, 1) * c%jacinv(i, 1, 1, 1)
             c%G33(i, 1, 1, 1) = c%G33(i, 1, 1, 1) * c%jacinv(i, 1, 1, 1)
          end do
          !$omp end do
          !$omp do
          do i = 1, ntot
             c%G12(i, 1, 1, 1) = c%drdx(i, 1, 1, 1) * c%dsdx(i, 1, 1, 1) &
                  + c%drdy(i, 1, 1, 1) * c%dsdy(i, 1, 1, 1) &
                  + c%drdz(i, 1, 1, 1) * c%dsdz(i, 1, 1, 1)
 
             c%G13(i, 1, 1, 1) = c%drdx(i, 1, 1, 1) * c%dtdx(i, 1, 1, 1) &
                  + c%drdy(i, 1, 1, 1) * c%dtdy(i, 1, 1, 1) &
                  + c%drdz(i, 1, 1, 1) * c%dtdz(i, 1, 1, 1)
 
             c%G23(i, 1, 1, 1) = c%dsdx(i, 1, 1, 1) * c%dtdx(i, 1, 1, 1) &
                  + c%dsdy(i, 1, 1, 1) * c%dtdy(i, 1, 1, 1) &
                  + c%dsdz(i, 1, 1, 1) * c%dtdz(i, 1, 1, 1)
          end do
          !$omp end do
          !$omp do
          do i = 1, ntot
             c%G12(i, 1, 1, 1) = c%G12(i, 1, 1, 1) * c%jacinv(i, 1, 1, 1)
             c%G13(i, 1, 1, 1) = c%G13(i, 1, 1, 1) * c%jacinv(i, 1, 1, 1)
             c%G23(i, 1, 1, 1) = c%G23(i, 1, 1, 1) * c%jacinv(i, 1, 1, 1)
          end do
          !$omp end do
          !$omp do
          do e = 1, c%msh%nelv
             do concurrent(i = 1:lxyz)
                c%G11(i,1,1,e) = c%G11(i,1,1,e) * c%Xh%w3(i,1,1)
                c%G22(i,1,1,e) = c%G22(i,1,1,e) * c%Xh%w3(i,1,1)
                c%G12(i,1,1,e) = c%G12(i,1,1,e) * c%Xh%w3(i,1,1)
 
                c%G33(i,1,1,e) = c%G33(i,1,1,e) * c%Xh%w3(i,1,1)
                c%G13(i,1,1,e) = c%G13(i,1,1,e) * c%Xh%w3(i,1,1)
                c%G23(i,1,1,e) = c%G23(i,1,1,e) * c%Xh%w3(i,1,1)
             end do
          end do
          !$omp end do
          !$omp end parallel
       end if
    end if
 
  end subroutine coef_generate_geo
 
  subroutine coef_generate_geo_compressed(c)
    type(coef_t), intent(inout) :: c
    integer :: e, m, i, lxyz, m_max
    integer, allocatable :: c_inds_rev(:) ! reverse compression indices map
    real(kind=rp) :: ctol = 1.0e-7_rp
    real(kind=rp) :: diff = 0.0_rp
 
    ! First step, allocate full-size lookup structure for entire mesh
    allocate(c%compression_inds(c%msh%nelv))
    allocate(c_inds_rev(c%msh%nelv))
 
    ! Second step, loop over all elements, compute compression mapping
    ! First entry must be itself to get started
    m_max = 1
    c%compression_inds(1) = 1
    c_inds_rev = 0
    c_inds_rev(1) = 1
 
    ! Loop over elements, but skip first
    lxyz = c%Xh%lx * c%Xh%ly * c%Xh%lz
    do e = 2, c%msh%nelv
       ! Loop over possible compression candidates
       do m = 1, m_max
          diff = 0.0_rp
          ! Loop over quadrature points
          do i = 1, lxyz
             ! diff += abs( \| G(i,:,:,e) - G(i,:,:,reverse(m)) \|_l1 )
             diff = diff + abs(c%G11(i,1,1,e) - c%G11(i,1,1,c_inds_rev(m))) &
                         + 2.0*abs(c%G12(i,1,1,e) - c%G12(i,1,1,c_inds_rev(m))) &
                         + 2.0*abs(c%G13(i,1,1,e) - c%G13(i,1,1,c_inds_rev(m))) &
                         + abs(c%G22(i,1,1,e) - c%G22(i,1,1,c_inds_rev(m))) &
                         + 2.0*abs(c%G23(i,1,1,e) - c%G23(i,1,1,c_inds_rev(m))) &
                         + abs(c%G33(i,1,1,e) - c%G33(i,1,1,c_inds_rev(m)))
          end do
 
          ! match is found; mapping(e) is redundant
          if ( diff .le. ctol ) then
             c%compression_inds(e) = m
             exit
          end if
       end do
 
       ! never found a match
       if ( diff .gt. ctol ) then
          m_max = m_max + 1
          c%compression_inds(e) = m_max
          c_inds_rev(m_max) = e
       end if
    end do
 
    write(*,*)
    write(*,*) '------Mapping Compression-----'
    write(*,*) 'Compressed from ', c%msh%nelv, ' to ', m_max
 
    ! Third step, allocate and fill Gij_compressed objects
    allocate(c%G11_compressed(c%Xh%lx, c%Xh%ly, c%Xh%lz, m_max))
    allocate(c%G22_compressed(c%Xh%lx, c%Xh%ly, c%Xh%lz, m_max))
    allocate(c%G33_compressed(c%Xh%lx, c%Xh%ly, c%Xh%lz, m_max))
    allocate(c%G12_compressed(c%Xh%lx, c%Xh%ly, c%Xh%lz, m_max))
    allocate(c%G13_compressed(c%Xh%lx, c%Xh%ly, c%Xh%lz, m_max))
    allocate(c%G23_compressed(c%Xh%lx, c%Xh%ly, c%Xh%lz, m_max))
    do m = 1, m_max
       do i = 1, lxyz
          c%G11_compressed(i,1,1,m) = c%G11(i,1,1,c_inds_rev(m))
          c%G22_compressed(i,1,1,m) = c%G22(i,1,1,c_inds_rev(m))
          c%G33_compressed(i,1,1,m) = c%G33(i,1,1,c_inds_rev(m))
          c%G12_compressed(i,1,1,m) = c%G12(i,1,1,c_inds_rev(m))
          c%G13_compressed(i,1,1,m) = c%G13(i,1,1,c_inds_rev(m))
          c%G23_compressed(i,1,1,m) = c%G23(i,1,1,c_inds_rev(m))
       end do
    end do
 
    deallocate(c_inds_rev)
 
  end subroutine coef_generate_geo_compressed
 
  subroutine coef_generate_mass(c)
    type(coef_t), intent(inout) :: c
    integer :: e, i, lxyz, ntot
 
    lxyz = c%Xh%lx * c%Xh%ly * c%Xh%lz
    ntot = c%dof%size()
 
    if (neko_bcknd_device .eq. 1) then
       call device_coef_generate_mass(c%B_d, c%Binv_d, c%jac_d, c%Xh%w3_d, &
            lxyz, c%msh%nelv)
       ! copy to host only at initialization.
       if (.not. c%coef_metrics_initialized) then
          call device_memcpy(c%B, c%B_d, ntot, device_to_host, sync = .false.)
       end if
    else
       do concurrent(e = 1:c%msh%nelv)
          ! Here we need to handle things differently for axis symmetric elements
          do concurrent(i = 1:lxyz)
             c%B(i,1,1,e) = c%jac(i,1,1,e) * c%Xh%w3(i,1,1)
             c%Binv(i,1,1,e) = c%B(i,1,1,e)
          end do
       end do
    end if
 
    call c%gs_h%op(c%Binv, ntot, gs_op_add)
 
    if (neko_bcknd_device .eq. 1) then
       call device_invcol1(c%Binv_d, ntot)
       ! copy to host only at initialization.
       if (.not. c%coef_metrics_initialized) then
          call device_memcpy(c%Binv, c%Binv_d, ntot, &
               device_to_host, sync = .true.)
       end if
    else
       call invcol1(c%Binv, ntot)
    end if
 
    if (neko_bcknd_device .eq. 1) then
       c%volume = device_glsum(c%B_d, ntot)
    else
       c%volume = glsum(c%B, ntot)
    end if
 
  end subroutine coef_generate_mass
 
  pure function coef_get_normal(this, i, j, k, e, facet) result(normal)
    class(coef_t), intent(in) :: this
    integer, intent(in) :: i, j, k, e, facet
    real(kind=rp) :: normal(3)
 
    select case (facet)
    case (1, 2)
       normal(1) = this%nx(j, k, facet, e)
       normal(2) = this%ny(j, k, facet, e)
       normal(3) = this%nz(j, k, facet, e)
    case (3, 4)
       normal(1) = this%nx(i, k, facet, e)
       normal(2) = this%ny(i, k, facet, e)
       normal(3) = this%nz(i, k, facet, e)
    case (5, 6)
       normal(1) = this%nx(i, j, facet, e)
       normal(2) = this%ny(i, j, facet, e)
       normal(3) = this%nz(i, j, facet, e)
    end select
  end function coef_get_normal
 
  pure function coef_get_area(this, i, j, k, e, facet) result(area)
    class(coef_t), intent(in) :: this
    integer, intent(in) :: i, j, k, e, facet
    real(kind=rp) :: area
 
    select case (facet)
    case (1, 2)
       area = this%area(j, k, facet, e)
    case (3, 4)
       area = this%area(i, k, facet, e)
    case (5, 6)
       area = this%area(i, j, facet, e)
    end select
  end function coef_get_area
 
 
  subroutine coef_generate_area_and_normal(coef)
    type(coef_t), intent(inout) :: coef
    real(kind=rp), allocatable :: a(:,:,:,:)
    real(kind=rp), allocatable :: b(:,:,:,:)
    real(kind=rp), allocatable :: c(:,:,:,:)
    real(kind=rp), allocatable :: dot(:,:,:,:)
    integer :: n, m, e, i, j, k, lx
    real(kind=rp) :: weight, len
    n = coef%dof%size()
    lx = coef%Xh%lx
 
    if (neko_bcknd_device .eq. 1) then
 
       call device_coef_generate_area_and_normal( &
            coef%area_d, coef%nx_d, coef%ny_d, coef%nz_d, &
            coef%dxdr_d, coef%dydr_d, coef%dzdr_d, &
            coef%dxds_d, coef%dyds_d, coef%dzds_d, &
            coef%dxdt_d, coef%dydt_d, coef%dzdt_d, &
            coef%Xh%wx_d, coef%Xh%wy_d, coef%Xh%wz_d, &
            lx, coef%msh%nelv, neko_eps)
 
       ! Here, we always copy back to host.
       m = size(coef%area)
       call device_memcpy(coef%area, coef%area_d, m, &
            device_to_host, sync = .false.)
       call device_memcpy(coef%nx, coef%nx_d, m, &
            device_to_host, sync = .false.)
       call device_memcpy(coef%ny, coef%ny_d, m, &
            device_to_host, sync = .false.)
       call device_memcpy(coef%nz, coef%nz_d, &
            m, device_to_host, sync = .true.)
 
    else
 
       allocate(a(coef%Xh%lx, coef%Xh%lx, coef%Xh%lx, coef%msh%nelv))
       allocate(b(coef%Xh%lx, coef%Xh%lx, coef%Xh%lx, coef%msh%nelv))
       allocate(c(coef%Xh%lx, coef%Xh%lx, coef%Xh%lx, coef%msh%nelv))
       allocate(dot(coef%Xh%lx, coef%Xh%lx, coef%Xh%lx, coef%msh%nelv))
 
       !$omp parallel private (e, i, j, k, weight, len)
 
       ! ds x dt
       !$omp do simd
       do i = 1, n
          a(i, 1, 1, 1) = coef%dyds(i, 1, 1, 1) * coef%dzdt(i, 1, 1, 1) &
               - coef%dzds(i, 1, 1, 1) * coef%dydt(i, 1, 1, 1)
 
          b(i, 1, 1, 1) = coef%dzds(i, 1, 1, 1) * coef%dxdt(i, 1, 1, 1) &
               - coef%dxds(i, 1, 1, 1) * coef%dzdt(i, 1, 1, 1)
 
          c(i, 1, 1, 1) = coef%dxds(i, 1, 1, 1) * coef%dydt(i, 1, 1, 1) &
               - coef%dyds(i, 1, 1, 1) * coef%dxdt(i, 1, 1, 1)
       end do
       !$omp end do simd
       !$omp do simd
       do i = 1, n
          dot(i, 1, 1, 1) = a(i, 1, 1, 1) * a(i, 1, 1, 1) &
               + b(i, 1, 1, 1) * b(i, 1, 1, 1) &
               + c(i, 1, 1, 1) * c(i, 1, 1, 1)
       end do
       !$omp end do simd
       !$omp do
       do e = 1, coef%msh%nelv
          do concurrent(k = 1:coef%Xh%lx)
             do concurrent(j = 1:coef%Xh%lx)
                weight = coef%Xh%wy(j) * coef%Xh%wz(k)
                coef%area(j, k, 2, e) = sqrt(dot(lx, j, k, e)) * weight
                coef%area(j, k, 1, e) = sqrt(dot(1, j, k, e)) * weight
                coef%nx(j,k, 1, e) = -a(1, j, k, e)
                coef%nx(j,k, 2, e) = a(lx, j, k, e)
                coef%ny(j,k, 1, e) = -b(1, j, k, e)
                coef%ny(j,k, 2, e) = b(lx, j, k, e)
                coef%nz(j,k, 1, e) = -c(1, j, k, e)
                coef%nz(j,k, 2, e) = c(lx, j, k, e)
             end do
          end do
       end do
       !$omp end do
 
       ! dr x dt
       !$omp do simd
       do i = 1, n
          a(i, 1, 1, 1) = coef%dydr(i, 1, 1, 1) * coef%dzdt(i, 1, 1, 1) &
               - coef%dzdr(i, 1, 1, 1) * coef%dydt(i, 1, 1, 1)
 
          b(i, 1, 1, 1) = coef%dzdr(i, 1, 1, 1) * coef%dxdt(i, 1, 1, 1) &
               - coef%dxdr(i, 1, 1, 1) * coef%dzdt(i, 1, 1, 1)
 
          c(i, 1, 1, 1) = coef%dxdr(i, 1, 1, 1) * coef%dydt(i, 1, 1, 1) &
               - coef%dydr(i, 1, 1, 1) * coef%dxdt(i, 1, 1, 1)
       end do
       !$omp end do simd
       !$omp do simd
       do i = 1, n
          dot(i, 1, 1, 1) = a(i, 1, 1, 1) * a(i, 1, 1, 1) &
               + b(i, 1, 1, 1) * b(i, 1, 1, 1) &
               + c(i, 1, 1, 1) * c(i, 1, 1, 1)
       end do
       !$omp end do simd
       !$omp do
       do e = 1, coef%msh%nelv
          do concurrent(k = 1:coef%Xh%lx)
             do concurrent(j = 1:coef%Xh%lx)
                weight = coef%Xh%wx(j) * coef%Xh%wz(k)
                coef%area(j, k, 3, e) = sqrt(dot(j, 1, k, e)) * weight
                coef%area(j, k, 4, e) = sqrt(dot(j, lx, k, e)) * weight
                coef%nx(j,k, 3, e) = a(j, 1, k, e)
                coef%nx(j,k, 4, e) = -a(j, lx, k, e)
                coef%ny(j,k, 3, e) = b(j, 1, k, e)
                coef%ny(j,k, 4, e) = -b(j, lx, k, e)
                coef%nz(j,k, 3, e) = c(j, 1, k, e)
                coef%nz(j,k, 4, e) = -c(j, lx, k, e)
             end do
          end do
       end do
       !$omp end do
       ! dr x ds
       !$omp do simd
       do i = 1, n
          a(i, 1, 1, 1) = coef%dydr(i, 1, 1, 1) * coef%dzds(i, 1, 1, 1) &
               - coef%dzdr(i, 1, 1, 1) * coef%dyds(i, 1, 1, 1)
 
          b(i, 1, 1, 1) = coef%dzdr(i, 1, 1, 1) * coef%dxds(i, 1, 1, 1) &
               - coef%dxdr(i, 1, 1, 1) * coef%dzds(i, 1, 1, 1)
 
          c(i, 1, 1, 1) = coef%dxdr(i, 1, 1, 1) * coef%dyds(i, 1, 1, 1) &
               - coef%dydr(i, 1, 1, 1) * coef%dxds(i, 1, 1, 1)
       end do
       !$omp end do simd
       !$omp do simd
       do i = 1, n
          dot(i, 1, 1, 1) = a(i, 1, 1, 1) * a(i, 1, 1, 1) &
               + b(i, 1, 1, 1) * b(i, 1, 1, 1) &
               + c(i, 1, 1, 1) * c(i, 1, 1, 1)
       end do
       !$omp end do simd
       !$omp do
       do e = 1, coef%msh%nelv
          do concurrent(k = 1:coef%Xh%lx)
             do concurrent(j = 1:coef%Xh%lx)
                weight = coef%Xh%wx(j) * coef%Xh%wy(k)
                coef%area(j, k, 5, e) = sqrt(dot(j, k, 1, e)) * weight
                coef%area(j, k, 6, e) = sqrt(dot(j, k, lx, e)) * weight
                coef%nx(j,k, 5, e) = -a(j, k, 1, e)
                coef%nx(j,k, 6, e) = a(j, k, lx, e)
                coef%ny(j,k, 5, e) = -b(j, k, 1, e)
                coef%ny(j,k, 6, e) = b(j, k, lx, e)
                coef%nz(j,k, 5, e) = -c(j, k, 1, e)
                coef%nz(j,k, 6, e) = c(j, k, lx, e)
             end do
          end do
       end do
       !$omp end do
       ! Normalize
       !$omp do
       do j = 1, size(coef%nz)
          len = sqrt(coef%nx(j,1,1,1)**2 + &
               coef%ny(j,1,1,1)**2 + coef%nz(j,1,1,1)**2)
          if (len .gt. neko_eps) then
             coef%nx(j,1,1,1) = coef%nx(j,1,1,1) / len
             coef%ny(j,1,1,1) = coef%ny(j,1,1,1) / len
             coef%nz(j,1,1,1) = coef%nz(j,1,1,1) / len
          end if
       end do
       !$omp end do
       !$omp end parallel
 
       deallocate(dot)
       deallocate(c)
       deallocate(b)
       deallocate(a)
 
    end if
 
  end subroutine coef_generate_area_and_normal
 
  subroutine coef_generate_cyclic_bc(this)
    class(coef_t), intent(inout) :: this
    real(kind=rp) :: un(3), len, d
    integer :: lx, ly, lz, np, np_glb, ierr
    integer :: i, j, k, pf, pe, n, nc, ncyc
 
    if (.not. this%cyclic) return
 
    np = this%msh%periodic%size
    call mpi_allreduce(np, np_glb, 1, &
         mpi_integer, mpi_sum, neko_comm, ierr)
 
    if (np_glb .eq. 0) then
       call neko_error("There are no periodic boundaries. " // &
            "Switch cyclic off in the case file.")
    end if
 
    if (np .eq. 0) return
 
    lx = this%Xh%lx
    ly = this%Xh%ly
    lz = this%Xh%lz
    ncyc = this%cyc_msk(0) - 1
    nc = 1
    do n = 1, np
       pf = this%msh%periodic%facet_el(n)%x(1)
       pe = this%msh%periodic%facet_el(n)%x(2)
       do k = 1, lz
          do j = 1, ly
             do i = 1, lx
                if (index_is_on_facet(i, j, k, lx, ly, lz, pf)) then
                   un = this%get_normal(i, j, k, pe, pf)
                   len = sqrt(un(1) * un(1) + un(2) * un(2))
                   if (len .gt. neko_eps) then
                      d = this%dof%y(i, j, k, pe) * un(1) &
                           - this%dof%x(i, j, k, pe) * un(2)
 
                      this%cyc_msk(nc) = linear_index(i, j, k, pe, lx, ly, lz)
                      this%R11(nc) = un(1) / len * sign(1.0_rp, d)
                      this%R12(nc) = un(2) / len * sign(1.0_rp, d)
                      nc = nc + 1
                   else
                      call neko_error("x and y components of surface " // &
                           "normals are zero. Cyclic rotations must be " // &
                           "around z-axis.")
                   end if
                end if
             end do
          end do
       end do
    end do
 
    if (nc - 1 /= ncyc) then
       call neko_error("The number of cyclic GLL points were " // &
            "not estimated correctly.")
    end if
 
    if (neko_bcknd_device .eq. 1) then
       call device_memcpy(this%cyc_msk, this%cyc_msk_d, ncyc+1, &
            host_to_device, sync = .false.)
       call device_memcpy(this%R11, this%R11_d, ncyc, &
            host_to_device, sync = .false.)
       call device_memcpy(this%R12, this%R12_d, ncyc, &
            host_to_device, sync = .false.)
    end if
 
  end subroutine coef_generate_cyclic_bc
 
 
  subroutine coef_recompute_metrics(this)
    class(coef_t), intent(inout) :: this
 
    call coef_generate_dxyzdrst(this)
    call coef_generate_geo(this)
    call coef_generate_area_and_normal(this)
    call coef_generate_mass(this)
    if (this%cyclic) then
       call coef_generate_cyclic_bc(this)
    end if
  end subroutine coef_recompute_metrics
 
 
  subroutine coef_enable_lagged_mass(this)
    class(coef_t), intent(inout), target :: this
    integer :: n
 
    ! Return if already allocated distinctly
    if (.not. associated(this%Blag, this%B)) return
 
    n = this%Xh%lx * this%Xh%ly * this%Xh%lz * this%msh%nelv
 
 
    nullify(this%Blag)
    nullify(this%Blaglag)
 
    allocate(this%Blag(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
    allocate(this%Blaglag(this%Xh%lx, this%Xh%ly, this%Xh%lz, this%msh%nelv))
 
    this%Blag = this%B
    this%Blaglag = this%B
 
    if (neko_bcknd_device .eq. 1) then
 
       this%Blag_d = c_null_ptr
       this%Blaglag_d = c_null_ptr
 
       call device_map(this%Blag, this%Blag_d, n)
       call device_map(this%Blaglag, this%Blaglag_d, n)
 
       call device_memcpy(this%Blag, this%Blag_d, n, &
            host_to_device, sync = .false.)
       call device_memcpy(this%Blaglag, this%Blaglag_d, n, &
            host_to_device, sync = .true.)
    end if
 
  end subroutine coef_enable_lagged_mass
 
 
  subroutine coef_update_lagged_mass(this)
    class(coef_t), intent(inout), target :: this
    integer :: n
 
    ! If this%Blag does not have separate memory, we don't need to update it.
    if (associated(this%Blag, this%B)) return
    n = this%Xh%lx * this%Xh%ly * this%Xh%lz * this%msh%nelv
    if (neko_bcknd_device .eq. 1) then
       call device_copy(this%Blaglag_d, this%Blag_d, n)
       call device_copy(this%Blag_d, this%B_d, n)
    else
       this%Blaglag = this%Blag
       this%Blag = this%B
    end if
 
  end subroutine coef_update_lagged_mass
 
end module coefs