! 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 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
?? include 'header/lbm_deriveMacros.inc'
module mus_Vreman_module
  ! treelm modules
  use env_module,                    only: rk
  use tem_float_module,              only: operator(.flt.)

  ! 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_Vreman_3D
  public :: mus_turbVisc_Vreman_2D

contains

  ! ************************************************************************** !
  !> Calculate eddy viscosity with Vreman model for 3D stencil
  !! Fortran implementation of this model:
  !! http://www.vremanresearch.nl/Vreman_Subgridmodel_Fortran.txt
  !! \todo add reference and formula
  subroutine mus_turbVisc_Vreman_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
    ! --------------------------------------------------------------------------
    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
    integer :: iElem
    !> gradient of velocity
    real(kind=rk) :: gradU(3,3)
    ! --------------------------------------------------------------------------
    ! viscosity coeff
    visc_coeff = (turbConfig%coeff%C_v * dxL )**2.0_rk

    do iElem = 1, nSolve

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

      ! 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
        turbVisc(iElem) = 0.0_rk
      else
        turbVisc(iElem) = visc_coeff*sqrt(bbeta/gradU_sqr) / dtL
      end if
    end do

  end subroutine mus_turbVisc_Vreman_3D
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> Calculate eddy viscosity with Vreman model for 2D stencil
  !! model
  !! \todo add reference and formula
  subroutine mus_turbVisc_Vreman_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) :: b_11, b_12, b_22
    real(kind=rk) :: gradU_sqr, visc_coeff
    real(kind=rk) :: bbeta
    !> gradient of velocity
    real(kind=rk) :: gradU(2,2)
    ! --------------------------------------------------------------------------
    ! viscosity coeff
    visc_coeff = (turbConfig%coeff%C_v * dxL )**2.0_rk

    do iElem = 1, nSolve

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

      ! 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
        turbVisc(iElem) = 0.0_rk
      else
        turbVisc(iElem) = visc_coeff*sqrt(bbeta/gradU_sqr) / dtL
      end if
    end do

  end subroutine mus_turbVisc_Vreman_2D
  ! ************************************************************************** !

end module mus_Vreman_module
