! See copyright notice in the COPYRIGHT file.
! This utest program tests if getStrainRate acoustic diffusive produce the same results
program mus_strainRate_test
  use env_module,         only: rk, eps
  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
  use mus_derivedQuantities_module2, only: getShearRateTensor_acoustic, &
    &                                      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) :: Sa(3,3), Sd(3,3)
  real(kind=rk) :: omega, rho0
  integer :: iVal, jVal, iOmega

  call tem_start( 'strain rate utest', 'utest', general )
  error = .false.
  tolerance = eps * 2500._rk
  write(*,*) 'tolerance = ', tolerance

  ! Initialize stencil
  write(*,*) "Initialize stencil D3Q19"
  call mus_define_d3q19( layout = layout, nElems = 1 )

  rho0  = 1.0_rk

  ! Initialize pdf
  write(*,*) "Initialize pdf"
  allocate( Pdf( layout%fStencil%QQ ) )
  do iVal = 1, layout%fStencil%QQ
    PDF( iVal ) = real(iVal,rk) * 0.1_rk
  enddo

  ! loop over omega range: 1.05 -> 1.95
  do iOmega = 1, 19

    omega = 1.0_rk + 0.05_rk * iOmega
    write(*,"(A,F5.2)") 'omega = ', omega

    ! calculate strain rate
    write(*,*) "Calculate diffusive strain rate"
    Sd = getShearRateTensor_diffusive( PDF, omega, layout )
    write(*,*) "Calculate acoustic  strain rate"
    Sa = getShearRateTensor_acoustic( PDF, omega, layout, rho0 )

    ! check Sd and Sa
    write(*,*) "Calculate strain difference between acoustic and diffusive"
    max_error = maxval( abs( Sd - Sa ) )
    write(*, "(A, ES20.10)") 'max error = ', max_error
    if ( max_error > tolerance ) then
      error = .true.
      do jVal = 1, 3
        do iVal = 1, 3
          write(*, "(2A5, 2A20)") 'iVal', 'jVal', 'Sa', 'Sd'
          write(*, "(2I5, 2ES20.10)") iVal, jVal, Sa(iVal, jVal), Sd(iVal, jVal)
        end do
      end do
    end if
    write(*,"(A)") ' '

  end do ! iOmega

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

end program mus_strainRate_test
!******************************************************************************
