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, &
       neko_bcknd_sx, neko_bcknd_xsmm
  use num_types, only : rp, sp, dp
  use dofmap, only : dofmap_t
  use space, only : space_t
  use math, only : rone, invcol1, addcol3, subcol3, copy, &
       chsign, rzero, invers2, glsum, glmax, neko_eps, neko_eps_sp, &
       eig_sym2, eig_sym3
  use logger, only : neko_log, log_size
  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
 
  real(kind=rp), public, parameter :: neko_metric_cond_sp = 1.0e4_rp
 
  real(kind=rp), public, parameter :: neko_metric_perturb_max = 1.0e-5_rp
 
  integer, public, parameter :: coef_full = 0
  integer, public, parameter :: coef_operator = 1
 
  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) :: metric_cond = 0.0_rp
     real(kind=rp) :: metric_scaled_cond = 0.0_rp
     integer :: metric_degenerate = 0
     real(kind=rp) :: metric_perturb = 0.0_rp
     logical :: metric_sp_safe = .false.
     real(kind=rp), allocatable :: g11_compr(:,:,:,:)
     real(kind=rp), allocatable :: g22_compr(:,:,:,:)
     real(kind=rp), allocatable :: g33_compr(:,:,:,:)
     real(kind=rp), allocatable :: g12_compr(:,:,:,:)
     real(kind=rp), allocatable :: g13_compr(:,:,:,:)
     real(kind=rp), allocatable :: g23_compr(:,:,:,:)
     integer, allocatable :: compression_inds(:)
     logical :: geo_compression = .false.
     real(kind=rp) :: geo_compression_tol = 0.0_rp
 
     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.
 
     integer :: scope = coef_full
 
 
     real(kind=rp) :: volume = 0.0_rp
 
     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) :: g11_compr_d = c_null_ptr
     type(c_ptr) :: g22_compr_d = c_null_ptr
     type(c_ptr) :: g33_compr_d = c_null_ptr
     type(c_ptr) :: g12_compr_d = c_null_ptr
     type(c_ptr) :: g13_compr_d = c_null_ptr
     type(c_ptr) :: g23_compr_d = c_null_ptr
     type(c_ptr) :: compression_inds_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
 
     integer :: metrics_version = 0
 
   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) :: require_facets => coef_require_facets
     procedure, pass(this) :: generate_cyclic_bc => coef_generate_cyclic_bc
     procedure, pass(this) :: recompute_metrics => coef_recompute_metrics
     procedure, pass(this) :: metric_condition => coef_metric_condition
     procedure, pass(this) :: enable_geo_compression => &
          coef_enable_geo_compression
     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, scope)
    class(coef_t), intent(inout), target :: this
    type(gs_t), intent(inout), target :: gs_h
    integer, intent(in), optional :: scope
    integer :: n, m, ncyc
 
    call this%free()
 
    ! After free(), so that it survives into the initialized state
    if (present(scope)) then
       if (scope .ne. coef_full .and. scope .ne. coef_operator) then
          call neko_error('Unknown coefficient scope')
       end if
       this%scope = scope
    end if
 
    call neko_log%section('Coefficients')
 
    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))
 
    if (this%scope .eq. coef_full) then
       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))
    end if
 
    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
 
       if (this%scope .eq. coef_full) then
          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
 
    end if
 
    call coef_generate_dxyzdrst(this)
 
    call coef_generate_geo(this)
 
    ! Both are dead weight under COEF_OPERATOR: nothing reads the facet
    ! metrics, and the conditioning diagnostic is two per-point eigenvalue
    ! solves plus three reductions spent on something no one consumes yet.
    if (this%scope .eq. coef_full) then
       call coef_metric_condition(this)
 
       call coef_generate_area_and_normal(this)
    end if
 
    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
 
    call coef_release_scratch(this)
 
    call neko_log%end_section()
 
  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
 
    call coef_free_geo_compr(this)
    this%geo_compression = .false.
    this%geo_compression_tol = 0.0_rp
 
    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)
 
    this%coef_metrics_initialized = .false.
    this%scope = coef_full
 
  end subroutine coef_free
 
  subroutine coef_release_scratch(this)
    class(coef_t), intent(inout), target :: this
 
    if (this%scope .eq. coef_full) return
 
    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%drdy)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%drdy, this%drdy_d)
       deallocate(this%drdy)
    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%dsdx)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%dsdx, this%dsdx_d)
       deallocate(this%dsdx)
    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%dsdz)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%dsdz, this%dsdz_d)
       deallocate(this%dsdz)
    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%dtdy)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%dtdy, this%dtdy_d)
       deallocate(this%dtdy)
    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%Binv)) then
       if (neko_bcknd_device .eq. 1) call device_unmap(this%Binv, this%Binv_d)
       deallocate(this%Binv)
    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
 
  end subroutine coef_release_scratch
 
  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%x, y => c%dof%y%x, z => c%dof%z%x, &
         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%x_d, c%dof%y%x_d, c%dof%z%x_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_metric_condition(this)
    class(coef_t), intent(inout) :: this
    real(kind=dp) :: e1, e2, e3, scal
    real(kind=dp) :: c1, c2, c3, c12, c13, c23
    real(kind=rp) :: kmax, kcmax, tmp(1)
    integer :: i, j, k, e, n, ndeg, ndeg_glb, ierr
    character(len=LOG_SIZE) :: log_buf
 
    n = this%dof%size()
 
    ! Post-initialization the host copy is stale on device builds
    if (neko_bcknd_device .eq. 1 .and. this%coef_metrics_initialized) then
       call device_memcpy(this%G11, this%G11_d, n, device_to_host, &
            sync = .false.)
       call device_memcpy(this%G22, this%G22_d, n, device_to_host, &
            sync = .false.)
       call device_memcpy(this%G33, this%G33_d, n, device_to_host, &
            sync = .false.)
       call device_memcpy(this%G12, this%G12_d, n, device_to_host, &
            sync = .false.)
       call device_memcpy(this%G13, this%G13_d, n, device_to_host, &
            sync = .false.)
       call device_memcpy(this%G23, this%G23_d, n, device_to_host, &
            sync = .true.)
    end if
 
    kmax = 0.0_rp
    kcmax = 0.0_rp
    ndeg = 0
 
    do e = 1, this%msh%nelv
       do k = 1, this%Xh%lz
          do j = 1, this%Xh%ly
             do i = 1, this%Xh%lx
 
                ! Scale out the magnitude before the eigenvalue solve; the
                ! condition number is invariant under it and w3*J spans a
                ! wide range within an element
                scal = max(abs(real(this%G11(i,j,k,e), dp)), &
                           abs(real(this%G22(i,j,k,e), dp)))
                scal = max(scal, abs(real(this%G33(i,j,k,e), dp)))
                scal = max(scal, abs(real(this%G12(i,j,k,e), dp)))
                scal = max(scal, abs(real(this%G13(i,j,k,e), dp)))
                scal = max(scal, abs(real(this%G23(i,j,k,e), dp)))
 
                if (scal .le. 0.0_dp) then
                   ndeg = ndeg + 1
                   cycle
                end if
 
                if (this%msh%gdim .eq. 2) then
                   call eig_sym2(real(this%G11(i,j,k,e), dp) / scal, &
                        real(this%G22(i,j,k,e), dp) / scal, &
                        real(this%G12(i,j,k,e), dp) / scal, e1, e3)
                else
                   call eig_sym3(real(this%G11(i,j,k,e), dp) / scal, &
                        real(this%G22(i,j,k,e), dp) / scal, &
                        real(this%G33(i,j,k,e), dp) / scal, &
                        real(this%G12(i,j,k,e), dp) / scal, &
                        real(this%G13(i,j,k,e), dp) / scal, &
                        real(this%G23(i,j,k,e), dp) / scal, e1, e2, e3)
                end if
 
                if (e3 .le. 0.0_dp) then
                   ndeg = ndeg + 1
                else
                   kmax = max(kmax, real(e1 / e3, rp))
 
                   ! Jacobi scale the metric: C_ij = G_ij/sqrt(G_ii G_jj).
                   ! kappa(C) is what bounds the effect of rounding G to a
                   ! lower precision, and unlike kappa(G) it is independent
                   ! of the element aspect ratio. The point is positive
                   ! definite here, so the diagonal is strictly positive and
                   ! the square roots are safe. An orthogonal element gives
                   ! the identity, which eig_sym3 returns exactly through
                   ! its isotropic branch.
                   c12 = real(this%G12(i,j,k,e), dp) &
                        / sqrt(real(this%G11(i,j,k,e), dp) &
                        * real(this%G22(i,j,k,e), dp))
 
                   if (this%msh%gdim .eq. 2) then
                      call eig_sym2(1.0_dp, 1.0_dp, c12, c1, c3)
                   else
                      c13 = real(this%G13(i,j,k,e), dp) &
                           / sqrt(real(this%G11(i,j,k,e), dp) &
                           * real(this%G33(i,j,k,e), dp))
                      c23 = real(this%G23(i,j,k,e), dp) &
                           / sqrt(real(this%G22(i,j,k,e), dp) &
                           * real(this%G33(i,j,k,e), dp))
                      call eig_sym3(1.0_dp, 1.0_dp, 1.0_dp, &
                           c12, c13, c23, c1, c2, c3)
                   end if
 
                   if (c3 .gt. 0.0_dp) then
                      kcmax = max(kcmax, real(c1 / c3, rp))
                   else
                      ! C is a positive diagonal congruence of a metric
                      ! already found definite, so this is unreachable
                      ! analytically. Fail safe rather than underreport.
                      kcmax = huge(0.0_rp)
                   end if
                end if
 
             end do
          end do
       end do
    end do
 
    tmp(1) = kmax
    this%metric_cond = glmax(tmp, 1)
 
    tmp(1) = kcmax
    this%metric_scaled_cond = glmax(tmp, 1)
 
    call mpi_allreduce(ndeg, ndeg_glb, 1, mpi_integer, mpi_sum, &
         neko_comm, ierr)
    this%metric_degenerate = ndeg_glb
 
    ! The two condition numbers enter through different unit roundoffs, so
    ! the estimate is additive: storage rounding of G contributes
    ! eps_sp*kappa(C) and accumulation in rp contributes eps_rp*kappa(G).
    ! On a double build the second term all but vanishes and skew decides; on
    ! an rp = sp build it dominates and aspect ratio decides. No test on the
    ! build is needed, which is the point of writing it this way.
    this%metric_perturb = real(neko_eps_sp, rp) * this%metric_scaled_cond &
         + neko_eps * this%metric_cond
 
    this%metric_sp_safe = (this%metric_degenerate .eq. 0) .and. &
         (this%metric_perturb .le. neko_metric_perturb_max)
 
    write(log_buf, '(A,ES12.5)') 'Metric condition : ', this%metric_cond
    call neko_log%message(log_buf)
    write(log_buf, '(A,ES12.5)') 'Metric skew cond : ', &
         this%metric_scaled_cond
    call neko_log%message(log_buf)
    write(log_buf, '(A,ES12.5)') 'Metric perturb   : ', this%metric_perturb
    call neko_log%message(log_buf)
    write(log_buf, '(A,L1)') 'Metric fp32 safe : ', this%metric_sp_safe
    call neko_log%message(log_buf)
 
    ! Three independent questions, reported independently: can the solve
    ! break down, how far wrong is the operator, and is the mesh valid at
    ! all. They are not grades of one thing, so none suppresses another -- a
    ! wall resolved boundary layer in a single precision build trips the
    ! first two at once, and both are worth saying, because the definiteness
    ! limit is a possibility while the perturbation is a certainty.
 
    ! Past the arithmetic limit the local operator can lose positive
    ! definiteness and a Krylov solve can break down rather than degrade,
    ! which is only reachable when rp is sp. Reported rather than fatal: the
    ! threshold keeps three orders of margin below 1/eps_sp, so crossing it
    ! makes breakdown possible, not certain, and the run may well be fine.
    if (rp .eq. sp .and. this%metric_cond .gt. neko_metric_cond_sp) then
       write(log_buf, '(A,ES12.5)') &
            'Metric too ill conditioned for single precision, limit ', &
            neko_metric_cond_sp
       call neko_log%warning(log_buf)
       call neko_log%message('Geometric factors may lose positive ' // &
            'definiteness, consider a double precision build')
    end if
 
    ! How much the single precision factors actually perturb the operator,
    ! regardless of whether the definiteness limit above was also crossed.
    ! On a single precision build that error is being incurred now, so it
    ! warrants a warning; on a double precision build it describes a storage
    ! path that may not be in use, so it is recorded without one. Which of
    ! the two terms dominates is the actionable part: skew is a meshing
    ! problem, aspect ratio is a precision problem, and they have different
    ! remedies.
    if (this%metric_perturb .gt. neko_metric_perturb_max) then
       if (rp .eq. sp) then
          write(log_buf, '(A,ES12.5,A,ES12.5)') &
               'Single precision metric error ', this%metric_perturb, &
               ', tolerance ', neko_metric_perturb_max
          call neko_log%warning(log_buf)
       else
          call neko_log%message('Single precision storage of the ' // &
               'geometric factors would exceed the error tolerance')
       end if
 
       if (real(neko_eps_sp, rp) * this%metric_scaled_cond .ge. &
            neko_eps * this%metric_cond) then
          call neko_log%message('Dominated by element skew')
       else
          call neko_log%message('Dominated by element aspect ratio')
       end if
    end if
 
    ! A metric that is not positive definite is a degenerate or inverted
    ! element, which is a mesh problem in any precision. Reported rather than
    ! fatal, since this is a diagnostic and the rest of the setup may still
    ! want to run
    if (this%metric_degenerate .gt. 0) then
       write(log_buf, '(A,I0)') &
            'Non positive definite metric at points: ', this%metric_degenerate
       call neko_log%error(log_buf)
    end if
 
  end subroutine coef_metric_condition
 
  subroutine coef_enable_geo_compression(this, tol)
    class(coef_t), intent(inout) :: this
    real(kind=rp), intent(in), optional :: tol
 
    if (.not. this%coef_metrics_initialized) then
       call neko_error('Geometric factors must be generated before they ' // &
            'are compressed')
    end if
 
    if (present(tol)) then
       if (tol .lt. 0.0_rp) then
          call neko_error('Geometric factor compression tolerance must ' // &
               'be non-negative')
       end if
    end if
 
    if (neko_bcknd_xsmm .eq. 1) then
       if (this%scope .eq. coef_full) then
          call neko_log%section('Compressed G')
          call neko_log%message('Not built for the XSMM backend')
          call neko_log%end_section()
       end if
       return
    end if
 
    this%geo_compression_tol = 1.0e-7_rp
    if (present(tol)) then
       this%geo_compression_tol = tol
    end if
 
    this%geo_compression = .true.
    call coef_generate_geo_compr(this, this%scope .eq. coef_full)
 
  end subroutine coef_enable_geo_compression
 
  subroutine coef_generate_geo_compr(c, verbose)
    type(coef_t), intent(inout) :: c
    logical, intent(in) :: verbose
    integer :: e, m, i, lxyz, m_max, ierr, n
    integer, allocatable :: c_inds_rev(:) ! reverse compression indices map
    integer :: stats(2), glb_stats(2)
    real(kind=rp) :: ctol, diff
    character(len=LOG_SIZE) :: log_buf
 
    call coef_free_geo_compr(c)
 
    ! On a device backend the host factors are only current after
    ! initialization, while recompute_metrics() updates the device ones
    if (neko_bcknd_device .eq. 1) then
       n = c%dof%size()
       call device_memcpy(c%G11, c%G11_d, n, device_to_host, sync = .false.)
       call device_memcpy(c%G22, c%G22_d, n, device_to_host, sync = .false.)
       call device_memcpy(c%G33, c%G33_d, n, device_to_host, sync = .false.)
       call device_memcpy(c%G12, c%G12_d, n, device_to_host, sync = .false.)
       call device_memcpy(c%G13, c%G13_d, n, device_to_host, sync = .false.)
       call device_memcpy(c%G23, c%G23_d, n, device_to_host, sync = .true.)
    end if
 
    ctol = c%geo_compression_tol
    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 = 0
    if (c%msh%nelv .gt. 0) then
       m_max = 1
       c%compression_inds(1) = 1
       c_inds_rev = 0
       c_inds_rev(1) = 1
    end if
 
    ! 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_rp*abs(c%G12(i,1,1,e) - c%G12(i,1,1,c_inds_rev(m))) &
                  + 2.0_rp*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_rp*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, written so that a NaN never matches either
       if (.not. (diff .le. ctol)) then
          m_max = m_max + 1
          c%compression_inds(e) = m_max
          c_inds_rev(m_max) = e
       end if
    end do
 
    ! Third step, allocate and fill Gij_compr objects
    allocate(c%G11_compr(c%Xh%lx, c%Xh%ly, c%Xh%lz, m_max))
    allocate(c%G22_compr(c%Xh%lx, c%Xh%ly, c%Xh%lz, m_max))
    allocate(c%G33_compr(c%Xh%lx, c%Xh%ly, c%Xh%lz, m_max))
    allocate(c%G12_compr(c%Xh%lx, c%Xh%ly, c%Xh%lz, m_max))
    allocate(c%G13_compr(c%Xh%lx, c%Xh%ly, c%Xh%lz, m_max))
    allocate(c%G23_compr(c%Xh%lx, c%Xh%ly, c%Xh%lz, m_max))
    do m = 1, m_max
       do i = 1, lxyz
          c%G11_compr(i,1,1,m) = c%G11(i,1,1,c_inds_rev(m))
          c%G22_compr(i,1,1,m) = c%G22(i,1,1,c_inds_rev(m))
          c%G33_compr(i,1,1,m) = c%G33(i,1,1,c_inds_rev(m))
          c%G12_compr(i,1,1,m) = c%G12(i,1,1,c_inds_rev(m))
          c%G13_compr(i,1,1,m) = c%G13(i,1,1,c_inds_rev(m))
          c%G23_compr(i,1,1,m) = c%G23(i,1,1,c_inds_rev(m))
       end do
    end do
 
    deallocate(c_inds_rev)
 
    if (neko_bcknd_device .eq. 1) then
       n = lxyz * m_max
       call device_map(c%G11_compr, c%G11_compr_d, n)
       call device_map(c%G22_compr, c%G22_compr_d, n)
       call device_map(c%G33_compr, c%G33_compr_d, n)
       call device_map(c%G12_compr, c%G12_compr_d, n)
       call device_map(c%G13_compr, c%G13_compr_d, n)
       call device_map(c%G23_compr, c%G23_compr_d, n)
       call device_map(c%compression_inds, c%compression_inds_d, c%msh%nelv)
       call device_memcpy(c%G11_compr, c%G11_compr_d, n, host_to_device, &
            sync = .false.)
       call device_memcpy(c%G22_compr, c%G22_compr_d, n, host_to_device, &
            sync = .false.)
       call device_memcpy(c%G33_compr, c%G33_compr_d, n, host_to_device, &
            sync = .false.)
       call device_memcpy(c%G12_compr, c%G12_compr_d, n, host_to_device, &
            sync = .false.)
       call device_memcpy(c%G13_compr, c%G13_compr_d, n, host_to_device, &
            sync = .false.)
       call device_memcpy(c%G23_compr, c%G23_compr_d, n, host_to_device, &
            sync = .false.)
       call device_memcpy(c%compression_inds, c%compression_inds_d, &
            c%msh%nelv, host_to_device, sync = .true.)
    end if
 
    if (verbose) then
       stats(1) = c%msh%nelv
       stats(2) = m_max
       call mpi_allreduce(stats, glb_stats, 2, mpi_integer, mpi_sum, &
            neko_comm, ierr)
 
       call neko_log%section('Compressed G')
       write(log_buf, '(A,I0)') 'Copies     : ', glb_stats(2)
       call neko_log%message(log_buf)
       write(log_buf, '(A,I0)') 'Elements   : ', glb_stats(1)
       call neko_log%message(log_buf)
       write(log_buf, '(A,ES12.5)') 'Tolerance  : ', ctol
       call neko_log%message(log_buf)
       call neko_log%end_section()
    end if
 
  end subroutine coef_generate_geo_compr
 
  subroutine coef_free_geo_compr(c)
    type(coef_t), intent(inout) :: c
 
    if (allocated(c%compression_inds)) then
       if (neko_bcknd_device .eq. 1) then
          call device_unmap(c%compression_inds, c%compression_inds_d)
       end if
       deallocate(c%compression_inds)
    end if
 
    if (allocated(c%G11_compr)) then
       if (neko_bcknd_device .eq. 1) then
          call device_unmap(c%G11_compr, c%G11_compr_d)
       end if
       deallocate(c%G11_compr)
    end if
 
    if (allocated(c%G22_compr)) then
       if (neko_bcknd_device .eq. 1) then
          call device_unmap(c%G22_compr, c%G22_compr_d)
       end if
       deallocate(c%G22_compr)
    end if
 
    if (allocated(c%G33_compr)) then
       if (neko_bcknd_device .eq. 1) then
          call device_unmap(c%G33_compr, c%G33_compr_d)
       end if
       deallocate(c%G33_compr)
    end if
 
    if (allocated(c%G12_compr)) then
       if (neko_bcknd_device .eq. 1) then
          call device_unmap(c%G12_compr, c%G12_compr_d)
       end if
       deallocate(c%G12_compr)
    end if
 
    if (allocated(c%G13_compr)) then
       if (neko_bcknd_device .eq. 1) then
          call device_unmap(c%G13_compr, c%G13_compr_d)
       end if
       deallocate(c%G13_compr)
    end if
 
    if (allocated(c%G23_compr)) then
       if (neko_bcknd_device .eq. 1) then
          call device_unmap(c%G23_compr, c%G23_compr_d)
       end if
       deallocate(c%G23_compr)
    end if
 
 
  end subroutine coef_free_geo_compr
 
  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
          ! Under COEF_OPERATOR the Binv transfer below is skipped, so this
          ! becomes the last one of the routine and has to synchronise.
          call device_memcpy(c%B, c%B_d, ntot, device_to_host, &
               sync = (c%scope .eq. coef_operator))
       end if
    else
       !$omp parallel do private(e, i)
       do e = 1, c%msh%nelv
          ! Here we need to handle things differently for axis symmetric elements
          do 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
       !$omp end parallel do
    end if
 
    ! Neither Binv nor the volume is read when only applying an operator,
    ! and assembling Binv costs a gather-scatter round on top of that.
    if (c%scope .ne. coef_full) return
 
    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
 
  subroutine coef_require_facets(this, who)
    class(coef_t), intent(in) :: this
    character(len=*), intent(in) :: who
 
    if (this%scope .ne. coef_full) then
       call neko_error(who // ' needs the facet areas and normals, which ' // &
            'a COEF_OPERATOR coef does not build')
    end if
 
  end subroutine coef_require_facets
 
  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
 
    ! Builds the rotation matrices from get_normal()
    if (this%scope .ne. coef_full) then
       call neko_error('Cyclic boundaries need the facet normals, ' // &
            'which COEF_OPERATOR does not build')
    end if
 
    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%x(i, j, k, pe) * un(1) &
                           - this%dof%x%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
 
    if (this%scope .ne. coef_full) then
       call neko_error('Rebuilding the geometry needs the derivative ' // &
            'arrays, which COEF_OPERATOR releases')
    end if
 
    call coef_generate_dxyzdrst(this)
    call coef_generate_geo(this)
    if (this%geo_compression) then
       call coef_generate_geo_compr(this, .false.)
    end if
    call coef_generate_area_and_normal(this)
    call coef_generate_mass(this)
    if (this%cyclic) then
       call coef_generate_cyclic_bc(this)
    end if
    this%metrics_version = this%metrics_version + 1
  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