! 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
?? include 'header/lbm_macros.inc'
module mus_gradData_module
  ! 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_stencil_module,  only: tem_stencilHeader_type
  use tem_param_module,    only: qN00, q100, q0N0, q010, q00N, q001
  use tem_construction_module, only: tem_levelDesc_type

  implicit none
  private

  public :: mus_gradData_type
  public :: mus_init_gradData

  !> Contains information required to compute gradient like position
  !! of six direct neigbors in the state array and coeff for the gradient
  !! operation
  type mus_gradData_type
    !> 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(:,:)
  end type mus_gradData_type

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)
    !---------------------------------------------------------------------------
    !> Gradient type
    type(mus_gradData_type), intent(out) :: me
    !> neighbor connectivity array
    integer, intent(in) :: neigh(:)
    !> levelDesc to access communication buffers of state array
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> stencil header
    type(tem_stencilHeader_type), intent(in) :: stencil
    !> 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
    real(kind=rk) :: coeffForGrad
    integer :: nghElem
    !---------------------------------------------------------------------------
    write(logUnit(3),*) 'Allocate gradient data ...'
    nDims = stencil%nDims
    ! allocate gradient data type
    allocate(me%neighPos( nSolve, nDims, 2 ) )
    me%neighPos = -1
    allocate(me%FDcoeff( nSolve, nDims ) )
    me%FDcoeff = -1.0_rk

    write(logUnit(1),*) 'Filling gradient data ...'
    ! fill gradient data with face neighbor information for turbulence model
!write(dbgUnit(1),*) 'Filling gradient data ...'
    ! Number of faces per element
    nFaces = nDims*2

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

      iFace = 0

      ! Loop over every link
      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)
        coeffForGrad = 0.5_rk

        ! 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 = ?NghElemIDX?(neigh, iDir, stencil, iElem, nScalars, nSize)

        ! 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
          me%FDcoeff(iElem,1) = coeffForGrad
        case (q100) !+x
          iFace = iFace+1
          me%neighPos(iElem,1,2) = nghElem
          me%FDcoeff(iElem,1) = coeffForGrad
        case (q0N0) !-y
          iFace = iFace+1
          me%neighPos(iElem,2,1) = nghElem
          me%FDcoeff(iElem,2) = coeffForGrad
        case (q010) !+y
          iFace = iFace+1
          me%neighPos(iElem,2,2) = nghElem
          me%FDcoeff(iElem,2) = coeffForGrad
        case (q00N) !-z
          iFace = iFace+1
          me%neighPos(iElem,3,1) = nghElem
          me%FDcoeff(iElem,3) = coeffForGrad
        case (q001) !+z
          iFace = iFace+1
          me%neighPos(iElem,3,2) = nghElem
          me%FDcoeff(iElem,3) = coeffForGrad
        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))
!      call tem_abort('Verify fill gradData')
  end subroutine mus_init_gradData
  ! ************************************************************************* !

end module mus_gradData_module
