LCOV - code coverage report
Current view: | top level - dynamics/se/dycore - dof_mod.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 59 | 143 | 41.3 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 4 | 13 | 30.8 %

              Line data    Source code
           1             : module dof_mod
           2             :   use shr_kind_mod,   only: r8=>shr_kind_r8, i8=>shr_kind_i8
           3             :   use dimensions_mod, only: np, npsq, nelem, nelemd
           4             :   use quadrature_mod, only: quadrature_t
           5             :   use element_mod,    only: element_t,index_t
           6             :   use spmd_utils,     only: mpi_integer
           7             :   use parallel_mod,   only: parallel_t
           8             :   use edge_mod,       only: initedgebuffer,freeedgebuffer,            &
           9             :                longedgevpack, longedgevunpackmin
          10             :   use edgetype_mod,   only: longedgebuffer_t
          11             :   use bndry_mod,      only: bndry_exchange
          12             : implicit none
          13             : private
          14             :   ! public data
          15             :   ! public subroutines
          16             :   public :: global_dof
          17             :   public :: UniquePoints
          18             :   public :: PutUniquePoints
          19             :   public :: UniqueNcolsP
          20             :   public :: UniqueCoords
          21             :   public :: CreateUniqueIndex
          22             :   public :: SetElemOffset
          23             :   public :: CreateMetaData
          24             :
          25             :   interface UniquePoints
          26             :      module procedure UniquePoints2D
          27             :      module procedure UniquePoints3D
          28             :      module procedure UniquePoints4D
          29             :   end interface
          30             :   interface PutUniquePoints
          31             :      module procedure PutUniquePoints2D
          32             :      module procedure PutUniquePoints3D
          33             :      module procedure PutUniquePoints4D
          34             :   end interface
          35             :
          36             :
          37             : contains
          38             :
          39       21600 :   subroutine genLocalDof(ig,npts,ldof)
          40             :
          41             :     integer, intent(in) :: ig
          42             :     integer, intent(in) :: npts
          43             :     integer, intent(inout) :: ldof(:,:)
          44             :
          45             :     integer :: i,j,npts2
          46             :
          47             :
          48       21600 :      npts2=npts*npts
          49      108000 :      do j=1,npts
          50      453600 :        do i=1,npts
          51      432000 :            ldof(i,j) = (ig-1)*npts2 + (j-1)*npts + i
          52             :        enddo
          53             :      enddo
          54             :
          55       21600 :   end subroutine genLocalDOF
          56             :
          57             : ! ===========================================
          58             : ! global_dof
          59             : !
          60             : ! Compute the global degree of freedom for each element...
          61             : ! ===========================================
          62             :
          63        1536 :   subroutine global_dof(par,elem)
          64             :
          65             :     type (parallel_t),intent(in) :: par
          66             :     type (element_t)             :: elem(:)
          67             :
          68             :     type (LongEdgeBuffer_t)    :: edge
          69             :
          70             :     real(kind=r8)  da                     ! area element
          71             :
          72             :     type (quadrature_t) :: gp
          73             :
          74        3072 :     integer :: ldofP(np,np,nelemd)
          75             :
          76             :     integer ii
          77             :     integer i,j,ig,ie
          78             :     integer kptr
          79             :     integer iptr
          80             :
          81             :     ! ===================
          82             :     ! begin code
          83             :     ! ===================
          84        1536 :     call initEdgeBuffer(edge,1)
          85             :
          86             :     ! =================================================
          87             :     ! mass matrix on the velocity grid
          88             :     ! =================================================
          89             :
          90             :
          91       12336 :     do ie=1,nelemd
          92       10800 :        ig = elem(ie)%vertex%number
          93       10800 :        call genLocalDOF(ig,np,ldofP(:,:,ie))
          94             :
          95       10800 :        kptr=0
          96       12336 :        call LongEdgeVpack(edge,ldofP(:,:,ie),1,kptr,elem(ie)%desc)
          97             :     end do
          98             :
          99             :     ! ==============================
         100             :     ! Insert boundary exchange here
         101             :     ! ==============================
         102             :
         103        1536 :     call bndry_exchange(par,edge)
         104             :
         105       12336 :     do ie=1,nelemd
         106             :        ! we should unpack directly into elem(ie)%gdofV, but we dont have
         107             :        ! a VunpackMIN that takes integer*8.  gdofV integer*8 means
         108             :        ! more than 2G grid points.
         109       10800 :        kptr=0
         110       10800 :        call LongEdgeVunpackMIN(edge,ldofP(:,:,ie),1,kptr,elem(ie)%desc)
         111      228336 :        elem(ie)%gdofP(:,:)=ldofP(:,:,ie)
         112             :     end do
         113             : !$OMP BARRIER
         114        1536 :     call FreeEdgeBuffer(edge)
         115             :
         116        1536 :   end subroutine global_dof
         117             :
         118             :
         119           0 :   subroutine UniquePoints2D(idxUnique,src,dest)
         120             :     type (index_t) :: idxUnique
         121             :     real (kind=r8) :: src(:,:)
         122             :     real (kind=r8) :: dest(:)
         123             :
         124             :     integer :: i,j,ii
         125             :
         126             :
         127           0 :     do ii=1,idxUnique%NumUniquePts
         128           0 :        i=idxUnique%ia(ii)
         129           0 :        j=idxUnique%ja(ii)
         130           0 :        dest(ii)=src(i,j)
         131             :     enddo
         132             :
         133           0 :   end subroutine UniquePoints2D
         134             :
         135             : ! putUniquePoints first zeros out the destination array, then fills the unique points of the
         136             : ! array with values from src.  A boundary communication should then be called to fill in the
         137             : ! redundent points of the array
         138             :
         139           0 :   subroutine putUniquePoints2D(idxUnique,src,dest)
         140             :     type (index_t) :: idxUnique
         141             :     real (kind=r8),intent(in) :: src(:)
         142             :     real (kind=r8),intent(out) :: dest(:,:)
         143             :
         144             :     integer :: i,j,ii
         145             :
         146           0 :     dest=0.0D0
         147           0 :     do ii=1,idxUnique%NumUniquePts
         148           0 :        i=idxUnique%ia(ii)
         149           0 :        j=idxUnique%ja(ii)
         150           0 :        dest(i,j)=src(ii)
         151             :     enddo
         152             :
         153           0 :   end subroutine putUniquePoints2D
         154             :
         155           0 :   subroutine UniqueNcolsP(elem,idxUnique,cid)
         156             :     use element_mod, only : GetColumnIdP, element_t
         157             :     type (element_t), intent(in) :: elem
         158             :     type (index_t), intent(in) :: idxUnique
         159             :     integer,intent(out) :: cid(:)
         160             :     integer :: i,j,ii
         161             :
         162             :
         163           0 :     do ii=1,idxUnique%NumUniquePts
         164           0 :        i=idxUnique%ia(ii)
         165           0 :        j=idxUnique%ja(ii)
         166           0 :        cid(ii)=GetColumnIdP(elem,i,j)
         167             :     enddo
         168             :
         169           0 :   end subroutine UniqueNcolsP
         170             :
         171             :
         172           0 :   subroutine UniqueCoords(idxUnique,src,lat,lon)
         173             :
         174             :     use coordinate_systems_mod, only  : spherical_polar_t
         175             :     type (index_t), intent(in) :: idxUnique
         176             :
         177             :     type (spherical_polar_t) :: src(:,:)
         178             :     real (kind=r8), intent(out) :: lat(:)
         179             :     real (kind=r8), intent(out) :: lon(:)
         180             :
         181             :     integer :: i,j,ii
         182             :
         183           0 :     do ii=1,idxUnique%NumUniquePts
         184           0 :        i=idxUnique%ia(ii)
         185           0 :        j=idxUnique%ja(ii)
         186           0 :        lat(ii)=src(i,j)%lat
         187           0 :        lon(ii)=src(i,j)%lon
         188             :     enddo
         189             :
         190           0 :   end subroutine UniqueCoords
         191             :
         192           0 :   subroutine UniquePoints3D(idxUnique,nlyr,src,dest)
         193             :     type (index_t) :: idxUnique
         194             :     integer :: nlyr
         195             :     real (kind=r8) :: src(:,:,:)
         196             :     real (kind=r8) :: dest(:,:)
         197             :
         198             :     integer :: i,j,k,ii
         199             :
         200           0 :     do ii=1,idxUnique%NumUniquePts
         201           0 :        i=idxUnique%ia(ii)
         202           0 :        j=idxUnique%ja(ii)
         203           0 :        do k=1,nlyr
         204           0 :           dest(ii,k)=src(i,j,k)
         205             :        enddo
         206             :     enddo
         207             :
         208           0 :   end subroutine UniquePoints3D
         209           0 :   subroutine UniquePoints4D(idxUnique,d3,d4,src,dest)
         210             :     type (index_t) :: idxUnique
         211             :     integer :: d3,d4
         212             :     real (kind=r8) :: src(:,:,:,:)
         213             :     real (kind=r8) :: dest(:,:,:)
         214             :
         215             :     integer :: i,j,k,n,ii
         216             :
         217           0 :     do n=1,d4
         218           0 :        do k=1,d3
         219           0 :           do ii=1,idxUnique%NumUniquePts
         220           0 :              i=idxUnique%ia(ii)
         221           0 :              j=idxUnique%ja(ii)
         222           0 :              dest(ii,k,n)=src(i,j,k,n)
         223             :           enddo
         224             :        end do
         225             :     enddo
         226             :
         227           0 :   end subroutine UniquePoints4D
         228             :
         229             : ! putUniquePoints first zeros out the destination array, then fills the unique points of the
         230             : ! array with values from src.  A boundary communication should then be called to fill in the
         231             : ! redundent points of the array
         232             :
         233           0 :   subroutine putUniquePoints3D(idxUnique,nlyr,src,dest)
         234             :     type (index_t) :: idxUnique
         235             :     integer :: nlyr
         236             :     real (kind=r8),intent(in) :: src(:,:)
         237             :     real (kind=r8),intent(out) :: dest(:,:,:)
         238             :
         239             :     integer :: i,j,k,ii
         240             :
         241           0 :     dest=0.0D0
         242           0 :     do k=1,nlyr
         243           0 :        do ii=1,idxUnique%NumUniquePts
         244           0 :           i=idxUnique%ia(ii)
         245           0 :           j=idxUnique%ja(ii)
         246           0 :           dest(i,j,k)=src(ii,k)
         247             :        enddo
         248             :     enddo
         249             :
         250           0 :   end subroutine putUniquePoints3D
         251             :
         252           0 :   subroutine putUniquePoints4D(idxUnique,d3,d4,src,dest)
         253             :     type (index_t) :: idxUnique
         254             :     integer :: d3,d4
         255             :     real (kind=r8),intent(in) :: src(:,:,:)
         256             :     real (kind=r8),intent(out) :: dest(:,:,:,:)
         257             :
         258             :     integer :: i,j,k,n,ii
         259             :
         260           0 :     dest=0.0D0
         261           0 :     do n=1,d4
         262           0 :        do k=1,d3
         263           0 :           do ii=1,idxunique%NumUniquePts
         264           0 :              i=idxUnique%ia(ii)
         265           0 :              j=idxUnique%ja(ii)
         266           0 :              dest(i,j,k,n)=src(ii,k,n)
         267             :           enddo
         268             :        enddo
         269             :     end do
         270           0 :   end subroutine putUniquePoints4D
         271             :
         272        1536 :   subroutine SetElemOffset(par,elem,GlobalUniqueColsP)
         273             :     use spmd_utils, only : mpi_sum
         274             :
         275             :     type (parallel_t)    :: par
         276             :     type (element_t)     :: elem(:)
         277             :     integer, intent(out) :: GlobalUniqueColsP
         278             :
         279        1536 :     integer, allocatable :: numElemP(:),numElem2P(:)
         280             :     integer, allocatable :: numElemV(:),numElem2V(:)
         281        1536 :     integer, allocatable :: gOffset(:)
         282             :
         283             :     integer            :: ie, ig, nprocs, ierr
         284             :     logical, parameter :: Debug = .FALSE.
         285             :
         286        1536 :     nprocs = par%nprocs
         287        4608 :     allocate(numElemP(nelem))
         288        3072 :     allocate(numElem2P(nelem))
         289        3072 :     allocate(gOffset(nelem))
         290    24887808 :     numElemP=0;numElem2P=0;gOffset=0
         291             :
         292       12336 :     do ie = 1, nelemd
         293       10800 :       ig = elem(ie)%GlobalId
         294       12336 :       numElemP(ig) = elem(ie)%idxP%NumUniquePts
         295             :     end do
         296        1536 :     call MPI_Allreduce(numElemP,numElem2P,nelem,MPI_INTEGER,MPI_SUM,par%comm,ierr)
         297             :
         298        1536 :     gOffset(1)=1
         299     8294400 :     do ig = 2, nelem
         300     8294400 :       gOffset(ig) = gOffset(ig-1)+numElem2P(ig-1)
         301             :     end do
         302       12336 :     do ie = 1, nelemd
         303       10800 :       ig = elem(ie)%GlobalId
         304       12336 :       elem(ie)%idxP%UniquePtOffset=gOffset(ig)
         305             :     end do
         306        1536 :     GlobalUniqueColsP = gOffset(nelem)+numElem2P(nelem)-1
         307             :
         308        1536 :     deallocate(numElemP)
         309        1536 :     deallocate(numElem2P)
         310        1536 :     deallocate(gOffset)
         311        1536 :   end subroutine SetElemOffset
         312             :
         313       10800 :   subroutine CreateUniqueIndex(ig,gdof,idx)
         314             :
         315             :     integer :: ig
         316             :     type (index_t) :: idx
         317             :     integer(i8) :: gdof(:,:)
         318             :
         319       10800 :     integer, allocatable :: ldof(:,:)
         320             :     integer :: i,j,ii,npts
         321             :
         322             :
         323       10800 :     npts = size(gdof,dim=1)
         324       43200 :     allocate(ldof(npts,npts))
         325             :     ! ====================
         326             :     ! Form the local DOF
         327             :     ! ====================
         328       10800 :     call genLocalDOF(ig,npts,ldof)
         329             :
         330       10800 :     ii=1
         331             :
         332       54000 :     do j=1,npts
         333      226800 :        do i=1,npts
         334             :           ! ==========================
         335             :           ! check for point ownership
         336             :           ! ==========================
         337      216000 :           if(gdof(i,j) .eq. ldof(i,j)) then
         338       97204 :              idx%ia(ii) = i
         339       97204 :              idx%ja(ii) = j
         340       97204 :              ii=ii+1
         341             :           endif
         342             :        enddo
         343             :     enddo
         344             :
         345       10800 :     idx%NumUniquePts=ii-1
         346       10800 :     deallocate(ldof)
         347             :
         348       10800 :   end subroutine CreateUniqueIndex
         349             :
         350             :
         351           0 :   subroutine CreateMetaData(par,elem,subelement_corners, fdofp)
         352             :     type (parallel_t),          intent(in)  :: par
         353             :     type (element_t), target                :: elem(:)
         354             :
         355             :     integer,          optional, intent(out) :: subelement_corners((np-1)*(np-1)*nelemd,4)
         356             :     integer,          optional              :: fdofp(np,np,nelemd)
         357             :
         358             :     type (index_t),   pointer               :: idx
         359             :     type (LongEdgeBuffer_t)                 :: edge
         360             :     integer                                 :: i, j, ii, ie, base
         361             :     integer(i8),      pointer               :: gdof(:,:)
         362           0 :     integer                                 :: fdofp_local(np,np,nelemd)
         363             :
         364           0 :     call initEdgeBuffer(edge,1)
         365           0 :     fdofp_local=0
         366             :
         367           0 :     do ie=1,nelemd
         368           0 :        idx => elem(ie)%idxP
         369           0 :        do ii=1,idx%NumUniquePts
         370           0 :           i=idx%ia(ii)
         371           0 :           j=idx%ja(ii)
         372             :
         373           0 :           fdofp_local(i,j,ie) = -(idx%UniquePtoffset+ii-1)
         374             :        end do
         375           0 :        call LongEdgeVpack(edge,fdofp_local(:,:,ie),1,0,elem(ie)%desc)
         376             :     end do
         377           0 :     call bndry_exchange(par,edge)
         378           0 :     do ie=1,nelemd
         379           0 :        base = (ie-1)*(np-1)*(np-1)
         380           0 :        call LongEdgeVunpackMIN(edge,fdofp_local(:,:,ie),1,0,elem(ie)%desc)
         381           0 :        if(present(subelement_corners)) then
         382             :           ii=0
         383           0 :           do j=1,np-1
         384           0 :              do i=1,np-1
         385           0 :                 ii=ii+1
         386           0 :                 subelement_corners(base+ii,1) = -fdofp_local(i,j,ie)
         387           0 :                 subelement_corners(base+ii,2) = -fdofp_local(i,j+1,ie)
         388           0 :                 subelement_corners(base+ii,3) = -fdofp_local(i+1,j+1,ie)
         389           0 :                 subelement_corners(base+ii,4) = -fdofp_local(i+1,j,ie)
         390             :              end do
         391             :           end do
         392             :        end if
         393             :     end do
         394           0 :     if(present(fdofp)) then
         395           0 :        fdofp=-fdofp_local
         396             :     end if
         397             :
         398             :
         399             :
         400           0 :   end subroutine CreateMetaData
         401             :
         402             : end module dof_mod

Generated by: LCOV version 1.14