! Copyright (c) 2019 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.
! ****************************************************************************** !
!> This module contains function to compute eddy viscosity for
!! Wall-Adapting Local Eddy-Viscosity turbulence
!! model.
!! This implementation follows the LES described by Weickert et al.
!! Weickert, M., Teike, G., Schmidt, O., & Sommerfeld, M. (2010).
!! Investigation of the LES WALE turbulence model within the lattice Boltzmann
!! framework. Computers and Mathematics with Applications, 59(7), 2200–2214.
!! author: Kannan Masilamani
module mus_WALE_module
  ! treelm modules
  use env_module,                    only: rk, eps
  use tem_param_module,              only: div1_2, div1_3
  use tem_compileconf_module,        only: vlen

  ! musubi modules
  use mus_turbulence_module,         only: mus_turbulence_config_type
  use mus_gradData_module,           only: mus_gradData_type
  use mus_derivedQuantities_module2, only: getGradU

  implicit none
  private

  public :: mus_turbVisc_WALE_3D
  public :: mus_turbVisc_WALE_2D

contains

  ! ************************************************************************** !
  !> Calculate eddy viscosity with WALE (Wall-Adapting Local Eddy-viscosity)
  !! model
  !! \todo add reference and formula
  subroutine mus_turbVisc_WALE_3D(turbVisc, turbConfig, gradData, auxField, &
    &                             velPos, nSolve, nAuxScalars, dxL, dtL)
    ! --------------------------------------------------------------------------
    !> output: turbulent viscosity
    real(kind=rk), intent(out) :: turbVisc(:)
    !> Contains turbulenct coefficients
    type(mus_turbulence_config_type), intent(in) :: turbConfig
    !> gradient data
    type(mus_gradData_type), intent(in) :: gradData
    !> Auxiliary field variable array
    real(kind=rk), intent(in) :: auxField(:)
    !> position of velocity components in auxField
    integer, intent(in) :: velPos(3)
    !> Number of element to solve in this level
    integer, intent(in) :: nSolve
    !> number of scalars in auxField array
    integer, intent(in) :: nAuxScalars
    !> current level lattice element size
    real(kind=rk), intent(in) :: dxL
    !> current level lattice time step size
    real(kind=rk), intent(in) :: dtL
    ! --------------------------------------------------------------------------
    integer :: iElem
    real(kind=rk) :: SR(6), Sd(6), oneThird_trSd, Sd_sqr, SR_sqr, OP1, OP2
    real(kind=rk) :: visc_coeff
    !> gradient of velocity
    real(kind=rk):: gradU(3,3,vlen)
    real(kind=rk):: gradU_sqr(3,3,vlen)
    integer :: ndims
    integer :: leftngh(3), rightngh(3)
    real(kind=rk) :: leftvel(3,3), rightvel(3,3)
    integer :: nChunks, iChunks, nChunkElems, low_bound, elempos
    ! --------------------------------------------------------------------------
    ! viscosity coeff
    visc_coeff = (turbConfig%coeff%C_w * dxL )**2.0_rk
    ndims = 3
    nChunks = ceiling(real(nSolve, kind=rk) / real(vlen, kind=rk))

    do iChunks = 1, nChunks
      ! calculate the end  number of iElem loop
      nChunkElems = min(vlen, nSolve - ((iChunks - 1) * vlen))
      low_bound = (iChunks - 1) * vlen

      gradU(:,:,1:nChunkElems) = getGradU( auxField = auxField,    &
        &                   gradData           = gradData,    &
        &                   velPos             = velPos,      &
        &                   nAuxScalars        = nAuxScalars, &
        &                   nDims              = 3,           &
        &                   nSolve             = nChunkElems, &
        &                   elemOffset         = low_bound    )

      do iElem = 1, nChunkElems
        gradU_sqr(:,:,iElem) = matmul(gradU(:,:,iElem), gradU(:,:,iElem))
      end do !iElem

      do iElem = 1, nChunkElems
        elempos = low_bound + iElem
        ! traceless symmetric part of the square of the velocity gradient tensor
        ! Sd_ij = 1/2(du_k/dx_i du_j/dx_k + du_k/dx_j du_i/dx_k)
        !       - 1/3\delta_ij du_k/dx_l du_l/dx_k
        oneThird_trSd = (gradU_sqr(1,1,iElem) + gradU_sqr(2,2,iElem)  &
          &             + gradU_sqr(3,3,iElem))*div1_3
        Sd(1) = gradU_sqr(1,1,iElem) - oneThird_trSd                  !XX
        Sd(2) = gradU_sqr(2,2,iElem) - oneThird_trSd                  !YY
        Sd(3) = gradU_sqr(3,3,iElem) - oneThird_trSd                  !ZZ
        Sd(4) = 0.5_rk*(gradU_sqr(1,2,iElem)+gradU_sqr(2,1,iElem))    !XY
        Sd(5) = 0.5_rk*(gradU_sqr(2,3,iElem)+gradU_sqr(3,2,iElem))    !YZ
        Sd(6) = 0.5_rk*(gradU_sqr(1,3,iElem)+gradU_sqr(3,1,iElem))    !XZ

        ! double inner product of Sd: Sd_ij Sd_ij
        Sd_sqr = Sd(1)**2 + Sd(2)**2 + Sd(3)**2            &
          &    + 2.0_rk * ( Sd(4)**2 + Sd(5)**2 + Sd(6)**2 )

        ! symmetric strain rate tensors
        SR(1) = gradU(1,1,iElem)                               !XX
        SR(2) = gradU(2,2,iElem)                               !YY
        SR(3) = gradU(3,3,iElem)                               !ZZ
        SR(4) = (gradU(1,2,iElem)+gradU(2,1,iElem))*0.5_rk     !XY
        SR(5) = (gradU(2,3,iElem)+gradU(3,2,iElem))*0.5_rk     !YZ
        SR(6) = (gradU(1,3,iElem)+gradU(3,1,iElem))*0.5_rk     !XZ

        ! double inner product of tensor
        SR_sqr = SR(1)**2 + SR(2)**2 + SR(3)**2            &
          &    + 2.0_rk * ( SR(4)**2 + SR(5)**2 + SR(6)**2 )

        ! sub-grid scale kinetic energy
        ! k_sgs = (C_w^2 * dx /C_k)^2 (OP1/OP2)^2
        ! subgrid scale eddy viscosity
        ! nu_kgs = C_k dx sqrt(k_sgs) = (C_w * dx)^2 (OP1/OP2)
        ! OP1 = (Sd_ij Sd_ij)^(3/2)
        ! OP2 = (SR_ij SR_ij)^(5/2) + (Sd_ij Sd_ij)^(5/4)
        ! Add small fraction to denominator to avoid division by zero
        OP1 = Sd_sqr**1.5_rk
        OP2 = SR_sqr**2.5_rk + Sd_sqr**1.25_rk + eps
        ! turbulent viscosity
        turbVisc(elempos) = visc_coeff * (OP1/OP2) / dtL

      end do !iElem
    end do !iChunks

  end subroutine mus_turbVisc_WALE_3D
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> Calculate eddy viscosity with WALE (Wall-Adapting Local Eddy-viscosity)
  !! model
  !! \todo add reference and formula
  subroutine mus_turbVisc_WALE_2D(turbVisc, turbConfig, gradData, auxField, &
    &                             velPos, nSolve, nAuxScalars, dxL, dtL)
    ! --------------------------------------------------------------------------
    !> output: turbulent viscosity
    real(kind=rk), intent(out) :: turbVisc(:)
    !> Contains turbulenct coefficients
    type(mus_turbulence_config_type), intent(in) :: turbConfig
    !> gradient data
    type(mus_gradData_type), intent(in) :: gradData
    !> Auxiliary field variable array
    real(kind=rk), intent(in) :: auxField(:)
    !> position of velocity components in auxField
    integer, intent(in) :: velPos(3)
    !> Number of element to solve in this level
    integer, intent(in) :: nSolve
    !> number of scalars in auxField array
    integer, intent(in) :: nAuxScalars
    !> current level lattice element size
    real(kind=rk), intent(in) :: dxL
    !> current level lattice time step size
    real(kind=rk), intent(in) :: dtL
    ! --------------------------------------------------------------------------
    integer :: iElem
    real(kind=rk) :: gradU(2,2,vlen)
    real(kind=rk) :: gradU_sqr(2,2,vlen)
    real(kind=rk) :: SR(3), Sd(3), onehalf_trSd, Sd_sqr, SR_sqr, OP1, OP2
    integer :: ndims
    integer :: leftngh(2), rightngh(2)
    real(kind=rk) :: leftvel(2,2), rightvel(2,2)
    integer :: nChunks, iChunks, nChunkElems, low_bound, elempos
    real(kind=rk) :: visc_coeff
    !> gradient of velocity
    ! --------------------------------------------------------------------------
    visc_coeff = (turbConfig%coeff%C_w * dxL )**2.0_rk
    ndims = 2

    do iChunks = 1, nChunks
      ! calculate the end  number of iElem loop
      nChunkElems = min(vlen, nSolve - ((iChunks - 1) * vlen))
      low_bound = (iChunks - 1) * vlen

      gradU(:,:,1:nSolve) = getGradU( auxField = auxField,    &
        &                   gradData           = gradData,    &
        &                   velPos             = velPos,      &
        &                   nAuxScalars        = nAuxScalars, &
        &                   nDims              = 2,           &
        &                   nSolve             = nChunkElems, &
        &                   elemOffset         = low_bound    )

      do iElem = 1, nChunkElems
      ! square of velocity gradient. gradU . gradU
        gradU_sqr(:,:,iElem) = matmul(gradU(:,:,iElem), gradU(:,:,iElem))
      end do !iElem
      
      do iElem = 1, nChunkElems
        elempos = low_bound + iElem
        ! traceless symmetric part of the square of the velocity gradient tensor
        ! Sd_ij = 1/2(du_k/dx_i du_j/dx_k + du_k/dx_j du_i/dx_k)
        !       - 1/3\delta_ij du_k/dx_l du_l/dx_k
        onehalf_trSd = (gradU_sqr(1,1,iElem) + gradU_sqr(2,2,iElem))*div1_2
        Sd(1) = gradU_sqr(1,1,iElem) - onehalf_trSd            !XX
        Sd(2) = gradU_sqr(2,2,iElem) - onehalf_trSd            !YY
        Sd(3) = 0.5_rk*(gradU_sqr(1,2,iElem)+gradU_sqr(2,1,iElem))    !XY

        ! double inner product of Sd: Sd_ij Sd_ij
        Sd_sqr = Sd(1)**2 + Sd(2)**2 + 2.0_rk*Sd(3)**2

        ! symmetric strain rate tensors
        SR(1) = gradU(1,1,iElem)
        SR(2) = gradU(2,2,iElem)
        SR(3) = (gradU(1,2,iElem)+gradU(2,1,iElem))*0.5_rk

        ! double inner product of tensor
        SR_sqr = SR(1)**2 + SR(2)**2 + 2.0_rk*SR(3)**2

        ! sub-grid scale kinetic energy
        ! k_sgs = (C_w^2 * dx /C_k)^2 (OP1/OP2)^2
        ! subgrid scale eddy viscosity
        ! nu_kgs = C_k dx sqrt(k_sgs) = (C_w * dx)^2 (OP1/OP2)
        ! OP1 = (Sd_ij Sd_ij)^(3/2)
        ! OP2 = (SR_ij SR_ij)^(5/2) + (Sd_ij Sd_ij)^(5/4)
        ! Add small fraction to denominator to avoid division by zero
        OP1 = Sd_sqr**1.5_rk
        OP2 = SR_sqr**2.5_rk + Sd_sqr**1.25_rk + eps

        ! turbulent viscosity
        turbVisc(elempos) = visc_coeff * (OP1/OP2) / dtL

      end do !iElem
    end do !iChunks

  end subroutine mus_turbVisc_WALE_2D
  ! ************************************************************************** !

end module mus_WALE_module
