! See copyright notice in the COPYRIGHT file.
! This unit test check equilibirium definition of d3q27
program mus_feq_d3q27_test
  use env_module,         only: rk, eps
  use tem_general_module, only: tem_start, tem_finalize, tem_general_type
  use tem_stencil_module, only: d3q27_cxDir
  use mus_d3q27_module,   only: feq_d3q27

  implicit none

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

  type( tem_general_type ) :: general
  real(kind=rk) :: rho, u,v,w

  call tem_start( 'D3Q27 f equilibirium test', 'utest', general )
  error = .false.
  tolerance = eps * 2500._rk
  write(*,"(A,F20.14)") 'tolerance = ', tolerance

  rho = 1.0_rk
  u   = 0.37_rk
  v   = -0.79_rk
  w   = 0.61_rk
  if ( check_conservation(rho, u,v,w) > tolerance ) then
    error = .true.
  end if

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

contains

  function check_conservation( rho, u, v, w ) result( max_error )
    real(kind=rk), intent(in) :: rho, u,v,w
    integer :: ii, i,j,k
    real(kind=rk) :: feq(27)
    real(kind=rk) :: post_u, post_v, post_w, post_rho
    real(kind=rk) :: error_rho, error_u, error_v, error_w, max_error

      ! write(*,"(3A3, A7)") 'i', 'j', 'k', 'feq'
    do ii = 1, 27
      i = d3q27_cxDir(1,ii)
      j = d3q27_cxDir(2,ii)
      k = d3q27_cxDir(3,ii)
      feq(ii) = feq_d3q27( i, j, k, rho, u, v, w )
      ! write(*,"(3I3, F7.3)") i, j, k, feq_d3q27( i, j, k, rho, u, v, w )
    end do

    post_rho = sum( feq )
    post_u = 0.0_rk
    post_v = 0.0_rk
    post_w = 0.0_rk
    do ii = 1, 27
      i = d3q27_cxDir(1,ii)
      j = d3q27_cxDir(2,ii)
      k = d3q27_cxDir(3,ii)
      post_u = post_u + feq(ii) * i
      post_v = post_v + feq(ii) * j
      post_w = post_w + feq(ii) * k
    end do

    error_rho = abs(post_rho - rho)
    error_u   = abs(post_u   - u)
    error_v   = abs(post_v   - v)
    error_w   = abs(post_w   - w)
    max_error = max( error_rho, error_u, error_v, error_w )

    write(*,"(3A10)") 'pre', 'post', 'error'
    write(*,"(3F10.5)") rho, post_rho, error_rho
    write(*,"(3F10.5)")   u, post_u,   error_u
    write(*,"(3F10.5)")   v, post_v,   error_v
    write(*,"(3F10.5)")   w, post_w,   error_w

  end function check_conservation

end program mus_feq_d3q27_test
