! Copyright (c) 2019-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2020 Peter Vitt <peter.vitt2@uni-siegen.de>
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice,
! this list of conditions and the following disclaimer.
!
! 2. 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.
!
! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “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 UNIVERSITY OF SIEGEN 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.
! ****************************************************************************** !
!> author: Kannan Masilamani
!! author: Jiaxing Qi
!! This module provides the MUSUBI specific functions for calculating
!! macroscopic quantities from the state variables.
!!
!! The depending common interface between MUSUBI and ATELES is defined in the
!! tem_derived_module. The functionality for accessing a variable from the state
!! and evaluating a lua function are also provided in the tem_derived module.
!!
!! Do not use get_Element or get_Point routines to update the state !
!!
! Copyright (c) 2011-2013 Manuel Hasert <m.hasert@grs-sim.de>
! Copyright (c) 2011 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2011 Konstantin Kleinheinz <k.kleinheinz@grs-sim.de>
! Copyright (c) 2011-2012 Simon Zimny <s.zimny@grs-sim.de>
! Copyright (c) 2012, 2014-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2012 Kartik Jain <kartik.jain@uni-siegen.de>
! Copyright (c) 2013-2015, 2019 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice,
! this list of conditions and the following disclaimer.
!
! 2. 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.
!
! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “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 UNIVERSITY OF SIEGEN 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.
! Copyright (c) 2013 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2013-2014 Nikhil Anand <nikhil.anand@uni-siegen.de>
! Copyright (c) 2014, 2016 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2015, 2018, 2020 Peter Vitt <peter.vitt2@uni-siegen.de>
! Copyright (c) 2016 Verena Krupp <verena.krupp@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice, this
! list of conditions and the following disclaimer.
!
! 2. 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.
!
! 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 HOLDER 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.
!--------------------------------------------
!    A O S - Array of structures layout new
!-------------------------------------------
! Access to get_point value output
! Access to get_element value output
module mus_derQuanPS_module
  use iso_c_binding, only: c_loc, c_ptr, c_f_pointer

  ! include treelm modules
  use tem_param_module,         only: div1_2, div1_3, div1_54, div1_9, div3_4, &
    &                                 sqrt3, cs2inv, cs2, t2cs2inv, t2cs4inv,  &
    &                                 cs4inv, q000
  use env_module,               only: rk, long_k, labelLen
  use tem_float_module,         only: operator(.feq.), operator(.fge.), &
    &                                 operator(.fle.)
  use tem_variable_module,      only: tem_variable_type
  use tem_stencil_module,       only: tem_stencilHeader_type
  use tem_topology_module,      only: tem_levelOf
  use tem_time_module,          only: tem_time_type
  use treelmesh_module,         only: treelmesh_type
  use tem_logging_module,       only: logUnit
  use tem_varSys_module,        only: tem_varSys_type, tem_varSys_op_type,     &
    &                                 tem_varSys_append_derVar,                &
    &                                 tem_varSys_append_auxFieldVar,           &
    &                                 tem_varSys_proc_point,                   &
    &                                 tem_varSys_proc_element,                 &
    &                                 tem_varSys_proc_setParams,               &
    &                                 tem_varSys_proc_getParams,               &
    &                                 tem_varSys_proc_setupIndices,            &
    &                                 tem_varSys_proc_getValOfIndex,           &
    &                                 tem_varSys_getPoint_dummy,               &
    &                                 tem_varSys_getElement_dummy,             &
    &                                 tem_varSys_setupIndices_dummy,           &
    &                                 tem_varSys_getValOfIndex_dummy,          &
    &                                 tem_varSys_setParams_dummy,              &
    &                                 tem_varSys_getParams_dummy
  use tem_aux_module,           only: tem_abort
  use tem_operation_var_module, only: tem_evalMag_forElement,     &
    &                                 tem_evalMag_forPoint,       &
    &                                 tem_evalMag_fromIndex,      &
    &                                 tem_opVar_setupIndices,     &
    &                                 tem_get_new_varSys_data_ptr
  use tem_property_module,      only: prp_hasBnd, prp_hasQval
  use tem_tools_module,         only: tem_PositionInSorted
  use tem_debug_module,         only: dbgUnit
  use tem_grow_array_module,    only: grw_labelarray_type, append

  ! include musubi modules
  use mus_source_type_module,        only: mus_source_op_type
  use mus_pdf_module,                only: pdf_data_type
  use mus_scheme_header_module,      only: mus_scheme_header_type
  use mus_scheme_layout_module,      only: mus_scheme_layout_type
  use mus_scheme_type_module,        only: mus_scheme_type
  use mus_derivedQuantities_module2, only: secondMom, getShearRate
  use mus_varSys_module,             only: mus_varSys_data_type,             &
    &                                      mus_varSys_solverData_type,       &
    &                                      mus_get_new_solver_ptr,           &
    &                                      mus_deriveVar_forPoint,           &
    &                                      mus_generic_varFromPDF_fromIndex, &
    &                                      mus_generic_fromPDF_forElement,   &
    &                                      mus_derive_fromPDF
  use mus_operation_var_module,      only: mus_opVar_setupIndices
  use mus_derVarPos_module,          only: mus_derVarPos_type
  use mus_derQuan_module,            only: derivePressure,        &
    &                                      derivePressure_fromIndex
  use mus_physics_module,            only: mus_convertFac_type

  implicit none

  private

  public :: mus_append_derVar_lbmPS
  public :: deriveEquilPS_FromMacro
  public :: deriveEquilPS_fromAux
  public :: deriveAuxPS_fromState

  ! source variable
  public :: derive_injectionPS
  public :: derive_equalInjectionPS

  ! source update
  public :: applySrc_injectionPS
  public :: applySrc_equalInjectionPS

contains

  ! **************************************************************************** !
  !> subroutine to add derive variables for weakly compressible LBM
  !! (schemekind = 'passive_scalar') to the varsys.
  !! for passive scalar contains only one derive variable:
  !! density
  subroutine mus_append_derVar_lbmPS( varSys, solverData, fldLabel, derVarName )
    ! ---------------------------------------------------------------------------
    !> global variable system
    type(tem_varSys_type), intent(inout)  :: varSys

    !> Contains pointer to solver data types
    type(mus_varSys_solverData_type), target, intent(in) :: solverData

    !> array of field label prefix. Size=nFields
    character(len=*), intent(in)              :: fldLabel

    !> array of derive physical variables
    type(grw_labelarray_type), intent(inout) :: derVarName
    ! ---------------------------------------------------------------------------
    ! number of derive variables
    integer :: addedPos
    logical :: wasAdded
    character(len=labelLen), allocatable ::  input_varname(:)
    character(len=labelLen)  ::  varName
    procedure(tem_varSys_proc_point), pointer :: get_point => NULL()
    procedure(tem_varSys_proc_element), pointer :: get_element => NULL()
    procedure(tem_varSys_proc_setParams), pointer :: set_params => null()
    procedure(tem_varSys_proc_getParams), pointer :: get_params => null()
    procedure(tem_varSys_proc_setupIndices), pointer :: &
      &                                      setup_indices => null()
    procedure(tem_varSys_proc_getValOfIndex), pointer :: &
      &                                       get_valOfIndex => null()
    character(len=labelLen) :: derVarName_loc
    ! ---------------------------------------------------------------------------
    nullify(get_point, get_element, set_params, get_params, setup_indices, &
      &     get_valOfIndex)

    derVarName_loc = 'pressure'

    call append(derVarName, derVarName_loc)

    ! set pointers for pressure variable
    get_element => derivePressure
    get_point => mus_deriveVar_forPoint
    setup_indices => mus_opVar_setupIndices
    get_valOfIndex => derivePressure_fromIndex
    set_params => tem_varSys_setParams_dummy
    get_params => tem_varSys_getParams_dummy

    ! update variable names with field label
    varname = trim(fldLabel)//trim(adjustl(derVarName_loc))
    allocate(input_varname(1))
    input_varname(1) = trim(fldLabel)//'density'

    ! append variable to varSys
    call tem_varSys_append_derVar(                            &
      &  me             = varSys,                             &
      &  varName        = trim(varname),                      &
      &  nComponents    = 1,                                  &
      &  input_varname  = input_varname,                      &
      &  method_data    = mus_get_new_solver_ptr(solverData), &
      &  get_point      = get_point,                          &
      &  get_element    = get_element,                        &
      &  set_params     = set_params,                         &
      &  get_params     = get_params,                         &
      &  setup_indices  = setup_indices,                      &
      &  get_valOfIndex = get_valOfIndex,                     &
      &  pos            = addedPos,                           &
      &  wasAdded       = wasAdded                            )

    if (wasAdded) then
      write(logUnit(10),*) 'Appended variable: '//trim(varname)
    else if (addedpos < 1) then
      write(logUnit(1),*) 'Error: variable '//trim(varname)// &
        &                 ' is not added to variable system'
      call tem_abort()
    end if
    deallocate(input_varname)

  end subroutine mus_append_derVar_lbmPS
  ! **************************************************************************** !


! ****************************************************************************** !
!       Subroutines with common interface for the function pointers            !
! ****************************************************************************** !

! ****************************************************************************** !
  !> This routine computes equilbrium from density and velocity
  !! This must comply with mus_variable_module%derive_FromMacro
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[derive_FromMacro]] in derived/[[mus_derVarPos_module]].f90 in order to be
  !! callable via [[mus_derVarPos_type:equilFromMacro]] function pointer.
  subroutine deriveEquilPS_FromMacro( density, velocity, iField, nElems, &
    &                                 varSys, layout, res                )
    ! -------------------------------------------------------------------- !
    !> Array of density.
    !! Single species: dens_1, dens_2 .. dens_n
    !! multi-species: dens_1_sp1, dens_1_sp2, dens_2_sp1, dens_2_sp2 ...
    !!                dens_n_sp1, dens_n_sp2
    real(kind=rk), intent(in) :: density(:)

    !> Array of velocity.
    !! Size: dimension 1: n*nFields. dimension 2: 3 (nComp)
    !! 1st dimension arrangement for multi-species is same as density
    real(kind=rk), intent(in) :: velocity(:, :)

    !> Current field
    integer, intent(in) :: iField

    !> number of elements
    integer, intent(in) :: nElems

    !> variable system which is required to access fieldProp
    !! information via variable method data c_ptr
    type(tem_varSys_type), intent(in) :: varSys

    !> scheme layout contains stencil definition and lattice weights
    type(mus_scheme_layout_type), intent(in) :: layout

    !> Output of this routine
    !! Dimension: n*nComponents of res
    real(kind=rk), intent(out) :: res(:)
    ! -------------------------------------------------------------------- !
    real(kind=rk) :: fEq(layout%fStencil%QQ), vel(3), ucx
    integer :: QQ, iElem, iDir
    ! ---------------------------------------------------------------------------
    QQ = layout%fStencil%QQ
    do iElem = 1, nElems
      vel = velocity(:,iElem)

      do iDir = 1, QQ
        ucx = dot_product( layout%fStencil%cxDirRK(:, iDir), vel )

        ! calculate equilibrium
        fEq( iDir ) = layout%weight( iDir ) * density(iElem) &
          & * ( 1._rk +  cs2inv * ucx )
      enddo

      res( (iElem-1)*QQ+1: iElem*QQ ) = fEq
    end do
  end subroutine deriveEquilPS_FromMacro
! ****************************************************************************** !


! **************************************************************************** !
  !> This routine computes auxField from state array
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[derive_auxFromState]] in derived/[[mus_derVarPos_module]].f90 in order to
  !! be callable via [[mus_derVarPos_type:auxFieldFromState]] function pointer.
  subroutine deriveAuxPS_fromState( derVarPos, state, neigh, iField, nElems, &
    &                               nSize, iLevel, stencil, varSys, auxField )
    ! -------------------------------------------------------------------- !
    !> Position of derive variable in variable system
    class(mus_derVarPos_type), intent(in) :: derVarPos
    !> Array of state
    !! n * layout%stencil(1)%QQ * nFields
    real(kind=rk), intent(in) :: state(:)

    !> connectivity vector
    integer, intent(in) :: neigh(:)

    !> Current field
    integer, intent(in) :: iField

    !> number of elements
    integer, intent(in) :: nElems

    !> number of elements in state array
    integer, intent(in) :: nSize

    !> current level
    integer, intent(in) :: iLevel

    !> stencil header contains discrete velocity vectors
    type(tem_stencilHeader_type), intent(in) :: stencil

    !> variable system which is required to access fieldProp
    !! information via variable method data c_ptr
    type(tem_varSys_type), intent(in) :: varSys

    !> Output of this routine
    !! Size: nElems*nAuxScalars
    real(kind=rk), intent(inout) :: auxField(:)
    ! -------------------------------------------------------------------- !
    integer :: iElem, iDir
    real(kind=rk) :: pdf( stencil%QQ )
    ! ------------------------------------------------------------------------ !
    !NEC$ ivdep
    do iElem = 1, nElems
      !NEC$ shortloop
      do iDir = 1, stencil%QQ
        pdf(iDir) = state(                                                     &
          & ( ielem-1)* varsys%nscalars+idir+( 1-1)* stencil%qq)
      end do

      ! element offset is not required because passive scalar has only
      ! one aux scalar
      ! density
      auxField(iElem) = sum(pdf)
    end do

  end subroutine deriveAuxPS_fromState
! **************************************************************************** !

  ! ************************************************************************** !
  !> This routine computes equilbrium from auxField
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[derive_equilFromAux]] in derived/[[mus_derVarPos_module]].f90 in order to
  !! be callable via [[mus_derVarPos_type:equilFromAux]] function pointer.
  subroutine deriveEquilPS_fromAux( derVarPos, auxField, iField, nElems, &
    &                               varSys, layout, fEq                  )
    ! -------------------------------------------------------------------- !
    !> Position of derive variable in variable system
    class(mus_derVarPos_type), intent(in) :: derVarPos
    !> Array of auxField.
    !! Single species: dens_1, vel_1, dens_2, vel_2, .. dens_n, vel_n
    !! multi-species: dens_1_sp1, vel_1_spc1, dens_1_sp2, vel_1_spc2,
    !!                dens_2_sp1, vel_2_spc2, dens_2_sp2, vel_2_spc2 ...
    !!                dens_n_sp1, vel_n_sp1, dens_n_sp2, vel_n_spc2
    !! Access: (iElem-1)*nAuxScalars + auxField_varPos
    real(kind=rk), intent(in) :: auxField(:)

    !> Current field
    integer, intent(in) :: iField

    !> number of elements
    integer, intent(in) :: nElems

    !> variable system which is required to access fieldProp
    !! information via variable method data c_ptr
    type(tem_varSys_type), intent(in) :: varSys

    !> scheme layout contains stencil definition and lattice weights
    type(mus_scheme_layout_type), intent(in) :: layout

    !> Output of this routine
    !! Dimension: n*QQ of res
    real(kind=rk), intent(out) :: fEq(:)
    ! -------------------------------------------------------------------- !
    write(logUnit(1),*) 'ERROR: Equilibrium calculation requires transport '
    write(logUnit(1),*) 'velocity and it is not provided to this routine'
    write(logUnit(1),*) 'Solution: use deriveEquilPS_fromMacro'
    call tem_abort()
  end subroutine deriveEquilPS_fromAux
  ! ************************************************************************** !

! ****************************************************************************** !
   !> Derive injection variable defined as a source term for lbm passive scalar.
   !! It evaluates spacetime function defined in lua file for injection variable
   !! and convert it to state value which is to be added to the state
  recursive subroutine derive_injectionPS(fun, varsys, elempos, time, tree, &
    &                                     nElems, nDofs, res                )
    ! -------------------------------------------------------------------- !
    !> Description of the method to obtain the variables, here some preset
    !! values might be stored, like the space time function to use or the
    !! required variables.
    class(tem_varSys_op_type), intent(in) :: fun

    !> The variable system to obtain the variable from.
    type(tem_varSys_type), intent(in) :: varSys

    !> Position of the TreeID of the element to get the variable for in the
    !! global treeID list.
    integer, intent(in) :: elempos(:)

    !> Point in time at which to evaluate the variable.
    type(tem_time_type), intent(in)  :: time

    !> global treelm mesh info
    type(treelmesh_type), intent(in) :: tree

    !> Number of values to obtain for this variable (vectorized access).
    integer, intent(in) :: nElems

    !> Number of degrees of freedom within an element.
    integer, intent(in) :: nDofs

    !> Resulting values for the requested variable.
    !!
    !! Linearized array dimension:
    !! (n requested entries) x (nComponents of this variable)
    !! x (nDegrees of freedom)
    !! Access: (iElem-1)*fun%nComponents*nDofs +
    !!         (iDof-1)*fun%nComponents + iComp
    real(kind=rk), intent(out) :: res(:)
    ! -------------------------------------------------------------------- !
    type(mus_varSys_data_type), pointer :: fPtr
    type(mus_scheme_type), pointer :: scheme
    real(kind=rk) :: luaRes(nElems)
    real(kind=rk) :: velocity(nElems*3)
    real(kind=rk) :: local_vel(3)
    real(kind=rk), allocatable :: linkwise_luaRes(:)
    integer :: trans_velPos, data_varPos
    integer :: iElem, iDir
    integer :: nCompPDF, nInputStates
    ! -------------------------------------------------------------------- !
    ! convert c pointer to solver type fortran pointer
    call c_f_pointer( fun%method_data, fPtr )
    scheme => fPtr%solverData%scheme

    ! number of pdf states this source depends on
    ! last input is spacetime function so it is neglected
    nInputStates = fun%nInputs - 1
    ! constant parameter
    nCompPDF = scheme%layout%fStencil%QQ
    allocate(linkwise_luaRes(nCompPDF))

    ! last input var is spacetime function
    data_varPos = fun%input_varPos(fun%nInputs)
    ! Get injection which might be defined to vary in space and time
    call varSys%method%val(data_varPos)%get_element( varSys  = varSys,  &
      &                                              elemPos = elemPos, &
      &                                              time    = time,    &
      &                                              tree    = tree,    &
      &                                              nElems  = nElems,  &
      &                                              ndofs   = 1,       &
      &                                              res     = luaRes   )

    ! Get transport velocity field from stFun
    trans_velPos = scheme%transVar%method(1)%data_varPos
    call varSys%method%val(trans_velPos)%get_element( varSys  = varSys,  &
      &                                               elemPos = elemPos, &
      &                                               time    = time,    &
      &                                               tree    = tree,    &
      &                                               nElems  = nElems,  &
      &                                               nDofs   = 1,       &
      &                                               res     = velocity )

    do iElem = 1, nElems

      ! velocity from dependent scheme
      local_vel = velocity((iElem-1)*3+1 : iElem*3)

      do iDir = 1, nCompPDf
        linkwise_luaRes(iDir) =                                              &
          & luaRes(iElem) /  6._rk * (1._rk + 3._rk * (                      &
          & dot_product(                                                     &
          & scheme%layout%fStencil%cxDirRK(:,iDir), local_vel(:)) ) )
      end do

      ! Update souce depends on nInputStates
      ! if nInputStates = 1, it is field source else it is global source
      res( (iElem-1)*fun%nComponents + 1 : iElem*fun%nComponents) &
        & = linkwise_luaRes

    end do !iElem

  end subroutine derive_injectionPS
! ****************************************************************************** !


! ****************************************************************************** !
   !> Derive injection variable defined as a source term for lbm passive scalar.
   !! It evaluates spacetime function defined in lua file for injection variable
   !! and convert it to state value which is to be added to the state
  recursive subroutine derive_equalInjectionPS(fun, varsys, elempos, time, &
    &                                          tree, nElems, nDofs, res    )
    ! -------------------------------------------------------------------- !
    !> Description of the method to obtain the variables, here some preset
    !! values might be stored, like the space time function to use or the
    !! required variables.
    class(tem_varSys_op_type), intent(in) :: fun

    !> The variable system to obtain the variable from.
    type(tem_varSys_type), intent(in) :: varSys

    !> Position of the TreeID of the element to get the variable for in the
    !! global treeID list.
    integer, intent(in) :: elempos(:)

    !> Point in time at which to evaluate the variable.
    type(tem_time_type), intent(in)  :: time

    !> global treelm mesh info
    type(treelmesh_type), intent(in) :: tree

    !> Number of values to obtain for this variable (vectorized access).
    integer, intent(in) :: nElems

    !> Number of degrees of freedom within an element.
    integer, intent(in) :: nDofs

    !> Resulting values for the requested variable.
    !!
    !! Linearized array dimension:
    !! (n requested entries) x (nComponents of this variable)
    !! x (nDegrees of freedom)
    !! Access: (iElem-1)*fun%nComponents*nDofs +
    !!         (iDof-1)*fun%nComponents + iComp
    real(kind=rk), intent(out) :: res(:)
    ! -------------------------------------------------------------------- !
    type(mus_varSys_data_type), pointer :: fPtr
    real(kind=rk) :: luaRes(nElems)
    real(kind=rk), allocatable :: linkwise_luaRes(:)
    integer :: iElem, iDir, stFun_pos
    integer :: iField, nCompPDF, nInputStates
    real(kind=rk) :: srcTerm(fun%nComponents)
    ! -------------------------------------------------------------------- !
    ! convert c pointer to solver type fortran pointer
    call c_f_pointer( fun%method_data, fPtr )

    ! number of pdf states this source depends on
    ! last input is spacetime function so it is neglected
    nInputStates = fun%nInputs - 1
    ! constant parameter
    nCompPDF = fPtr%solverData%scheme%layout%fStencil%QQ
    allocate(linkwise_luaRes(nCompPDF))

    ! last input var is spacetime function
    stFun_pos = fun%input_varPos(fun%nInputs)

    ! Get injection which might be defined to vary in space and time
    call varSys%method%val(stfun_pos)%get_element( varSys  = varSys,  &
      &                                            elemPos = elemPos, &
      &                                            time    = time,    &
      &                                            tree    = tree,    &
      &                                            ndofs   = 1,       &
      &                                            nElems  = nElems,  &
      &                                            res     = luaRes   )

    ! loop for all elements
    do iElem = 1, nElems

      ! compute the linkwise equal distribution
      do iDir = 1, nCompPDf
        linkwise_luaRes(iDir) = luaRes(iElem) /  real(nCompPDF, kind=rk)
      end do

      ! Update souce depends on nInputStates
      ! if nInputStates = 1, it is field source else it is global source
      do iField = 1, nInputStates
        ! loop over all directions
        do iDir = 1, nCompPDF
          srcTerm((iField-1)*nCompPDF + iDir) = linkwise_luaRes(iDir)
        end do
      end do
    end do !iElem

  end subroutine derive_equalInjectionPS
! ****************************************************************************** !

! ****************************************************************************** !
  !> Update state with source variable "injection"
  !! Similar to derive routine but it updates the state whereas derive
  !! is used for tracking
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[proc_apply_source]] in derived/[[mus_source_type_module]].f90 in order to
  !! be callable via [[mus_source_op_type:applySrc]] function pointer.
  subroutine applySrc_injectionPS( fun, inState, outState, neigh, auxField,    &
    &                              nPdfSize, iLevel, varSys, time, phyConvFac, &
    &                              derVarPos                                   )
    ! -------------------------------------------------------------------- !
    !> Description of method to apply source terms
    class(mus_source_op_type), intent(in) :: fun

    !> input  pdf vector
    real(kind=rk), intent(in) :: inState(:)

    !> output pdf vector
    real(kind=rk), intent(inout) :: outState(:)

    !> connectivity Array corresponding to state vector
    integer,intent(in) :: neigh(:)

    !> auxField array
    real(kind=rk), intent(in) :: auxField(:)

    !> number of elements in state Array
    integer, intent(in) :: nPdfSize

    !> current level
    integer, intent(in) :: iLevel

    !> variable system
    type(tem_varSys_type), intent(in) :: varSys

    !> Point in time at which to evaluate the variable.
    type(tem_time_type), intent(in)  :: time

    !> Physics conversion factor for current level
    type(mus_convertFac_type), intent(in) :: phyConvFac

    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------- !
    type(mus_varSys_data_type), pointer :: fPtr
    type(mus_scheme_type), pointer :: scheme
    real(kind=rk) :: luaRes(fun%elemLvl(iLevel)%nElems)
    real(kind=rk) :: transVel(fun%elemLvl(iLevel)%nElems*3)
    real(kind=rk) :: local_vel(3)
    real(kind=rk), allocatable :: linkwise_luaRes(:)
    integer :: iElem, nElems, iDir, posInTotal
    integer :: iField, depField, nScalars, QQ, nInputStates
    integer :: trans_varPos
    ! ---------------------------------------------------------------------------
    ! convert c pointer to solver type fortran pointer
    call c_f_pointer( varSys%method%val( fun%srcTerm_varPos )%method_data, &
      &               fPtr )
    scheme => fPtr%solverData%scheme

    ! Number of elements to apply source terms
    nElems = fun%elemLvl(iLevel)%nElems

    ! number of pdf states this source depends on
    ! last input is spacetime function so it is neglected
    nInputStates = varSys%method%val(fun%srcTerm_varPos)%nInputs - 1
    ! constant parameter
    QQ = scheme%layout%fStencil%QQ
    nScalars = varSys%nScalars
    allocate(linkwise_luaRes(QQ))

    ! Get injection which might be defined to vary in space and time
    call varSys%method%val(fun%data_varPos)%get_valOfIndex( &
      & varSys  = varSys,                                   &
      & time    = time,                                     &
      & iLevel  = iLevel,                                   &
      & idx     = fun%elemLvl(iLevel)%idx(1:nElems),        &
      & nVals   = nElems,                                   &
      & res     = luaRes                                    )

    ! Get transport velocity field from stFun
    ! Setup indices for transport variable is set for nSolve elements
    ! in total list so to get only velocity on active source space,
    ! use idx of source on transVar%index
    trans_varPos = scheme%transVar%method(1)%data_varPos
    call varSys%method%val(trans_varPos)%get_valOfIndex(   &
      & varSys  = varSys,                                  &
      & time    = time,                                    &
      & iLevel  = iLevel,                                  &
      & idx     = scheme%transVar%method(1)%pntIndex       &
      &           %indexLvl(iLevel)                        &
      &           %val(fun%elemLvl(iLevel)%idx(1:nElems)), &
      & nVals   = nElems,                                  &
      & res     = transVel                                 )

    ! convert physical to lattice
    transVel = transVel / fPtr%solverData%physics%fac(iLevel)%vel

    ! now update the state vector for all elements
    do iElem = 1, nElems

      ! position of source element in levelwise state array
      posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem)

      ! velocity from dependent scheme
      local_vel = transVel((iElem-1)*3+1 : iElem*3)

      do iDir = 1, QQ
        linkwise_luaRes(iDir) =                                              &
          & luaRes(iElem) /  6._rk * (1._rk + 3._rk * (                      &
          & dot_product(                                                     &
          & fPtr%solverData%scheme%layout%fStencil%cxDirRK(:,iDir),          &
          & local_vel(:)) ) )
      end do

      ! Update souce depends on nInputStates
      ! if nInputStates = 1, it is field source else it is global source
      do iField = 1, nInputStates
        depField = varSys%method%val(fun%srcTerm_varPos)%input_varPos(iField)
        ! loop over all directions
        do iDir = 1, QQ
          ! add the linkwise lua res for the current element
          outState( (posintotal-1)*nscalars+idir+(depfield-1)*qq )&
            & = outState(                                                      &
            & (posintotal-1)*nscalars+idir+(depfield-1)*qq )      &
            & + linkwise_luaRes(iDir)
        end do
      end do
    end do !iElem
  end subroutine applySrc_injectionPS
! ****************************************************************************** !


! ****************************************************************************** !
  !> Update state with source variable "equalInjection"
  !! Simuilar to derive routine but it updates the state whereas derive
  !! is used for tracking
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[proc_apply_source]] in derived/[[mus_source_type_module]].f90 in order to
  !! be callable via [[mus_source_op_type:applySrc]] function pointer.
  subroutine applySrc_equalInjectionPS( fun, inState, outState, neigh,      &
    &                                   auxField, nPdfSize, iLevel, varSys, &
    &                                   time, phyConvFac, derVarPos         )
    ! -------------------------------------------------------------------- !
    !> Description of method to apply source terms
    class(mus_source_op_type), intent(in) :: fun

    !> input  pdf vector
    real(kind=rk), intent(in) :: inState(:)

    !> output pdf vector
    real(kind=rk), intent(inout) :: outState(:)

    !> connectivity Array corresponding to state vector
    integer,intent(in) :: neigh(:)

    !> auxField array
    real(kind=rk), intent(in) :: auxField(:)

    !> number of elements in state Array
    integer, intent(in) :: nPdfSize

    !> current level
    integer, intent(in) :: iLevel

    !> variable system
    type(tem_varSys_type), intent(in) :: varSys

    !> Point in time at which to evaluate the variable.
    type(tem_time_type), intent(in)  :: time

    !> Physics conversion factor for current level
    type(mus_convertFac_type), intent(in) :: phyConvFac

    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------- !
    type(mus_varSys_data_type), pointer :: fPtr
    real(kind=rk) :: luaRes(fun%elemLvl(iLevel)%nElems)
    real(kind=rk), allocatable :: linkwise_luaRes(:)
    integer :: iElem, nElems, iDir, posInTotal
    integer :: iField, depField, nScalars, QQ, nInputStates
    ! ---------------------------------------------------------------------------
    ! convert c pointer to solver type fortran pointer
    call c_f_pointer( varSys%method%val( fun%srcTerm_varPos )%method_data, &
      &               fPtr )

    nElems = fun%elemLvl(iLevel)%nElems
    ! number of pdf states this source depends on
    ! last input is spacetime function so it is neglected
    nInputStates = varSys%method%val(fun%srcTerm_varPos)%nInputs - 1
    ! constant parameter
    QQ = fPtr%solverData%scheme%layout%fStencil%QQ
    nScalars = varSys%nScalars
    allocate(linkwise_luaRes(QQ))

    ! Get injection which might be defined to vary in space and time
    call varSys%method%val(fun%data_varPos)%get_valOfIndex( &
      & varSys  = varSys,                                   &
      & time    = time,                                     &
      & iLevel  = iLevel,                                   &
      & idx     = fun%elemLvl(iLevel)%idx(1:nElems),        &
      & nVals   = nElems,                                   &
      & res     = luaRes                                    )

    ! now update the state vector for all elements
    do iElem = 1, nElems

      ! position of source element in levelwise state array
      posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem)

      ! compute the linkwise equal distribution
      do iDir = 1, QQ
        linkwise_luaRes(iDir) = luaRes(iElem) /  real(QQ, kind=rk)
      end do

      ! Update souce depends on nInputStates
      ! if nInputStates = 1, it is field source else it is global source
      do iField = 1, nInputStates
        depField = varSys%method%val(fun%srcTerm_varPos)%input_varPos(iField)
        ! loop over all directions
        do iDir = 1, QQ
          ! add the linkwise lua res for the current element
          outState( (posintotal-1)*nscalars+idir+(depfield-1)*qq ) &
            & = outState(                                                  &
            & (posintotal-1)*nscalars+idir+(depfield-1)*qq )       &
            & + linkwise_luaRes(iDir)
        end do
      end do
    end do !iElem
  end subroutine applySrc_equalInjectionPS
! ****************************************************************************** !


end module mus_derQuanPS_module
