src/util/math/rcregions_norm_check.f90

Go to the documentation of this file.

module rcregions_norm_check
  
  ! モジュール引用 ; Use statement
  !
  use dc_types, only: DP  ! 倍精度実数型. Double precision. 
  
  ! 宣言文 ; Declaration statement
  !
  implicit none
  private

  ! 公開手続き
  ! Public procedure
  !
  public calc_norm1, calc_norm2, calc_norm3

  ! 公開変数
  ! Public variables
  !
  integer, parameter :: RC_REGIONS_NUM = 10

contains

!
!
!
function calc_norm1(true_field, num_field) result(normVal)
  ! 宣言文 ; Declaration statements
  !
  real(DP), intent(in) :: true_field(:,:,:,:)
  real(DP), intent(in) :: num_field(:,:,:,:)
  real(DP) normVal

  ! 作業変数
  ! Work varaibles
  !
  real(DP) fieldT_sum
  real(DP) fieldError_sum
  integer :: rcID, i,j
  integer :: EMax
  
  ! 実行文 ; Executable statements
  !

  EMax = size(true_field, 2)
  fieldT_sum = 0.0d0
  fieldError_sum = 0.0d0

  do rcID=1, RC_REGIONS_NUM
    do j=1,EMax-1
      do i=2,EMax
        fieldError_sum = fieldError_sum + abs_Val(true_field(rcID,i,j,:) - num_field(rcID, i,j, :))
        fieldT_sum = fieldT_sum + abs_Val(true_field(rcID,i,j,:))
        
      end do
    end do
  end do
  
  fieldError_sum = fieldError_sum + abs_Val(true_field(1,1,1,:) - num_field(1,1,1, :))
  fieldT_sum = fieldT_sum + abs_Val(true_field(1,1,1,:))
    
  fieldError_sum = fieldError_sum + abs_Val(true_field(6,EMax,EMax,:) - num_field(6,EMax,EMax, :))
  fieldT_sum = fieldT_sum + abs_Val(true_field(6,EMax,EMax,:))
  
  normVal = fieldError_sum / fieldT_sum

end function calc_norm1

!
function calc_norm2(true_field, num_field) result(normVal)
  ! 宣言文 ; Declaration statements
  !
  real(DP), intent(in) :: true_field(:,:,:,:)
  real(DP), intent(in) :: num_field(:,:,:,:)
  real(DP) normVal

  ! 作業変数
  ! Work varaibles
  !
  real(DP) fieldT_sum
  real(DP) fieldError_sum
  integer :: rcID, i,j
  integer :: EMax


  ! 実行文 ; Executable statements
  !

  EMax = size(true_field, 2)
  fieldT_sum = 0.0d0
  fieldError_sum = 0.0d0

  do rcID=1, RC_REGIONS_NUM
    do j=1,EMax-1
      do i=2,EMax
        fieldError_sum = fieldError_sum + abs_Val(true_field(rcID,i,j,:) - num_field(rcID, i,j, :))**2
        fieldT_sum = fieldT_sum + abs_Val(true_field(rcID,i,j,:))**2
      end do
    end do
  end do

  fieldError_sum = fieldError_sum + abs_Val(true_field(1,1,1,:) - num_field(1,1,1, :))**2
  fieldT_sum = fieldT_sum + abs_Val(true_field(1,1,1,:))**2
  
  fieldError_sum = fieldError_sum + abs_Val(true_field(6,EMax,EMax,:) - num_field(6,EMax,EMax, :))**2
  fieldT_sum = fieldT_sum + abs_Val(true_field(6,EMax,EMax,:))**2
  
  normVal = dsqrt(fieldError_sum) / dsqrt(fieldT_sum)
    
end function calc_norm2

!
function calc_norm3(true_field, num_field, output_info_flag ) result(normVal)
  ! 宣言文 ; Declaration statements
  !
  real(DP), intent(in) :: true_field(:,:,:,:)
  real(DP), intent(in) :: num_field(:,:,:,:)
  logical, intent(in), optional :: output_info_flag
  real(DP) normVal

  ! 作業変数
  ! Work varaibles
  !
  real(DP) :: fieldTMax
  real(DP) :: fieldErrorMax
  real(DP) absfieldError, absfieldT 
  integer :: rcID, i,j

  integer :: max_rcID, max_i, max_j
  integer :: EMax

  ! 実行文 ; Executable statements
  !

  EMax = size(true_field, 2)
  fieldTMax = 0.0d0
  fieldErrorMax = 0.0d0
  max_rcID = 1; max_i = 1; max_j = 1
 
  do rcID=1, RC_REGIONS_NUM
    do j=1,EMax
      do i=1,EMax
        
        absfieldError = abs_Val(true_field(rcID,i,j,:) - num_field(rcID,i,j,:))
        absfieldT = abs_Val(true_field(rcID,i,j,:))
        
        if( absfieldError >= fieldErrorMax ) then
          fieldErrorMax = absfieldError; 
          max_rcID=rcID; max_i=i; max_j=j;
        end if
        
        if( absfieldT > fieldTMax ) then
          fieldTMax = absfieldT
         ! write(*,*) rcID, i,j, absfieldError, absfieldT
        end if
        
      end do
    end do
  end do

  normVal = fieldErrorMax / fieldTMax

  if ( present(output_info_flag) .and. output_info_flag ) then
     write(*,*) 'norm3 maximum error info:', max_rcID, max_i, max_j 
     write(*,*) 'fieldTMax:', fieldTMax, ', true:', true_field(max_rcID,max_i,max_j,:), ', num:',  num_field(max_rcID,max_i,max_j,:)
  end if
end function  calc_norm3

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 rcregions_norm_check