scale_mesh_cubedspheredom2d.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Mesh / Cubed-sphere 2D domain
!!
!! @par Description
!!      Manage mesh data of cubed-sphere 2D domain for element-based methods
!!
!! @author Yuta Kawai, Team SCALE
!<
#include "scaleFElib.h"
module scale_mesh_cubedspheredom2d
 
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_io
  use scale_precision
 
  use scale_mesh_base2d, only: &
    meshbase2d, meshbase2d_init, meshbase2d_final,    &
    meshbase2d_setgeometricinfo,                      &
    meshbase2d_dimtype_num,                           &
    meshbase2d_dimtypeid_x, meshbase2d_dimtypeid_y,   &
    meshbase2d_dimtypeid_xy, meshbase2d_dimtypeid_xyt
 
  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 cubed-sphere 2D computational domain
  type, extends(meshbase2d), public :: meshcubedspheredom2d
    integer :: negx !< Number of elements in X-direction in a panel of cubed-sphere mesh
    integer :: negy !< Number of elements in Y-direction in a panel of cubed-sphere mesh
    
    real(rp), public :: xmin_gl !< Minimum X-coordinate in a panel of cubed-sphere mesh
    real(rp), public :: xmax_gl !< Maximum X-coordinate in a panel of cubed-sphere mesh
    real(rp), public :: ymin_gl !< Minimum Y-coordinate in a panel of cubed-sphere mesh
    real(rp), public :: ymax_gl !< Maximum Y-coordinate in a panel of cubed-sphere mesh
 
    integer, allocatable :: rcdomijp2lcmeshid(:,:,:) !< Mapping from (i,j,panel) to local mesh ID
 
    real(rp) :: rplanet     !< Planetary radius
    integer :: halosize_1d  !< Halo size in 1D direction
  contains
    procedure :: init => meshcubedspheredom2d_init
    procedure :: final => meshcubedspheredom2d_final
    procedure :: generate => meshcubedspheredom2d_generate
    procedure :: assigndomid => meshcubedspheredom2d_assigndomid
  end type meshcubedspheredom2d
 
  public :: meshcubedspheredom2d_check_division_params
  public :: meshcubedspheredom2d_setuplocaldom
  
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  
  private :: meshcubedspheredom2d_calc_normal
  private :: meshcubedspheredom2d_coord_conv
  private :: fill_halo_metric
 
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
 
contains
  subroutine meshcubedspheredom2d_init(this, &
    NeGX, NeGY, RPlanet,                     &
    refElem, NLocalMeshPerPrc,               &
    nproc, myrank                            )
    
    use scale_const, only: &
      pi => const_pi
    implicit none
 
    class(meshcubedspheredom2d), intent(inout) :: this
    integer, intent(in) :: NeGX
    integer, intent(in) :: NeGY
    real(RP), intent(in) :: RPlanet
    type(quadrilateralelement), intent(in), target :: refElem
    integer, intent(in) :: NLocalMeshPerPrc
    integer, intent(in), optional :: nproc
    integer, intent(in), optional :: myrank
    !-----------------------------------------------------------------------------
 
    this%NeGX = negx
    this%NeGY = negy
 
    this%xmin_gl = - 0.25_rp * pi 
    this%xmax_gl = + 0.25_rp * pi 
    this%ymin_gl = - 0.25_rp * pi 
    this%ymax_gl = + 0.25_rp * pi 
    this%RPlanet = rplanet
    this%dom_vol = 4.0_rp * pi * rplanet**2
 
 
    call meshbase2d_init( this, refelem, nlocalmeshperprc, &
      nproc, myrank )
 
    call this%SetDimInfo( meshbase2d_dimtypeid_x, "x", "1", "X-coordinate" )
    call this%SetDimInfo( meshbase2d_dimtypeid_y, "y", "1", "Y-coordinate" )
    call this%SetDimInfo( meshbase2d_dimtypeid_xy, "xy", "1", "XY-coordinate" )
    call this%SetDimInfo( meshbase2d_dimtypeid_xyt, "xyt", "1", "XY-coordinate" )
  
    return
  end subroutine meshcubedspheredom2d_init
 
  subroutine meshcubedspheredom2d_final( this )
    implicit none
    class(meshcubedspheredom2d), intent(inout) :: this
    !-----------------------------------------------------------------------------
  
    if (this%isGenerated) then
      if ( allocated(this%rcdomIJP2LCMeshID) ) then
        !$acc exit data delete(this%rcdomIJP2LCMeshID)
        deallocate( this%rcdomIJP2LCMeshID )
      end if
    end if
 
    call meshbase2d_final( this )
 
    return
  end subroutine meshcubedspheredom2d_final
 
  subroutine meshcubedspheredom2d_generate( this )
    implicit none
 
    class(meshcubedspheredom2d), intent(inout), target :: this
 
    integer :: n
    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 :: NprcX_lc, NprcY_lc
    integer :: tileID
   
    !-----------------------------------------------------------------------------
 
    call meshcubedspheredom2d_check_division_params( &
      nprcx_lc, nprcy_lc,                            &
      this%PRC_NUM, this%LOCAL_MESH_NUM_global,      &
      .true. )
 
    !--- Construct the connectivity of patches  (only master node)
 
    call this%AssignDomID( & 
      nprcx_lc, nprcy_lc,          & ! (in)
      tileid_table, panelid_table, & ! (out)
      pi_table, pj_table )           ! (out)
 
    do n=1, this%LOCAL_MESH_NUM
      mesh => this%lcmesh_list(n)
      tileid = tileid_table(n, mesh%PRC_myrank+1)
 
      call meshcubedspheredom2d_setuplocaldom( mesh, &
        tileid,  panelid_table(tileid),                                       &
        pi_table(tileid), pj_table(tileid), nprcx_lc, nprcy_lc,               &
        this%xmin_gl, this%xmax_gl, this%ymin_gl, this%ymax_gl, this%RPlanet, &
        this%NeGX/nprcx_lc, this%NeGY/nprcy_lc )
 
      !---
      ! 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 meshcubedspheredom2d_generate
 
  subroutine meshcubedspheredom2d_check_division_params( &
    NprcX_lc, NprcY_lc,                                  &
    PRC_NUM, LOCAL_MESH_NUM_global,                      &
    call_prc_abort )
    
    use scale_prc, only: prc_abort
    implicit none
 
    integer, intent(in) :: prc_num
    integer, intent(in) :: local_mesh_num_global
    integer, intent(out) :: nprcx_lc
    integer, intent(out) :: nprcy_lc
    logical, intent(in), optional :: call_prc_abort
 
    integer :: tile_num_per_panel
    logical :: call_prc_abort_
    !-----------------------------------------------------------------------------
 
    if (present(call_prc_abort)) then
      call_prc_abort_ = call_prc_abort
    else
      call_prc_abort_ = .false.
    end if
 
    if ( mod(local_mesh_num_global, 6) /= 0 ) then
      log_error("MeshCubedSphereDom2D_division_params",*) "The total number of local mesh must be a multiple of 6. Check!"
      if (call_prc_abort_) call prc_abort
    end if
 
    tile_num_per_panel = local_mesh_num_global / 6
 
    if ( prc_num <= 6 ) then
      if ( ( prc_num == 1 ) .or. &
          ( prc_num <= 6 .and. (mod(prc_num,2)==0 .or. mod(prc_num,3)==0)) ) then
        nprcx_lc = 1; nprcy_lc = 1
      else
        log_error("MeshCubedSphereDom2D_division_params",*) "The number of proceses is inappropriate. Check!"
        if (call_prc_abort_) call prc_abort
      end if
    else
      if ( mod(prc_num,6) == 0 ) then
        nprcx_lc = int(sqrt(dble(tile_num_per_panel)))
        nprcy_lc = tile_num_per_panel / nprcx_lc
        if ( nprcx_lc /= nprcy_lc ) then
          log_error("MeshCubedSphereDom2D_division_params",*) "The number of proceses is inappropriate. Check!"
          if (call_prc_abort_) call prc_abort  
        end if
      else
        log_error("MeshCubedSphereDom2D_division_params",*) "The number of proceses is inappropriate. Check!"
        if (call_prc_abort_) call prc_abort        
      end if
    end if
 
    return
  end subroutine meshcubedspheredom2d_check_division_params
 
  subroutine meshcubedspheredom2d_setuplocaldom( lcmesh,   &
    tileID, panelID,                                       &
    i, j, NprcX, NprcY,                                    &
    dom_xmin, dom_xmax, dom_ymin, dom_ymax, planet_radius, &
    NeX, NeY )
 
    use scale_meshutil_cubedsphere2d, only: &
      meshutilcubedsphere2d_genconnectivity,   &
      meshutilcubedsphere2d_genrectdomain,     &
      meshutilcubedsphere2d_buildinteriormap,  &
      meshutilcubedsphere2d_genpatchboundarymap
 
    use scale_cubedsphere_coord_cnv, only: &
      cubedspherecoordcnv_getmetric,      &
      cubedspherecoordcnv_cs2lonlatpos
    
    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
    real(rp), intent(in) :: planet_radius
    integer, intent(in) :: nex, ney
    
    class(elementbase2d), pointer :: elem
    real(rp) :: delx, dely
    integer :: ke
 
    real(rp), allocatable :: gam(:,:)
    !-----------------------------------------------------------------------------
 
    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 meshutilcubedsphere2d_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, meshcubedspheredom2d_coord_conv, meshcubedspheredom2d_calc_normal )
    
    call cubedspherecoordcnv_getmetric( &
      lcmesh%pos_en(:,:,1), lcmesh%pos_en(:,:,2), elem%Np * lcmesh%Ne, planet_radius, & ! (in)
      lcmesh%G_ij, lcmesh%GIJ, lcmesh%Gsqrt(:,lcmesh%NeS:lcmesh%NeE)                  ) ! (out)
    !$acc update host(lcmesh%G_ij, lcmesh%GIJ)
 
    allocate( gam(elem%Np,lcmesh%Ne) )
    gam(:,:) = 1.0_rp
    !$acc data copyin(gam)
 
    call cubedspherecoordcnv_cs2lonlatpos( &
      lcmesh%panelID, lcmesh%pos_en(:,:,1), lcmesh%pos_en(:,:,2), gam(:,:), & ! (in)
      lcmesh%Ne * lcmesh%refElem2D%Np,                                      & ! (in)
      lcmesh%lon(:,:), lcmesh%lat(:,:)                                      ) ! (out)
    !$acc update host(lcmesh%lon, lcmesh%lat)
 
    !$acc end data
    
    !---
 
    call meshutilcubedsphere2d_genconnectivity( lcmesh%EToE, lcmesh%EToF, & ! (out)
      lcmesh%EToV, lcmesh%Ne, elem%Nfaces )                                 ! (in)
    !$acc update device(lcmesh%EToE, lcmesh%EToF)
 
    !---
    call meshutilcubedsphere2d_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 meshutilcubedsphere2d_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)
    
    !--
    call fill_halo_metric( lcmesh%Gsqrt, &
      lcmesh%VMapM, lcmesh%VMapP, lcmesh, elem )
    !$acc update host(lcmesh%Gsqrt)
 
    return
  end subroutine meshcubedspheredom2d_setuplocaldom
 
  !- private ------------------------------
 
  subroutine meshcubedspheredom2d_assigndomid( this, &
    NprcX_lc, NprcY_lc,                               &
    tileID_table, panelID_table,                      &
    pi_table, pj_table )
  
    use scale_meshutil_cubedsphere2d, only: &       
      meshutilcubedsphere2d_buildglobalmap
    
    implicit none
 
    class(meshcubedspheredom2d), target, intent(inout) :: this    
    integer, intent(in) :: nprcx_lc
    integer, intent(in) :: nprcy_lc
    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, ps_lc
    integer :: ilc_count, jlc_count, plc_count
    integer :: ilc, jlc, plc
    integer :: npanel_lc
    
    type(localmesh2d), pointer :: lcmesh
    !-----------------------------------------------------------------------------
 
    call meshutilcubedsphere2d_buildglobalmap( &
      panelid_table, pi_table, pj_table,                                            & ! (out)
      this%tileID_globalMap, this%tileFaceID_globalMap, this%tilePanelID_globalMap, & ! (out)
      this%LOCAL_MESH_NUM_global )                                                    ! (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
          ps_lc = panelid_table(tileid); plc_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
        if(ps_lc < panelid_table(tileid)) plc_count = plc_count + 1
      end if 
    end do
    end do
    !$acc update device(this%tileID_global2localMap, this%PRCRank_globalMap)
 
    allocate( this%rcdomIJP2LCMeshID(ilc_count,jlc_count,plc_count) )
    do plc=1, plc_count
    do jlc=1, jlc_count
    do ilc=1, ilc_count
      this%rcdomIJP2LCMeshID(ilc,jlc,plc) = ilc + (jlc - 1)*ilc_count + (plc-1)*ilc_count*jlc_count
    end do
    end do
    end do
    !$acc enter data copyin(this%rcdomIJP2LCMeshID)
    
    return
  end subroutine meshcubedspheredom2d_assigndomid
 
  subroutine meshcubedspheredom2d_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 meshcubedspheredom2d_coord_conv  
 
  subroutine meshcubedspheredom2d_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 meshcubedspheredom2d_calc_normal 
 
  !--
 
!OCL SERIAL
  subroutine fill_halo_metric( Gsqrt, vmapM, vmapP, lmesh, elem )
    implicit none
    class(localmesh2d), intent(in) :: lmesh
    class(elementbase2d), intent(in) :: elem
    integer, intent(in) :: vmapm(elem%nfptot*lmesh%ne)
    integer, intent(in) :: vmapp(elem%nfptot*lmesh%ne)
    real(rp), intent(inout) :: gsqrt(elem%np*lmesh%nea)
 
    integer :: i, im, ip
    integer :: npxne
    !------------------------------------------------
 
    npxne = elem%Np * lmesh%Ne
    
    !$omp parallel do private(i, iM, iP)
    !$acc parallel loop present(Gsqrt, vmapM, vmapP)
    do i=1, elem%NfpTot*lmesh%Ne
      im = vmapm(i); ip = vmapp(i)
      if ( ip > npxne ) then
        gsqrt(ip) = gsqrt(im)
      end if
    end do  
    return
  end subroutine fill_halo_metric
 
end module scale_mesh_cubedspheredom2d