scale_meshfield_analysis_numerror_base.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Data / Statistics / numerical error
!!
!! @par Description
!!           This module provides a base class useful for evaluating numerical errors
!!
!! @author Yuta Kawai, Team SCALE
!!
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_meshfield_analysis_numerror_base
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prc, only: &
    prc_ismaster, prc_abort
  use scale_time_manager, only: &
    time_nstep
 
  use scale_element_base, only: elementbase
  use scale_element_line, only: lineelement
  use scale_element_quadrilateral, only: quadrilateralelement
  use scale_element_hexahedral, only: hexahedralelement
  use scale_localmesh_base, only: localmeshbase
  use scale_localmeshfield_base, only: localmeshfieldbase
  use scale_mesh_base, only: meshbase
  use scale_meshfield_base, only: meshfieldbase
  
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  ! 
 
  !> Derived type saving information for numerical error analysis
  type, abstract, public :: meshfieldanalysisnumerrorinfobase
  end type meshfieldanalysisnumerrorinfobase
 
  !> Base type for numerical error analysis of mesh field
  type, public :: meshfieldanalysisnumerrorbase
    integer :: var_num             !< Number of variables to be analyzed
    integer :: log_fid             !< File ID for log output
    integer :: log_step_interval   !< Interval of time step for log output
    integer :: log_rstep           !< Remaining step for log output
    logical :: output_error_first  !< Flag for outputting error at the first time
 
    class(meshbase), pointer :: mesh !< Pointer to an object of type MeshBase
 
    integer :: polyordererrorcheck         !< Polynomial order when evaluating numerical errors
    integer :: intrp_np                    !< Number of interpolation points
    integer :: ndim                        !< Number of dimensions
    real(rp), allocatable :: intrpmat(:,:)   !< Interpolation matrix for numerical error evaluation
    real(rp), allocatable :: intw_intrp(:)   !< Weights of integration for interpolation points
    real(rp), allocatable :: epos_intrp(:,:) !< Positions of interpolation points
 
    class(meshfieldanalysisnumerrorinfobase), pointer :: info !< Pointer to information for numerical error analysis
  contains
    procedure :: init_base => meshfield_analysis_numerror_base_init 
    generic :: init => init_base
    procedure :: final => meshfield_analysis_numerror_base_final
    procedure :: regist => meshfield_analysis_numerror_base_regist
    procedure :: evaluate_base => meshfield_analysis_numerror_base_evaluate
    !
    procedure :: evaluate_error_lc => meshfield_analysis_numerror_base_evaluate_error_lc
    procedure :: evaluate_covariance_lc => meshfield_analysis_numerror_base_evaluate_covariance_lc
  end type meshfieldanalysisnumerrorbase
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedures & parameters & variables
  !
  abstract interface 
    subroutine evaluate_error_interface( this_, tsec, &
      num_error_l1_lc_, num_error_l2_lc_, num_error_linf_lc_, &
      numsol_mean, exactsol_mean )
      import rp
      import meshfieldanalysisnumerrorbase
      class(meshfieldanalysisnumerrorbase), intent(in) :: this_
      real(RP), intent(in) :: tsec
      real(RP), intent(inout) :: num_error_l1_lc_(this_%var_num)
      real(RP), intent(inout) :: num_error_l2_lc_(this_%var_num)
      real(RP), intent(inout) :: num_error_linf_lc_(this_%var_num)    
      real(RP), intent(inout) :: numsol_mean(this_%var_num) 
      real(RP), intent(inout) :: exactsol_mean(this_%var_num) 
    end subroutine evaluate_error_interface
 
    subroutine calc_covariance_interface( this_, tsec, &
      cov_numsol_numsol_lc, cov_numsol_exactsol_lc, cov_exactsol_exactsol_lc, &
      numsol_mean, exactsol_mean                      )
      import rp
      import meshfieldanalysisnumerrorbase
      class(meshfieldanalysisnumerrorbase), intent(in) :: this_
      real(RP), intent(in) :: tsec
      real(RP), intent(inout) :: cov_numsol_numsol_lc(this_%var_num)
      real(RP), intent(inout) :: cov_numsol_exactsol_lc(this_%var_num)
      real(RP), intent(inout) :: cov_exactsol_exactsol_lc(this_%var_num)
      real(RP), intent(in) :: numsol_mean(this_%var_num) 
      real(RP), intent(in) :: exactsol_mean(this_%var_num)   
    end subroutine calc_covariance_interface
  end interface
 
  private :: numerror_do_step
 
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
 
contains
 
  !-----------------------------------------------------------------------------
  !> Initialize an object for numerical error analysis
!OCL SERIAL
  subroutine meshfield_analysis_numerror_base_init( this,                      &
    porder_error_check, ndim, np, intrp_np, log_fname_base, log_step_interval, &
    mesh, numerror_analysis_info )
    implicit none
    class(meshfieldanalysisnumerrorbase), intent(inout) :: this
    integer, intent(in) :: porder_error_check      !< Polynomial order when evaluating numerical errors
    integer, intent(in) :: ndim                    !< Number of dimensions
    integer, intent(in) :: np                      !< Number of nodes per element
    integer, intent(in) :: intrp_np                !< Number of interpolation points for numerical error evaluation
    character(len=*), intent(in) :: log_fname_base !< Base name of log file for numerical error analysis
    integer, intent(in) :: log_step_interval       !< Interval of time step for log output
    class(meshbase), intent(in), target :: mesh !< Mesh for numerical error analysis
    class(meshfieldanalysisnumerrorinfobase), intent(in), target :: numerror_analysis_info !< Information for numerical error analysis
 
    character(len=H_MID) :: fname
    integer :: ierr
    !---------------------------------------------------------------------------
 
    this%var_num = 0
 
    this%log_step_interval = log_step_interval
    this%log_rstep         = log_step_interval
 
    this%output_error_first = .true. 
 
    !--
    if ( prc_ismaster ) then
      this%log_fid = io_get_available_fid()
      fname = trim(log_fname_base)//".peall"
      open( unit   = this%log_fid,                 &
        file   = fname,                            &
        form   = 'formatted',                      &
        iostat = ierr                              )
      if ( ierr /= 0 ) then
        log_error('MeshField_NumErrorAnalysis_Init',*) 'File open error! :', trim(fname)
        call prc_abort
      endif 
    end if
 
    !--
    this%PolyOrderErrorCheck = porder_error_check
    this%intrp_np = intrp_np
    this%ndim = ndim
    allocate( this%IntrpMat(this%intrp_np,np) )
    allocate( this%intw_intrp(this%intrp_np) )
    allocate( this%epos_intrp(this%intrp_np,this%ndim) )
 
    this%info => numerror_analysis_info
    this%mesh => mesh
    return
  end subroutine meshfield_analysis_numerror_base_init
 
  !-----------------------------------------------------------------------------
  !> Finalize an object for numerical error analysis
!OCL SERIAL
  subroutine meshfield_analysis_numerror_base_final( this )
    implicit none
    class(meshfieldanalysisnumerrorbase), intent(inout) :: this
    !---------------------------------------------------------------------------
    
    deallocate( this%IntrpMat, this%intw_intrp, this%epos_intrp )
 
    close( this%log_fid )
    this%log_fid = -1
 
    this%info => null()
    this%mesh => null()
    return
  end subroutine meshfield_analysis_numerror_base_final
 
  !> Register a variable for numerical error analysis
!OCL SERIAL
  subroutine meshfield_analysis_numerror_base_regist( this, varname, unit, varid )
    implicit none
    class(meshfieldanalysisnumerrorbase), intent(inout) :: this
    character(len=*), intent(in) :: varname !< Name of variable to be analyzed
    character(len=*), intent(in) :: unit    !< Unit of variable to be analyzed
    integer, intent(out) :: varid           !< ID of variable to be analyzed
    !---------------------------------------------------------------------------
 
    this%var_num = this%var_num + 1
    varid = this%var_num
 
    if ( prc_ismaster ) then
      write(this%log_fid,'(A25)',advance='no') 'L1_error   ('//trim(varname)//')'
      write(this%log_fid,'(A25)',advance='no') 'L2_error   ('//trim(varname)//')'
      write(this%log_fid,'(A25)',advance='no') 'Linf_error ('//trim(varname)//')'
      write(this%log_fid,'(A25)',advance='no') 'Ediss      ('//trim(varname)//')'
      write(this%log_fid,'(A25)',advance='no') 'Edisp      ('//trim(varname)//')'
    end if
    return
  end subroutine meshfield_analysis_numerror_base_regist
 
  !> Evaluate numerical errors
!OCL SERIAL
  subroutine meshfield_analysis_numerror_base_evaluate( &
    this, tstep, tsec, dom_vol, evaluate_error, calc_covariance )
    use mpi
    use scale_prc, only: &
      prc_local_comm_world
    implicit none
    class(meshfieldanalysisnumerrorbase), intent(inout) :: this
    integer, intent(in) :: tstep
    real(RP) :: tsec
    real(RP), intent(in) :: dom_vol
    procedure(evaluate_error_interface) :: evaluate_error
    procedure(calc_covariance_interface) :: calc_covariance
 
    logical :: do_check_numerror
 
    integer, parameter :: SUM_BUF_L1_ID = 1
    integer, parameter :: SUM_BUF_L2_ID = 2
    integer, parameter :: SUM_BUF_MEAN_NUMSOL_ID   = 3
    integer, parameter :: SUM_BUF_MEAN_EXACTSOL_ID = 4
    integer, parameter :: SUM_BUF_NUM   = 4
 
    real(RP) :: num_error_linf_lc(this%var_num)
    real(RP) :: sum_buf_lc(this%var_num,SUM_BUF_NUM)
    real(RP) :: num_error_linf(this%var_num)
    real(RP) :: sum_buf(this%var_num,SUM_BUF_NUM)
 
    real(RP) :: covariance_lc(this%var_num,3)
    real(RP) :: covariance(this%var_num,3)
    real(RP) :: Ediss(this%var_num), Edisp(this%var_num)
 
    integer :: iv
    integer :: ierr
    !---------------------------------------------------------------------------
 
    if ( this%var_num == 0 ) return
    
    if ( this%output_error_first ) then
      if ( prc_ismaster ) write(this%log_fid,*)
      this%output_error_first = .false.
    end if
 
    do_check_numerror = numerror_do_step( this, tstep )
    if ( .not. do_check_numerror  ) then
      !if ( PRC_ismaster .and. TIME_NOWSTEP == TIME_NSTEP ) close( NUMERROR_LOG_FID )
      return
    end if
 
    !---
    num_error_linf_lc(:) = 0.0_rp
    sum_buf_lc(:,:)      = 0.0_rp
 
    call evaluate_error( this, tsec, &
      sum_buf_lc(:,sum_buf_l1_id), sum_buf_lc(:,sum_buf_l2_id), num_error_linf_lc(:), &
      sum_buf_lc(:,sum_buf_mean_numsol_id), sum_buf_lc(:,sum_buf_mean_exactsol_id)    )
 
    call mpi_allreduce( sum_buf_lc(:,:), sum_buf(:,:), &
      this%var_num * sum_buf_num,                      &
      mpi_double_precision,   &
      mpi_sum,                &
      prc_local_comm_world,   &
      ierr                    )
    call mpi_allreduce( num_error_linf_lc(:), num_error_linf(:), &
      this%var_num,           &
      mpi_double_precision,   &
      mpi_max,                &
      prc_local_comm_world,   &
      ierr                    )
 
    sum_buf(:,sum_buf_mean_numsol_id) = sum_buf(:,sum_buf_mean_numsol_id) / dom_vol
    sum_buf(:,sum_buf_mean_exactsol_id) = sum_buf(:,sum_buf_mean_exactsol_id) / dom_vol
  
    !--
    covariance_lc(:,:) = 0.0_rp
 
    call calc_covariance( this, tsec, &
      covariance_lc(:,1), covariance_lc(:,2), covariance_lc(:,3),            &
      sum_buf(:,sum_buf_mean_numsol_id), sum_buf(:,sum_buf_mean_exactsol_id) )
    
    call mpi_allreduce( covariance_lc(:,:), covariance(:,:), &
      this%var_num * 3,                                      &
      mpi_double_precision,   &
      mpi_sum,                &
      prc_local_comm_world,   &
      ierr                    )
 
    if ( prc_ismaster ) then
      
      ediss(:) = ( sqrt(covariance(:,1) / dom_vol) -  sqrt(covariance(:,3) / dom_vol) )**2      &
               + ( sum_buf(:,sum_buf_mean_numsol_id) - sum_buf(:,sum_buf_mean_exactsol_id) )**2
      
      edisp(:) = 2.0_rp * (   sqrt( covariance(:,1)  *  covariance(:,3) ) &
                            - covariance(:,2)                             ) / dom_vol
 
      write(this%log_fid,'(A,ES18.8)',advance='no') 'tsec=', tsec
      do iv=1, this%var_num
       write(this%log_fid,'(A,ES18.8)',advance='no') ' ', sum_buf(iv,sum_buf_l1_id) / dom_vol
       write(this%log_fid,'(A,ES18.8)',advance='no') ' ', sqrt(sum_buf(iv,sum_buf_l2_id) / dom_vol)
       write(this%log_fid,'(A,ES18.8)',advance='no') ' ', num_error_linf(iv)
       write(this%log_fid,'(A,ES18.8)',advance='no') ' ', ediss(iv)
       write(this%log_fid,'(A,ES18.8)',advance='no') ' ', edisp(iv)
      end do
      write(this%log_fid,*)
 
!      if ( tstep > TIME_NSTEP ) close( NUMERROR_LOG_FID )
    end if
 
    return
  end subroutine meshfield_analysis_numerror_base_evaluate
 
  !> Evaluate numerical errors at local mesh level
!OCL SERIAL
  subroutine meshfield_analysis_numerror_base_evaluate_error_lc( base, &
    num_error_l1_lc, num_error_l2_lc, num_error_linf_lc, &
    numsol_mean_lc, exactsol_mean_lc,                    &
    q, qexact, qexact_intrp, lcmesh, elem                )
    implicit none
    class(meshfieldanalysisnumerrorbase), intent(in) :: base
    class(localmeshbase), intent(in) :: lcmesh
    class(elementbase), intent(in) :: elem
    real(RP), intent(inout) :: num_error_l1_lc(base%var_num)
    real(RP), intent(inout) :: num_error_l2_lc(base%var_num)
    real(RP), intent(inout) :: num_error_linf_lc(base%var_num)
    real(RP), intent(inout) :: numsol_mean_lc(base%var_num) 
    real(RP), intent(inout) :: exactsol_mean_lc(base%var_num) 
    real(RP), intent(in) :: q(elem%Np,lcmesh%Ne,base%var_num)
    real(RP), intent(in) :: qexact(elem%Np,lcmesh%Ne,base%var_num)
    real(RP), intent(in) :: qexact_intrp(base%intrp_np,lcmesh%Ne,base%var_num)
    
    integer :: iv
    integer :: ke
    
    real(RP) :: linf_max_tmp(lcmesh%Ne,base%var_num)
    real(RP) :: q_intrp(base%intrp_np)
    real(RP) :: JGsqrtxIntw_intrp(base%intrp_np)
    !---------------------------------------------------------------------------
 
    do iv=1, base%var_num
      !$omp parallel do private(ke, q_intrp, JGsqrtxIntw_intrp) reduction(+: num_error_l1_lc, num_error_l2_lc, numsol_mean_lc, exactsol_mean_lc )
      do ke=lcmesh%NeS, lcmesh%NeE
        q_intrp(:) = matmul( base%IntrpMat, q(:,ke,iv) )
        jgsqrtxintw_intrp(:) = matmul( base%IntrpMat, lcmesh%J(:,ke) * lcmesh%Gsqrt(:,ke) )
        jgsqrtxintw_intrp(:) = jgsqrtxintw_intrp(:) * base%intw_intrp(:)
 
        num_error_l1_lc(iv) = num_error_l1_lc(iv) + sum( jgsqrtxintw_intrp(:) * abs( q_intrp(:) - qexact_intrp(:,ke,iv) ) )
        num_error_l2_lc(iv) = num_error_l2_lc(iv) + sum( jgsqrtxintw_intrp(:) * ( q_intrp(:) - qexact_intrp(:,ke,iv) )**2 )
        linf_max_tmp(ke,iv) = maxval(abs(q(:,ke,iv) - qexact(:,ke,iv)))
 
        numsol_mean_lc(iv) = numsol_mean_lc(iv) + sum( jgsqrtxintw_intrp(:) * q_intrp(:) )
        exactsol_mean_lc(iv) = exactsol_mean_lc(iv) + sum( jgsqrtxintw_intrp(:) * qexact_intrp(:,ke,iv) )
      end do
      
      num_error_linf_lc(iv) = max( num_error_linf_lc(iv), maxval( linf_max_tmp(:,iv) ) )
    end do
 
    return
  end subroutine meshfield_analysis_numerror_base_evaluate_error_lc
 
  !> Evaluate covariance of numerical and exact solutions at local mesh level
!OCL SERIAL
  subroutine meshfield_analysis_numerror_base_evaluate_covariance_lc( base, &
    cov_numsol_numsol_lc, cov_numsol_exactsol_lc, cov_exactsol_exactsol_lc, &
    q, qexact, qexact_intrp, numsol_mean, exactsol_mean,                    &
    lcmesh, elem                                                            )
  
    implicit none
    class(meshfieldanalysisnumerrorbase), intent(in) :: base
    class(localmeshbase), intent(in) :: lcmesh
    class(elementbase), intent(in) :: elem
    real(RP), intent(inout) :: cov_numsol_numsol_lc(base%var_num)
    real(RP), intent(inout) :: cov_numsol_exactsol_lc(base%var_num)
    real(RP), intent(inout) :: cov_exactsol_exactsol_lc(base%var_num)
    real(RP), intent(in) :: q(elem%Np,lcmesh%Ne,base%var_num)
    real(RP), intent(in) :: qexact(elem%Np,lcmesh%Ne,base%var_num)
    real(RP), intent(in) :: qexact_intrp(base%intrp_np,lcmesh%Ne,base%var_num)
    real(RP), intent(in) :: numsol_mean(base%var_num) 
    real(RP), intent(in) :: exactsol_mean(base%var_num) 
 
    integer :: iv
    integer :: ke
    
    real(RP) :: dq_intrp(base%intrp_np)
    real(RP) :: JGsqrtxIntw_intrp(base%intrp_np)
    real(RP) :: dq_exact_intrp(base%intrp_np)
    !---------------------------------------------------------------------------
 
    do iv=1, base%var_num
      !$omp parallel do private(ke, dq_intrp, dq_exact_intrp, JGsqrtxIntw_intrp) reduction(+: cov_numsol_numsol_lc, cov_numsol_exactsol_lc, cov_exactsol_exactsol_lc )
      do ke=lcmesh%NeS, lcmesh%NeE
        dq_intrp(:) = matmul( base%IntrpMat, q(:,ke,iv) ) - numsol_mean(iv)
        dq_exact_intrp(:) = qexact_intrp(:,ke,iv)- exactsol_mean(iv)
        jgsqrtxintw_intrp(:) = matmul( base%IntrpMat, lcmesh%J(:,ke) * lcmesh%Gsqrt(:,ke) )
        jgsqrtxintw_intrp(:) = jgsqrtxintw_intrp(:) * base%intw_intrp(:)
 
        cov_numsol_numsol_lc(iv)   = cov_numsol_numsol_lc(iv) + sum( jgsqrtxintw_intrp(:) * dq_intrp(:) * dq_intrp(:) )
        cov_numsol_exactsol_lc(iv) = cov_numsol_exactsol_lc(iv) + sum( jgsqrtxintw_intrp(:) * dq_intrp(:) * dq_exact_intrp(:) )
        cov_exactsol_exactsol_lc(iv) = cov_exactsol_exactsol_lc(iv) + sum( jgsqrtxintw_intrp(:) * dq_exact_intrp(:) * dq_exact_intrp(:) )
      end do
    end do
 
    return
  end subroutine meshfield_analysis_numerror_base_evaluate_covariance_lc
 
 
!--- private
 
  !> Get a flag whether to evaluate numerical errors at the current time step
!OCL SERIAL
  function numerror_do_step( this, step ) result(do_flag)
    implicit none
 
    type(meshfieldanalysisnumerrorbase), intent(inout) :: this
    integer, intent(in) :: step
    logical :: do_flag    
    !------------------------------------------------------------------------
 
    if ( step == 1 .or. time_nstep < step ) then
      do_flag = .true.
      return
    end if
 
    this%log_rstep = this%log_rstep - 1
    if ( this%log_rstep == 0 ) then
      do_flag = .true.
      this%log_rstep  = this%log_step_interval
    else
      do_flag = .false. 
    end if
 
    return
  end function numerror_do_step
 
end module scale_meshfield_analysis_numerror_base