! See copyright notice in the COPYRIGHT file.
! This utest program tests the getStrainRate_diffusive routine.
! It first initializes array of PDF, then calculate strain rate from it.
! Finally compare calculated with reference values.
program mus_getStrainRate_diffusive_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_diffusive

  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(:)
  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_res(6), error_rel
  integer :: nDims = 3
  integer :: QQ
  integer :: iDir

  call tem_start( 'getStrainRate diffusive 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'
  ! call mus_init_stencil( stencil = layout%stencil, nStencils = 1 )
  if ( nDims == 2 ) then
    call mus_define_d2q9( layout = layout, nElems = 1 )
  else
    call mus_define_d3q19( layout = layout, nElems = 1 )
  end if
  ! call mus_finish_stencil( layout = layout )
  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 Diffusive
  write(*,*) 'Calculate strain rate by diffusive'
  S = getShearRateTensor_diffusive( pdf, omega, layout )
  strain_res(1) = S(1,1)
  strain_res(2) = S(2,2)
  strain_res(3) = S(3,3)
  strain_res(4) = S(1,2)
  strain_res(5) = S(2,3)
  strain_res(6) = S(1,3)

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

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

  ! Evaluate strain ref and strain res
  error_rel = abs( (strain_ref(4) - strain_res(4)) / strain_ref(4) )
  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_res(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_diffusive_test
!******************************************************************************
