scale_mesh_base1d.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Mesh / Base 1D
!!
!! @par Description
!!      Base module to manage 1D meshes for element-based methods
!!
!! @author Yuta Kawai, Team SCALE
!<
#include "scaleFElib.h"
module scale_mesh_base1d
 
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
 
  use scale_localmesh_1d, only: &
    localmesh1d, localmesh1d_init, localmesh1d_final
  
  use scale_mesh_base, only: meshbase, &
    meshbase_init, meshbase_final, &
    meshbase_setgeometricinfo
  
  use scale_element_base, only: elementbase1d
  use scale_element_line, only: lineelement
 
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  ! 
  !> Derived type to manage a computational mesh (base type for 1D domain)
  type, abstract, public, extends(meshbase) :: meshbase1d
    type(localmesh1d), allocatable :: lcmesh_list(:) !< Array of local meshes in each MPI process
    class(elementbase1d), pointer :: refelem1d       !< Pointer to an object with a reference element for 1D
    
    integer, public :: neg
    integer, public :: nprc
 
    real(rp), public :: xmin_gl !< Minimum x-coordinate of the global mesh
    real(rp), public :: xmax_gl !< Maximum x-coordinate of the global mesh
 
    real(rp), public, allocatable :: fx(:)
  contains
    procedure(meshbase1d_generate), deferred :: generate 
    procedure :: getlocalmesh => meshbase1d_get_localmesh
  end type meshbase1d
 
  interface 
    subroutine meshbase1d_generate(this)
      import meshbase1d
      class(meshbase1d), intent(inout), target :: this
    end subroutine meshbase1d_generate
  end interface
 
  public :: meshbase1d_init, meshbase1d_final
  public :: meshbase1d_setgeometricinfo, meshbase1d_assigndomid, meshbase1d_setuplocaldom
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  integer, public :: meshbase1d_dimtype_num   = 2
  integer, public :: meshbase1d_dimtypeid_x   = 1
  integer, public :: meshbase1d_dimtypeid_xt  = 2
  
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
 
contains
  !> Initialize an object to manage a 1D computational mesh
!OCL SERIAL
  subroutine meshbase1d_init( this,         &
    NeG,                                    &
    dom_xmin, dom_xmax,                     &
    refElem, NLocalMeshPerPrc,              &
    nprocs, myrank, FX                      )
    
    implicit none
 
    class(meshbase1d), intent(inout) :: this
    integer, intent(in) :: neg
    real(rp), intent(in) :: dom_xmin
    real(rp), intent(in) :: dom_xmax   
    class(elementbase1d), intent(in), target :: refelem
    integer, intent(in) :: nlocalmeshperprc
    integer, intent(in), optional :: nprocs
    integer, intent(in), optional :: myrank
    real(rp), intent(in), optional :: fx(neg+1)
 
    integer :: n
    integer :: k
    real(rp) :: dx
    !-----------------------------------------------------------------------------
    
    this%NeG = neg
 
    this%xmin_gl       = dom_xmin
    this%xmax_gl       = dom_xmax
 
    !- Fx
    allocate( this%FX(neg+1) )
    if ( present(fx) ) then
      this%FX(:) = fx(:)
    else
      this%FX(1    ) = dom_xmin
      this%FX(neg+1) = dom_xmax
      dx = (dom_xmax - dom_xmin) / dble(neg)
      do k=2, neg
        this%FX(k) = this%FX(k-1) + dx
      end do
    end if
    !$acc enter data create(this%FX)
 
    this%refElem1D => refelem
    call meshbase_init( this,          &
      meshbase1d_dimtype_num, refelem, &
      nlocalmeshperprc, 2,             &
      nprocs                           )
 
    this%Nprc = this%PRC_NUM
         
    allocate( this%lcmesh_list(this%LOCAL_MESH_NUM) )
    do n=1, this%LOCAL_MESH_NUM
      call localmesh1d_init( this%lcmesh_list(n), n, refelem, myrank )
    end do
 
    call this%SetDimInfo( meshbase1d_dimtypeid_x, "x", "m", "X-coordinate" )
    call this%SetDimInfo( meshbase1d_dimtypeid_xt, "xt", "m", "X-coordinate" )
 
    return
  end subroutine meshbase1d_init
 
  !> Finalize an object to manage a 1D computational mesh
!OCL SERIAL
  subroutine meshbase1d_final( this )
    implicit none    
    class(meshbase1d), intent(inout) :: this
 
    integer :: n
    !-----------------------------------------------------------------------------
  
    if ( allocated ( this%lcmesh_list ) ) then 
      do n=1, this%LOCAL_MESH_NUM
        call localmesh1d_final( this%lcmesh_list(n), this%isGenerated )
      end do
  
      deallocate( this%lcmesh_list )
    end if
    
    call meshbase_final(this)
 
    return
  end subroutine meshbase1d_final
  
!OCL SERIAL
  subroutine meshbase1d_get_localmesh( this, id, ptr_lcmesh )
    use scale_localmesh_base, only: localmeshbase
    implicit none
 
    class(meshbase1d), target, intent(in) :: this
    integer, intent(in) :: id
    class(localmeshbase), pointer, intent(out) :: ptr_lcmesh
    !-------------------------------------------------------------
 
    ptr_lcmesh => this%lcmesh_list(id)
    return
  end subroutine meshbase1d_get_localmesh
 
!OCL SERIAL
  subroutine meshbase1d_setgeometricinfo( lcmesh )
    implicit none
    
    type(localmesh1d), intent(inout) :: lcmesh
 
    class(elementbase1d), pointer :: refelem
    integer :: ke
    integer :: f
    integer :: i
 
    integer :: node_ids(lcmesh%refelem%nv)
    real(rp) :: vx(lcmesh%refelem%nv)
    real(rp) :: xr(lcmesh%refelem%np)
    real(dp) :: escale(1,1,lcmesh%refelem%np)
    integer :: fmask(lcmesh%refelem%nfptot)
    integer :: fid(lcmesh%refelem1d%nfp,lcmesh%refelem1d%nfaces)
 
  !-----------------------------------------------------------------------------
 
    call meshbase_setgeometricinfo(lcmesh, 1)
    refelem => lcmesh%refElem1D
    
    fmask(:) = reshape(refelem%Fmask, shape(fmask))
    do f=1, refelem%Nfaces
    do i=1, refelem%Nfp
       fid(i,f) = i + (f-1)*refelem%Nfp
    end do
    end do
  
    do ke=1, lcmesh%Ne
      node_ids(:) = lcmesh%EToV(ke,:)
      vx(:) = lcmesh%pos_ev(node_ids(:),1)
      lcmesh%pos_en(:,ke,1) = vx(1) + 0.5_rp*(refelem%x1(:) + 1.0_rp)*(vx(2) - vx(1))
      
      xr(:) = 0.5_rp*(vx(2) - vx(1)) !matmul(refElem%Dr,mesh%x(:,n))
 
      lcmesh%J(:,ke) = xr
      lcmesh%Escale(:,ke,1,1) =   1.0_rp/lcmesh%J(:,ke)
 
      !* Face
 
      !
      !mesh%fx(:,n) = mesh%x(fmask(:),n)
      !mesh%fy(:,n) = mesh%y(fmask(:),n)
 
      ! Calculate normal vectors
      lcmesh%normal_fn(fid(:,1),ke,1) = - 1.0_rp
      lcmesh%normal_fn(fid(:,2),ke,1) = + 1.0_rp
      lcmesh%sJ(:,ke) = 1.0_rp
      lcmesh%Fscale(:,ke) = lcmesh%sJ(:,ke)/lcmesh%J(fmask(:),ke)
      lcmesh%Gsqrt(:,ke) = 1.0_rp
    end do
    !$acc update device(lcmesh%pos_en, lcmesh%normal_fn, lcmesh%sJ, lcmesh%J, lcmesh%Escale, lcmesh%Fscale, lcmesh%Gsqrt)
 
    return
  end subroutine meshbase1d_setgeometricinfo
 
!OCL SERIAL
  subroutine meshbase1d_assigndomid( this, &
    tileID_table, panelID_table,        &
    pi_table )
  
    use scale_meshutil_1d, only: &       
      meshutil1d_buildglobalmap
    implicit none
    
    class(meshbase1d), 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 :: n
    integer :: p
    integer :: tileid
    !-----------------------------------------------------------------------------
    
    call meshutil1d_buildglobalmap( &
      panelid_table, pi_table,                                                      & ! (out)
      this%tileID_globalMap, this%tileFaceID_globalMap, this%tilePanelID_globalMap, & ! (out)
      this%LOCAL_MESH_NUM_global )                                                    ! (in)
 
    do p=1, this%PRC_NUM
    do n=1, this%LOCAL_MESH_NUM
      tileid = n + (p-1)*this%LOCAL_MESH_NUM
 
      tileid_table(n,p)                   = tileid
      this%tileID_global2localMap(tileid) = n
      this%PRCRank_globalMap(tileid)      = p - 1
    end do
    end do
    !$acc update device(this%tileID_globalMap, this%tileFaceID_globalMap,                  &
    !$acc   this%tilePanelID_globalMap, this%tileID_global2localMap, this%PRCRank_globalMap)
 
    return
  end subroutine meshbase1d_assigndomid
 
!OCL SERIAL
  subroutine meshbase1d_setuplocaldom( lcmesh,  &
    tileID, panelID,                            &
    i, Nprc,                                    &
    dom_xmin, dom_xmax,                         &
    Ne, FX )
 
    use scale_meshutil_1d, only: &
      meshutil1d_genconnectivity,   &
      meshutil1d_genlinedomain,     &
      meshutil1d_buildinteriormap,  &
      meshutil1d_genpatchboundarymap
 
    use scale_localmesh_base, only: bctype_interior
    implicit none
      
    type(localmesh1d), intent(inout) :: lcmesh
    integer, intent(in) :: tileid
    integer, intent(in) :: panelid
    integer, intent(in) :: i
    integer, intent(in) :: nprc
    real(rp) :: dom_xmin, dom_xmax
    integer, intent(in) ::ne
    real(rp), intent(in) :: fx(ne*nprc+1)
    
    class(elementbase1d), pointer :: elem
    real(rp) :: delx
    real(rp) :: fx_lc(ne+1)    
    !-----------------------------------------------------------------------------
 
    elem => lcmesh%refElem1D
    lcmesh%tileID = tileid
    lcmesh%panelID = panelid
    !$acc update device(lcmesh%tileID, lcmesh%panelID)
    
    !--
 
    lcmesh%Ne  = ne
    lcmesh%Nv  = ne + 1
    lcmesh%NeS = 1
    lcmesh%NeE = lcmesh%Ne
    lcmesh%NeA = lcmesh%Ne + 2
    !$acc update device(lcmesh%Ne, lcmesh%Nv, lcmesh%NeS, lcmesh%NeE, lcmesh%NeA)
 
    !delx = (dom_xmax - dom_xmin)/dble(Nprc)
    fx_lc(:) = fx((i-1)*ne+1:i*ne+1)
    lcmesh%xmin = fx_lc(1)
    lcmesh%xmax = fx_lc(ne+1)
    !$acc update device(lcmesh%xmin, lcmesh%xmax)
 
    allocate( lcmesh%pos_ev(lcmesh%Nv,1) )
    allocate( lcmesh%EToV(lcmesh%Ne,2) )
    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 meshutil1d_genlinedomain( lcmesh%pos_ev, lcmesh%EToV,   & ! (out)
        lcmesh%Ne, lcmesh%xmin, lcmesh%xmax, fx=fx_lc   )          ! (in)
    !$acc update device(lcmesh%pos_ev, lcmesh%EToV)
 
    !---
    call meshbase1d_setgeometricinfo( lcmesh )
    
    !---
 
    call meshutil1d_genconnectivity( lcmesh%EToE, lcmesh%EToF, & ! (out)
        & lcmesh%EToV, lcmesh%Ne, elem%Nfaces )                  ! (in)
    !$acc update device(lcmesh%EToE, lcmesh%EToF)
 
    !---
    call meshutil1d_buildinteriormap( lcmesh%VmapM, lcmesh%VMapP, lcmesh%MapM, lcmesh%MapP,  & ! (out)
      & lcmesh%pos_en, lcmesh%pos_ev, lcmesh%EToE, lcmesh%EtoF, lcmesh%EtoV,                 & ! (in)
      & elem%Fmask, lcmesh%Ne, elem%Np, elem%Nfp, elem%Nfaces, lcmesh%Nv )                     ! (in)
 
    call meshutil1d_genpatchboundarymap( lcmesh%VMapB, lcmesh%MapB, lcmesh%VMapP, & ! (inout)
      & lcmesh%pos_en, lcmesh%xmin, lcmesh%xmax,                                  & ! (in)
      & elem%Fmask, lcmesh%Ne, elem%Np, elem%Nfp, elem%Nfaces, lcmesh%Nv )          ! (in)
    !$acc update device(lcmesh%VMapM, lcmesh%VMapP, lcmesh%MapM, lcmesh%MapP)
    !$acc enter data copyin(lcmesh%VMapB, lcmesh%MapB)
 
    return
  end subroutine meshbase1d_setuplocaldom
 
end module scale_mesh_base1d