scale_meshfield_analysis_numerror.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Data / Statistics / numerical error
!!
!! @par Description
!!           This module provides classes to evaluate numerical error for 1D, 2D, and 3D fields
!!
!! @author Yuta Kawai, Team SCALE
!!
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_meshfield_analysis_numerror
  !-----------------------------------------------------------------------------
  !
  !++ 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,                &
    elementbase1d, elementbase2d, elementbase3d
  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_mesh_base1d, only: meshbase1d
  use scale_mesh_base2d, only: meshbase2d
  use scale_mesh_base3d, only: meshbase3d
  use scale_localmesh_1d, only: localmesh1d
  use scale_localmesh_2d, only: localmesh2d
  use scale_localmesh_3d, only: localmesh3d
  use scale_meshfield_base, only: &
    meshfieldbase, meshfield1d, meshfield2d, meshfield3d
    
  use scale_meshfield_analysis_numerror_base, only: &
    meshfieldanalysisnumerrorbase,   &
    meshfieldanalysisnumerrorinfobase
 
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  ! 
 
  !- 1D
 
  !> Derived type for numerical error analysis of 1D field
  type, public, extends(meshfieldanalysisnumerrorbase) :: meshfieldanalysisnumerror1d
    class(meshbase1d), pointer :: mesh1d               !< Pointer to an object of 1D mesh
    procedure(set_data_lc_1d), pointer :: set_data_lc
  contains
    procedure :: init_1d => meshfield_analysis_numerror_1d_init
    generic :: init => init_1d
    procedure :: evaluate => meshfield_analysis_numerror_1d_evaluate
  end type meshfieldanalysisnumerror1d
  interface 
    subroutine set_data_lc_1d( this_, q, qexact, qexact_intrp, lcmesh, elem1D, intrp_epos, tsec )
      import rp
      import elementbase1d
      import localmesh1d
      import meshfieldanalysisnumerror1d
      class(meshfieldanalysisnumerror1d), intent(in) :: this_
      class(localmesh1d), intent(in) :: lcmesh
      class(elementbase1d), intent(in) :: elem1D
      real(RP), intent(out) :: q(elem1D%Np,lcmesh%Ne,this_%var_num)
      real(RP), intent(out) :: qexact(elem1D%Np,lcmesh%Ne,this_%var_num)
      real(RP), intent(out) :: qexact_intrp(this_%intrp_np,lcmesh%Ne,this_%var_num)
      real(RP), intent(in) :: intrp_epos(this_%intrp_np,this_%ndim)
      real(RP), intent(in) :: tsec
    end subroutine set_data_lc_1d
  end interface
 
  !- 2D
  !> Derived type for numerical error analysis of 2D field
  type, public, extends(meshfieldanalysisnumerrorbase) :: meshfieldanalysisnumerror2d
    class(meshbase2d), pointer :: mesh2d               !< Pointer to an object of 2D mesh
    procedure(set_data_lc_2d), pointer :: set_data_lc
  contains
    procedure :: init_2d => meshfield_analysis_numerror_2d_init
    generic :: init => init_2d
    procedure :: evaluate => meshfield_analysis_numerror_2d_evaluate
  end type meshfieldanalysisnumerror2d
  interface 
    subroutine set_data_lc_2d( this_, q, qexact, qexact_intrp, lcmesh, elem2D, intrp_epos, tsec )
      import rp
      import elementbase2d
      import localmesh2d
      import meshfieldanalysisnumerror2d
      class(meshfieldanalysisnumerror2d), intent(in) :: this_
      class(localmesh2d), intent(in) :: lcmesh
      class(elementbase2d), intent(in) :: elem2D
      real(RP), intent(out) :: q(elem2D%Np,lcmesh%Ne,this_%var_num)
      real(RP), intent(out) :: qexact(elem2D%Np,lcmesh%Ne,this_%var_num)
      real(RP), intent(out) :: qexact_intrp(this_%intrp_np,lcmesh%Ne,this_%var_num)
      real(RP), intent(in) :: intrp_epos(this_%intrp_np,this_%ndim)
      real(RP), intent(in) :: tsec
    end subroutine set_data_lc_2d
  end interface
 
  !- 3D
  !> Derived type for numerical error analysis of 3D field
  type, public, extends(meshfieldanalysisnumerrorbase) :: meshfieldanalysisnumerror3d
    class(meshbase3d), pointer :: mesh3d               !< Pointer to an object of 3D mesh
    procedure(set_data_lc_3d), pointer :: set_data_lc
  contains
    procedure :: init_3d => meshfield_analysis_numerror_3d_init
    generic :: init => init_3d
    procedure :: evaluate => meshfield_analysis_numerror_3d_evaluate
  end type meshfieldanalysisnumerror3d
  interface 
    subroutine set_data_lc_3d( this_, q, qexact, qexact_intrp, lcmesh, elem3D, intrp_epos, tsec )
      import rp
      import elementbase3d
      import localmesh3d
      import meshfieldanalysisnumerror3d
      class(meshfieldanalysisnumerror3d), intent(in) :: this_
      class(localmesh3d), intent(in) :: lcmesh
      class(elementbase3d), intent(in) :: elem3D
      real(RP), intent(out) :: q(elem3D%Np,lcmesh%Ne,this_%var_num)
      real(RP), intent(out) :: qexact(elem3D%Np,lcmesh%Ne,this_%var_num)
      real(RP), intent(out) :: qexact_intrp(this_%intrp_np,lcmesh%Ne,this_%var_num)
      real(RP), intent(in) :: intrp_epos(this_%intrp_np,this_%ndim)
      real(RP), intent(in) :: tsec
    end subroutine set_data_lc_3d
  end interface
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
 
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
 
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
 
contains
 
!-- 1D --
 
  !-----------------------------------------------------------------------------
  !> Initialize a object to evaluate numerical error for 1D field
!OCL SERIAL
  subroutine meshfield_analysis_numerror_1d_init( this,    &
    porder_error_check, log_fname_base, log_step_interval, &
    mesh, refElem1D,                                       &
    set_data_lc, numerror_analysis_info )
    implicit none
    class(meshfieldanalysisnumerror1d), intent(inout) :: this
    integer, intent(in) :: porder_error_check
    character(len=*), intent(in) :: log_fname_base
    integer, intent(in) :: log_step_interval
    class(meshbase1d), intent(in), target :: mesh
    type(lineelement), intent(in) :: refElem1D
    procedure(set_data_lc_1d) :: set_data_lc      
    class(meshfieldanalysisnumerrorinfobase), intent(in), target :: numerror_analysis_info
    !---------------------------------------------------------------------------
 
    call this%MeshFieldAnalysisNumerrorBase%Init( &
      porder_error_check, 1, refelem1d%Np, porder_error_check**1,     &
      log_fname_base, log_step_interval, mesh, numerror_analysis_info )
 
    this%IntrpMat(:,:) = refelem1d%GenIntGaussLegendreIntrpMat( &
      this%PolyOrderErrorCheck,                                 &  ! (in)
      this%intw_intrp, this%epos_intrp(:,1)                     )  ! (out)
 
    this%mesh1D => mesh
    this%set_data_lc => set_data_lc
    return
  end subroutine meshfield_analysis_numerror_1d_init
  
!OCL SERIAL
  subroutine meshfield_analysis_numerror_1d_evaluate( this, &
    tstep, tsec )
    implicit none
    class(meshfieldanalysisnumerror1d), intent(inout) :: this
    integer, intent(in) :: tstep
    real(RP), intent(in) :: tsec
    !---------------------------------------------------------------------------
 
    call this%Evaluate_base( tstep, tsec, this%mesh1D%dom_vol, evaluate_error_core, calc_covariance_core )
    return
  end subroutine meshfield_analysis_numerror_1d_evaluate
 
!-- 2D --
 
  !-----------------------------------------------------------------------------
  !> Initialize a object to evaluate numerical error for 2D field
!OCL SERIAL
  subroutine meshfield_analysis_numerror_2d_init( this, &
    porder_error_check, log_fname_base, log_step_interval, &
    mesh, refElem2D,                                       &
    set_data_lc, numerror_analysis_info )
    implicit none
    class(meshfieldanalysisnumerror2d), intent(inout) :: this
    integer, intent(in) :: porder_error_check
    character(len=*), intent(in) :: log_fname_base
    integer, intent(in) :: log_step_interval
    class(meshbase2d), intent(in), target :: mesh
    type(quadrilateralelement), intent(in) :: refElem2D
    class(meshfieldanalysisnumerrorinfobase), intent(in), target :: numerror_analysis_info
    procedure(set_data_lc_2d) :: set_data_lc
    !---------------------------------------------------------------------------
 
    call this%MeshFieldAnalysisNumerrorBase%Init( &
      porder_error_check, 2, refelem2d%Np, porder_error_check**2,     &
      log_fname_base, log_step_interval, mesh, numerror_analysis_info )
 
    this%IntrpMat(:,:) = refelem2d%GenIntGaussLegendreIntrpMat( &
      this%PolyOrderErrorCheck,                                   &  ! (in)
      this%intw_intrp, this%epos_intrp(:,1), this%epos_intrp(:,2) )  ! (out)
 
    this%mesh2D => mesh
    this%set_data_lc => set_data_lc
    return
  end subroutine meshfield_analysis_numerror_2d_init
 
!OCL SERIAL
  subroutine meshfield_analysis_numerror_2d_evaluate( this, &
    tstep, tsec )
    implicit none
    class(meshfieldanalysisnumerror2d), intent(inout) :: this
    integer, intent(in) :: tstep
    real(RP), intent(in) :: tsec
    !---------------------------------------------------------------------------
 
    call this%Evaluate_base( tstep, tsec, this%mesh2D%dom_vol, evaluate_error_core, calc_covariance_core )
    return
  end subroutine meshfield_analysis_numerror_2d_evaluate
 
!-- 3D --
 
  !-----------------------------------------------------------------------------
  !> Initialize a object to evaluate numerical error for 3D field
!OCL SERIAL
  subroutine meshfield_analysis_numerror_3d_init( this,    &
    porder_error_check, log_fname_base, log_step_interval, &
    mesh, refElem3D,                                       &
    set_data_lc, numerror_analysis_info )
    implicit none
    class(meshfieldanalysisnumerror3d), intent(inout) :: this
    integer, intent(in) :: porder_error_check
    character(len=*), intent(in) :: log_fname_base
    integer, intent(in) :: log_step_interval
    class(meshbase3d), intent(in), target :: mesh
    type(hexahedralelement), intent(in) :: refElem3D
    class(meshfieldanalysisnumerrorinfobase), intent(in), target :: numerror_analysis_info
    procedure(set_data_lc_3d) :: set_data_lc
    !---------------------------------------------------------------------------
 
    call this%MeshFieldAnalysisNumerrorBase%Init( &
      porder_error_check, 3, refelem3d%Np, porder_error_check**3,     &
      log_fname_base, log_step_interval, mesh, numerror_analysis_info )
 
 
    this%IntrpMat(:,:) = refelem3d%GenIntGaussLegendreIntrpMat( &
      this%PolyOrderErrorCheck,                                                         &  ! (in)
      this%intw_intrp, this%epos_intrp(:,1), this%epos_intrp(:,2), this%epos_intrp(:,3) )  ! (out)
 
    this%mesh3D => mesh
    this%set_data_lc => set_data_lc
    return
  end subroutine meshfield_analysis_numerror_3d_init
 
!OCL SERIAL
  subroutine meshfield_analysis_numerror_3d_evaluate( &
    this, tstep, tsec )
    implicit none
    class(meshfieldanalysisnumerror3d), intent(inout) :: this
    integer, intent(in) :: tstep
    real(RP), intent(in) :: tsec
    procedure(set_data_lc_3d) :: set_data_lc
    !---------------------------------------------------------------------------
 
    call this%Evaluate_base( tstep, tsec, this%mesh3D%dom_vol, evaluate_error_core, calc_covariance_core )
    return
  end subroutine meshfield_analysis_numerror_3d_evaluate
 
!--- private ---------------------
 
!OCL SERIAL
  subroutine evaluate_error_core( base, tsec, &
    num_error_l1_lc, num_error_l2_lc, num_error_linf_lc, &
    numsol_mean_lc, exactsol_mean_lc                     )
    implicit none
    class(meshfieldanalysisnumerrorbase), intent(in) :: base
    real(RP), intent(in) :: tsec
    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), allocatable :: q(:,:,:)
    real(RP), allocatable :: qexact(:,:,:)
    real(RP), allocatable :: qexact_intrp(:,:,:)
 
    integer :: lcdomid      
    integer :: Np, Ne
    class(elementbase), pointer :: elem
    class(localmeshbase), pointer :: lcmesh
    !---------------------------------------------------------------------------
 
    do lcdomid=1, base%mesh%LOCAL_MESH_NUM
      call prepare_lc_data( base, lcdomid, tsec,      & ! (in)
        q, qexact, qexact_intrp, lcmesh, elem, np, ne ) ! (out)
 
      call base%Evaluate_error_lc( &
        num_error_l1_lc, num_error_l2_lc, num_error_linf_lc, &
        numsol_mean_lc, exactsol_mean_lc,                    &
        q, qexact, qexact_intrp, lcmesh, elem                )
 
      deallocate( q, qexact, qexact_intrp )
    end do
    return
  end subroutine evaluate_error_core
 
!OCL SERIAL
  subroutine calc_covariance_core( base, tsec,&
    cov_numsol_numsol_lc, cov_numsol_exactsol_lc, cov_exactsol_exactsol_lc, &
    numsol_mean, exactsol_mean                      )
    implicit none
    class(meshfieldanalysisnumerrorbase), intent(in) :: base
    real(RP), intent(in) :: tsec
    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) :: numsol_mean(base%var_num) 
    real(RP), intent(in) :: exactsol_mean(base%var_num) 
 
    real(RP), allocatable :: q(:,:,:)
    real(RP), allocatable :: qexact(:,:,:)
    real(RP), allocatable :: qexact_intrp(:,:,:)
 
    integer :: lcdomid      
    integer :: Np, Ne
    class(localmeshbase), pointer :: lcmesh
    class(elementbase), pointer :: elem
    !---------------------------------------------------------------------------
 
    do lcdomid=1, base%mesh%LOCAL_MESH_NUM
      call prepare_lc_data( base, lcdomid, tsec,      & ! (in)
        q, qexact, qexact_intrp, lcmesh, elem, np, ne ) ! (out)
 
      call base%Evaluate_covariance_lc( &
        cov_numsol_numsol_lc, cov_numsol_exactsol_lc, cov_exactsol_exactsol_lc, &
        q, qexact, qexact_intrp, numsol_mean, exactsol_mean, lcmesh, elem       )
      
      deallocate( q, qexact, qexact_intrp )
    end do
    return
  end subroutine calc_covariance_core    
 
!OCL SERIAL
  subroutine prepare_lc_data( base, lcdomid, tsec, &
    q, qexact, qexact_intrp, lcmesh, elem, Np, Ne  )
    implicit none
    class(meshfieldanalysisnumerrorbase), intent(in) :: base
    integer, intent(in) :: lcdomid
    real(RP), intent(in) :: tsec
    real(RP), allocatable, intent(out) :: q(:,:,:)
    real(RP), allocatable, intent(out) :: qexact(:,:,:)
    real(RP), allocatable, intent(out) :: qexact_intrp(:,:,:)
    class(localmeshbase), pointer, intent(out) :: lcmesh
    class(elementbase),   pointer, intent(out) :: elem
    integer, intent(out) :: Np, Ne
 
    class(localmesh1d), pointer :: lcmesh1D
    class(localmesh2d), pointer :: lcmesh2D
    class(localmesh3d), pointer :: lcmesh3D
    !-------------------------------------------------------------------
 
    !- Decide Np/Ne and pointers
    select type(base)
    class is (meshfieldanalysisnumerror1d)
      lcmesh1d => base%mesh1D%lcmesh_list(lcdomid)
      np = lcmesh1d%refElem1D%Np;  ne = lcmesh1d%Ne
      lcmesh => lcmesh1d;         elem => lcmesh1d%refElem1D
    class is (meshfieldanalysisnumerror2d)
      lcmesh2d => base%mesh2D%lcmesh_list(lcdomid)
      np = lcmesh2d%refElem2D%Np;  ne = lcmesh2d%Ne
      lcmesh => lcmesh2d;         elem => lcmesh2d%refElem2D
    class is (meshfieldanalysisnumerror3d)
      lcmesh3d => base%mesh3D%lcmesh_list(lcdomid)
      np = lcmesh3d%refElem3D%Np;  ne = lcmesh3d%Ne
      lcmesh => lcmesh3d;         elem => lcmesh3d%refElem3D
    end select
 
    allocate( q(np,ne,base%var_num) )
    allocate( qexact(np,ne,base%var_num) )
    allocate( qexact_intrp(base%intrp_np,ne,base%var_num) )
 
    !- Fill arrays via set_data_lc
    select type(base)
    class is (meshfieldanalysisnumerror1d)
      call base%set_data_lc( q, qexact, qexact_intrp,       &
        lcmesh1d, lcmesh1d%refElem1D, base%epos_intrp, tsec )
    class is (meshfieldanalysisnumerror2d)
      call base%set_data_lc( q, qexact, qexact_intrp,       &
        lcmesh2d, lcmesh2d%refElem2D, base%epos_intrp, tsec )
    class is (meshfieldanalysisnumerror3d)
      call base%set_data_lc( q, qexact, qexact_intrp,       &
        lcmesh3d, lcmesh3d%refElem3D, base%epos_intrp, tsec )
    end select
    return
  end subroutine prepare_lc_data
 
end module scale_meshfield_analysis_numerror