! See copyright notice in the COPYRIGHT file.
! ****************************************************************************** !
!> 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, labelLen
  use tem_param_module,     only: div1_2, div1_3
  use tem_tools_module,     only: tem_horizontalSpacer
  use tem_aux_module,       only: tem_abort
  use tem_logging_module,   only: logUnit
  use tem_float_module,     only: operator(.flt.)
  
  ! musubi modules
  use mus_turbulence_module, only: mus_turbulence_config_type

  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
  pure function mus_turbVisc_WALE_3D(turbConfig, gradU, dxL, dtL) result(res)
    ! --------------------------------------------------------------------------
    !> Contains turbulenct coefficients
    type(mus_turbulence_config_type), intent(in) :: turbConfig
    !> gradient of velocity
    real(kind=rk), intent(in) :: gradU(:,:)
    !> current level lattice element size
    real(kind=rk), intent(in) :: dxL
    !> current level lattice time step size
    real(kind=rk), intent(in) :: dtL
    !> output: turbulent viscosity
    real(kind=rk) :: res
    ! --------------------------------------------------------------------------
    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
    ! --------------------------------------------------------------------------
    ! viscosity coeff
    visc_coeff = (turbConfig%coeff%C_w * dxL )**2.0_rk
    
    ! 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*( Sd(4)**2 + Sd(5)**2 + Sd(6)**2 )

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

    ! double inner product of tensor
    SR_sqr = SR(1)**2 + SR(2)**2 + SR(3)**2    &
      &    + 2*(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)
    ! OP1 = (Sd_ij Sd_ij)^3
    ! OP2 = ( (SR_ij SR_ij)^(5/2) + (Sd_ij Sd_ij)^(5/4) )^2
    OP1 = Sd_sqr**3.0_rk
    OP2 = (SR_sqr**2.5_rk + Sd_sqr**1.25_rk)**2.0_rk
     
    ! subgrid scale eddy viscosity 
    ! nu_kgs = C_k dx sqrt(k_sgs) = (C_w * dx)^2 sqrt(OP1/OP2)
    ! if denominator is two small set S_ij S_ij = 0
    if (OP2 .flt. 1e-12_rk ) then
      res = 0.0_rk
    else  
      ! turbulent viscosity
      res = visc_coeff * sqrt(OP1/OP2) / dtL
    end if  

  end function mus_turbVisc_WALE_3D
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> Calculate eddy viscosity with WALE (Wall-Adapting Local Eddy-viscosity) 
  !! model
  !! \todo add reference and formula
  pure function mus_turbVisc_WALE_2D(turbConfig, gradU, dxL, dtL) result(res)
    ! --------------------------------------------------------------------------
    !> Contains turbulenct coefficients
    type(mus_turbulence_config_type), intent(in) :: turbConfig
    !> gradient of velocity
    real(kind=rk), intent(in) :: gradU(:,:)
    !> current level lattice element size
    real(kind=rk), intent(in) :: dxL
    !> current level lattice time step size
    real(kind=rk), intent(in) :: dtL
    !> output: turbulent viscosity
    real(kind=rk) :: res
    ! --------------------------------------------------------------------------
    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
    ! --------------------------------------------------------------------------
    visc_coeff = (turbConfig%coeff%C_w * dxL )**2.0_rk
    
    ! 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)
    ! OP1 = (Sd_ij Sd_ij)^3
    ! OP2 = ( (SR_ij SR_ij)^(5/2) + (Sd_ij Sd_ij)^(5/4) )^2
    OP1 = Sd_sqr**3.0_rk
    OP2 = (SR_sqr**2.5_rk + Sd_sqr**1.25_rk)**2.0_rk
     
    ! subgrid scale eddy viscosity 
    ! nu_kgs = C_k dx sqrt(k_sgs) = (C_w * dx)^2 sqrt(OP1/OP2)
    ! if denominator is two small set S_ij S_ij = 0
    if (OP2 .flt. 1e-12_rk ) then 
      !OP2 = (Sd_sqr**1.25_rk)**2.0_rk
      res = 0.0_rk
    else  
      ! turbulent viscosity
      res = visc_coeff * sqrt(OP1/OP2) / dtL
    end if  

  end function mus_turbVisc_WALE_2D
  ! ************************************************************************** !

end module mus_WALE_module
