## Source Code

! Copyright (c) 2019-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2019 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.
! **************************************************************************** !
!> This module contains data types, function and routines for gradient
!! computation.
!!
!! author: Kannan Masilamani
! 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.
! include treelm modules
use env_module,          only: rk
use tem_debug_module,    only: dbgUnit
use tem_logging_module,  only: logUnit
use tem_aux_module,      only: tem_abort
use tem_param_module,    only: qN00, q100, q0N0, q010, q00N, q001
use tem_construction_module, only: tem_levelDesc_type

implicit none
private

!> Contains information required to compute gradient like position
!! of six direct neigbors in the state array and coeff for the gradient
!! operation
!> Stores position of 6 direct face neighbors in the state array for
!! each element
!! Size: nSize, stencil%nDims, 2
!! last index refers to left and right side
!! 1 is left/negative side and 2 is right/positive side
integer, allocatable :: neighPos(:,:,:)

!> coeff to calculate 1st order derivative
!! use forward difference for element with boundary (coeff = 1.0)
!! and central difference for inner fluid element (coeff=0.5)
!! size: nSize, stencil%nDims
real(kind=rk), allocatable :: FDcoeff(:,:)

contains

! ************************************************************************* !
!> This routine initialize gradData with direct neighbors in state and
!! finite difference coefficients.
subroutine mus_init_gradData(me, neigh, levelDesc, stencil, nSize, nSolve, &
&                          nScalars)
!---------------------------------------------------------------------------
!> neighbor connectivity array
integer, intent(in) :: neigh(:)
!> levelDesc to access communication buffers of state array
type(tem_levelDesc_type), intent(in) :: levelDesc
!> Number of elements in state array
integer, intent(in) :: nSize
!> Number of elements solved in compute kernel i.e. excluding halos
integer, intent(in) :: nSolve
!> number of scalars in state array
integer, intent(in) :: nScalars
!---------------------------------------------------------------------------
integer :: nDims
integer :: iElem, iDir, iFace, iMeshDir
integer :: nFaces
integer :: nghElem
!---------------------------------------------------------------------------
nDims = stencil%nDims
allocate(me%neighPos( nSolve, nDims, 2 ) )
me%neighPos = -1
allocate(me%FDcoeff( nSolve, nDims ) )
me%FDcoeff = -1.0_rk

! fill gradient data with face neighbor information for turbulence model
! Number of faces per element
nFaces = nDims*2

elemLoop: do iElem = 1, nSolve
!write(dbgUnit(1),*) 'iElem: ', iElem, ' treeID: ', levelDesc%total(iElem)

iFace = 0

neighLoop: do iDir  = 1, stencil%QQN
!write(dbgUnit(1),*) 'iDir: ', iDir

! coeff to calculate 1st order derivative
! use forward difference for element with boundary (coeff = 1.0)
! and central difference for inner fluid element (coeff=0.5)

! Find neighor using neigh array because for boundary, neigh array
! points to current element so no need to check for existence of the
! neighbor element
nghElem = int((neigh(( stencil%cxdirinv( idir)-1)* nsize+ ielem)-1)/ nscalars)+1

! if neighbor element is my element then my neighbor is boundary
! use Forward difference to approximate gradient
if (nghElem == iElem) coeffForGrad = 1.0_rk

! map neighbor direction to treelm stencil to identify the face
iMeshDir = stencil%map(iDir)
select case(iMeshDir)
case (qN00) !-x
iFace = iFace+1
me%neighPos(iElem,1,1) = nghElem
case (q100) !+x
iFace = iFace+1
me%neighPos(iElem,1,2) = nghElem
case (q0N0) !-y
iFace = iFace+1
me%neighPos(iElem,2,1) = nghElem
case (q010) !+y
iFace = iFace+1
me%neighPos(iElem,2,2) = nghElem
case (q00N) !-z
iFace = iFace+1
me%neighPos(iElem,3,1) = nghElem
case (q001) !+z
iFace = iFace+1
me%neighPos(iElem,3,2) = nghElem
end select
!write(dbgUnit(1),*) 'iFace: ', iFace, 'iMeshDir: ', iMeshDir, 'coeff:', coeffForGrad
!write(dbgUnit(1),*) 'neighIDX: ', nghElem, ' neighID: ', levelDesc%total(nghElem)
! exit iDir loop
if (iFace == nFaces) exit neighLoop
end do neighLoop
end do elemLoop
!      write(logUnit(1),*) 'Done with fill grad data.'
!      write(dbgUnit(1),*) 'Done with fill grad data.'
!      flush(dbgUnit(1))