! Copyright (c) 2016-2017 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2016 Verena Krupp <verena.krupp@uni-siegen.de>
! Copyright (c) 2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2017 Raphael Haupt <raphael.haupt@uni-siegen.de>
! Copyright (c) 2020 Peter Vitt <peter.vitt2@uni-siegen.de>
! Copyright (c) 2021-2022 Gregorio Gerardo Spinelli <gregoriogerardo.spinelli@dlr.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
! ****************************************************************************** !
!> This module provides the routine for applying operators. Currently it is
!! only implemented for 3D and needs to be extended to 2d
module mus_operation_var_module
  use, intrinsic :: iso_c_binding,   only: c_f_pointer, c_ptr
  use env_module,                    only: rk, long_k, labellen

  use tem_logging_module,            only: logUnit
  use tem_debug_module,              only: dbgUnit
  use tem_varSys_module,             only: tem_varSys_type,                    &
    &                                      tem_varSys_op_type,                 &
    &                                      tem_varSys_proc_element,            &
    &                                      tem_varSys_solverData_evalElem_type
  use tem_time_module,               only: tem_time_type
  use treelmesh_module,              only: treelmesh_type
  use tem_topology_module,           only: tem_levelOf
  use tem_grow_array_module,         only: grw_intArray_type, append, truncate
  use tem_operation_var_module,      only: tem_opVar_fill_inputIndex, &
    &                                      tem_varSys_op_Data_type

  ! include musubi modules
  use mus_varSys_module,             only: mus_varSys_data_type,       &
    &                                      mus_varSys_solverData_type, &
    &                                      mus_get_new_solver_ptr

  implicit none

  private

  public :: mus_opVar_setupIndices
  public :: mus_set_opVar_getElement
  public :: mus_opVar_gradU_forElement
  public :: mus_opVar_vorticity_forElement
  public :: mus_opVar_QCriterion_forElement

contains

  ! ************************************************************************** !
  !> Routine to store musubi varSys Data in operation variable solver_bundle.
  !! Unline Ateles, Musubi operations does not require any special treatment so
  !! it uses to generic routines in treelm
  subroutine mus_set_opVar_getElement( solData_evalElem, fun )
    ! --------------------------------------------------------------------------
    !> Description on how to set the element retrieval function for stfuns.
    class(tem_varSys_solverData_evalElem_type), intent(in) :: solData_evalElem

    !> 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.
    type(tem_varSys_op_type), intent(inout) :: fun
    ! --------------------------------------------------------------------------
    type(tem_varSys_op_Data_type), pointer :: fptr
    type(mus_varSys_solverData_type), pointer :: fSDptr
    ! --------------------------------------------------------------------------

    write(logunit(10),*) "Setting different solver_bundle and " &
      & // " get_element routine for variable at position ",    &
      & fun%myPos
    call C_F_Pointer(fun%method_data, fptr)
    call c_f_pointer(solData_evalElem%solver_bundle, fSDptr)
    fptr%solver_bundle = mus_get_new_solver_ptr( fSDptr )

  end subroutine mus_set_opVar_getElement
  ! ************************************************************************** !

  ! ************************************************************************** !
  recursive subroutine mus_opVar_setupIndices( fun, varSys, point, offset_bit,&
    &                                iLevel, tree, nPnts, idx        )
    ! -------------------------`-----------------------------------------------!
    !> 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

    !> List of space coordinate points to store as growing array in
    !! method_data
    real(kind=rk), intent(in) :: point(:,:)

    !> Offset bit encoded as character for every point.
    !!
    !! Offset integer coord(3) is converted into a character with
    !! offset_bit = achar( (coord(1)+1) + (coord(2)+1)*4 + (coord(3)+1)*16 )
    !! Backward transformation form character to 3 integer:
    !! coord(1) = mod(ichar(offset_bit),4) - 1
    !! coord(2) = mod(ichar(offset_bit),16)/4 - 1
    !! coord(3) = ichar(offset_bit)/16 - 1
    !!
    !! If not present default is to center i.e offset_bit = achar(1+4+16)
    character, optional, intent(in) :: offset_bit(:)

    !> Level to which input points belong to
    integer, intent(in) :: iLevel

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

    !> Number of points to add in method_data of this variable
    integer, intent(in) :: nPnts

    !> Index of points in the growing array and variable val array.
    !! Size: nPoints
    !!
    !! This must be stored in boundary or source depends on who
    !! calls this routine.
    !! This index is required to return a value using getValOfIndex.
    integer, intent(out) :: idx(:)
    ! -------------------------------------------------------------------------!
    type(mus_varSys_data_type), pointer :: fPtr
    integer :: iPnt, iDep, nVals_prev
    type(grw_intArray_type), allocatable :: inputIndex_loc(:)
    integer, allocatable :: idxPerPnt(:)
    ! --------------------------------------------------------------------------
    write(dbgUnit(4),*) 'setup indices for the points of derived variable ', &
      &                 trim(varSys%varname%val(fun%myPos))

    call C_F_POINTER( fun%method_Data, fPtr )

    ! allcoate the index array for all inpits
    if (.not. allocated(fPtr%opData%input_pntIndex)) then
      allocate( fPtr%OpData%input_pntIndex(fun%nInputs) )
    end if

    ! allocate temporary inputIndex with size of nInputs and initialize
    ! growing array with length nPnts
    allocate(inputIndex_loc(fun%nInputs))

    ! store which is the last entry in the indexLvl to contiguous fill the index
    ! array string from this position
    ! all input variables get the same points, we just take the nVals entry
    ! from the first input variable
    nVals_prev = fPtr%OpData%input_pntIndex(1)%indexLvl(iLevel)%nVals

    ! Now fill in the index arrays for the inputs
    call tem_opVar_fill_inputIndex( fun        = fun,           &
      &                             varSys     = varSys,        &
      &                             point      = point,         &
      &                             offset_bit = offset_bit,    &
      &                             iLevel     = iLevel,        &
      &                             tree       = tree,          &
      &                             nPnts      = nPnts,         &
      &                             inputIndex = inputIndex_loc )

    ! fill the index array of the derived variable, it starts with the first
    ! entry in this call = nVals_prev and is continguous until nVals_prev+nVals
    allocate(idxPerPnt(fun%nInputs))
    idx = 0
    do iPnt = 1, nPnts
      do iDep = 1, fun%nInputs
        idxPerPnt(iDep) = inputIndex_loc(iDep)%val(iPnt)
      end do
      ! set index only when any of dependent variable has valid index
      if (any(idxPerPnt > 0)) then
        do iDep = 1, fun%nInputs
          call append( me = fPtr%opData%input_pntIndex(iDep)%indexLvl(iLevel), &
            &         val = inputIndex_loc(iDep)%val(iPnt)                     )
        end do
        ! set index to last position in input_pntIndex of dep var 1 of
        ! indexLvl of iLevel
        idx(iPnt) = fPtr%opData%input_pntIndex(1)%indexLvl(iLevel)%nVals
      end if
    end do

    do iDep = 1, fun%nInputs
      call truncate (fPtr%opData%input_pntIndex(iDep)%indexLvl(iLevel) )
    end do

  end subroutine mus_opVar_setupIndices
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> This routine returns the velocity gradient from velocity in auxField
  !!
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_element]].
  recursive subroutine mus_opVar_gradU_forElement(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(:)
    ! -------------------------------------------------------------------- !
    integer :: statePos, iElem, iLevel
    type(mus_varSys_data_type), pointer :: fPtr
    real(kind=rk) :: gradU(3,3,1)
    integer :: velPos(3), elemOff, nDims
    ! -------------------------------------------------------------------- !
    call C_F_POINTER( fun%method_Data, fPtr )

    associate ( levelPointer => fPtr%solverData%geometry%levelPointer,   &
      &         gradData     => fPtr%solverData%scheme%gradData,         &
      &         auxField     => fPtr%solverData%scheme%auxField,         &
      &         derVarPos    => fPtr%solverData%scheme%derVarPos(1),     &
      &         stencil      => fPtr%solverData%scheme%layout%fStencil,  &
      &         Grad         => fPtr%solverData%scheme%Grad              )


      nDims = stencil%nDims
      velPos = varSys%method%val(derVarPos%velocity)%auxField_varPos(1:3)

      ! res is always AOS layout
      res = 0.0_rk
      do iElem = 1, nElems
        ! to acess Grad calculation pointers
        iLevel = tem_levelOf( tree%treeID( elemPos(iElem) ) )

        ! if state array is defined level wise then use levelPointer(pos)
        ! to access state array
        statePos = levelPointer( elemPos(iElem) )

        ! compute velocity gradient from velocity field stored in auxField
        gradU = 0.0_rk

        gradU(:nDims,:nDims,:1) = Grad%U_ptr(                        &
          &                   auxField    = auxField(iLevel)%val(:), &
          &                   gradData    = gradData(iLevel),        &
          &                   velPos      = velPos,                  &
          &                   nAuxScalars = varSys%nAuxScalars,      &
          &                   nDims       = nDims,                   &
          &                   nSolve      = 1,                       &
          &                   elemOffset  = statePos - 1             )

        ! store gradU in res array
        elemOff= (iElem-1)*9
        res( elemOff + 1 ) = gradU(1,1,1) !dudx
        res( elemOff + 2 ) = gradU(1,2,1) !dudy
        res( elemOff + 3 ) = gradU(1,3,1) !dudz
        res( elemOff + 4 ) = gradU(2,1,1) !dvdx
        res( elemOff + 5 ) = gradU(2,2,1) !dvdy
        res( elemOff + 6 ) = gradU(2,3,1) !dvdz
        res( elemOff + 7 ) = gradU(3,1,1) !dwdx
        res( elemOff + 8 ) = gradU(3,2,1) !dwdy
        res( elemOff + 9 ) = gradU(3,3,1) !dwdz
      end do !iElem

    end associate

  end subroutine mus_opVar_gradU_forElement
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> This routine computes vorticity from curl of velocity in auxField.
  !! $$vorticity = \nabla \times u$$
  !!
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_element]].
  recursive subroutine mus_opVar_vorticity_forElement(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(:)
    ! -------------------------------------------------------------------- !
    integer :: statePos, iElem, iLevel
    type(mus_varSys_data_type), pointer :: fPtr
    real(kind=rk) :: gradU(3,3,1)
    real(kind=rk) :: vorticity(3)
    integer :: velPos(3), elemOff,nDims
    ! -------------------------------------------------------------------- !
    call C_F_POINTER( fun%method_Data, fPtr )

    associate ( levelPointer => fPtr%solverData%geometry%levelPointer,   &
      &         gradData     => fPtr%solverData%scheme%gradData,         &
      &         auxField     => fPtr%solverData%scheme%auxField,         &
      &         derVarPos    => fPtr%solverData%scheme%derVarPos(1),     &
      &         stencil      => fPtr%solverData%scheme%layout%fStencil,  &
      &         Grad         => fPtr%solverData%scheme%Grad              )

      nDims = stencil%nDims
      velPos = varSys%method%val(derVarPos%velocity)%auxField_varPos(1:3)

      ! res is always AOS layout
      res = 0.0_rk
      vorticity = 0.0_rk
      do iElem = 1, nElems

        iLevel = tem_levelOf( tree%treeID( elemPos(iElem) ) )

        ! if state array is defined level wise then use levelPointer(pos)
        ! to access state array
        statePos = levelPointer( elemPos(iElem) )

        ! compute velocity gradient from velocity field stored in auxField
        gradU(:nDims,:nDims,:) = Grad%U_ptr(                     &
          &               auxField    = auxField(iLevel)%val(:), &
          &               gradData    = gradData(iLevel),        &
          &               velPos      = velPos,                  &
          &               nAuxScalars = varSys%nAuxScalars,      &
          &               nDims       = nDims,                   &
          &               nSolve      = 1,                       &
          &               elemOffset  = statePos - 1             )

        ! compute vorticity
        if (nDims > 1) then
          vorticity(3) = gradU(2,1,1) - gradU(1,2,1) !dvdx - dudy
          if (nDims > 2) then
            vorticity(1) = gradU(3,2,1) - gradU(2,3,1) !dwdy - dvdz
            vorticity(2) = gradU(1,3,1) - gradU(3,1,1) !dxdz - dwdx
          end if
        end if

        ! store gradU in res array
        elemOff= (iElem-1)*3
        res( elemOff + 1 : elemOff + 3 ) = vorticity(:)
      end do !iElem

    end associate

  end subroutine mus_opVar_vorticity_forElement
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> This routine computes Q-criterion from velocity in auxField.
  !! $$Q = 1/2 ( (tr(\nabla u))^2 - tr(\nabla u \cdot \nabla u) )
  !! = 1/2 |\Omega|^2 - |S|^2 $$,
  !! where $$\Omega$$ and S are asymmetric (vorticity tensor) and
  !! symmetric (rate of strain) part of velocity gradient.
  !! i.e $$\Omega = 1/2(du_i/dx_j - du_j/dx_i)$$ and
  !! $$S=1/2(du_i/dx_j + du_j/dx_i)$$.
  !!
  !! Ref: http://www.ctr.maths.lu.se/media/thesis/2013/FMA820.pdf
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_element]].
  recursive subroutine mus_opVar_QCriterion_forElement(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(:)
    ! -------------------------------------------------------------------- !
    integer :: statePos, iElem, iLevel, iDir
    type(mus_varSys_data_type), pointer :: fPtr
    real(kind=rk):: gradU(3,3,1), gradU_sqr(3,3)
    real(kind=rk) :: tr_gradU, tr_gradU_sqr
    integer :: velPos(3), nDims
    ! -------------------------------------------------------------------- !
    call C_F_POINTER( fun%method_Data, fPtr )

    associate ( levelPointer => fPtr%solverData%geometry%levelPointer,  &
      &         gradData     => fPtr%solverData%scheme%gradData,        &
      &         auxField     => fPtr%solverData%scheme%auxField,        &
      &         derVarPos    => fPtr%solverData%scheme%derVarPos(1),    &
      &         stencil      => fPtr%solverData%scheme%layout%fStencil, &
      &         Grad         => fPtr%solverData%scheme%Grad             )

      nDims = stencil%nDims
      velPos = varSys%method%val(derVarPos%velocity)%auxField_varPos(1:3)

      ! res is always AOS layout
      res = 0.0_rk
      do iElem = 1, nElems

        iLevel = tem_levelOf( tree%treeID( elemPos(iElem) ) )

        ! if state array is defined level wise then use levelPointer(pos)
        ! to access state array
        statePos = levelPointer( elemPos(iElem) )
        ! compute velocity gradient from velocity field stored in auxField
        gradU(:nDims,:nDims,:) = Grad%U_ptr(              &
          &               auxField    = auxField(iLevel)%val(:), &
          &               gradData    = gradData(iLevel),        &
          &               velPos      = velPos,                  &
          &               nAuxScalars = varSys%nAuxScalars,      &
          &               nDims       = nDims,                   &
          &               nSolve      = 1,                       &
          &               elemOffset  = statePos - 1             )


        ! square of velocity gradient. gradU . gradU
        gradU_sqr(:nDims, :nDims) = matmul( gradU(:nDims,:nDims,1), &
          &                                 gradU(:nDims,:nDims,1)  )
        ! trace of gradU_sqr
        tr_gradU_sqr = gradU_sqr(1,1)
        ! trace of gradU
        tr_gradU = gradU(1,1,1)
        do iDir=2,nDims
          tr_gradU_sqr = tr_gradU_sqr + gradU_sqr(iDir,iDir)
          tr_gradU = tr_gradU + gradU(iDir, iDir, 1)
        end do

        ! compute and store Q-criterion
        res( iElem ) = 0.5_rk * (tr_gradU*tr_gradU - tr_gradU_sqr)

      end do !iElem

    end associate

  end subroutine mus_opVar_QCriterion_forElement
  ! ************************************************************************** !

end module mus_operation_var_module
! **************************************************************************** !
