scale_mesh_base3d.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Mesh / Base 3D
!!
!! @par Description
!!      Base module to manage 3D meshes for element-based methods
!!
!! @author Yuta Kawai, Team SCALE
!<
#include "scaleFElib.h"
module scale_mesh_base3d
 
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
 
  use scale_localmesh_3d, only: &
    localmesh3d, localmesh3d_init, localmesh3d_final
  use scale_localmesh_2d, only: &
    localmesh2d
  use scale_mesh_base, only: &
    meshbase, meshbase_init, meshbase_final
  use scale_mesh_base2d, only: &
    meshbase2d
  use scale_element_base, only: elementbase3d
 
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  ! 
  !> Derived type to manage a computational mesh (base type for 3D domain)
  type, abstract, public, extends(meshbase) :: meshbase3d
    type(localmesh3d), allocatable :: lcmesh_list(:)
    type(elementbase3d), pointer :: refelem3d
  contains
    procedure(meshbase3d_generate), deferred :: generate 
    procedure(meshbase3d_getmesh2d), deferred  :: getmesh2d
    procedure(meshbase3d_set_geometric_with_vcoord), deferred :: set_geometric_with_vcoord
    procedure :: getlocalmesh => meshbase3d_get_localmesh
  end type meshbase3d
 
  interface 
    subroutine meshbase3d_generate(this)
      import meshbase3d
      class(meshbase3d), intent(inout), target :: this
    end subroutine meshbase3d_generate
 
    subroutine meshbase3d_getmesh2d(this, ptr_mesh2D)
      import meshbase3d
      import meshbase2d
      class(meshbase3d), intent(in), target :: this
      class(meshbase2d), pointer, intent(out) :: ptr_mesh2D
    end subroutine meshbase3d_getmesh2d
 
    subroutine meshbase3d_set_geometric_with_vcoord(this, lcdomID, GsqrtV_lc, zlev_lc, G13_lc, G23_lc)
      import rp
      import meshbase3d
      class(meshbase3d), intent(inout), target :: this
      integer, intent(in) :: lcdomID
      real(RP), intent(in) :: GsqrtV_lc(this%refElem3D%Np,this%lcmesh_list(lcdomID)%NeA)
      real(RP), intent(in) :: zlev_lc(this%refElem3D%Np,this%lcmesh_list(lcdomID)%NeA)
      real(RP), intent(in) :: G13_lc(this%refElem3D%Np,this%lcmesh_list(lcdomID)%NeA)
      real(RP), intent(in) :: G23_lc(this%refElem3D%Np,this%lcmesh_list(lcdomID)%NeA)
    end subroutine meshbase3d_set_geometric_with_vcoord
  end interface
 
  public :: meshbase3d_init, meshbase3d_final
  public :: meshbase3d_setgeometricinfo
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  integer, public :: meshbase3d_dimtype_num    = 8
  integer, public :: meshbase3d_dimtypeid_x    = 1
  integer, public :: meshbase3d_dimtypeid_y    = 2
  integer, public :: meshbase3d_dimtypeid_z    = 3
  integer, public :: meshbase3d_dimtypeid_zt   = 4
  integer, public :: meshbase3d_dimtypeid_xy   = 5
  integer, public :: meshbase3d_dimtypeid_xyt  = 6
  integer, public :: meshbase3d_dimtypeid_xyz  = 7
  integer, public :: meshbase3d_dimtypeid_xyzt = 8
  
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
 
contains
!OCL SERIAL
  subroutine meshbase3d_init(this,        &
    refElem, NLocalMeshPerPrc, NsideTile, &
    nproc, myrank )
    
    implicit none
 
    class(meshbase3d), intent(inout) :: this
    class(elementbase3d), intent(in), target :: refelem
    integer, intent(in) :: nlocalmeshperprc
    integer, intent(in) :: nsidetile
    integer, intent(in), optional :: nproc
    integer, intent(in), optional :: myrank
 
    integer :: n
    !-----------------------------------------------------------------------------
    
    this%refElem3D => refelem
    call meshbase_init( this,            &
      meshbase3d_dimtype_num, refelem,   &
      nlocalmeshperprc, nsidetile, nproc )
 
    allocate( this%lcmesh_list(this%LOCAL_MESH_NUM) )
    do n=1, this%LOCAL_MESH_NUM
      call localmesh3d_init( this%lcmesh_list(n), n, refelem, myrank )
    end do
 
    call this%SetDimInfo( meshbase3d_dimtypeid_x, "x", "m", "X-coordinate" )
    call this%SetDimInfo( meshbase3d_dimtypeid_y, "y", "m", "Y-coordinate" )
    call this%SetDimInfo( meshbase3d_dimtypeid_z, "z", "m", "Z-coordinate" )
    call this%SetDimInfo( meshbase3d_dimtypeid_zt, "z", "m", "Z-coordinate" )
    call this%SetDimInfo( meshbase3d_dimtypeid_xy, "xy", "m", "XY-coordinate" )
    call this%SetDimInfo( meshbase3d_dimtypeid_xyt, "xyt", "m", "XY-coordinate" )
    call this%SetDimInfo( meshbase3d_dimtypeid_xyz, "xyz", "m", "XYZ-coordinate" )
    call this%SetDimInfo( meshbase3d_dimtypeid_xyzt, "xyzt", "m", "XYZ-coordinate" )
 
    return
  end subroutine meshbase3d_init
 
!OCL SERIAL
  subroutine meshbase3d_final( this )
    implicit none
    class(meshbase3d), intent(inout) :: this
 
    integer :: n
    !-----------------------------------------------------------------------------
    
    if ( allocated ( this%lcmesh_list ) ) then 
      do n=1, this%LOCAL_MESH_NUM
        call localmesh3d_final( this%lcmesh_list(n), this%isGenerated )
      end do
  
      deallocate( this%lcmesh_list )
    end if
    
    call meshbase_final(this)
 
    return
  end subroutine meshbase3d_final
  
!OCL SERIAL
  subroutine meshbase3d_get_localmesh( this, id, ptr_lcmesh )
    use scale_localmesh_base, only: localmeshbase
    implicit none
 
    class(meshbase3d), target, intent(in) :: this
    integer, intent(in) :: id
    class(localmeshbase), pointer, intent(out) :: ptr_lcmesh
    !-------------------------------------------------------------
 
    ptr_lcmesh => this%lcmesh_list(id)
    return
  end subroutine meshbase3d_get_localmesh
 
!OCL SERIAL
  subroutine meshbase3d_setgeometricinfo( lcmesh, coord_conv, calc_normal )
    use scale_mesh_base, only: meshbase_setgeometricinfo
    implicit none
    
    type(localmesh3d), intent(inout) :: lcmesh
    interface
      subroutine coord_conv( x, y, z, xX, xY, xZ, yX, yY, yZ, zX, zY, zZ, &
        vx, vy, vz, elem )
        import elementbase3d
        import rp
        type(elementbase3d), intent(in) :: elem
        real(rp), intent(out) :: x(elem%np), y(elem%np), z(elem%np)
        real(rp), intent(out) :: xx(elem%np), xy(elem%np), xz(elem%np)
        real(rp), intent(out) :: yx(elem%np), yy(elem%np), yz(elem%np)
        real(rp), intent(out) :: zx(elem%np), zy(elem%np), zz(elem%np)
        real(rp), intent(in) :: vx(elem%nv), vy(elem%nv), vz(elem%nv)        
      end subroutine coord_conv
      subroutine calc_normal( normal_fn, &
        Escale_f, fid_h, fid_v, elem )
        import elementbase3d
        import rp
        type(elementbase3d), intent(in) :: elem
        real(rp), intent(out) :: normal_fn(elem%nfptot,3)
        integer, intent(in) :: fid_h(elem%nfp_h,elem%nfaces_h)
        integer, intent(in) :: fid_v(elem%nfp_v,elem%nfaces_v)        
        real(rp), intent(in) :: escale_f(elem%nfptot,3,3)
      end subroutine calc_normal
    end interface
 
    class(elementbase3d), pointer :: refelem
    integer :: ke, ke2d
    integer :: f
    integer :: i, j
    integer :: d
    integer :: fmask(lcmesh%refelem%nfptot)
    integer :: fid_h(lcmesh%refelem3d%nfp_h,lcmesh%refelem3d%nfaces_h)
    integer :: fid_v(lcmesh%refelem3d%nfp_v,lcmesh%refelem3d%nfaces_v)
    real(rp) :: escale_f(lcmesh%refelem%nfptot,3,3)
 
    integer :: node_ids(lcmesh%refelem%nv)
    real(rp) :: vx(lcmesh%refelem%nv), vy(lcmesh%refelem%nv), vz(lcmesh%refelem%nv)    
    real(rp) :: xx(lcmesh%refelem%np), xy(lcmesh%refelem%np), xz(lcmesh%refelem%np)
    real(rp) :: yx(lcmesh%refelem%np), yy(lcmesh%refelem%np), yz(lcmesh%refelem%np)
    real(rp) :: zx(lcmesh%refelem%np), zy(lcmesh%refelem%np), zz(lcmesh%refelem%np)
    !-----------------------------------------------------------------------------
 
    refelem => lcmesh%refElem3D
 
    call meshbase_setgeometricinfo( lcmesh, 3 )
 
    allocate( lcmesh%zS(refelem%Np,lcmesh%Ne) )
    allocate( lcmesh%Sz(refelem%Np,lcmesh%Ne) )
    allocate( lcmesh%zlev(refelem%Np,lcmesh%Ne) )
    allocate( lcmesh%gam(refelem%Np,lcmesh%NeA) )
    allocate( lcmesh%GsqrtH(refelem%Nfp_v,lcmesh%Ne2D) )
    allocate( lcmesh%G_ij(refelem%Nfp_v,lcmesh%Ne2D,2,2) )
    allocate( lcmesh%GIJ (refelem%Nfp_v,lcmesh%Ne2D,2,2) )
    allocate( lcmesh%GI3 (refelem%Np,lcmesh%NeA,2) )
    allocate( lcmesh%lon2D(refelem%Nfp_v,lcmesh%Ne2D) )
    allocate( lcmesh%lat2D(refelem%Nfp_v,lcmesh%Ne2D) )
    !$acc enter data create( lcmesh%zS, lcmesh%Sz, lcmesh%zlev, lcmesh%gam, &
    !$acc   lcmesh%GsqrtH, lcmesh%G_ij, lcmesh%GIJ, lcmesh%GI3,             &
    !$acc   lcmesh%lon2D, lcmesh%lat2D )
 
    do f=1, refelem%Nfaces_h
    do i=1, refelem%Nfp_h
      fid_h(i,f) = i + (f-1)*refelem%Nfp_h
      fmask(fid_h(i,f)) = refelem%Fmask_h(i,f)
    end do
    end do
    do f=1, refelem%Nfaces_v
    do i=1, refelem%Nfp_v
      fid_v(i,f) = i + refelem%Nfaces_h*refelem%Nfp_h + (f-1)*refelem%Nfp_v
      fmask(fid_v(i,f)) = refelem%Fmask_v(i,f)
    end do
    end do
 
    !$omp parallel private( ke, node_ids, vx, vy, vz,    &
    !$omp xX, xY, xZ, yX, yY, yZ, zX, zY, zZ,            &
    !$omp i, j, Escale_f, d                              )
 
    !$omp do
    do ke=1, lcmesh%Ne
      node_ids(:) = lcmesh%EToV(ke,:)
      vx(:) = lcmesh%pos_ev(node_ids(:),1)
      vy(:) = lcmesh%pos_ev(node_ids(:),2)
      vz(:) = lcmesh%pos_ev(node_ids(:),3)
      call coord_conv( &
        lcmesh%pos_en(:,ke,1), lcmesh%pos_en(:,ke,2), lcmesh%pos_en(:,ke,3), & ! (out)
        xx, xy, xz, yx, yy, yz, zx, zy, zz,                                  & ! (out)
        vx, vy, vz, refelem )                                                  ! (in)
 
      lcmesh%J(:,ke) =   xx(:)*(yy(:)*zz(:) - zy(:)*yz) & 
                       - yx(:)*(xy(:)*zz(:) - zy(:)*xz) &
                       + zx(:)*(xy(:)*yz(:) - yy(:)*xz)
 
      lcmesh%Escale(:,ke,1,1) =   (yy(:)*zz(:) - zy(:)*yz(:))/lcmesh%J(:,ke)
      lcmesh%Escale(:,ke,1,2) = - (xy(:)*zz(:) - zy(:)*xz(:))/lcmesh%J(:,ke)
      lcmesh%Escale(:,ke,1,3) =   (xy(:)*yz(:) - yy(:)*xz(:))/lcmesh%J(:,ke)
      
      lcmesh%Escale(:,ke,2,1) = - (yx(:)*zz(:) - zx(:)*yz(:))/lcmesh%J(:,ke)
      lcmesh%Escale(:,ke,2,2) =   (xx(:)*zz(:) - zx(:)*xz(:))/lcmesh%J(:,ke)
      lcmesh%Escale(:,ke,2,3) = - (xx(:)*yz(:) - yx(:)*xz(:))/lcmesh%J(:,ke)
      
      lcmesh%Escale(:,ke,3,1) =   (yx(:)*zy(:) - zx(:)*yy(:))/lcmesh%J(:,ke)
      lcmesh%Escale(:,ke,3,2) = - (xx(:)*zy(:) - zx(:)*xy(:))/lcmesh%J(:,ke)
      lcmesh%Escale(:,ke,3,3) =   (xx(:)*yy(:) - yx(:)*xy(:))/lcmesh%J(:,ke)
      
      !* Face
 
      !
      !mesh%fx(:,n) = mesh%x(fmask(:),n)
      !mesh%fy(:,n) = mesh%y(fmask(:),n)
 
      ! Calculate normal vectors
      do j=1, 3
      do i=1, 3
        escale_f(:,i,j) = lcmesh%Escale(fmask(:),ke,i,j)
      end do
      end do
      call calc_normal( lcmesh%normal_fn(:,ke,:), & ! (out)
         escale_f, fid_h, fid_v, refelem )          ! (in)
 
      lcmesh%sJ(:,ke) = sqrt( &
        lcmesh%normal_fn(:,ke,1)**2 + lcmesh%normal_fn(:,ke,2)**2 + lcmesh%normal_fn(:,ke,3)**2 )
      do d=1, 3
        lcmesh%normal_fn(:,ke,d) = lcmesh%normal_fn(:,ke,d)/lcmesh%sJ(:,ke)
      end do
      lcmesh%sJ(:,ke) = lcmesh%sJ(:,ke)*lcmesh%J(fmask(:),ke)
 
      lcmesh%Fscale(:,ke) = lcmesh%sJ(:,ke)/lcmesh%J(fmask(:),ke)
      lcmesh%zlev(:,ke) = lcmesh%pos_en(:,ke,3)
    end do
    !$omp end do
    !$acc update device(lcmesh%pos_en, lcmesh%normal_fn, lcmesh%sJ, lcmesh%J, lcmesh%Escale, lcmesh%Fscale, lcmesh%zlev)
 
    !$omp workshare
    lcmesh%Gsqrt (:,:)   = 1.0_rp
    lcmesh%GsqrtH(:,:)   = 1.0_rp
    lcmesh%GIJ   (:,:,1,1) = 1.0_rp
    lcmesh%GIJ   (:,:,2,1) = 0.0_rp
    lcmesh%GIJ   (:,:,1,2) = 0.0_rp
    lcmesh%GIJ   (:,:,2,2) = 1.0_rp
    lcmesh%G_ij  (:,:,1,1) = 1.0_rp
    lcmesh%G_ij  (:,:,2,1) = 0.0_rp
    lcmesh%G_ij  (:,:,1,2) = 0.0_rp
    lcmesh%G_ij  (:,:,2,2) = 1.0_rp
    lcmesh%GI3   (:,:,1)   = 0.0_rp
    lcmesh%GI3   (:,:,2)   = 0.0_rp
    lcmesh%gam   (:,:)     = 1.0_rp
    !$omp end workshare
    !$acc update device(lcmesh%Gsqrt, lcmesh%GsqrtH, lcmesh%GIJ, lcmesh%G_ij, lcmesh%GI3, lcmesh%gam)
    
    !$omp end parallel
 
    return
  end subroutine meshbase3d_setgeometricinfo
  
end module scale_mesh_base3d