scale_localmesh_3d.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Mesh / Local 3D
!!
!! @par Description
!!      Module to manage 3D local mesh for element-based methods
!!
!! @author Yuta Kawai, Team SCALE
!<
#include "scaleFElib.h"
module scale_localmesh_3d
 
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_localmesh_base, only: &
    localmeshbase, localmeshbase_init, localmeshbase_final
  use scale_localmesh_2d, only: localmesh2d
  use scale_element_base, only: elementbase, elementbase3d
 
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  !
  
  !> Derived type to manage a local 3D computational domain 
  type, extends(localmeshbase), public :: localmesh3d
 
    type(elementbase3d), pointer :: refelem3d !< Pointer to a object of 3D reference element
    real(rp) :: xmin  !< Minimum x-coordinate value of the local computational domain
    real(rp) :: xmax  !< Maximum x-coordinate value of the local computational domain
    real(rp) :: ymin  !< Minimum y-coordinate value of the local computational domain
    real(rp) :: ymax  !< Maximum y-coordinate value of the local computational domain
    real(rp) :: zmin  !< Minimum z-coordinate value of the local computational domain
    real(rp) :: zmax  !< Maximum z-coordinate value of the local computational domain
 
    integer :: nex    !< Number of finite elements for the x-direction in the local computational domain
    integer :: ney    !< Number of finite elements for the y-direction in the local computational domain
    integer :: nez    !< Number of finite elements for the z-direction in the local computational domain
    integer :: ne2d   !< NeX * NeY
    integer :: ne2da  !< NeX * NeY + halo size
 
    real(rp), allocatable :: sz(:,:)
    real(rp), allocatable :: zs(:,:)
    real(rp), allocatable :: gi3(:,:,:)    !< The contravariant component of metric tensor with vertical general coordinate 
    real(rp), allocatable :: gsqrth(:,:)   !< The Jacobian of horizontal transformation in the computational coordinate
    real(rp), allocatable :: zlev(:,:)     !< z-coordinates (actual level)
    real(rp), allocatable :: gam(:,:)      !< Factor for approximation with spherical shell domain (= r/a)
 
    class(localmesh2d), pointer :: lcmesh2d !< Pointer to an object to manage 2D local computational domain
    real(rp), allocatable :: lon2d(:,:)     !< Longitude coordinate with 2D mesh
    real(rp), allocatable :: lat2d(:,:)     !< Latitude coordinate with 2D mesh
 
    integer, allocatable :: emap3dto2d(:)   !< Array to map 3D element ID to 2D element ID
  contains
    procedure :: setlocalmesh2d => localmesh3d_setlocalmesh2d 
    procedure :: getvmapz1d => localmesh3d_getvmapz1d 
    procedure :: getvmapz3d => localmesh3d_getvmapz3d
  end type localmesh3d
 
  public :: localmesh3d_init, localmesh3d_final
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
 
contains
!> Initialize an object to manage a 3D local computational domain
!OCL SERIAL
  subroutine localmesh3d_init( this, &
    lcdomID, refElem, myrank )
    implicit none
 
    class(localmesh3d), intent(inout) :: this
    integer, intent(in) :: lcdomid
    class(elementbase3d), intent(in), target :: refelem
    integer, intent(in), optional :: myrank
    !-------------------------------------------------
 
    this%refElem3D    => refelem
    nullify( this%lcmesh2D )
 
    call localmeshbase_init( this, lcdomid, refelem, 3, myrank )
    !$acc enter data copyin(this)
    return
  end subroutine localmesh3d_init
 
!> Finalize an object to manage a 3D local computational domain
!OCL SERIAL
  subroutine localmesh3d_final( this, is_generated )
    implicit none
    type(localmesh3d), intent(inout) :: this
    logical, intent(in) :: is_generated
    !-------------------------------------------------
 
    call localmeshbase_final( this, is_generated )
    if (is_generated) then
      !$acc exit data delete(this%G_ij, this%GIJ)
      deallocate( this%G_ij, this%GIJ )
      !$acc exit data delete( this%Sz, this%zS, this%GI3, this%GsqrtH, this%zlev, this%gam, this%lon2D, this%lat2D, this%EMap3Dto2D )
      deallocate( this%zS, this%Sz )
      deallocate( this%GI3, this%GsqrtH )
      deallocate( this%zlev )
      deallocate( this%gam )
      deallocate( this%lon2D, this%lat2D )
      deallocate( this%EMap3Dto2D )
    end if
    !$acc exit data delete(this)
 
    return
  end subroutine localmesh3d_final
 
!> Set a pointer to an object to mange 2D local computational domain
!OCL SERIAL
  subroutine localmesh3d_setlocalmesh2d( this, lcmesh2D )
    implicit none
    class(localmesh3d), intent(inout) :: this
    class(localmesh2d), intent(in), target :: lcmesh2d
    !-------------------------------------------------
 
    this%lcmesh2D => lcmesh2d
    return
  end subroutine localmesh3d_setlocalmesh2d
 
!> Get arrays for vertical mapping node index with vertical element boundary to that with DG data
!OCL SERIAL
  subroutine localmesh3d_getvmapz1d( this, vmapM, vmapP )
 
    implicit none
    class(localmesh3d), intent(in), target :: this
    integer, intent(out) :: vmapm(this%refelem3d%nfptot,this%nez)
    integer, intent(out) :: vmapp(this%refelem3d%nfptot,this%nez)    
 
    integer :: ke_z
    integer :: f
    integer :: vs, ve
 
    class(elementbase3d), pointer :: elem
    !------------------------------
 
    elem => this%refElem3D
 
    do ke_z=1, this%NeZ
      do f=1, elem%Nfaces_h
        vs = 1 + (f-1)*elem%Nfp_h
        ve = vs + elem%Nfp_h - 1
        vmapm(vs:ve,ke_z) = elem%Fmask_h(:,f) + (ke_z-1)*elem%Np
      end do
      do f=1, elem%Nfaces_v
        vs = elem%Nfp_h*elem%Nfaces_h + 1 + (f-1)*elem%Nfp_v
        ve = vs + elem%Nfp_v - 1
        vmapm(vs:ve,ke_z) = elem%Fmask_v(:,f) + (ke_z-1)*elem%Np
      end do
      vmapp(:,ke_z) = vmapm(:,ke_z)
    end do
 
    do ke_z=1, this%NeZ
      vs = elem%Nfp_h*elem%Nfaces_h + 1
      ve = vs + elem%Nfp_v - 1
      if (ke_z > 1) &
        vmapp(vs:ve,ke_z) = elem%Fmask_v(:,2) + (ke_z-2)*elem%Np
 
      vs = elem%Nfp_h*elem%Nfaces_h + elem%Nfp_v + 1
      ve = vs + elem%Nfp_v - 1
      if (ke_z < this%NeZ) &
        vmapp(vs:ve,ke_z) = elem%Fmask_v(:,1) + ke_z*elem%Np
    end do
 
    return
  end subroutine localmesh3d_getvmapz1d
 
!> Get arrays for vertical mapping node index with vertical element boundary to that with DG data
!OCL SERIAL
  subroutine localmesh3d_getvmapz3d( this, vmapM, vmapP )
 
    implicit none
    class(localmesh3d), intent(in), target :: this
    integer, intent(out) :: vmapm(this%refelem3d%nfptot,this%ne)
    integer, intent(out) :: vmapp(this%refelem3d%nfptot,this%ne)    
 
    integer :: ke, ke_xy, ke_z
    integer :: ke_neigh
    integer :: f
    integer :: vs, ve
 
    class(elementbase3d), pointer :: elem
    !------------------------------
 
    elem => this%refElem3D
 
    !$omp parallel do collapse(2) private(ke,ke_xy,ke_z,ke_neigh,f,vs,ve)
    do ke_z=1, this%NeZ
    do ke_xy=1, this%NeX * this%NeY
      ke = ke_xy + (ke_z-1) * this%NeX * this%NeY
 
      do f=1, elem%Nfaces_h
        vs = 1 + (f-1)*elem%Nfp_h
        ve = vs + elem%Nfp_h - 1
        vmapm(vs:ve,ke) = elem%Fmask_h(:,f) + (ke-1)*elem%Np
      end do
      do f=1, elem%Nfaces_v
        vs = elem%Nfp_h*elem%Nfaces_h + 1 + (f-1)*elem%Nfp_v
        ve = vs + elem%Nfp_v - 1
        vmapm(vs:ve,ke) = elem%Fmask_v(:,f) + (ke-1)*elem%Np
      end do
      vmapp(:,ke) = vmapm(:,ke)
 
      !--
      vs = elem%Nfp_h*elem%Nfaces_h + 1
      ve = vs + elem%Nfp_v - 1
      if (ke_z > 1) then
        ke_neigh = ke_xy + (ke_z-2) * this%NeX * this%NeY
        vmapp(vs:ve,ke) = elem%Fmask_v(:,2) + (ke_neigh-1)*elem%Np
      end if
      vs = elem%Nfp_h*elem%Nfaces_h + elem%Nfp_v + 1
      ve = vs + elem%Nfp_v - 1
      if (ke_z < this%NeZ) then
        ke_neigh = ke_xy + ke_z * this%NeX * this%NeY
        vmapp(vs:ve,ke) = elem%Fmask_v(:,1) + (ke_neigh-1)*elem%Np
      end if
    end do
    end do
 
    return
  end subroutine localmesh3d_getvmapz3d
 
end module scale_localmesh_3d