! See copyright notice in the COPYRIGHT file.
! This utest program tests the getStrainRate_acoustic routine.
! It first initializes array of PDF, then calculate strain rate from it.
program mus_getStrainRate_acoustic_test
  use env_module,         only: rk
  use tem_general_module, only: tem_start, tem_finalize, tem_general_type

  use mus_scheme_layout_module, only: mus_scheme_layout_type, &
    &                                 mus_define_d3q19, &
    &                                 mus_define_d2q9
  use mus_derivedQuantities_module2, only: getShearRateTensor_acoustic, &
    &                                      getShearRateTensor_diffusive, &
    &                                      getEquilibriumIncomp, convPrePost, &
    &                                      secondMom

  implicit none

  logical :: error
  real(kind=rk) :: tolerance, max_error

  type( tem_general_type ) :: general
  type( mus_scheme_layout_type ) :: layout
  real(kind=rk), allocatable :: pdf(:), fEq(:), fNeq(:)
  real(kind=rk) :: rho0  = 1.0_rk
  real(kind=rk) :: omega = 1.5_rk
  real(kind=rk) :: nu
  real(kind=rk) :: S(3,3), strain_ref(6), strain_exp(6), error_rel
  real(kind=rk) :: strain_acs(6), strain_dff(6)
  integer :: nDims = 3
  real(kind=rk) :: rho, vel(3)

  integer :: iDir, QQ

  call tem_start( 'getStrainRate acoustic utest', 'utest', general )
  error = .false.
  tolerance = 1E-05
  max_error = 0.0_rk
  write(*,"(A, F14.9)") 'tolerance = ', tolerance

  write(*,"(A)")    'Flow parameters:'
  write(*,"(A, F10.6)") ' rho0  = ', rho0
  write(*,"(A, F10.6)") ' omega = ', omega
  nu = (1.0_rk / omega - 0.5_rk ) / 3.0_rk
  write(*,"(A, F10.6)") ' nu    = ', nu

  write(*,*) 'Initialize stencil'
  if ( nDims == 2 ) then
    call mus_define_d2q9( layout = layout, nElems = 1 )
  else
    call mus_define_d3q19( layout = layout, nElems = 1 )
  end if
  QQ = layout%fStencil%QQ

  ! initialize PDF
  write(*,*) 'Initialize strain_ref'
  strain_ref = [ -0.667092462303E-05,  -0.294694068991E-07, &
    &             0.335019701549E-05,  -0.315414531135E-02, &
    &             0.312250225676E-16,   0.312250225676E-16 ]
  write(*,*) 'Initialize PDF'
  allocate( pdf(QQ) )
  pdf(:) =  [ 0.552684362566E-01,   0.557908803218E-01,                       &
    &         0.557916264472E-01,   0.563198180271E-01,   0.557931137865E-01, &
    &         0.557916264472E-01,   0.278954401609E-01,   0.278954401609E-01, &
    &         0.278965568932E-01,   0.278965568932E-01,   0.276342181283E-01, &
    &         0.281599090135E-01,   0.276342181283E-01,   0.281599090135E-01, &
    &         0.274578450095E-01,   0.278060770311E-01,   0.283377645600E-01, &
    &         0.279850754323E-01,   0.334749758683E+00 ]

  ! Calculate strain rate acoustic explicitly inside this routine
  write(*,*) 'Calculate strain rate acoustic explicitly in utest'
  allocate( fEq(QQ) )
  allocate( fNeq(QQ) )
  rho = sum( PDF )
  vel = 0._rk
  do iDir = 1, QQ
    Vel(:) = Vel(:) + PDF(iDir) * dble(layout%fStencil%cxDir(:,iDir))
  end do
  vel = vel / rho0
  fEq(:) = getEquilibriumIncomp( rho, vel, layout, rho0 )
  fNeq   = ( PDF - fEq ) * convPrePost( omega )
  strain_exp(:) = secondMom( layout%fStencil%cxcx, fNeq, QQ )
  strain_exp(:) =  -1.5_rk * omega / rho0 * strain_exp

  ! print out PDF and fNeq values
  write(*, "(A)") "print out pdf and fNeq in utest routine"
  write(*, "(A4, 2A20)") "iDir", "PDF", "fNeq"
  do iDir = 1, QQ
    write(*, "(I4, 2F20.15)") iDir, PDF(iDir), fNeq(iDir)
  end do
  write(*,*) ''

  ! calculate strain rate acoustic by function call
  write(*,*) 'Calculate strain rate acoustic by function call'
  S = 0._rk
  strain_acs = 0._rk
  S = getShearRateTensor_acoustic( pdf, omega, layout, rho0  )
  strain_acs(1) = S(1,1)
  strain_acs(2) = S(2,2)
  strain_acs(3) = S(3,3)
  strain_acs(4) = S(1,2)
  strain_acs(5) = S(2,3)
  strain_acs(6) = S(1,3)

  ! calculate strain rate acoustic by function call
  write(*,*) 'Calculate strain rate diffusive by function call'
  S = 0._rk
  strain_dff = 0._rk
  S = getShearRateTensor_diffusive( pdf, omega, layout )
  strain_dff(1) = S(1,1)
  strain_dff(2) = S(2,2)
  strain_dff(3) = S(3,3)
  strain_dff(4) = S(1,2)
  strain_dff(5) = S(2,3)
  strain_dff(6) = S(1,3)

  ! Compare strain_ref and strain_res
  write(*,"(4A20)") "Compare strain_ref", "strain_exp", "strain_acs", "strain_dff"
  do iDir = 1, 6
    write(*,"(4F20.15)") strain_ref(iDir), strain_exp(iDir), strain_acs(iDir), strain_dff(iDir)
  end do
  write(*,*) ''

  ! Evaluate strain ref and strain res
  write(*,*) " Evaluating error"
  error_rel = abs( (strain_ref(4) - strain_acs(4)) / strain_ref(4) )
  write(*,"(A,F20.15)") "      error: ", error_rel
  write(*,"(A,F20.15)") "  tolerance: ", tolerance
  if ( error_rel > tolerance ) then
    error = .true.
    write(*, "(A)")    "Max Error exceeds tolerance!!!"
    write(*, "(3A20)") "strain_ref", "strain_res", "difference"
    write(*, "(3F20.15)") strain_ref(4), strain_acs(4), error_rel
  else
    write(*, "(A)") "Max Error within tolerance!"
    error = .false.
  end if

  call tem_finalize(general)
  if (.not. error) then
    write(*,'(A)') 'PASSED'
  else
    write(*,'(A)') 'FAILED'
  end if

end program mus_getStrainRate_acoustic_test
!******************************************************************************
