scale_mesh_hierarchy_3d.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Mesh / 3D domain
!!
!! @par Description
!!      Manage mesh hierarchy of 3D domain for element-based methods
!!
!! @author Yuta Kawai, Team SCALE
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_mesh_hierarchy_3d
 
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  !
  use scale_precision
  use scale_io
  use scale_prc, only: prc_abort
 
  use scale_element_base, only: elementbase3d
  use scale_element_hexahedral, only: hexahedralelement
 
  use scale_localmesh_3d, only: localmesh3d
  use scale_mesh_base3d, only: meshbase3d
  use scale_mesh_cubedom3d, only: meshcubedom3d
 
  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 MeshBase3D
  type :: meshptr3d
    class(MeshBase3D), pointer :: ptr => null()
  end type meshptr3d
 
  !> Derived type to manage 3D local mesh data for multigrid
  type, extends(meshhierarchylocalmgdatabase), public :: meshhierarchylocalmgdata3d
  contains
    procedure :: init => meshhierarchylocalmgdata3d_init
    procedure :: final => meshhierarchylocalmgdata3d_final
  end type meshhierarchylocalmgdata3d
 
  !> Derived type to represent mesh hierarchy level in 3D domain
  type, extends(meshhierarchylevelbase), public :: meshhierarchylevel3d
    type(meshptr3d), pointer :: fine_mesh => null()
    type(meshptr3d), pointer :: coarse_mesh => null()
    type(meshhierarchylocalmgdata3d), allocatable :: mg_local(:)
 
    real(rp), allocatable :: pmat1d_f2c(:,:)
    real(rp), allocatable :: pmat1d_c2f(:,:)
  contains
    procedure :: init => meshhierarchylevel3d_init
    procedure :: final => meshhierarchylevel3d_final
  end type meshhierarchylevel3d
 
  !> Derived type for mesh hierarchy in 3D domain
  type, extends(meshhierarchybase), public :: meshhierarchy3d
    type(meshhierarchylevel3d), allocatable :: p_level(:)
    type(meshhierarchylevel3d), allocatable :: h_level(:)
 
    type(meshptr3d), allocatable :: p_mesh_list(:)
    type(meshptr3d), allocatable :: h_mesh_list(:)
 
    type(hexahedralelement), allocatable :: elem3d_list(:)
  contains
    procedure :: init => meshhierarchy3d_init
    procedure :: final => meshhierarchy3d_final
  end type meshhierarchy3d
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
contains
  !> Initialize an object for mesh hierarchy in 3D domain
!OCL SERIAL
  subroutine meshhierarchy3d_init( this, &
    parent_mesh,                                  &
    porder_list, p_LEVEL_NUM,                     &
    NeGX_list, NeGY_list, NeGZ_list, h_LEVEL_NUM  )
    implicit none
    class(meshhierarchy3d), intent(inout), target :: this
    class(meshbase3d), 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, intent(in) :: NeGZ_list(h_LEVEL_NUM)         !< Number of elements in z-direction for h-mesh levels
 
    integer :: poly_lev
    integer :: h_lev
 
    integer :: ldom_id
 
    type(meshhierarchylevel3d), pointer :: level_ptr
    class(meshbase3d), pointer :: mesh3D_ptr
 
    class(meshptr3d), pointer :: fine_mesh_ptr
    class(meshptr3d), pointer :: coarse_mesh_ptr
 
    integer :: NeGZ
    !-------------------------------------------------------------
 
    call meshhierarchybase_init( this, p_level_num, h_level_num )
 
    !--
    if ( parent_mesh%refElem3D%PolyOrder_h /= parent_mesh%refElem3D%PolyOrder_v ) then
      log_info('MeshHierarchy3D_Init',*) 'Currently, PolyOrder_h should equal PolyOrder_v. Check!'
      call prc_abort
    end if
    if ( porder_list(p_level_num) /= 1 .and. h_level_num > 0 ) then
      log_info('MeshHierarchy3D_Init',*) 'Currently, only p=1 is supported for the coarsest p-mesh level. Check!'
      call prc_abort
    end if
    if ( h_level_num > 0 ) then
      negz = negz_list(1)
      do h_lev=2, h_level_num
        if ( negz_list(h_lev) /= negz ) then
          log_info('MeshHierarchy3D_Init',*) 'Currently, the number of elements should be same in the vertical direction. Check!'
          call prc_abort
        end if
      end do
    end if
 
    !- Setup p-mesh hierarchy
 
    allocate( this%elem3D_list(p_level_num) )
    do poly_lev=1, p_level_num
      call this%elem3D_list(poly_lev)%Init( porder_list(poly_lev), 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 (meshcubedom3d)
        call construct_cubedom3d_mesh( this, this%p_mesh_list(poly_lev)%ptr,    &
          parent_mesh%NeGX, parent_mesh%NeGY, parent_mesh%NeGZ, parent_mesh%FZ, &
          parent_mesh, this%elem3D_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
 
    if ( this%NUM_hMG_LEVEL > 0 ) then
      allocate( this%h_mesh_list(this%NUM_hMG_LEVEL) )
 
      this%h_mesh_list(mesh_hierarchy_hmg_finest_level)%ptr => parent_mesh
      do h_lev=mesh_hierarchy_hmg_finest_level + 1, this%NUM_hMG_LEVEL
        select type(parent_mesh)
        type is (meshcubedom3d)
          call construct_cubedom3d_mesh( this, this%h_mesh_list(h_lev)%ptr,       &
            negx_list(h_lev), negy_list(h_lev), negz_list(h_lev), parent_mesh%FZ, &
            parent_mesh, this%elem3D_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
    end if
 
    return
  end subroutine meshhierarchy3d_init
 
  !> Finalize an object for mesh hierarchy in 3D domain
!OCL SERIAL
  subroutine meshhierarchy3d_final(this)
    implicit none
    class(meshhierarchy3d), intent(inout) :: this
 
    class(meshbase3d), pointer :: mesh3D_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
      mesh3d_ptr => this%p_mesh_list(poly_lev)%ptr
      select type(mesh3d_ptr)
      type is (meshcubedom3d)
        call mesh3d_ptr%Final()
      end select
    end do
    deallocate( this%p_mesh_list )
 
    ! h-mesh hierarchy finalization
    if ( this%NUM_hMG_LEVEL > 0 ) then
      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
        mesh3d_ptr => this%h_mesh_list(h_lev)%ptr
        select type(mesh3d_ptr)
        type is (meshcubedom3d)
          call mesh3d_ptr%Final()
        end select
      end do
      deallocate( this%h_mesh_list )
    end if
    
    !-
    call meshhierarchybase_final( this )
    return
  end subroutine meshhierarchy3d_final
 
!-- private --------------------------------------------------------------
 
  !> Initialize an object for mesh hierarchy level in 3D domain
!OCL SERIAL
  subroutine meshhierarchylevel3d_init( this, &
    level_id, mesh3D, mesh_list, LEVEL_NUM,   &
    hierarchy_type )
    use scale_mesh_hierarchy_base, only: &
      meshhierarchy_construct_pmg_mat1d
    implicit none
    class(meshhierarchylevel3d), intent(inout) :: this
    integer, intent(in) :: level_id
    class(meshbase3d), intent(in) :: mesh3D
    integer, intent(in) :: LEVEL_NUM
    class(meshptr3d), 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(mesh3d%refElem3D%Nnode_h1D,this%coarse_mesh%ptr%refElem3D%Nnode_h1D) )
      call meshhierarchy_construct_pmg_mat1d( this%pMat1D_c2f,   &
        this%coarse_mesh%ptr%refElem3D%Nnode_h1D, mesh3d%refElem3D%Nnode_h1D )
 
      allocate( this%pMat1D_f2c(this%coarse_mesh%ptr%refElem3D%Nnode_h1D,mesh3d%refElem3D%Nnode_h1D) )
      call meshhierarchy_construct_pmg_mat1d( this%pMat1D_f2c,   &
        mesh3d%refElem3D%Nnode_h1D, this%coarse_mesh%ptr%refElem3D%Nnode_h1D )
    end if
 
    !- h-hierarchy
    if ( hierarchy_type == mesh_hierarchy_type_hmg &
         .and. level_id < level_num                ) then
 
      allocate( this%mg_local( mesh3d%LOCAL_MESH_NUM ) )
      
      do ldom_id=1, mesh3d%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%refElem3D           )
      end do
    end if
 
    return
  end subroutine meshhierarchylevel3d_init
 
  !> Finalize an object for mesh hierarchy level in 3D domain
!OCL SERIAL
  subroutine meshhierarchylevel3d_final( this )
    implicit none
    class(meshhierarchylevel3d), intent(inout) :: this
 
    integer :: ldom_id
    !-------------------------------------------------------------
 
    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 meshhierarchylevel3d_final
 
!> Initialize an object to manage 3D local mesh data for multigrid
!OCL SERIAL
  subroutine meshhierarchylocalmgdata3d_init( this, &
    lcmesh3D, coarse_lcmesh_list, elem3D )
    use scale_mesh_hierarchy_base, only: meshhierarchylocalmgdatabase_init
    implicit none
    class(meshhierarchylocalmgdata3d), intent(inout) :: this
    class(localmesh3d), intent(in) :: lcmesh3D
    class(localmesh3d), intent(in), target :: coarse_lcmesh_list(:)
    class(elementbase3d), intent(in) :: elem3D
 
    integer :: i, j, k
    integer :: ke
 
    integer :: i_c, j_c, k_c
    integer :: ke_c
 
    integer :: ke2i(lcmesh3D%Ne)
    integer :: ke2j(lcmesh3D%Ne)
    integer :: ke2k(lcmesh3D%Ne)
 
    integer :: lcdomID_c
    integer :: lcTileID_c
    class(localmesh3d), pointer :: lcmesh_c
 
    real(RP) :: vx_c(elem3D%Nv)
    real(RP) :: vy_c(elem3D%Nv)
    real(RP) :: vz_c(elem3D%Nv)
    integer :: i_EtoV(elem3D%Nv)
 
    integer :: p
    real(RP) :: r_c, s_c, t_c
 
    integer :: l
    !-------------------------------------------------------------  
 
    call meshhierarchylocalmgdatabase_init( this, &
      lcmesh3d )
 
    ! Current implementation assumes identical tileID & lcdomID in both coarsened meshes          
    this%CoarseLocalMesh_tileIDlist(1) = lcmesh3d%tileID
    this%CoarseLocalMesh_lcdomIDlist(1) = 1
 
    ! Preparation
    do k=1, lcmesh3d%NeZ
    do j=1, lcmesh3d%NeY
    do i=1, lcmesh3d%NeX
      ke = i + (j-1)*lcmesh3d%NeX + (k-1)*lcmesh3d%NeX*lcmesh3d%NeY
      ke2i(ke) = i; ke2j(ke) = j; ke2k(ke) = k
    end do
    end do
    end do
 
    ! Set the relation between coarse and fine grid,
    ! and construct interpolation operator
 
    allocate( this%If2c_emap(4,lcmesh3d%Ne/4) )
 
    do ke=lcmesh3d%NeS, lcmesh3d%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; k_c = ke2k(ke)
      ke_c = i_c + (j_c-1)*lcmesh_c%NeX + (k_c-1)*lcmesh_c%NeX*lcmesh_c%NeY
      this%Ic2f_emap(ke) = ke_c
 
      i = ke2i(ke) - 2*(i_c-1)
      j = ke2j(ke) - 2*(j_c-1)
!      k = ke2k(ke) - 2*(k_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)
      vz_c(:) = lcmesh_c%pos_ev(i_etov(:),3)
      do p=1, elem3d%Np
        r_c = -1.0_rp + 2.0_rp * ( lcmesh3d%pos_en(p,ke,1) - vx_c(1) ) / ( vx_c(2) - vx_c(1) )
        s_c = -1.0_rp + 2.0_rp * ( lcmesh3d%pos_en(p,ke,2) - vy_c(1) ) / ( vy_c(4) - vy_c(1) )
        t_c = -1.0_rp + 2.0_rp * ( lcmesh3d%pos_en(p,ke,3) - vz_c(1) ) / ( vz_c(5) - vz_c(1) )
 
        this%Ic2f(p,:,ke) = 0.125_rp * &
            (/ ( 1.0_rp - r_c ) * ( 1.0_rp - s_c ) * ( 1.0_rp - t_c ), ( 1.0_rp + r_c ) * ( 1.0_rp - s_c ) * ( 1.0_rp - t_c ),  &
               ( 1.0_rp - r_c ) * ( 1.0_rp + s_c ) * ( 1.0_rp - t_c ), ( 1.0_rp + r_c ) * ( 1.0_rp + s_c ) * ( 1.0_rp - t_c ),  &
               ( 1.0_rp - r_c ) * ( 1.0_rp - s_c ) * ( 1.0_rp + t_c ), ( 1.0_rp + r_c ) * ( 1.0_rp - s_c ) * ( 1.0_rp + t_c ),  &
               ( 1.0_rp - r_c ) * ( 1.0_rp + s_c ) * ( 1.0_rp + t_c ), ( 1.0_rp + r_c ) * ( 1.0_rp + s_c ) * ( 1.0_rp + t_c )  /)
      end do
    end do  
    return
  end subroutine meshhierarchylocalmgdata3d_init
 
!> Finalize an object to manage 3D local mesh data for multigrid
!OCL SERIAL
  subroutine meshhierarchylocalmgdata3d_final( this )
    use scale_mesh_hierarchy_base, only: meshhierarchylocalmgdatabase_final
    implicit none
    class(meshhierarchylocalmgdata3d), intent(inout) :: this
    !-------------------------------------------------------------  
    call meshhierarchylocalmgdatabase_final( this )
    return
  end subroutine meshhierarchylocalmgdata3d_final
 
!-
!OCL SERIAL
  subroutine construct_cubedom3d_mesh( this, child_mesh_base_ptr, &
    NeGX, NeGY, NeGZ, FZ, &
    parent_mesh, elem3D   )
    implicit none
    class(meshhierarchy3d), intent(inout) :: this
    class(meshbase3d), intent(out), pointer :: child_mesh_base_ptr
    integer, intent(in) :: NeGX
    integer, intent(in) :: NeGY
    integer, intent(in) :: NeGZ
    real(RP), intent(in) :: FZ(NeGZ+1)
    type(meshcubedom3d), intent(in) :: parent_mesh
    type(hexahedralelement), intent(in) :: elem3D
 
    type(meshcubedom3d), pointer :: child_mesh_ptr
    !-------------------------------------------------------------
 
    allocate( child_mesh_ptr )
    call child_mesh_ptr%Init( negx, negy, negz,                                            &
      parent_mesh%xmin_gl, parent_mesh%xmax_gl, parent_mesh%ymin_gl, parent_mesh%ymax_gl,  &
      parent_mesh%zmin_gl, parent_mesh%zmax_gl,                                            &
      parent_mesh%isPeriodicX, parent_mesh%isPeriodicY, parent_mesh%isPeriodicZ,           &
      elem3d, 1, parent_mesh%NprcX, parent_mesh%NprcY,                                     &
      fz=fz )
 
    call child_mesh_ptr%Generate()
 
    child_mesh_base_ptr => child_mesh_ptr
    return
  end subroutine construct_cubedom3d_mesh
end module scale_mesh_hierarchy_3d