! 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
?? include 'header/lbm_deriveMacros.inc'
module mus_WALE_module
  ! treelm modules
  use env_module,                    only: rk, eps
  use tem_param_module,              only: div1_2, div1_3

  ! 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) :: gradU_sqr(3,3)
    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)
    ! --------------------------------------------------------------------------
    ! viscosity coeff
    visc_coeff = (turbConfig%coeff%C_w * dxL )**2.0_rk

    do iElem = 1, nSolve

      gradU = getGradU(iElem, auxField, gradData, velPos, nAuxScalars, 3)

      ! square of velocity gradient. gradU . gradU
      gradU_sqr = matmul(gradU, gradU)

      ! 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) + gradU_sqr(2,2) + gradU_sqr(3,3))*div1_3
      Sd(1) = gradU_sqr(1,1) - oneThird_trSd            !XX
      Sd(2) = gradU_sqr(2,2) - oneThird_trSd            !YY
      Sd(3) = gradU_sqr(3,3) - oneThird_trSd            !ZZ
      Sd(4) = 0.5_rk*(gradU_sqr(1,2)+gradU_sqr(2,1))    !XY
      Sd(5) = 0.5_rk*(gradU_sqr(2,3)+gradU_sqr(3,2))    !YZ
      Sd(6) = 0.5_rk*(gradU_sqr(1,3)+gradU_sqr(3,1))    !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)                         !XX
      SR(2) = gradU(2,2)                         !YY
      SR(3) = gradU(3,3)                         !ZZ
      SR(4) = (gradU(1,2)+gradU(2,1))*0.5_rk     !XY
      SR(5) = (gradU(2,3)+gradU(3,2))*0.5_rk     !YZ
      SR(6) = (gradU(1,3)+gradU(3,1))*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(iElem) = visc_coeff * (OP1/OP2) / dtL

    end do

  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_sqr(2,2)
    real(kind=rk) :: SR(3), Sd(3), onehalf_trSd, Sd_sqr, SR_sqr, OP1, OP2
    real(kind=rk) :: visc_coeff
    !> gradient of velocity
    real(kind=rk) :: gradU(2,2)
    ! --------------------------------------------------------------------------
    visc_coeff = (turbConfig%coeff%C_w * dxL )**2.0_rk

    do iElem = 1, nSolve

      gradU = getGradU(iElem, auxField, gradData, velPos, nAuxScalars, 2)

      ! square of velocity gradient. gradU . gradU
      gradU_sqr = matmul(gradU, gradU)

      ! 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) + gradU_sqr(2,2))*div1_2
      Sd(1) = gradU_sqr(1,1) - onehalf_trSd            !XX
      Sd(2) = gradU_sqr(2,2) - onehalf_trSd            !YY
      Sd(3) = 0.5_rk*(gradU_sqr(1,2)+gradU_sqr(2,1))    !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)
      SR(2) = gradU(2,2)
      SR(3) = (gradU(1,2)+gradU(2,1))*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(iElem) = visc_coeff * (OP1/OP2) / dtL

    end do

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

end module mus_WALE_module
