atl_modg_nerplanck_kernel_module.f90 Source File


This file depends on

sourcefile~~atl_modg_nerplanck_kernel_module.f90~~EfferentGraph sourcefile~atl_modg_nerplanck_kernel_module.f90 atl_modg_nerplanck_kernel_module.f90 sourcefile~atl_cube_elem_module.f90 atl_cube_elem_module.f90 sourcefile~atl_modg_nerplanck_kernel_module.f90->sourcefile~atl_cube_elem_module.f90 sourcefile~atl_equation_module.f90 atl_equation_module.f90 sourcefile~atl_modg_nerplanck_kernel_module.f90->sourcefile~atl_equation_module.f90 sourcefile~atl_kerneldata_module.f90 atl_kerneldata_module.f90 sourcefile~atl_modg_nerplanck_kernel_module.f90->sourcefile~atl_kerneldata_module.f90 sourcefile~atl_nerplanck_flux_module.f90 atl_nerplanck_flux_module.f90 sourcefile~atl_modg_nerplanck_kernel_module.f90->sourcefile~atl_nerplanck_flux_module.f90 sourcefile~atl_scheme_module.f90 atl_scheme_module.f90 sourcefile~atl_modg_nerplanck_kernel_module.f90->sourcefile~atl_scheme_module.f90 sourcefile~ply_modg_basis_module.f90 ply_modg_basis_module.f90 sourcefile~atl_modg_nerplanck_kernel_module.f90->sourcefile~ply_modg_basis_module.f90 sourcefile~atl_eqn_acoustic_module.f90 atl_eqn_acoustic_module.f90 sourcefile~atl_equation_module.f90->sourcefile~atl_eqn_acoustic_module.f90 sourcefile~atl_eqn_advection_1d_module.f90 atl_eqn_advection_1d_module.f90 sourcefile~atl_equation_module.f90->sourcefile~atl_eqn_advection_1d_module.f90 sourcefile~atl_eqn_bbm_module.f90 atl_eqn_bbm_module.f90 sourcefile~atl_equation_module.f90->sourcefile~atl_eqn_bbm_module.f90 sourcefile~atl_eqn_euler_module.f90 atl_eqn_euler_module.f90 sourcefile~atl_equation_module.f90->sourcefile~atl_eqn_euler_module.f90 sourcefile~atl_eqn_heat_module.f90 atl_eqn_heat_module.f90 sourcefile~atl_equation_module.f90->sourcefile~atl_eqn_heat_module.f90 sourcefile~atl_eqn_lineareuler_module.f90 atl_eqn_linearEuler_module.f90 sourcefile~atl_equation_module.f90->sourcefile~atl_eqn_lineareuler_module.f90 sourcefile~atl_eqn_maxwell_module.f90 atl_eqn_maxwell_module.f90 sourcefile~atl_equation_module.f90->sourcefile~atl_eqn_maxwell_module.f90 sourcefile~atl_eqn_nerplanck_module.f90 atl_eqn_nerplanck_module.f90 sourcefile~atl_equation_module.f90->sourcefile~atl_eqn_nerplanck_module.f90 sourcefile~atl_eqn_nvrstk_module.f90 atl_eqn_nvrstk_module.f90 sourcefile~atl_equation_module.f90->sourcefile~atl_eqn_nvrstk_module.f90 sourcefile~atl_materialfun_module.f90 atl_materialFun_module.f90 sourcefile~atl_equation_module.f90->sourcefile~atl_materialfun_module.f90 sourcefile~ply_dof_module.f90 ply_dof_module.f90 sourcefile~atl_kerneldata_module.f90->sourcefile~ply_dof_module.f90 sourcefile~atl_scheme_module.f90->sourcefile~ply_modg_basis_module.f90 sourcefile~atl_modg_1d_scheme_module.f90 atl_modg_1d_scheme_module.f90 sourcefile~atl_scheme_module.f90->sourcefile~atl_modg_1d_scheme_module.f90 sourcefile~atl_modg_2d_scheme_module.f90 atl_modg_2d_scheme_module.f90 sourcefile~atl_scheme_module.f90->sourcefile~atl_modg_2d_scheme_module.f90 sourcefile~atl_modg_scheme_module.f90 atl_modg_scheme_module.f90 sourcefile~atl_scheme_module.f90->sourcefile~atl_modg_scheme_module.f90 sourcefile~atl_stabilization_module.f90 atl_stabilization_module.f90 sourcefile~atl_scheme_module.f90->sourcefile~atl_stabilization_module.f90 sourcefile~ply_modg_basis_module.f90->sourcefile~ply_dof_module.f90 sourcefile~ply_space_integration_module.f90 ply_space_integration_module.f90 sourcefile~ply_modg_basis_module.f90->sourcefile~ply_space_integration_module.f90 sourcefile~atl_eqn_acoustic_module.f90->sourcefile~atl_materialfun_module.f90 sourcefile~atl_eqn_heat_module.f90->sourcefile~atl_materialfun_module.f90 sourcefile~atl_eqn_lineareuler_module.f90->sourcefile~atl_materialfun_module.f90 sourcefile~atl_eqn_nvrstk_module.f90->sourcefile~atl_eqn_euler_module.f90 sourcefile~atl_modg_1d_scheme_module.f90->sourcefile~ply_dof_module.f90 sourcefile~atl_modg_2d_scheme_module.f90->sourcefile~ply_dof_module.f90 sourcefile~atl_modg_scheme_module.f90->sourcefile~ply_dof_module.f90 sourcefile~atl_cons_positivity_preserv_module.f90 atl_cons_positivity_preserv_module.f90 sourcefile~atl_stabilization_module.f90->sourcefile~atl_cons_positivity_preserv_module.f90 sourcefile~atl_covolume_module.f90 atl_covolume_module.f90 sourcefile~atl_stabilization_module.f90->sourcefile~atl_covolume_module.f90 sourcefile~atl_positivity_preserv_module.f90 atl_positivity_preserv_module.f90 sourcefile~atl_stabilization_module.f90->sourcefile~atl_positivity_preserv_module.f90 sourcefile~atl_spectral_viscosity_module.f90 atl_spectral_viscosity_module.f90 sourcefile~atl_stabilization_module.f90->sourcefile~atl_spectral_viscosity_module.f90 sourcefile~atl_covolume_module.f90->sourcefile~atl_spectral_viscosity_module.f90

Source Code

! Copyright (c) 2012 Jan Hueckelheim <j.hueckelheim@grs-sim.de>
! Copyright (c) 2012-2015 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2012 Jens Zudrop <j.zudrop@grs-sim.de>
! Copyright (c) 2012 Melven Zoellner <yameta@freenet.de>
! Copyright (c) 2013-2018, 2020 Peter Vitt <peter.vitt2@uni-siegen.de>
! Copyright (c) 2014 Verena Krupp <verena.krupp@uni-siegen.de>
! Copyright (c) 2014, 2016 Nikhil Anand <nikhil.anand@uni-siegen.de>
! Copyright (c) 2016 Tobias Girresser <tobias.girresser@student.uni-siegen.de>
! Copyright (c) 2017 Daniel PetrĂ³ <daniel.petro@student.uni-siegen.de>
!
! Permission to use, copy, modify, and distribute this software for any
! purpose with or without fee is hereby granted, provided that the above
! copyright notice and this permission notice appear in all copies.
!
! THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHORS DISCLAIM ALL WARRANTIES
! WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
! MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR
! ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
! WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
! ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
! OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
! **************************************************************************** !

! Copyright (c) 2014,2016-2017 Peter Vitt <peter.vitt2@uni-siegen.de>
! Copyright (c) 2014 Harald Klimach <harald.klimach@uni-siegen.de>
!
! Parts of this file were written by Peter Vitt and Harald Klimach for
! University of Siegen.
!
! Permission to use, copy, modify, and distribute this software for any
! purpose with or without fee is hereby granted, provided that the above
! copyright notice and this permission notice appear in all copies.
!
! THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHORS DISCLAIM ALL WARRANTIES
! WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
! MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR
! ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
! WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
! ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
! OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
! **************************************************************************** !
!
! Return the position of a given ansatz function combination in the
! linearized list of modal coefficients for Q-Tensor product polynomials.
! You must provide
! * Ansatzfunction index in x direction. Index starts with 1.
! * Ansatzfunction index in y direction. Index starts with 1.
! * Ansatzfunction index in z direction. Index starts with 1.
! * The maximal polynomial degree per spatial direction.
! * The variable to store the position of the modal coefficient in the list of
!   modal coefficients in.


! Return the position of a given ansatz function combination in the
! linearized list of modal coefficients for Q-Tensor product polynomials.
! You must provide
! * Ansatzfunction index in x direction. Index starts with 1.
! * Ansatzfunction index in y direction. Index starts with 1.
! * The maximal polynomial degree per spatial direction.
! * The variable to store the position of the modal coefficient in the list of
!   modal coefficients in.


! Return the position of a given ansatz function combination in the
! linearized list of modal coefficients for Q-Tensor product polynomials.
! You must provide
! * Ansatzfunction index in x direction. Index starts with 1.
! * The variable to store the position of the modal coefficient in the list of
!   modal coefficients in.


! Return the position of a given ansatz function combination in the
! linearized list of modal coefficients for P-Tensor product polynomials.
! You must provide
! * Ansatzfunction index in x direction. Index starts with 1.
! * Ansatzfunction index in y direction. Index starts with 1.
! * Ansatzfunction index in z direction. Index starts with 1.
! * The maximal polynomial degree per spatial direction.
! * The variable to store the position of the modal coefficient in the list of
!   modal coefficients in.


! Return the position of a given ansatz function combination in the
! linearized list of modal coefficients for P-Tensor product polynomials.
! You must provide
! * Ansatzfunction index in x direction. Index starts with 1.
! * Ansatzfunction index in y direction. Index starts with 1.
! * The maximal polynomial degree per spatial direction.
! * The variable to store the position of the modal coefficient in the list of
!   modal coefficients in.


! Return the position of a given ansatz function combination in the
! linearized list of modal coefficients for P-Tensor product polynomials.
! You must provide
! * Ansatzfunction index in x direction. Index starts with 1.
! * The variable to store the position of the modal coefficient in the list of
!   modal coefficients in.


! Return the number of degrees of freedom for Q polynomial space
! Your must provide:
! * The maximal polynomial degree per spatial direction
! * The variable to store the number of degrees of freedom for a Q tensor
!   product polynomial


! Return the number of degrees of freedom for broken polynomial space
! Your must provide:
! * The maximal polynomial degree per spatial direction (for P Tensor product
!   polynomials this assumed to be the same for each spatial direction).
! * The variable to store the number of degrees of freedom for a P tensor
!   product polynomial


! Return the number of degrees of freedom for Q polynomial space
! You must provide:
! * The maximal polynomial degree per spatial direction
! * The variable to store the number of degrees of freedom for a Q tensor
!   product polynomial


! Return the number of degrees of freedom for broken polynomial space
! You must provide:
! * The maximal polynomial degree per spatial direction (for P Tensor product
!   polynomials this assumed to be the same for each spatial direction).
! * A variable to store the number of degrees of freedom for a P tensor product
!   polynomial


! Return the number of degrees of freedom for Q polynomial space
! You must provide:
! * The maximal polynomial degree per spatial direction
! * The variable to store the number of degrees of freedom for a Q tensor
!   product polynomial


! Return the number of degrees of freedom for broken polynomial space
! You must provide:
! * The maximal polynomial degree per spatial direction (for P Tensor product
!   polynomials this assumed to be the same for each spatial direction).
! * The variable to store the number of degrees of freedom for a P tensor
!   product polynomial

! The x, y and z ansatz degrees are turned into the degrees of the next
! ansatz function in the layered P list
! You must provide:
! * Ansatz function index in x direction. First ansatz function has index 1.
! * Ansatz function index in y direction. First ansatz function has index 1.
! * Ansatz function index in z direction. First ansatz function has index 1.

! The x and y ansatz degrees are turned into the degrees of the next
! ansatz function in the layered P list
! You must provide:
! * Ansatz function index in x direction. First ansatz function has index 1.
! * Ansatz function index in y direction. First ansatz function has index 1.

! The x ansatz degree is turned into the degree of the next
! ansatz function in the layered P list
! You must provide:
! * Ansatz function index in x direction. First ansatz function has index 1.

! The x, y and z ansatz degrees are turned into the degrees of the next
! ansatz function in the linearized Q tensor
! You must provide:
! * Ansatz function index in x direction. First ansatz function has index 1.
! * Ansatz function index in y direction. First ansatz function has index 1.
! * Ansatz function index in z direction. First ansatz function has index 1.
! * Maximal polynomial degree

! The x and y ansatz degrees are turned into the degrees of the next
! ansatz function in the linearized Q tensor
! You must provide:
! * Ansatz function index in x direction. First ansatz function has index 1.
! * Ansatz function index in y direction. First ansatz function has index 1.
! * Maximal polynomial degree

! The x ansatz degree is turned into the degree of the next
! ansatz function in the linearized Q tensor
! You must provide:
! * Ansatz function index in x direction. First ansatz function has index 1.
!> Module for routines and datatypes of MOdal Discontinuous Galerkin (MODG)
!! scheme. This scheme is a spectral scheme for linear, purley hyperbolic
!! partial differential equation systems.
module atl_modg_nerplanck_kernel_module
  use env_module,               only: rk, long_k
  use tem_element_module,       only: eT_fluid
  use tem_element_module,       only: eT_fluid

  use atl_equation_module,       only: atl_equations_type
  use atl_kerneldata_module,     only: atl_kerneldata_type, atl_statedata_type
  use atl_cube_elem_module,      only: atl_cube_elem_type
  use atl_nerplanck_flux_module, only: atl_nerplanck_numflux_solve,      &
    &                                  atl_nerplanck_physflux_solve,     &
    &                                  atl_nerplanck_numflux_preprocess, &
    &                                  atl_nerplanck_physflux_preprocess
  use atl_scheme_module,         only: atl_modg_scheme_type
  use ply_modg_basis_module,     only: ply_faceValLeftBndAns, &
    &                                  ply_faceValLeftBndTest, &
    &                                  ply_faceValRightBndTest, &
    &                                  ply_scalProdDualLeg,     &
    &                                  ply_scalProdDualLegDiff

  implicit none
  private

  public :: atl_preprocess_modg_nerplanck_kernel

  !> interface for preprocessing the data for the kernel
  interface atl_preprocess_modg_nerplanck_kernel
    module procedure preprocess_modg_nerplanck_kernel
  end interface atl_preprocess_modg_nerplanck_kernel

contains

  subroutine preprocess_modg_nerplanck_kernel( equation, kerneldata, statedata, &
                                   & mesh, scheme)
    ! --------------------------------------------------------------------------
    !> The mesh you are working with.
    type(atl_cube_elem_type) :: mesh
    !> The equation you solve.
    type(atl_equations_type) :: equation
    !> The data of the kernel.
    type(atl_kerneldata_type) :: kerneldata
    !> THe state if the equation
    type(atl_statedata_type) :: statedata
    !> Parameters of the modal dg scheme
    type(atl_modg_scheme_type) :: scheme
    ! --------------------------------------------------------------------------

    ! Initialize the workspace for the modg kernel with zeros. The following routines
    ! will add the influence of their projections to this array.
    kerneldata%state_der = 0.0_rk

    ! Loop over all the elements and calculate the projection of
    ! the cell internal physical flux. Therefore, we do the following for each cell:
    ! First, we loop over the modal dofs and create a modal representation of the
    ! physical flux.
    ! Second, we loop over the test functions and project the modal representation
    ! of the physical flux onto them.

    ! add physical flux calculation.
    call modg_nerplanck_physFlux_pre( mesh, equation, statedata, scheme )

    ! Loop over all the faces and calculate the projection of the
    ! numerical flux. Therefore, we do the following for each faces:
    ! First, we loop over all modal dofs and calculate the modal representation
    ! of the flux on this face.
    ! Second, we loop over all test functions and project the modal representation
    ! of the numerical flux for this face onto the test functions.

    ! ... numerical flux across faces in x direction
    call modg_nerplanck_x_numFlux_pre( mesh, equation, statedata, scheme )

    ! ... numerical flux across faces in y direction
    call modg_nerplanck_y_numFlux_pre( mesh, equation, statedata, scheme )

    ! ... numerical flux across faces in z direction
    call modg_nerplanck_z_numFlux_pre( mesh, equation, statedata, scheme )

    ! Finally, we loop over all the cells and apply the cell local inverse mass matrix.
    ! All of the previous steps and this one include the correct Determinant of the
    ! Jacobians of the mappings from reference element/face to the physical element/face.
    ! It also includes the elemental timestepping of the time_integration_module.
    call modg_invMassMatrix_pre( mesh, scheme )

  end subroutine preprocess_modg_nerplanck_kernel


  !> Numerical flux calculation for Nernst-Planck equation across the faces in X direction.
  subroutine modg_nerplanck_x_numFlux_pre( mesh, equation, statedata, scheme )
    ! --------------------------------------------------------------------------
    !> The mesh you are working with.
    type(atl_cube_elem_type) :: mesh
    !> The equation you solve.
    type(atl_equations_type) :: equation
    !> THe state if the equation
    type(atl_statedata_type) :: statedata
    !> Parameters of the modal dg scheme
    type(atl_modg_scheme_type), intent(in) :: scheme
    ! --------------------------------------------------------------------------!
    ! Loop vars
    integer :: xAnsFunc, yAnsFunc, zAnsFunc, xTestFunc, yTestFunc, zTestFunc
    ! Positions for the modal coefficients
    integer :: modalCoeff, testPos, ansPos
    ! Modal coefficients for elements left and right of the face
    real(kind=rk) :: leftModalCoeff(4), rightModalCoeff(4), flux(3)
    ! Modal representation of the flux on the face
    real(kind=rk), allocatable :: numFluxBuffer(:,:)
    ! The values of the x ansatz and test functions at the face from left and right
    ! limit.
    real(kind=rk) :: faceValLeft, faceValRight
    ! Loop var for the faces.
    integer :: xside
    ! Element positions of the left and right element of the face.
    integer(kind=long_k) :: left_neighbor, right_neighbor
    ! Speed of light, i.e. maximum information propagation speed.
    real(kind=rk) :: diffusivity
    ! Scalar products between y and z ansatz and test functions.
    real(kind=rk) :: yScalProd, zScalProd
    ! Determinant of the jacobian of the mapping from reference face to
    ! the current face.
    real(kind=rk) :: jacobiDet
    ! The minimum indices for the relevant ansatz functions.
    integer :: zAnsFuncMin, yAnsFuncMin
    ! --------------------------------------------------------------------------

    allocate(numFluxBuffer((scheme%maxPolyDegree+1)**3,3))

    ! Calculate the speed of light
    diffusivity = equation%nerplanck%D

    ! Build the determinant of the Jacobian of the mapping from reference
    ! face to the physical face. Since we are considering cubes it is
    ! the same everywhere.
    jacobiDet = (mesh%length/2.0_rk)**2.0_rk

    ! Loop over all fluid the faces in x direction
    XFaceLoop: do xside = 1, size(mesh%faces%faces(1)%computeFace%leftPos)
      ! Get the fluid neighbors for this face.
      left_neighbor = mesh%faces%faces(1)%computeFace%leftPos(xside)
      right_neighbor = mesh%faces%faces(1)%computeFace%rightPos(xside)

      ! Loop over all the ansatz functions. Here, we use simple tensor products
      ! of one-dimensional ansatz functions.
      do xAnsFunc = 1, scheme%maxPolyDegree+1
        faceValRight = ply_faceValLeftBndAns(xAnsFunc)
        do yAnsFunc = 1, (scheme%maxPolyDegree+1)
          do zAnsFunc = 1, (scheme%maxPolyDegree+1)
            ! Get the modal coefficients of the left and right ansatz functions and
            ! calculate the modal flux with them.
  modalcoeff = xansfunc                                      &
    &      + ( ( yansfunc-1)                             &
    &      + (zansfunc-1)*(scheme%maxpolydegree+1))*(scheme%maxpolydegree+1)
            leftModalCoeff(:) = statedata%state(left_neighbor, modalCoeff, :)
            rightModalCoeff(:) = statedata%state(right_neighbor, modalCoeff, :) * faceValRight
            call atl_nerplanck_numflux_preprocess( &
              & left        = leftModalCoeff,      &
              & right       = rightModalCoeff,     &
              & diffusivity = diffusivity,         &
              & flux        = flux                 )
            ! Buffer the modal values of the numerical flux
            numFluxBuffer( modalCoeff, :) = flux
          end do
        end do
      end do

      ! Loop over all the test functions and project the numerical flux to them.
      do xTestFunc = 1,min(scheme%maxPolyDegree+1,2)
        faceValLeft = ply_faceValRightBndTest(xTestFunc)
        faceValRight = ply_faceValLeftBndTest(xTestFunc)
        do yTestFunc = 1,scheme%maxPolyDegree+1
          do zTestFunc = 1,scheme%maxPolyDegree+1

            ! position of the test functions in the kerneldata
  testpos = xtestfunc                                      &
    &      + ( ( ytestfunc-1)                             &
    &      + (ztestfunc-1)*(scheme%maxpolydegree+1))*(scheme%maxpolydegree+1)

            ! get the relevant indices of ansatz functions for the projection
            do xAnsFunc = 1, scheme%maxPolyDegree+1
              yAnsFuncMin = yTestFunc-2
              if( yAnsFuncMin < 1 ) then
                yAnsFuncMin = yTestFunc
              end if
              do yAnsFunc = yAnsFuncMin,yTestFunc,2
                zAnsFuncMin = zTestFunc-2
                if( zAnsFuncMin < 1 ) then
                  zAnsFuncMin = zTestFunc
                end if
                do zAnsFunc = zAnsFuncMin,zTestFunc,2

                  ! the position of the modal coefficeint of this ansatz functions
  anspos = xansfunc                                      &
    &      + ( ( yansfunc-1)                             &
    &      + (zansfunc-1)*(scheme%maxpolydegree+1))*(scheme%maxpolydegree+1)

                  ! calculate the projection of the ansatz and test function
                  yScalProd = ply_scalProdDualLeg(yAnsFunc, yTestFunc)
                  zScalProd = ply_scalProdDualLeg(zAnsFunc, zTestFunc)

!NA!                  ! buffer the result in kernel data and take care of the outer surface unit normal
!NA!                  ! ... for the element left of the face
!NA!                  outerNormalLeft = 1.0_rk
!NA!                  projectedFlux = outerNormalLeft * numFluxBuffer(ansPos,:)
!NA!                  kernelData%state_der(left_neighbor, testPos, 1, 2:4) &
!NA!                            & = kernelData%state_der(left_neighbor, testPos, 1, 2:4) &
!NA!                            & - faceValLeft * yScalProd * zScalProd &
!NA!                            & * projectedFlux * jacobiDet
!NA!                  ! ... for the element right of the face
!NA!                  outerNormalRight = -1.0_rk
!NA!                  projectedFlux = outerNormalRight * numFluxBuffer(ansPos,:)
!NA!                  kernelData%state_der(right_neighbor, testPos,1,2:4) &
!NA!                            & = kernelData%state_der(right_neighbor, testPos, 1, 2:4) &
!NA!                            & - faceValRight * yScalProd * zScalProd &
!NA!                            & * projectedFlux * jacobiDet

                end do
              end do
            end do

          end do
        end do
      end do

    end do XFaceLoop

    !GhostFaceLoop: do xside = 1, mesh%faces(1)%ghostfromfiner%leftTreeId%nVals


  end subroutine modg_nerplanck_x_numFlux_pre

  !> Numerical flux calculation for Nernst-Planck equation across the faces in Y direction.
  subroutine modg_nerplanck_y_numFlux_pre( mesh, equation, statedata, scheme )
    ! --------------------------------------------------------------------------
    !> The mesh you are working with.
    type(atl_cube_elem_type) :: mesh
    !> The equation you solve.
    type(atl_equations_type) :: equation
    !> THe state if the equation
    type(atl_statedata_type) :: statedata
    !> Parameters of the modal dg scheme
    type(atl_modg_scheme_type), intent(in) :: scheme
    ! --------------------------------------------------------------------------!
    ! Loop vars
    integer :: xAnsFunc, yAnsFunc, zAnsFunc, xTestFunc, yTestFunc, zTestFunc
    ! Positions for the modal coefficients
    integer :: modalCoeff, testPos, ansPos
    ! Modal coefficients for elements left and right of the face
    real(kind=rk) :: leftModalCoeff(4), rightModalCoeff(4), flux(3)
    ! Modal representation of the flux on the face
    real(kind=rk), allocatable :: numFluxBuffer(:,:)
    ! The values of the x ansatz and test functions at the face from left and right
    ! limit.
    real(kind=rk) :: faceValLeft, faceValRight
    ! Loop var for the faces.
    integer :: yside
    ! Element positions of the left and right element of the face.
    integer(kind=long_k) :: left_neighbor, right_neighbor
    ! Speed of light, i.e. maximum information propagation speed.
    real(kind=rk) :: diffusivity
    ! Scalar products between x and z ansatz and test functions.
    real(kind=rk) :: xScalProd, zScalProd
    ! Determinant of the jacobian of the mapping from reference face to
    ! the current face.
    real(kind=rk) :: jacobiDet
    ! The minimum indices for the relevant ansatz functions.
    integer :: zAnsFuncMin, xAnsFuncMin
    ! Rotation for the flux
    integer :: varRotation(4)
    ! --------------------------------------------------------------------------

    allocate(numFluxBuffer((scheme%maxPolyDegree+1)**3,3))

    ! Calculate the speed of light
    diffusivity = equation%nerplanck%D

    ! Build the determinant of the Jacobian of the mapping from reference
    ! face to the physical face. Since we are considering cubes it is
    ! the same everywhere.
    jacobiDet = (mesh%length/2.0_rk)**2.0_rk

    ! Loop over fluid the faces in y direction
    varRotation = equation%varRotation(2)%varTransformIndices(:)
    yFaceLoop: do yside = 1,size(mesh%faces%faces(2)%computeFace%leftPos)
      ! Get the fluid neighbors for this face.
      left_neighbor = mesh%faces%faces(2)%computeFace%leftPos(yside)
      right_neighbor = mesh%faces%faces(2)%computeFace%rightPos(yside)


      ! VK: old face descriptor
      ! Get the neighbors for this face.
      !left_neighbor  = mesh%side%ydir%left_elemid(yside)
      !right_neighbor = mesh%side%ydir%right_elemid(yside)

      ! Loop over all the ansatz functions. Here, we use simple tensor products
      ! of one-dimensional ansatz functions.
      do yAnsFunc = 1, scheme%maxPolyDegree+1
        faceValRight = ply_faceValLeftBndAns(yAnsFunc)
        do xAnsFunc = 1, (scheme%maxPolyDegree+1)
          do zAnsFunc = 1, (scheme%maxPolyDegree+1)
            ! Get the modal coefficients of the left and right ansatz functions and
            ! calculate the modal flux with them.
  modalcoeff = xansfunc                                      &
    &      + ( ( yansfunc-1)                             &
    &      + (zansfunc-1)*(scheme%maxpolydegree+1))*(scheme%maxpolydegree+1)
            ! ... keep in mind that we have to rotate the input and output of the flux
            leftModalCoeff(:) = statedata%state(left_neighbor, &
                                               & modalCoeff, varRotation)
            rightModalCoeff(:) = statedata%state(right_neighbor, &
                                               &  modalCoeff, varRotation) * faceValRight
            call atl_nerplanck_numflux_preprocess( &
              & left        = leftModalCoeff,      &
              & right       = rightModalCoeff,     &
              & diffusivity = diffusivity,         &
              & flux        = flux                 )
            ! Buffer the modal values of the numerical flux
            numFluxBuffer( modalCoeff, :) = flux
          end do
        end do
      end do

      ! Loop over all the test functions and project the numerical flux to them.
      do yTestFunc = 1,min(scheme%maxPolyDegree+1,2)
        faceValLeft = ply_faceValRightBndTest(yTestFunc)
        faceValRight = ply_faceValLeftBndTest(yTestFunc)
        do xTestFunc = 1,scheme%maxPolyDegree+1
          do zTestFunc = 1,scheme%maxPolyDegree+1

            ! position of the test functions in the kerneldata
  testpos = xtestfunc                                      &
    &      + ( ( ytestfunc-1)                             &
    &      + (ztestfunc-1)*(scheme%maxpolydegree+1))*(scheme%maxpolydegree+1)

            ! get the relevant indices of ansatz functions for the projection
            do yAnsFunc = 1, scheme%maxPolyDegree+1
              xAnsFuncMin = xTestFunc-2
              if( xAnsFuncMin < 1 ) then
                xAnsFuncMin = xTestFunc
              end if
              do xAnsFunc = xAnsFuncMin,xTestFunc,2
                zAnsFuncMin = zTestFunc-2
                if( zAnsFuncMin < 1 ) then
                  zAnsFuncMin = zTestFunc
                end if
                do zAnsFunc = zAnsFuncMin,zTestFunc,2

                  ! the position of the modal coefficeint of this ansatz functions
  anspos = xansfunc                                      &
    &      + ( ( yansfunc-1)                             &
    &      + (zansfunc-1)*(scheme%maxpolydegree+1))*(scheme%maxpolydegree+1)

                  ! calculate the projection of the ansatz and test function
                  xScalProd = ply_scalProdDualLeg(xAnsFunc, xTestFunc)
                  zScalProd = ply_scalProdDualLeg(zAnsFunc, zTestFunc)

!NA!                  ! buffer the result in kernel data and take care of the outer surface unit normal
!NA!                  ! ... for the element left of the face
!NA!                  outerNormalLeft = 1.0_rk
!NA!                  projectedFlux = outerNormalLeft * numFluxBuffer(ansPos,:)
!NA!                  kernelData%state_der(left_neighbor, testPos, 1, 2:4) &
!NA!                            & = kernelData%state_der(left_neighbor, testPos, 1, 2:4) &
!NA!                            & - faceValLeft * xScalProd * zScalProd &
!NA!                            & * projectedFlux * jacobiDet
!NA!                  ! ... for the element right of the face
!NA!                  outerNormalRight = -1.0_rk
!NA!                  projectedFlux = outerNormalRight * numFluxBuffer(ansPos,:)
!NA!                  kernelData%state_der(right_neighbor, testPos,1,2:4) &
!NA!                            & = kernelData%state_der(right_neighbor, testPos, 1, 2:4) &
!NA!                            & - faceValRight * xScalProd * zScalProd &
!NA!                            & * projectedFlux * jacobiDet

                end do
              end do
            end do

          end do
        end do
      end do

    end do YFaceLoop


  end subroutine modg_nerplanck_y_numFlux_pre

  !> Numerical flux calculation for Nernst-Planck equation across the faces in Z direction.
  subroutine modg_nerplanck_z_numFlux_pre( mesh, equation, statedata, scheme )
    ! --------------------------------------------------------------------------
    !> The mesh you are working with.
    type(atl_cube_elem_type) :: mesh
    !> The equation you solve.
    type(atl_equations_type) :: equation
    !> THe state if the equation
    type(atl_statedata_type) :: statedata
    !> Parameters of the modal dg scheme
    type(atl_modg_scheme_type), intent(in) :: scheme
    ! --------------------------------------------------------------------------!
    ! Loop vars
    integer :: xAnsFunc, yAnsFunc, zAnsFunc, xTestFunc, yTestFunc, zTestFunc
    ! Positions for the modal coefficients
    integer :: modalCoeff, testPos, ansPos
    ! Modal coefficients for elements left and right of the face
    real(kind=rk) :: leftModalCoeff(4), rightModalCoeff(4), flux(3)
    ! Modal representation of the flux on the face
    real(kind=rk), allocatable :: numFluxBuffer(:,:)
    ! The values of the x ansatz and test functions at the face from left and right
    ! limit.
    real(kind=rk) :: faceValLeft, faceValRight
    ! Loop var for the faces.
    integer :: zside
    ! Element positions of the left and right element of the face.
    integer(kind=long_k) :: left_neighbor, right_neighbor
    ! Speed of light, i.e. maximum information propagation speed.
    real(kind=rk) :: diffusivity
    ! Scalar products between x and z ansatz and test functions.
    real(kind=rk) :: xScalProd, yScalProd
    ! Determinant of the jacobian of the mapping from reference face to
    ! the current face.
    real(kind=rk) :: jacobiDet
    ! The minimum indices for the relevant ansatz functions.
    integer :: xAnsFuncMin, yAnsFuncMin
    ! Rotation for the flux
    integer :: varRotation(4)
    ! --------------------------------------------------------------------------

    allocate(numFluxBuffer((scheme%maxPolyDegree+1)**3,3))

    ! Calculate the speed of light
    diffusivity = equation%nerplanck%D

    ! Build the determinant of the Jacobian of the mapping from reference
    ! face to the physical face. Since we are considering cubes it is
    ! the same everywhere.
    jacobiDet = (mesh%length/2.0_rk)**2.0_rk

    ! Loop over all fluid faces in z direction
    varRotation = equation%varRotation(3)%varTransformIndices(:)
    zFaceLoop: do zside = 1,size(mesh%faces%faces(3)%computeFace%leftPos)
      ! Get the fluid neighbors for this face.
      left_neighbor = mesh%faces%faces(3)%computeFace%leftPos(zside)
      right_neighbor = mesh%faces%faces(3)%computeFace%rightPos(zside)


      ! Get the neighbors for this face.
      !left_neighbor  = mesh%side%zdir%left_elemid(zside)
      !right_neighbor = mesh%side%zdir%right_elemid(zside)

      ! Loop over all the ansatz functions. Here, we use simple tensor products
      ! of one-dimensional ansatz functions.
      do zAnsFunc = 1, scheme%maxPolyDegree+1
        faceValRight = ply_faceValLeftBndAns(zAnsFunc)
        do xAnsFunc = 1, (scheme%maxPolyDegree+1)
          do yAnsFunc = 1, (scheme%maxPolyDegree+1)
            ! Get the modal coefficients of the left and right ansatz functions and
            ! calculate the modal flux with them.
  modalcoeff = xansfunc                                      &
    &      + ( ( yansfunc-1)                             &
    &      + (zansfunc-1)*(scheme%maxpolydegree+1))*(scheme%maxpolydegree+1)
            ! ... keep in mind that we have to rotate the input and output of the flux
            leftModalCoeff(:) = statedata%state(left_neighbor, &
                                               & modalCoeff, varRotation)
            rightModalCoeff(:) = statedata%state(right_neighbor, &
                                               &  modalCoeff, varRotation) * faceValRight
            call atl_nerplanck_numflux_preprocess( &
              & left        = leftModalCoeff,      &
              & right       = rightModalCoeff,     &
              & diffusivity = diffusivity,         &
              & flux        = flux                 )

            ! Buffer the modal values of the numerical flux
            numFluxBuffer( modalCoeff, :) = flux
          end do
        end do
      end do

      ! Loop over all the test functions and project the numerical flux to them.
      do zTestFunc = 1,min(scheme%maxPolyDegree+1,2)
        faceValLeft = ply_faceValRightBndTest(zTestFunc)
        faceValRight = ply_faceValLeftBndTest(zTestFunc)
        do xTestFunc = 1,scheme%maxPolyDegree+1
          do yTestFunc = 1,scheme%maxPolyDegree+1

            ! position of the test functions in the kerneldata
  testpos = xtestfunc                                      &
    &      + ( ( ytestfunc-1)                             &
    &      + (ztestfunc-1)*(scheme%maxpolydegree+1))*(scheme%maxpolydegree+1)

            ! get the relevant indices of ansatz functions for the projection
            do zAnsFunc = 1, scheme%maxPolyDegree+1
              xAnsFuncMin = xTestFunc-2
              if( xAnsFuncMin < 1 ) then
                xAnsFuncMin = xTestFunc
              end if
              do xAnsFunc = xAnsFuncMin,xTestFunc,2
                yAnsFuncMin = yTestFunc-2
                if( yAnsFuncMin < 1 ) then
                  yAnsFuncMin = yTestFunc
                end if
                do yAnsFunc = yAnsFuncMin,yTestFunc,2

                  ! the position of the modal coefficeint of this ansatz functions
  anspos = xansfunc                                      &
    &      + ( ( yansfunc-1)                             &
    &      + (zansfunc-1)*(scheme%maxpolydegree+1))*(scheme%maxpolydegree+1)

                  ! calculate the projection of the ansatz and test function
                  xScalProd = ply_scalProdDualLeg(xAnsFunc, xTestFunc)
                  yScalProd = ply_scalProdDualLeg(yAnsFunc, yTestFunc)

!NA!                  ! buffer the result in kernel data and take care of the outer surface unit normal
!NA!                  ! ... for the element left of the face
!NA!                  outerNormalLeft = 1.0_rk
!NA!                  projectedFlux = outerNormalLeft * numFluxBuffer(ansPos,:)
!NA!                  kernelData%state_der(left_neighbor, testPos, 1, 2:4) &
!NA!                                  & = kernelData%state_der(left_neighbor, testPos, 1, 2:4) &
!NA!                                  & - faceValLeft * xScalProd * yScalProd &
!NA!                                  & * projectedFlux * jacobiDet
!NA!                  ! ... for the element right of the face
!NA!                  outerNormalRight = -1.0_rk
!NA!                  projectedFlux = outerNormalRight * numFluxBuffer(ansPos,:)
!NA!                  kernelData%state_der(right_neighbor, testPos,1, 2:4) &
!NA!                                  & = kernelData%state_der(right_neighbor, testPos, 1, 2:4) &
!NA!                                  & - faceValRight * xScalProd * yScalProd &
!NA!                                  & * projectedFlux * jacobiDet

                end do
              end do
            end do

          end do
        end do
      end do

    end do ZFaceLoop
  end subroutine modg_nerplanck_z_numFlux_pre


  !> Calculate the projection of the physical flux for the MODG scheme and
  !! Nernst-Planck equation.
  subroutine modg_nerplanck_physFlux_pre( mesh, equation, statedata, scheme )
    ! --------------------------------------------------------------------------
    !> The mesh you are working with.
    type(atl_cube_elem_type) :: mesh
    !> The equation you solve.
    type(atl_equations_type) :: equation
    !> THe state if the equation
    type(atl_statedata_type) :: statedata
    !> Parameters of the modal dg scheme
    type(atl_modg_scheme_type), intent(in) :: scheme
    ! --------------------------------------------------------------------------!
    ! Loop var for all the fluid elements
    integer :: iElem
    ! The modal coefficients of the current element in the loop.
    real(kind=rk), allocatable :: modalCoeffs(:,:)
    ! Loop vars for ansatz and test functions
    integer :: iAnsX, iAnsY, iAnsZ, iTestX, iTestY, iTestZ
    ! Position of ansatz and test functions in the linearized list of coefficients
    integer :: ansPos, testPos
    ! Determinant of the jacobian of the mapping of reference element to physical
    ! element.
    real(kind=rk) :: jacobiDet
    ! Scalar products for the volume integral
    real(kind=rk) :: scalProdX, scalProdY, scalProdZ
    ! Modal represenation of the physical flux in x, y and z direction
    real(kind=rk), allocatable :: physFlux(:,:)
    ! Buffer for the flux
    real(kind=rk) :: flux
    ! Material parameters of the Nernst-Planck equations
    real(kind=rk) :: diffusivity
    ! Lower loop bounds for the projections onto the test functions
    integer :: iAnsXMin, iAnsYMin, iAnsZMin
    ! --------------------------------------------------------------------------!

    allocate( modalCoeffs((scheme%maxPolyDegree+1)**3,4) )
    allocate( physFlux((scheme%maxPolyDegree+1)**3,1) )

    ! get the Nernst-Planck material parameters
    diffusivity = equation%nerplanck%D

    ! We have cubiic elements, so the determinant of the jacobian of the mapping
    ! from reference element to physical element is the same everywhere.
    ! Please notice, that the mapping of the element itself is usually
    ! (mesh%length/2.0)**3, but the derivative in the volume integral
    ! gives an additional prefactor of 2.0/mesh%length and therefore
    ! the following is the correct scaling factor for the volume integrals.
    jacobiDet = (mesh%length/2.0_rk)**2

    ! Loop over all the elements and calculate the projection
    elemLoop: do iElem = 1, mesh%descriptor%elem%nElems(eT_fluid)

      ! get the modal coefficients of the current cell (for all variables
      ! of the equation, therefore we use ":" for the third index).
      modalCoeffs = statedata%state(iElem,:,:)

      ! calculate the physical flux in modal representation, so we loop over all
      ! of them
      do iAnsX = 1, scheme%maxPolyDegree+1
        do iAnsY = 1, scheme%maxPolyDegree+1
          do iAnsZ = 1, scheme%maxPolyDegree+1
            ! position of the modal coefficient
  anspos = iansx                                      &
    &      + ( ( iansy-1)                             &
    &      + (iansz-1)*(scheme%maxpolydegree+1))*(scheme%maxpolydegree+1)
            call atl_nerplanck_physflux_preprocess(  &
              & state       = modalCoeffs(ansPos,:), &
              & diffusivity = diffusivity,           &
              & flux        = flux                   )
             physFlux(ansPos,:) = flux
          end do
        end do
      end do

      ! now, project onto all test functions
      do iTestX = 1, scheme%maxPolyDegree+1
        do iTestY = 1, scheme%maxPolyDegree+1
          do iTestZ = 1, scheme%maxPolyDegree+1
            ! the position of the current test functions
  testpos = itestx                                      &
    &      + ( ( itesty-1)                             &
    &      + (itestz-1)*(scheme%maxpolydegree+1))*(scheme%maxpolydegree+1)

            ! for x direction (x test function differentiated)
            ! get the relevant indices for the ansatz function
            if(iTestX > 1) then
              iAnsX = iTestX - 1
              scalProdX = ply_scalProdDualLegDiff(iAnsX, iTestX)
              iAnsYMin = iTestY-2
              if(iAnsYMin < 1) then
                iAnsYMin = iTestY
              end if
              do iAnsY = iAnsYMin, iTestY, 2
                scalProdY = ply_scalProdDualLeg(iAnsY, iTestY)
                iAnsZMin = iTestZ-2
                if(iAnsZMin < 1) then
                  iAnsZMin = iTestZ
                end if
                do iAnsZ = iAnsZMin, iTestZ, 2
                  scalProdZ = ply_scalProdDualLeg(iAnsZ, iTestZ)
                  ! now, we calculate the projection
  anspos = iansx                                      &
    &      + ( ( iansy-1)                             &
    &      + (iansz-1)*(scheme%maxpolydegree+1))*(scheme%maxpolydegree+1)
!NA!                  kernelData%state_der(iElem, testPos,1,1) &
!NA!                                  & = kernelData%state_der(iElem, testPos, 1, 1) &
!NA!                                  & + physFlux(ansPos,1) * jacobiDet * scalProdX &
!NA!                                  & * scalProdY * scalProdZ
                end do
              end do
            end if


            !  for y direction (y test function differentiated)
            if(iTestY > 1) then
              iAnsY = iTestY - 1
              scalProdY = ply_scalProdDualLegDiff(iAnsY, iTestY)
              iAnsXMin = iTestX-2
              if(iAnsXMin < 1) then
                iAnsXMin = iTestX
              end if
              do iAnsX = iAnsXMin, iTestX, 2
                scalProdX = ply_scalProdDualLeg(iAnsX, iTestX)
                iAnsZMin = iTestZ-2
                if(iAnsZMin < 1) then
                  iAnsZMin = iTestZ
                end if
                do iAnsZ = iAnsZMin, iTestZ, 2
                  scalProdZ = ply_scalProdDualLeg(iAnsZ, iTestZ)
                  ! now, we calculate the projection
  anspos = iansx                                      &
    &      + ( ( iansy-1)                             &
    &      + (iansz-1)*(scheme%maxpolydegree+1))*(scheme%maxpolydegree+1)
!NA!                  kernelData%state_der(iElem, testPos,1,1) &
!NA!                                  & = kernelData%state_der(iElem, testPos, 1, 1) &
!NA!                                  & + physFlux(ansPos,1) * jacobiDet * scalProdX &
!NA!                                  & * scalProdY * scalProdZ
                end do
              end do
            end if

            ! for z direction (z test function differentiated)
            if(iTestZ > 1) then
              iAnsZ = iTestZ - 1
              scalProdZ = ply_scalProdDualLegDiff(iAnsZ, iTestZ)
              iAnsYMin = iTestY-2
              if(iAnsYMin < 1) then
                iAnsYMin = iTestY
              end if
              do iAnsY = iAnsYMin, iTestY, 2
                scalProdY = ply_scalProdDualLeg(iAnsY, iTestY)
                iAnsXMin = iTestX-2
                if(iAnsXMin < 1) then
                  iAnsXMin = iTestX
                end if
                do iAnsX = iAnsXMin, iTestX, 2
                  scalProdX = ply_scalProdDualLeg(iAnsX, iTestX)
                  ! now, we calculate the projection
  anspos = iansx                                      &
    &      + ( ( iansy-1)                             &
    &      + (iansz-1)*(scheme%maxpolydegree+1))*(scheme%maxpolydegree+1)
!NA!                  kernelData%state_der(iElem, testPos,1,1) &
!NA!                                  & = kernelData%state_der(iElem, testPos, 1, 1) &
!NA!                                  & + physFlux(ansPos,1) * jacobiDet * scalProdX &
!NA!                                  & * scalProdY * scalProdZ
                end do
              end do
            end if

          end do
        end do
      end do

    end do elemLoop

  end subroutine modg_nerplanck_physFlux_pre


  !> Applies the inverse of the mass matrix for a 3D scheme.
  subroutine modg_invMassMatrix_pre( mesh, scheme )
    ! --------------------------------------------------------------------------
    !> The mesh you are working with.
    type(atl_cube_elem_type) :: mesh
    !> Parameters of the modal dg scheme
    type(atl_modg_scheme_type), intent(in) :: scheme
    ! --------------------------------------------------------------------------!
    ! Loop indices for ansatz functions
    integer :: iAnsX, iAnsY, iAnsZ
    ! Positions for the given ansatz functions
    integer :: ansLow, ansPrevLow, ansUp, ansPrevUp
    ! Loop var for the elements
    integer :: iElem
    ! Inverse of the determintant of the jacobian of the mapping from reference
    ! to physical element.
    real(kind=rk) :: inv_jacobiDet
    ! --------------------------------------------------------------------------!

    inv_jacobiDet = (2.0_rk/mesh%length)**3.0_rk

    ! Loop over all the elements and apply the inverse mass matrix
    elemLoop: do iElem = 1, mesh%descriptor%elem%nElems(eT_fluid)

      ! apply the 1D inverse of the mass matrix
      do iAnsX = 3, scheme%maxPolyDegree+1
        do iAnsY = 1, scheme%maxPolyDegree+1
          do iAnsZ = 1, scheme%maxPolyDegree+1
  anslow = iansx                                      &
    &      + ( ( iansy-1)                             &
    &      + (iansz-1)*(scheme%maxpolydegree+1))*(scheme%maxpolydegree+1)
  ansprevlow = iansx-2                                      &
    &      + ( ( iansy-1)                             &
    &      + (iansz-1)*(scheme%maxpolydegree+1))*(scheme%maxpolydegree+1)
!NA!            kerneldata%state_der(iElem, ansLow,1,:) = kerneldata%state_der(iElem, ansLow,1,:) &
!NA!                                                    & + kerneldata%state_der(iElem, ansPrevLow,1,:)
          end do
        end do
      end do
      do iAnsX = 1, scheme%maxPolyDegree+1
        do iAnsY = 1, scheme%maxPolyDegree+1
          do iAnsZ = 1, scheme%maxPolyDegree+1
  anslow = iansx                                      &
    &      + ( ( iansy-1)                             &
    &      + (iansz-1)*(scheme%maxpolydegree+1))*(scheme%maxpolydegree+1)
!NA!            kerneldata%state_der(iElem, ansLow,1,:) = kerneldata%state_der(iElem, ansLow,1,:) &
!NA!                                                    & * (2.0_rk*iAnsX-1.0_rk)/2.0_rk
          end do
        end do
      end do

      ! apply the 2D inverse of the mass matrix
      do iAnsY = 3, scheme%maxPolyDegree+1
        do iAnsZ = 1, scheme%maxPolyDegree+1
  anslow = 1                                      &
    &      + ( ( iansy-1)                             &
    &      + (iansz-1)*(scheme%maxpolydegree+1))*(scheme%maxpolydegree+1)
  ansup = scheme%maxpolydegree+1                                      &
    &      + ( ( iansy-1)                             &
    &      + (iansz-1)*(scheme%maxpolydegree+1))*(scheme%maxpolydegree+1)
  ansprevlow = 1                                      &
    &      + ( ( iansy-2-1)                             &
    &      + (iansz-1)*(scheme%maxpolydegree+1))*(scheme%maxpolydegree+1)
  ansprevup = scheme%maxpolydegree+1                                      &
    &      + ( ( iansy-2-1)                             &
    &      + (iansz-1)*(scheme%maxpolydegree+1))*(scheme%maxpolydegree+1)
!NA!          kerneldata%state_der(iElem, ansLow:ansUp,1,:) &
!NA!                                   & = kerneldata%state_der(iElem, ansLow:ansUp,1,:) &
!NA!                                   & + kerneldata%state_der(iElem, ansPrevLow:ansPrevUp,1,:)
        end do
      end do
      do iAnsY = 1, scheme%maxPolyDegree+1
        do iAnsZ = 1, scheme%maxPolyDegree+1
  anslow = 1                                      &
    &      + ( ( iansy-1)                             &
    &      + (iansz-1)*(scheme%maxpolydegree+1))*(scheme%maxpolydegree+1)
  ansup = scheme%maxpolydegree+1                                      &
    &      + ( ( iansy-1)                             &
    &      + (iansz-1)*(scheme%maxpolydegree+1))*(scheme%maxpolydegree+1)
!NA          kerneldata%state_der(iElem, ansLow:ansUp,1,:) &
!NA                                   & = kerneldata%state_der(iElem, ansLow:ansUp,1,:) &
!NA                                   & * (2.0_rk*iAnsY-1.0_rk)/2.0_rk
        end do
      end do

      ! apply the 3D inverse of the mass matrix
      do iAnsZ = 3, scheme%maxPolyDegree+1
  anslow = 1                                      &
    &      + ( ( 1-1)                             &
    &      + (iansz-1)*(scheme%maxpolydegree+1))*(scheme%maxpolydegree+1)
  ansup = scheme%maxpolydegree+1                                      &
    &      + ( ( scheme%maxpolydegree+1-1)                             &
    &      + (iansz-1)*(scheme%maxpolydegree+1))*(scheme%maxpolydegree+1)
  ansprevlow = 1                                      &
    &      + ( ( 1-1)                             &
    &      + (iansz-2-1)*(scheme%maxpolydegree+1))*(scheme%maxpolydegree+1)
  ansprevup = scheme%maxpolydegree+1                                      &
    &      + ( ( scheme%maxpolydegree+1-1)                             &
    &      + (iansz-2-1)*(scheme%maxpolydegree+1))*(scheme%maxpolydegree+1)
!NA!        kerneldata%state_der(iElem, ansLow:ansUp,1,:) &
!NA!                                   & = kerneldata%state_der(iElem, ansLow:ansUp,1,:) &
!NA!                                   & + kerneldata%state_der(iElem, ansPrevLow:ansPrevUp,1,:)
      end do
      do iAnsZ = 1, scheme%maxPolyDegree+1
  anslow = 1                                      &
    &      + ( ( 1-1)                             &
    &      + (iansz-1)*(scheme%maxpolydegree+1))*(scheme%maxpolydegree+1)
  ansup = scheme%maxpolydegree+1                                      &
    &      + ( ( scheme%maxpolydegree+1-1)                             &
    &      + (iansz-1)*(scheme%maxpolydegree+1))*(scheme%maxpolydegree+1)
!NA!        kerneldata%state_der(iElem, ansLow:ansUp,1,:) &
!NA!                                   & = kerneldata%state_der(iElem, ansLow:ansUp,1,:) &
!NA!                                   & * (2.0_rk*iAnsZ-1.0_rk)/2.0_rk
      end do
    end do elemLoop

  end subroutine modg_invMassMatrix_pre
end module atl_modg_nerplanck_kernel_module