! See copyright notice in the COPYRIGHT file.
! ****************************************************************************** !
!> This module contains function to compute eddy viscosity using
!! Vreman LES turbulence model.
!! Vreman, A. W. (2004). An eddy-viscosity subgrid-scale model for turbulent 
!! shear flow: Algebraic theory and applications. Physics of Fluids, 16(10), 
!! 3670–3681. 
!! model. 
!! author: Kannan Masilamani
module mus_Vreman_module
  ! treelm modules
  use env_module,           only: rk, labelLen
  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_Vreman_3D
  public :: mus_turbVisc_Vreman_2D

contains 

  ! ************************************************************************** !
  !> Calculate eddy viscosity with Vreman model 
  !! Fortran implementation of this model:
  !! http://www.vremanresearch.nl/Vreman_Subgridmodel_Fortran.txt
  !! \todo add reference and formula
  pure function mus_turbVisc_Vreman_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) :: b_11, b_12, b_13, b_22, b_23, b_33
    real(kind=rk) :: gradU_sqr, visc_coeff
    real(kind=rk) :: bbeta
    ! --------------------------------------------------------------------------
    ! viscosity coeff
    visc_coeff = (turbConfig%coeff%C_v * dxL )**2.0_rk

    ! beta_ij = gradU_ki gradU_kj 
    ! beta = gradU^T gradU
    b_11 = gradU(1,1)*gradU(1,1) + gradU(2,1)*gradU(2,1) + gradU(3,1)*gradU(3,1)
    b_12 = gradU(1,1)*gradU(1,2) + gradU(2,1)*gradU(2,2) + gradU(3,1)*gradU(3,2)
    b_13 = gradU(1,1)*gradU(1,3) + gradU(2,1)*gradU(2,3) + gradU(3,1)*gradU(3,3)
    b_22 = gradU(1,2)*gradU(1,2) + gradU(2,2)*gradU(2,2) + gradU(3,2)*gradU(3,2)
    b_23 = gradU(1,2)*gradU(1,3) + gradU(2,2)*gradU(2,3) + gradU(3,2)*gradU(3,3)
    b_33 = gradU(1,3)*gradU(1,3) + gradU(2,3)*gradU(2,3) + gradU(3,3)*gradU(3,3)
    
    ! double inner product of gradU
    gradU_sqr = gradU(1,1)**2 + gradU(1,2)**2 + gradU(1,3)**2 &
      &       + gradU(2,1)**2 + gradU(2,2)**2 + gradU(2,3)**2 &
      &       + gradU(3,1)**2 + gradU(3,2)**2 + gradU(3,3)**2

    ! numerator Bbeta
    bbeta=b_11*b_22-(b_12**2)+b_11*b_33-(b_13**2)+b_22*b_33-(b_23**2)

    if (bbeta .flt. 1e-12_rk) then
      res = 0.0_rk
    else
      res = visc_coeff*sqrt(bbeta/gradU_sqr) / dtL
    end if

  end function mus_turbVisc_Vreman_3D
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> Calculate eddy viscosity with Vreman model 
  !! model
  !! \todo add reference and formula
  pure function mus_turbVisc_Vreman_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) :: b_11, b_12, b_22
    real(kind=rk) :: gradU_sqr, visc_coeff
    real(kind=rk) :: bbeta
    ! --------------------------------------------------------------------------
    ! viscosity coeff
    visc_coeff = (turbConfig%coeff%C_v * dxL )**2.0_rk
    ! beta_ij = gradU_ki gradU_kj 
    ! beta = gradU^T gradU
    b_11 = gradU(1,1)*gradU(1,1) + gradU(2,1)*gradU(2,1)
    b_12 = gradU(1,1)*gradU(1,2) + gradU(2,1)*gradU(2,2)
    b_22 = gradU(1,2)*gradU(1,2) + gradU(2,2)*gradU(2,2)
    
    ! double inner product of gradU
    gradU_sqr = gradU(1,1)**2 + gradU(1,2)**2 &
      &       + gradU(2,1)**2 + gradU(2,2)**2

    ! numerator Bbeta
    bbeta=b_11*b_22-(b_12**2)

    if (bbeta .flt. 1e-12_rk) then
      res = 0.0_rk
    else
      res = visc_coeff*sqrt(bbeta/gradU_sqr) / dtL
    end if


  end function mus_turbVisc_Vreman_2D
  ! ************************************************************************** !

end module mus_Vreman_module
