!> \~japanese
!! \brief 数値計算に伴う誤差を評価するためのサブルーチンを提供するモジュール.
!!
!! \par Field_IcGrid2D_error_norm_check モジュールについて
!!
!! \~english
!! \brief This module provides some subroutines to evaluate the error for numerical methods.
!!
!! \par About Field_IcGrid2D_error_norm_check module
!!
!! \~
!!
!! <br><br>
!! \b Copyright (C) GFD Dennou Club, 2011-2012. All rights reserved. <br>
!! \b license ?? <br>
!! \author Yuta Kawai
!!
module Field_IcGrid2D_error_norm_check

  ! モジュール引用 ; Use statements
  !

  !
  !
  !
#ifdef LIB_MPI
  use mpi
#endif

  ! 種類型パラメータ
  ! Kind type parameter
  !
  use dc_types, only: DP


  ! 正二十面体格子データの管理
  ! Manager of icosahedral grid data
  !
  use IcGrid2D_FVM_Manager, only: &
    & IcGrid2D_FVM, &
    & get_EffSize_Min, get_EffSize_Max, &
    & get_LocalProcRectInfo

  !
  !
  !
  use LocalProcRect_Manager, only: &
    & LocalProcRectInfo, &
    & get_InteriorRegion_IndexInfo
   
  ! 正二十面体格子上の物理場データの管理
  ! Manager of the physical field data on icosahedral grid
  !
  use Field_IcGrid2D_Manager, only: &
    & Field_IcGrid2D, &
    & get_icgrid, get_field2D_rank
  
  ! 全矩形領域における厳密解に対する数値解の誤差チェック
  !
  !
  use rcregions_norm_check, only: &
       & calc_norm1, calc_norm2, calc_norm3

  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

  ! 公開手続き
  ! Public procedures
  !
  public :: numerical_error_norm1, numerical_error_norm2, numerical_error_norminfinity

contains

!
!> \~japanese
!! \brief 厳密解から数値解に対する誤差ノルム 1 を求める.
!!
!! @param[in] true_field ある場の厳密解を管理する構造型 Field_IcGrid2D の変数.
!! @param[in] num_field ある場の数値解を管理する構造型 Field_IcGrid2D の変数.
!! @return normVal 誤差ノルム 1.
!!
!! \~english
!! \brief Calculates the error norm 1 for the numerical solution using exact solution.
!!
!! @param[in] true_field The variable of derived type IcGrid Field which manages the exaact solution of a given field.
!! @param[in] num_field The variable of derived type IcGrid Field which manages the numerical solution of a given field.
!! @return normVal Error norm 1.
!!
function numerical_error_norm1( &
  & true_field, num_field &  ! (in)
  & ) result(normVal)

  ! 宣言文 ; Declaration statements
  !
  type(Field_IcGrid2D), intent(in) :: true_field
  type(Field_IcGrid2D), intent(in) :: num_field
  real(DP) normVal
  
  ! 作業変数
  ! Work variables
  !
  integer :: rcID, i, j
  integer :: glEMin, glEMax
  type(IcGrid2D_FVM), pointer :: icgrid
  integer :: iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax
  type(LocalProcRectInfo), pointer :: lcRCInfo
  
  real(DP) :: local_tmp1, local_tmp2
  real(DP) :: global_tmp1, global_tmp2
  integer :: my_rank, ierr
  
  integer :: jMaxAveBound, iMinAveBound  
  
  ! 実行文 ; Executable statements
  !
  icgrid => get_icgrid(true_field)
  glEMin = get_EffSize_Min(icgrid)
  glEMax = get_EffSize_Max(icgrid)
  
  lcRCInfo => get_LocalProcRectInfo(icgrid)
  call get_InteriorRegion_IndexInfo( &
    & lcRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out)
 
  local_tmp1 = 0d0; local_tmp2 = 0d0;

  jMaxAveBound = jEMax 
  iMinAveBound = iEMin 
  if ( jEMax == glEMax ) jMaxAveBound = jEMax - 1
  if ( iEMin == glEMin ) iMinAveBound = iEMin + 1
  
  do rcID=rcIDMin, rcIDMax
    do j=jEMin, jMaxAveBound
      do i=iMinAveBound, iEMax
        local_tmp1 = local_tmp1 + abs_Val(num_field%val(:,i,j,rcID) - true_field%val(:,i,j,rcID))
        local_tmp2 = local_tmp2 + abs_Val(true_field%val(:,i,j,rcID))
      end do
    end do
  end do
  
  if ( rcIDMax == 1 .and. iEMin == glEMin .and. jEMin == glEMin) then
    local_tmp1 = local_tmp1 + abs_Val(num_field%val(:,glEMin,glEMin,1) - true_field%val(:,glEMin,glEMin,1))
    local_tmp2 = local_tmp2 + abs_Val(true_field%val(:,glEMin,glEMin,1))
  end if

  if ( rcIDMin == 10 .and. iEMin == glEMax .and. jEMin == glEMax) then
    local_tmp1 = local_tmp1 + abs_Val(num_field%val(:,glEMax,glEMax,10) - true_field%val(:,glEMax,glEMax,10))
    local_tmp2 = local_tmp2 + abs_Val(true_field%val(:,glEMax,glEMax,10))
  end if
  
  if ( lcRCInfo%mpi_flag ) then

#ifdef LIB_MPI
    
    call MPI_Comm_rank(MPI_COMM_WORLD, my_rank, ierr)
       
    call MPI_AllReduce( local_tmp1, global_tmp1, 1, MPI_REAL8, &
      &                 MPI_SUM, MPI_COMM_WORLD, ierr )

    call MPI_AllReduce( local_tmp2, global_tmp2, 1, MPI_REAL8, &
      &                 MPI_SUM, MPI_COMM_WORLD, ierr )
      
#endif

  else
    my_rank = 0
    global_tmp1 = local_tmp1
    global_tmp2 = local_tmp2
  end if
  
  normVal = global_tmp1 / global_tmp2
  
end function numerical_error_norm1

!
!> \~japanese
!! \brief 厳密解から数値解に対する誤差ノルム 2 を求める.
!!
!! @param[in] true_field ある場の厳密解を管理する構造型 Field_IcGrid2D の変数.
!! @param[in] num_field ある場の数値解を管理する構造型 Field_IcGrid2D の変数.
!! @return normVal 誤差ノルム 2.
!!
!! \~english
!! \brief Calculates the error norm 2 for the numerical solution using exact solution.
!!
!! @param[in] true_field The variable of derived type IcGrid Field which manages the exaact solution of a given field.
!! @param[in] num_field The variable of derived type IcGrid Field which manages the numerical solution of a given field.
!! @return normVal Error norm 2.
!!
function numerical_error_norm2( &
  & true_field, num_field &  ! (in)
  & ) result(normVal)

  ! 宣言文 ; Declaration statements
  !
  type(Field_IcGrid2D), intent(in) :: true_field
  type(Field_IcGrid2D), intent(in) :: num_field
  real(DP) normVal
  
  ! 作業変数
  ! Work variables
  !
  integer :: rcID, i, j
  integer :: glEMin, glEMax
  type(IcGrid2D_FVM), pointer :: icgrid
  integer :: iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax
  type(LocalProcRectInfo), pointer :: lcRCInfo
  
  real(DP) :: local_tmp1, local_tmp2
  real(DP) :: global_tmp1, global_tmp2
  integer :: my_rank, ierr

  integer :: jMaxAveBound, iMinAveBound
  
  ! 実行文 ; Executable statements
  !
  icgrid => get_icgrid(true_field)
  glEMin = get_EffSize_Min(icgrid)
  glEMax = get_EffSize_Max(icgrid)
  
  lcRCInfo => get_LocalProcRectInfo(icgrid)
  call get_InteriorRegion_IndexInfo( &
    & lcRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out)
 
  local_tmp1 = 0d0; local_tmp2 = 0d0;

  jMaxAveBound = jEMax 
  iMinAveBound = iEMin 
  if ( jEMax == glEMax ) jMaxAveBound = jEMax - 1
  if ( iEMin == glEMin ) iMinAveBound = iEMin + 1
  
  do rcID=rcIDMin, rcIDMax
    do j=jEMin, jMaxAveBound
      do i=iMinAveBound, iEMax

        local_tmp1 = local_tmp1 + abs_Val(num_field%val(:,i,j,rcID) - true_field%val(:,i,j,rcID))**2
        local_tmp2 = local_tmp2 + abs_Val(true_field%val(:,i,j,rcID))**2
      end do
    end do
  end do
  
  if ( rcIDMin == 1 .and. iEMin == 1 .and. jEMin == 1) then
    local_tmp1 = local_tmp1 + abs_Val(num_field%val(:,1,1,1) - true_field%val(:,1,1,1))**2
    local_tmp2 = local_tmp2 + abs_Val(true_field%val(:,1,1,1))**2
  end if

  if ( rcIDMin == 10 .and. iEMin == glEMax .and. jEMin == glEMax) then
    local_tmp1 = local_tmp1 + abs_Val(num_field%val(:,glEMax,glEMax,10) - true_field%val(:,glEMax,glEMax,10))**2
    local_tmp2 = local_tmp2 + abs_Val(true_field%val(:,glEMax,glEMax,10))**2
  end if
  
  if ( lcRCInfo%mpi_flag ) then

#ifdef LIB_MPI
    !
    call MPI_Comm_rank(MPI_COMM_WORLD, my_rank, ierr)
    
    call MPI_AllReduce( local_tmp1, global_tmp1, 1, MPI_REAL8, &
      &                 MPI_SUM, MPI_COMM_WORLD, ierr )

    call MPI_AllReduce( local_tmp2, global_tmp2, 1, MPI_REAL8, &
      &                 MPI_SUM, MPI_COMM_WORLD, ierr )
      
#endif

  else
    my_rank = 0
    global_tmp1 = local_tmp1
    global_tmp2 = local_tmp2
  end if
        
  normVal = sqrt(global_tmp1) / sqrt(global_tmp2)

end function numerical_error_norm2

!
!> \~japanese
!! \brief 厳密解から数値解に対する誤差ノルム 3 を求める.
!!
!! @param[in] true_field ある場の厳密解を管理する構造型 Field_IcGrid2D の変数.
!! @param[in] num_field ある場の数値解を管理する構造型 Field_IcGrid2D の変数.
!! @param[in] maxmin_info_flag 厳密解と数値解の局所的なずれが最も大きい位置および最も小さい位置の情報を出力するかのフラグ.
!! デフォルトは false である.
!! @return normVal 誤差ノルム 3.
!!
!! \~english
!! \brief Calculates the error norm 3 for the numerical solution using exact solution.
!!
!! @param[in] true_field The variable of derived type IcGrid Field which manages the exaact solution of a given field.
!! @param[in] num_field The variable of derived type IcGrid Field which manages the numerical solution of a given field.
!! @param[in] maxmin_info_flag The flag indicating whether one needs the information about the location
!! where the local difference with exact and numerical solution is maximum or minimum. The default flag value is false.
!! @return normVal Error norm 3.
!!
function numerical_error_norminfinity( &
  & true_field, num_field, maxmin_info_flag &  ! (in)
  & ) result(normVal)

  ! 宣言文 ; Declaration statements
  !
  type(Field_IcGrid2D), intent(in) :: true_field
  type(Field_IcGrid2D), intent(in) :: num_field
  logical, intent(in), optional :: maxmin_info_flag
  real(DP) normVal
  
  ! 作業変数
  ! Work variables
  !
  logical :: info_flag
  
  integer :: rcID, i, j
  integer :: glEMin, glEMax
  type(IcGrid2D_FVM), pointer :: icgrid
  integer :: iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax
  type(LocalProcRectInfo), pointer :: lcRCInfo
  
  real(DP) absfieldError, absfieldT 
  integer :: max_rcID, max_i, max_j
  real(DP) :: local_fieldTMax(2,1),  local_fieldErrorMax(2,1)
  real(DP) :: global_fieldTMax(2,1), global_fieldErrorMax(2,1)

  integer :: my_rank, ierr
  
  ! 実行文 ; Executable statements
  !
  icgrid => get_icgrid(true_field)
  glEMin = get_EffSize_Min(icgrid)
  glEMax = get_EffSize_Max(icgrid)
  
  lcRCInfo => get_LocalProcRectInfo(icgrid)
  call get_InteriorRegion_IndexInfo( &
    & lcRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out)
 
  local_fieldTMax = 0.0d0
  local_fieldErrorMax = 0.0d0
  max_rcID = 1; max_i = 1; max_j = 1

  do rcID=rcIDMin, rcIDMax
    do j=jEMin,jEMax
      do i=iEMin,iEMax
        
        absfieldError = abs_Val(true_field%val(:, i,j,rcID) - num_field%val(:, i,j,rcID))
        absfieldT = abs_Val(true_field%val(:, i,j,rcID))
        
        if( absfieldError >= local_fieldErrorMax(1,1) ) then
          local_fieldErrorMax(1,1) = absfieldError; 
          max_rcID=rcID; max_i=i; max_j=j;
        end if
        
        if( absfieldT > local_fieldTMax(1,1) ) then
          local_fieldTMax(1,1) = absfieldT
        end if
        
      end do
    end do
  end do

  if ( lcRCInfo%mpi_flag ) then
#ifdef LIB_MPI
    !
    call MPI_Comm_rank(MPI_COMM_WORLD, my_rank, ierr)
    
    local_fieldErrorMax(2,1) = my_rank
    local_fieldTMax(2,1) = my_rank
    
    call MPI_AllReduce( local_fieldErrorMax, global_fieldErrorMax, 1, MPI_2DOUBLE_PRECISION, &
      &                 MPI_MAXLOC, MPI_COMM_WORLD, ierr )
      
    call MPI_AllReduce( local_fieldTMax, global_fieldTMax, 1, MPI_2DOUBLE_PRECISION, &
      &                 MPI_MAXLOC, MPI_COMM_WORLD, ierr )

#endif

  else
    my_rank = 0
    global_fieldErrorMax = local_fieldErrorMax
    global_fieldTMax = local_fieldTMax
  end if

  if ( present(maxmin_info_flag) .and. maxmin_info_flag ) then
    info_flag = .true.
  else
    info_flag = .false.
  end if
  
  if (       present(maxmin_info_flag) .and. maxmin_info_flag  &
    &  .and. my_rank == int(global_fieldErrorMax(2,1))         &
    & ) then
    
     write(*,*) '= norm3 maximum error info:'
     write(*,'(a,3(i6))') '* location(rcID, i, j):', max_rcID, max_i, max_j
     write(*,*) '* Max Val:', global_fieldTMax(1,1)
     write(*,*) '* True value:', true_field%val(:,max_i,max_j, max_rcID)
     write(*,*) '* Calculated value:', num_field%val(:,max_i,max_j, max_rcID)
  end if
    
  normVal = global_fieldErrorMax(1,1) / global_fieldTMax(1,1)

end function numerical_error_norminfinity

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 非公開手続き
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

function abs_Val(errorVal) result(absVal)

  real(DP), intent(in)  :: errorVal(:) 
  real(DP)absVal

  integer :: valDim
  valDim = size(errorVal, 1)

  if( valDim == 1 )then
    absVal = dabs(errorVal(1))
  else if( valDim == 3 ) then
    absVal = dsqrt(dot_product(errorVal, errorVal))
  else
    absVal = 0.0d0
    write(*,*) '!!!!!!!!!!!!!!!!!!!!!!!'
  end if

end function abs_Val

end module Field_IcGrid2D_error_norm_check
