scale_mesh_rectdom2d.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Mesh / Rectangle 2D domain
!!
!! @par Description
!!      Manage mesh data of rectangle 2D domain for element-based methods
!!
!! @author Yuta Kawai, Team SCALE
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_mesh_rectdom2d
 
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
 
  use scale_mesh_base2d, only: &
    meshbase2d, meshbase2d_init, meshbase2d_final, &
    meshbase2d_setgeometricinfo
 
  use scale_localmesh_2d, only: &
    localmesh2d
  use scale_element_base, only: elementbase2d
  use scale_element_quadrilateral, only: quadrilateralelement
 
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  ! 
  !> Derived type to manage a rectangular 2D computational domain
  type, extends(meshbase2d), public :: meshrectdom2d
    integer :: negx    !< Number of elements in X direction (global)
    integer :: negy    !< Number of elements in Y direction (global)
    
    integer :: nprcx   !< Number of processes in X direction for domain decomposition
    integer :: nprcy   !< Number of processes in Y direction for domain decomposition
 
    real(rp), public :: xmin_gl !< Minimum X coordinate of the global domain
    real(rp), public :: xmax_gl !< Maximum X coordinate of the global domain
    real(rp), public :: ymin_gl !< Minimum Y coordinate of the global domain
    real(rp), public :: ymax_gl !< Maximum Y coordinate of the global domain
 
    integer, allocatable :: rcdomij2lcmeshid(:,:)
 
    logical :: isperiodicx  !< Flag whether the domain is periodic in X direction
    logical :: isperiodicy  !< Flag whether the domain is periodic in Y direction
  contains
    procedure :: init => meshrectdom2d_init
    procedure :: final => meshrectdom2d_final
    procedure :: generate => meshrectdom2d_generate
    procedure :: assigndomid => meshrectdom2d_assigndomid
  end type meshrectdom2d
 
  public :: meshrectdom2d_coord_conv
  public :: meshrectdom2d_setuplocaldom
  
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
 
contains
  !> Initialize an object to manage a 2D rectangular computational mesh
!OCL SERIAL
  subroutine meshrectdom2d_init(this,       &
    NeGX, NeGY,                             &
    dom_xmin, dom_xmax, dom_ymin, dom_ymax, &
    isPeriodicX, isPeriodicY,               &
    refElem, NLocalMeshPerPrc,              &
    NprcX, NprcY,                           &
    nproc, myrank                           )
    
    implicit none
 
    class(meshrectdom2d), intent(inout) :: this
    integer, intent(in) :: NeGX          !< Number of elements in X direction (global)
    integer, intent(in) :: NeGY          !< Number of elements in Y direction (global)
    real(RP), intent(in) :: dom_xmin     !< Minimum X coordinate of the global domain
    real(RP), intent(in) :: dom_xmax     !< Maximum X coordinate of the global domain
    real(RP), intent(in) :: dom_ymin     !< Minimum Y coordinate of the global domain
    real(RP), intent(in) :: dom_ymax     !< Maximum Y coordinate of the global domain
    logical, intent(in) :: isPeriodicX   !< Flag whether the domain is periodic in X direction
    logical, intent(in) :: isPeriodicY   !< Flag whether the domain is periodic in Y direction
    type(quadrilateralelement), intent(in), target :: refElem !< Reference element for the 2D mesh
    integer, intent(in) :: NLocalMeshPerPrc !< Number of local meshes managed by each process
    integer, intent(in) :: NprcX            !< Number of processes in X direction for domain decomposition
    integer, intent(in) :: NprcY            !< Number of processes in Y direction for domain decomposition
    integer, intent(in), optional :: nproc  !< Total number of processes (if not provided, it will be determined from the parallel environment)
    integer, intent(in), optional :: myrank !< Rank of the current process (if not provided, it will be determined from the parallel environment)
 
    !-----------------------------------------------------------------------------
    
    this%NeGX = negx
    this%NeGY = negy
 
    this%xmin_gl       = dom_xmin
    this%xmax_gl       = dom_xmax
    this%ymin_gl       = dom_ymin
    this%ymax_gl       = dom_ymax
    this%dom_vol = (this%xmax_gl - this%xmin_gl) * (this%ymax_gl - this%ymin_gl)
 
    this%isPeriodicX = isperiodicx
    this%isPeriodicY = isperiodicy
 
    this%NprcX = nprcx
    this%NprcY = nprcy
 
    call meshbase2d_init( this, refelem, nlocalmeshperprc, &
      nproc, myrank )
 
    return
  end subroutine meshrectdom2d_init
 
  !> Finalize an object to manage a 2D rectangular computational mesh
!OCL SERIAL
  subroutine meshrectdom2d_final( this )
    implicit none
    class(meshrectdom2d), intent(inout) :: this
 
    integer :: n
    !-----------------------------------------------------------------------------
  
    if (this%isGenerated) then
      if ( allocated(this%rcdomIJ2LCMeshID) ) then
        !$acc exit data delete(this%rcdomIJ2LCMeshID)
        deallocate( this%rcdomIJ2LCMeshID )
      end if
    end if
 
    call meshbase2d_final( this )
 
    return
  end subroutine meshrectdom2d_final
  
  !> Generate meshes for the rectangular domain
!OCL SERIAL
  subroutine meshrectdom2d_generate( this )
    implicit none
    class(meshrectdom2d), intent(inout), target :: this
            
    integer :: n
    integer :: p
    type(localmesh2d), pointer :: mesh
 
    integer :: tileID_table(this%LOCAL_MESH_NUM, this%PRC_NUM)
    integer :: panelID_table(this%LOCAL_MESH_NUM*this%PRC_NUM)
    integer :: pi_table(this%LOCAL_MESH_NUM*this%PRC_NUM)
    integer :: pj_table(this%LOCAL_MESH_NUM*this%PRC_NUM)
 
    integer :: TILE_NUM_PER_PANEL
    real(RP) :: delx, dely
    integer :: tileID
    !-----------------------------------------------------------------------------
 
    tile_num_per_panel = this%LOCAL_MESH_NUM_global / 1
    
    
    !--- Construct the connectivity of patches  (only master node)
 
    call this%AssignDomID(  &
      tileid_table, panelid_table, & ! (out)
      pi_table, pj_table )           ! (out)
  
    !--- Setup local meshes managed by my process
 
    do n=1, this%LOCAL_MESH_NUM
      mesh => this%lcmesh_list(n)
      tileid = tileid_table(n, mesh%PRC_myrank+1)
 
      call meshrectdom2d_setuplocaldom( mesh, &
        tileid,  panelid_table(tileid),                             &
        pi_table(tileid), pj_table(tileid), this%NprcX, this%NprcY, &
        this%xmin_gl, this%xmax_gl, this%ymin_gl, this%ymax_gl,     &
        this%NeGX/this%NprcX, this%NeGY/this%NprcY                  )
 
      !---
      ! write(*,*) "** my_rank=", mesh%PRC_myrank
      ! write(*,*) " tileID:", mesh%tileID
      ! write(*,*) " pnlID:", mesh%panelID, "-- i,j (within a panel)=", pi_table(tileID), pj_table(tileID)
      ! write(*,*) " local mesh:", n, "( total", this%LOCAL_MESH_NUM, ")"
      ! write(*,*) " panel_connect:", this%tilePanelID_globalMap(:,mesh%tileID)
      ! write(*,*) " tile_connect:", this%tileID_globalMap(:,mesh%tileID)
      ! write(*,*) " face_connect:", this%tileFaceID_globalMap(:,mesh%tileID)
      ! write(*,*) " domain size"
      ! write(*,*) "   NeX, NeY:", mesh%NeX, mesh%NeY
      ! write(*,*) "   [X], [Y]:",  mesh%xmin, mesh%xmax, ":", mesh%ymin, mesh%ymax
    end do
 
    this%isGenerated = .true.
 
    return
  end subroutine meshrectdom2d_generate
 
  !- private ------------------------------------------------------
 
  subroutine meshrectdom2d_setuplocaldom( lcmesh, &
    tileID, panelID,                         &
    i, j, NprcX, NprcY,                      &
    dom_xmin, dom_xmax, dom_ymin, dom_ymax,  &
    NeX, NeY )
 
    use scale_meshutil_2d, only: &
      meshutil2d_genconnectivity,   &
      meshutil2d_genrectdomain,     &
      meshutil2d_buildinteriormap,  &
      meshutil2d_genpatchboundarymap
 
    use scale_localmesh_base, only: bctype_interior
    implicit none
      
    type(localmesh2d), intent(inout) :: lcmesh
    integer, intent(in) :: tileid
    integer, intent(in) :: panelid
    integer, intent(in) :: i, j
    integer, intent(in) :: nprcx, nprcy
    real(rp), intent(in) :: dom_xmin, dom_xmax
    real(rp), intent(in) :: dom_ymin, dom_ymax
    integer, intent(in) ::nex, ney
    
    class(elementbase2d), pointer :: elem
    real(rp) :: delx, dely
    
    !-----------------------------------------------------------------------------
 
    elem => lcmesh%refElem2D
    lcmesh%tileID = tileid
    lcmesh%panelID = panelid
    !$acc update device(lcmesh%tileID, lcmesh%panelID)
    
    !--
 
    lcmesh%Ne  = nex * ney
    lcmesh%Nv  = (nex + 1)*(ney + 1)
    lcmesh%NeS = 1
    lcmesh%NeE = lcmesh%Ne
    lcmesh%NeA = lcmesh%Ne + 2*(nex + ney)
    !$acc update device(lcmesh%Ne, lcmesh%Nv, lcmesh%NeS, lcmesh%NeE, lcmesh%NeA)
 
    lcmesh%NeX = nex
    lcmesh%NeY = ney
    !$acc update device(lcmesh%NeX, lcmesh%NeY)
 
    delx = (dom_xmax - dom_xmin)/dble(nprcx)
    dely = (dom_ymax - dom_ymin)/dble(nprcy)
 
    lcmesh%xmin = dom_xmin + (i-1)*delx
    lcmesh%xmax = dom_xmin +  i   *delx
    lcmesh%ymin = dom_ymin + (j-1)*dely
    lcmesh%ymax = dom_ymin +  j   *dely
    !$acc update device(lcmesh%xmin, lcmesh%xmax, lcmesh%ymin, lcmesh%ymax)
 
    allocate( lcmesh%pos_ev(lcmesh%Nv,2) )
    allocate( lcmesh%EToV(lcmesh%Ne,elem%Nv) )
    allocate( lcmesh%EToE(lcmesh%Ne,elem%Nfaces) )
    allocate( lcmesh%EToF(lcmesh%Ne,elem%Nfaces) )
    allocate( lcmesh%BCType(lcmesh%refElem%Nfaces,lcmesh%Ne) )
    allocate( lcmesh%VMapM(elem%NfpTot, lcmesh%Ne) )
    allocate( lcmesh%VMapP(elem%NfpTot, lcmesh%Ne) )
    allocate( lcmesh%MapM(elem%NfpTot, lcmesh%Ne) )
    allocate( lcmesh%MapP(elem%NfpTot, lcmesh%Ne) )
    !$acc enter data create(lcmesh%pos_ev, lcmesh%EToV, lcmesh%EToE, lcmesh%EToF, lcmesh%BCType, &
    !$acc   lcmesh%VMapM, lcmesh%VMapP, lcmesh%MapM, lcmesh%MapP)
 
    lcmesh%BCType(:,:) = bctype_interior
    !$acc update device(lcmesh%BCType)
 
    !----
 
    call meshutil2d_genrectdomain( lcmesh%pos_ev, lcmesh%EToV,   & ! (out)
      lcmesh%NeX, lcmesh%xmin, lcmesh%xmax,                      & ! (in)
      lcmesh%NeY, lcmesh%ymin, lcmesh%ymax )                       ! (in)
    !$acc update device(lcmesh%pos_ev, lcmesh%EToV)
    
    !---
    call meshbase2d_setgeometricinfo( lcmesh, meshrectdom2d_coord_conv, meshrectdom2d_calc_normal )
 
    !---
 
    call meshutil2d_genconnectivity( lcmesh%EToE, lcmesh%EToF, & ! (out)
      lcmesh%EToV, lcmesh%Ne, elem%Nfaces )                      ! (in)
    !$acc update device(lcmesh%EToE, lcmesh%EToF)
 
    !---
    call meshutil2d_buildinteriormap( lcmesh%VmapM, lcmesh%VMapP, lcmesh%MapM, lcmesh%MapP, &
      lcmesh%pos_en, lcmesh%pos_ev, lcmesh%EToE, lcmesh%EtoF, lcmesh%EtoV,                  &
      elem%Fmask, lcmesh%Ne, elem%Np, elem%Nfp, elem%Nfaces, lcmesh%Nv )
 
    call meshutil2d_genpatchboundarymap( lcmesh%VMapB, lcmesh%MapB, lcmesh%VMapP, &
      lcmesh%pos_en, lcmesh%xmin, lcmesh%xmax, lcmesh%ymin, lcmesh%ymax,          &
      elem%Fmask, lcmesh%Ne, elem%Np, elem%Nfp, elem%Nfaces, lcmesh%Nv)
    !$acc update device(lcmesh%VMapM, lcmesh%VMapP, lcmesh%MapM, lcmesh%MapP)
    !$acc enter data copyin(lcmesh%VMapB, lcmesh%MapB)
    
    return
  end subroutine meshrectdom2d_setuplocaldom
 
  subroutine meshrectdom2d_assigndomid( this, &
    tileID_table, panelID_table,              &
    pi_table, pj_table )
  
    use scale_meshutil_2d, only: &       
      meshutil2d_buildglobalmap
    use scale_io
    implicit none
 
    class(meshrectdom2d), target, intent(inout) :: this    
    integer, intent(out) :: tileid_table(this%local_mesh_num, this%prc_num)
    integer, intent(out) :: panelid_table(this%local_mesh_num*this%prc_num)
    integer, intent(out) :: pi_table(this%local_mesh_num*this%prc_num)
    integer, intent(out) :: pj_table(this%local_mesh_num*this%prc_num)
    
    integer :: n
    integer :: prc
    integer :: tileid
    integer :: is_lc, js_lc
    integer :: ilc_count, jlc_count
    integer :: ilc, jlc
    
    type(localmesh2d), pointer :: lcmesh    
    !-----------------------------------------------------------------------------
 
    call meshutil2d_buildglobalmap( &
      panelid_table, pi_table, pj_table,                                            & ! (out)
      this%tileID_globalMap, this%tileFaceID_globalMap, this%tilePanelID_globalMap, & ! (out)
      this%LOCAL_MESH_NUM_global, this%isPeriodicX, this%isPeriodicY,               & ! (in)
      this%NprcX, this%NprcY )                                                        ! (in)
    !$acc update device(this%tileID_globalMap, this%tileFaceID_globalMap, this%tilePanelID_globalMap)
 
    !----
    
    do prc=1, this%PRC_NUM
    do n=1, this%LOCAL_MESH_NUM
      tileid = n + (prc-1)*this%LOCAL_MESH_NUM
      lcmesh => this%lcmesh_list(n)      
      !-
      tileid_table(n,prc)                 = tileid
      this%tileID_global2localMap(tileid) = n
      this%PRCRank_globalMap(tileid)      = prc - 1
 
      !-
      if ( this%PRCRank_globalMap(tileid) == lcmesh%PRC_myrank ) then
        if (n==1) then
          is_lc = pi_table(tileid); ilc_count = 1
          js_lc = pj_table(tileid); jlc_count = 1
        end if
        if(is_lc < pi_table(tileid)) ilc_count = ilc_count + 1
        if(js_lc < pj_table(tileid)) jlc_count = jlc_count + 1
      end if 
    end do
    end do
    !$acc update device(this%tileID_global2localMap, this%PRCRank_globalMap)
 
    allocate( this%rcdomIJ2LCMeshID(ilc_count,jlc_count) )
    do jlc=1, jlc_count
    do ilc=1, ilc_count
      this%rcdomIJ2LCMeshID(ilc,jlc) = ilc + (jlc - 1)*ilc_count
    end do
    end do
    !$acc enter data copyin(this%rcdomIJ2LCMeshID)
 
    return
  end subroutine meshrectdom2d_assigndomid
  
  subroutine meshrectdom2d_coord_conv( x, y, xr, xs, yr, ys, &
    vx, vy, elem )
    implicit none
 
    type(elementbase2d), intent(in) :: elem
    real(rp), intent(out) :: x(elem%np), y(elem%np)
    real(rp), intent(out) :: xr(elem%np), xs(elem%np), yr(elem%np), ys(elem%np)
    real(rp), intent(in) :: vx(elem%nv), vy(elem%nv)
 
    !-------------------------------------------------
    
    x(:) = vx(1) + 0.5_rp*(elem%x1(:) + 1.0_rp)*(vx(2) - vx(1))
    y(:) = vy(1) + 0.5_rp*(elem%x2(:) + 1.0_rp)*(vy(3) - vy(1))
 
    xr(:) = 0.5_rp*(vx(2) - vx(1)) !matmul(refElem%Dx1,mesh%x1(:,n))
    xs(:) = 0.0_rp                 !matmul(refElem%Dx2,mesh%x1(:,n))
    yr(:) = 0.0_rp                 !matmul(refElem%Dx1,mesh%x2(:,n))
    ys(:) = 0.5_rp*(vy(3) - vy(1)) !matmul(refElem%Dx2,mesh%x2(:,n))
 
    return
  end subroutine meshrectdom2d_coord_conv  
 
  subroutine meshrectdom2d_calc_normal( normal_fn, &
    Escale_f, fid, elem )
    implicit none
 
    type(elementbase2d), intent(in) :: elem
    real(rp), intent(out) :: normal_fn(elem%nfptot,2)
    integer, intent(in) :: fid(elem%nfp,elem%nfaces)
    real(rp), intent(in) :: escale_f(elem%nfptot,2,2)
 
    integer :: d
    !-------------------------------------------------
 
    do d=1, 2
      normal_fn(fid(:,1),d) = - escale_f(fid(:,1),2,d)
      normal_fn(fid(:,2),d) = + escale_f(fid(:,2),1,d)
      normal_fn(fid(:,3),d) = + escale_f(fid(:,3),2,d)
      normal_fn(fid(:,4),d) = - escale_f(fid(:,4),1,d)  
    end do
 
    return
  end subroutine meshrectdom2d_calc_normal 
 
end module scale_mesh_rectdom2d