scale_mesh_hierarchy_2d.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Mesh / 2D domain
!!
!! @par Description
!!      Manage mesh hierarchy of 2D domain for element-based methods
!!
!! @author Yuta Kawai, Team SCALE
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_mesh_hierarchy_2d
 
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  !
  use scale_precision
  use scale_io
  use scale_prc, only: prc_abort
 
  use scale_element_base, only: elementbase2d
  use scale_element_quadrilateral, only: quadrilateralelement
 
  use scale_localmesh_2d, only: localmesh2d
  use scale_mesh_base2d, only: meshbase2d
  use scale_mesh_rectdom2d, only: meshrectdom2d
 
  use scale_mesh_hierarchy_base, only: &
    meshhierarchybase_init, &
    meshhierarchybase_final, &
    meshhierarchybase,       &
    meshhierarchylevelbase,  &
    meshhierarchylocalmgdatabase, &
    mesh_hierarchy_pmg_finest_level, &
    mesh_hierarchy_hmg_finest_level, &
    mesh_hierarchy_type_pmg,         &
    mesh_hierarchy_type_hmg
 
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  ! 
 
  !> Derived type to save a pointer to MeshBase2D
  type :: meshptr2d
    class(MeshBase2D), pointer :: ptr => null()
  end type meshptr2d
 
  !> Derived type to manage 2D local mesh data for multigrid
  type, extends(meshhierarchylocalmgdatabase), public :: meshhierarchylocalmgdata2d
  contains
    procedure :: init => meshhierarchylocalmgdata2d_init
    procedure :: final => meshhierarchylocalmgdata2d_final
  end type meshhierarchylocalmgdata2d
 
  !> Derived type to represent mesh hierarchy level in 2D domain
  type, extends(meshhierarchylevelbase), public :: meshhierarchylevel2d
    type(meshptr2d), pointer :: fine_mesh => null()
    type(meshptr2d), pointer :: coarse_mesh => null()
    type(meshhierarchylocalmgdata2d), allocatable :: mg_local(:)
 
    real(rp), allocatable :: pmat1d_f2c(:,:)
    real(rp), allocatable :: pmat1d_c2f(:,:)
  contains
    procedure :: init => meshhierarchylevel2d_init
    procedure :: final => meshhierarchylevel2d_final
  end type meshhierarchylevel2d
 
  !> Derived type for mesh hierarchy in 2D domain
  type, extends(meshhierarchybase), public :: meshhierarchy2d
    type(meshhierarchylevel2d), allocatable :: p_level(:)
    type(meshhierarchylevel2d), allocatable :: h_level(:)
 
    type(meshptr2d), allocatable :: p_mesh_list(:)
    type(meshptr2d), allocatable :: h_mesh_list(:)
 
    type(quadrilateralelement), allocatable :: elem2d_list(:)
  contains
    procedure :: init => meshhierarchy2d_init
    procedure :: final => meshhierarchy2d_final
  end type meshhierarchy2d
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
contains
  !> Initialize an object for mesh hierarchy in 2D domain
!OCL SERIAL
  subroutine meshhierarchy2d_init( this, &
    parent_mesh,                       &
    porder_list, p_LEVEL_NUM,          &
    NeGX_list, NeGY_list, h_LEVEL_NUM  )
    implicit none
    class(meshhierarchy2d), intent(inout), target :: this
    class(meshbase2d), intent(in), target :: parent_mesh  !< Pointer to the finest mesh
    integer, intent(in) :: p_LEVEL_NUM                    !< Number of p-mesh levels
    integer, intent(in) :: porder_list(p_LEVEL_NUM)       !< Polynomial order list for p-mesh levels
    integer, intent(in) :: h_LEVEL_NUM                    !< Number of h-mesh levels
    integer, intent(in) :: NeGX_list(h_LEVEL_NUM)         !< Number of elements in x-direction for h-mesh levels
    integer, intent(in) :: NeGY_list(h_LEVEL_NUM)         !< Number of elements in y-direction for h-mesh levels
    
    integer :: poly_lev
    integer :: h_lev
 
    integer :: ldom_id
 
    type(meshhierarchylevel2d), pointer :: level_ptr
    class(meshbase2d), pointer :: mesh2D_ptr
 
    class(meshptr2d), pointer :: fine_mesh_ptr
    class(meshptr2d), pointer :: coarse_mesh_ptr
    !-------------------------------------------------------------
 
    call meshhierarchybase_init( this, p_level_num, h_level_num )
 
    !- Setup p-mesh hierarchy
 
    if ( porder_list(p_level_num) /= 1 .and. h_level_num > 0 ) then
      log_info('MeshHierarchy2D_Init',*) 'Currently, only p=1 is supported for the coarsest p-mesh level when h_LEVEL_NUM > 0. Check!'
      call prc_abort
    end if
 
    !- Setup p-mesh hierarchy
 
    allocate( this%elem2D_list(p_level_num) )
    do poly_lev=1, p_level_num
      call this%elem2D_list(poly_lev)%Init( porder_list(poly_lev), .false. )
    end do
 
    allocate( this%p_mesh_list(this%NUM_pMG_LEVEL) )
 
    this%p_mesh_list(mesh_hierarchy_pmg_finest_level)%ptr => parent_mesh
    do poly_lev=mesh_hierarchy_pmg_finest_level + 1, this%NUM_pMG_LEVEL
      select type(parent_mesh)
      type is (meshrectdom2d)
        call construct_rectdom2d_mesh( this, this%p_mesh_list(poly_lev)%ptr, &
          parent_mesh%NeGX, parent_mesh%NeGY,     &
          parent_mesh, this%elem2D_list(poly_lev) )
      end select
    end do
 
    allocate( this%p_level(this%NUM_pMG_LEVEL) )
 
    do poly_lev=mesh_hierarchy_pmg_finest_level, this%NUM_pMG_LEVEL
      call this%p_level(poly_lev)%Init( poly_lev, this%p_mesh_list(poly_lev)%ptr, &
        this%p_mesh_list, this%NUM_pMG_LEVEL, mesh_hierarchy_type_pmg )
    end do
 
    !- Setup h-mesh hierarchy
 
    allocate( this%h_mesh_list(this%NUM_hMG_LEVEL) )
 
    this%h_mesh_list(mesh_hierarchy_hmg_finest_level)%ptr => this%p_mesh_list(this%NUM_pMG_LEVEL)%ptr
    do h_lev=mesh_hierarchy_hmg_finest_level + 1, this%NUM_hMG_LEVEL
      select type(parent_mesh)
      type is (meshrectdom2d)
        call construct_rectdom2d_mesh( this, this%h_mesh_list(h_lev)%ptr, &
          negx_list(h_lev), negy_list(h_lev),                             &
          parent_mesh, this%elem2D_list(p_level_num) )
      end select
    end do
 
    allocate( this%h_level(this%NUM_hMG_LEVEL) )
 
    do h_lev=mesh_hierarchy_hmg_finest_level, this%NUM_hMG_LEVEL
      call this%h_level(h_lev)%Init( h_lev, this%h_mesh_list(h_lev)%ptr, &
        this%h_mesh_list, this%NUM_hMG_LEVEL, mesh_hierarchy_type_hmg )
    end do
 
    return
  end subroutine meshhierarchy2d_init
 
  !> Finalize an object for mesh hierarchy in 2D domain
!OCL SERIAL
  subroutine meshhierarchy2d_final(this)
    implicit none
    class(meshhierarchy2d), intent(inout) :: this
 
    class(meshbase2d), pointer :: mesh2D_ptr
    integer :: poly_lev
    integer :: h_lev
    integer :: ldom_id
    !-------------------------------------------------------------
 
    ! p-mesh hierarchy finalization
    do poly_lev=mesh_hierarchy_pmg_finest_level, this%NUM_pMG_LEVEL
      call this%p_level(poly_lev)%Final()
    end do
    deallocate( this%p_level )
 
    do poly_lev=mesh_hierarchy_pmg_finest_level + 1, this%NUM_pMG_LEVEL
      mesh2d_ptr => this%p_mesh_list(poly_lev)%ptr
      select type(mesh2d_ptr)
      type is (meshrectdom2d)
        call mesh2d_ptr%Final()
      end select
    end do
    deallocate( this%p_mesh_list )
 
    ! h-mesh hierarchy finalization
    do h_lev=mesh_hierarchy_hmg_finest_level, this%NUM_hMG_LEVEL
      call this%h_level(h_lev)%Final()
    end do
    deallocate( this%h_level )
 
    do h_lev=mesh_hierarchy_hmg_finest_level + 1, this%NUM_hMG_LEVEL
      mesh2d_ptr => this%h_mesh_list(h_lev)%ptr
      select type(mesh2d_ptr)
      type is (meshrectdom2d)
        call mesh2d_ptr%Final()
      end select
    end do
    deallocate( this%h_mesh_list )
    
    !-
    call meshhierarchybase_final( this )
    return
  end subroutine meshhierarchy2d_final
 
!-- private --------------------------------------------------------------
 
  !> Initialize an object for mesh hierarchy level in 2D domain
!OCL SERIAL
  subroutine meshhierarchylevel2d_init( this, &
    level_id, mesh2D, mesh_list, LEVEL_NUM,   &
    hierarchy_type )
    use scale_mesh_hierarchy_base, only: &
      meshhierarchy_construct_pmg_mat1d
    implicit none
    class(meshhierarchylevel2d), intent(inout) :: this
    integer, intent(in) :: level_id
    class(meshbase2d), intent(in) :: mesh2D
    integer, intent(in) :: LEVEL_NUM
    class(meshptr2d), intent(in), target :: mesh_list(LEVEL_NUM)
    integer, intent(in) :: hierarchy_type
 
    integer :: ldom_id
    !-------------------------------------------------------------
 
    this%hierarchy_type = hierarchy_type
    this%level_id = level_id
 
    !- set mesh pointers with finer/coarser meshes
 
    if ( level_id > 1 ) then
      this%fine_mesh => mesh_list(level_id-1)
    else
      this%fine_mesh => null()
    end if
 
    if ( level_id < level_num ) then
      this%coarse_mesh => mesh_list(level_id+1)
    else
      this%coarse_mesh => null()
    end if
 
    !- p-hierarchy
    if ( hierarchy_type == mesh_hierarchy_type_pmg &
         .and. level_id < level_num                ) then
      
      allocate( this%pMat1D_c2f(mesh2d%refElem2D%Nfp,this%coarse_mesh%ptr%refElem2D%Nfp) )
      call meshhierarchy_construct_pmg_mat1d( this%pMat1D_c2f,   &
        this%coarse_mesh%ptr%refElem2D%Nfp, mesh2d%refElem2D%Nfp )
 
      allocate( this%pMat1D_f2c(this%coarse_mesh%ptr%refElem2D%Nfp,mesh2d%refElem2D%Nfp) )
      call meshhierarchy_construct_pmg_mat1d( this%pMat1D_f2c,   &
        mesh2d%refElem2D%Nfp, this%coarse_mesh%ptr%refElem2D%Nfp )
    end if
 
    !- h-hierarchy
    if ( hierarchy_type == mesh_hierarchy_type_hmg &
         .and. level_id < level_num                ) then
 
      allocate( this%mg_local( mesh2d%LOCAL_MESH_NUM ) )
 
      do ldom_id=1, mesh2d%LOCAL_MESH_NUM
        call this%mg_local(ldom_id)%Init( mesh_list(level_id)%ptr%lcmesh_list(ldom_id), &
          this%coarse_mesh%ptr%lcmesh_list, mesh_list(level_id)%ptr%refElem2D           )
      end do
    end if
 
    return
  end subroutine meshhierarchylevel2d_init
 
  !> Finalize an object for mesh hierarchy level in 2D domain
!OCL SERIAL
  subroutine meshhierarchylevel2d_final( this )
    implicit none
    class(meshhierarchylevel2d), intent(inout) :: this
 
    integer :: ldom_id
    !-------------------------------------------------------------
 
    !-
    if ( allocated(this%pMat1D_c2f) ) deallocate(this%pMat1D_c2f)
    if ( allocated(this%pMat1D_f2c) ) deallocate(this%pMat1D_f2c)
 
    !-
    if ( allocated(this%mg_local) ) then
      do ldom_id=1, size(this%mg_local)
        call this%mg_local(ldom_id)%Final()
      end do
      deallocate( this%mg_local )
    end if
 
    return
  end subroutine meshhierarchylevel2d_final
 
  !> Initialize an object for 2D local mesh data for multigrid
!OCL SERIAL
  subroutine meshhierarchylocalmgdata2d_init( this, &
    lcmesh2D, coarse_lcmesh_list, elem2D )
    use scale_mesh_hierarchy_base, only: meshhierarchylocalmgdatabase_init
    implicit none
    class(meshhierarchylocalmgdata2d), intent(inout) :: this
    class(localmesh2d), intent(in) :: lcmesh2D
    class(localmesh2d), intent(in), target :: coarse_lcmesh_list(:)
    class(elementbase2d), intent(in) :: elem2D
 
    integer :: i, j
    integer :: ke
 
    integer :: i_c, j_c
    integer :: ke_c
 
    integer :: ke2i(lcmesh2D%Ne)
    integer :: ke2j(lcmesh2D%Ne)
 
    integer :: lcdomID_c
    integer :: lcTileID_c
    class(localmesh2d), pointer :: lcmesh_c
 
    real(RP) :: vx_c(elem2D%Nv)
    real(RP) :: vy_c(elem2D%Nv)
    integer :: i_EtoV(elem2D%Nv)
 
    integer :: p
    real(RP) :: r_c, s_c
    !-------------------------------------------------------------  
 
    call meshhierarchylocalmgdatabase_init( this, &
      lcmesh2d )
 
    
    ! Current implementation assumes identical tileID & lcdomID in both coarsened meshes          
    this%CoarseLocalMesh_tileIDlist(1) = lcmesh2d%tileID
    this%CoarseLocalMesh_lcdomIDlist(1) = 1
 
    ! Preparation
    do j=1, lcmesh2d%NeY
    do i=1, lcmesh2d%NeX
      ke = i + (j-1)*lcmesh2d%NeX
      ke2i(ke) = i; ke2j(ke) = j
    end do
    end do
 
    ! Set the relation between coarse and fine grid,
    ! and construct interpolation operator
 
!    write(*,*) "  Constructing local MG data for h-mesh..."
 
    allocate( this%If2c_emap(4,lcmesh2d%Ne/4) )
 
    do ke=lcmesh2d%NeS, lcmesh2d%NeE
      lctileid_c = this%CoarseLocalMesh_tileIDlist(1)
      lcdomid_c = this%CoarseLocalMesh_lcdomIDlist(1)
 
      lcmesh_c => coarse_lcmesh_list(lcdomid_c)
 
      i_c = (ke2i(ke)+1)/2; j_c = (ke2j(ke)+1)/2
      ke_c = i_c + (j_c-1)*lcmesh2d%NeX/2
      this%Ic2f_emap(ke) = ke_c
 
      i = ke2i(ke) - 2*(i_c-1)
      j = ke2j(ke) - 2*(j_c-1)
      this%If2c_emap(i+(j-1)*2,ke_c) = ke
 
      i_etov(:) = lcmesh_c%EToV(ke_c,:)
      vx_c(:) = lcmesh_c%pos_ev(i_etov(:),1)
      vy_c(:) = lcmesh_c%pos_ev(i_etov(:),2)
 
      do p=1, elem2d%Np
        r_c = -1.0_rp + 2.0_rp * ( lcmesh2d%pos_en(p,ke,1) - vx_c(1) ) / ( vx_c(2) - vx_c(1) )
        s_c = -1.0_rp + 2.0_rp * ( lcmesh2d%pos_en(p,ke,2) - vy_c(1) ) / ( vy_c(3) - vy_c(1) )
 
        this%Ic2f(p,1:4,ke) = 0.25_rp * &
          (/ ( 1.0_rp - r_c ) * ( 1.0_rp - s_c ), ( 1.0_rp + r_c ) * ( 1.0_rp - s_c ),  &
             ( 1.0_rp - r_c ) * ( 1.0_rp + s_c ), ( 1.0_rp + r_c ) * ( 1.0_rp + s_c )  /)
      end do
    end do  
    return
  end subroutine meshhierarchylocalmgdata2d_init
 
  !> Finalize an object for 2D local mesh data for multigrid
!OCL SERIAL
  subroutine meshhierarchylocalmgdata2d_final( this )
    use scale_mesh_hierarchy_base, only: meshhierarchylocalmgdatabase_final
    implicit none
    class(meshhierarchylocalmgdata2d), intent(inout) :: this
    !-------------------------------------------------------------  
    call meshhierarchylocalmgdatabase_final( this )
    return
  end subroutine meshhierarchylocalmgdata2d_final
 
!-
!OCL SERIAL
  subroutine construct_rectdom2d_mesh( this, child_mesh_base_ptr, &
    NeGX, NeGY,         &
    parent_mesh, elem2D )
    implicit none
    class(meshhierarchy2d), intent(inout) :: this
    class(meshbase2d), intent(out), pointer :: child_mesh_base_ptr
    integer, intent(in) :: NeGX
    integer, intent(in) :: NeGY
    type(meshrectdom2d), intent(in) :: parent_mesh
    type(quadrilateralelement), intent(in) :: elem2D
 
    type(meshrectdom2d), pointer :: child_mesh_ptr
    !-------------------------------------------------------------
 
    allocate( child_mesh_ptr )
    call child_mesh_ptr%Init( negx, negy,                                                  &
      parent_mesh%xmin_gl, parent_mesh%xmax_gl, parent_mesh%ymin_gl, parent_mesh%ymax_gl,  &
      parent_mesh%isPeriodicX, parent_mesh%isPeriodicY,                                    &
      elem2d, 1, parent_mesh%NprcX, parent_mesh%NprcY )
 
    call child_mesh_ptr%Generate()
 
    child_mesh_base_ptr => child_mesh_ptr
    return
  end subroutine construct_rectdom2d_mesh
end module scale_mesh_hierarchy_2d