check_mesh_sph_1DeC_EM_file Subroutine

subroutine check_mesh_sph_1DeC_EM_file(correct, rc)

Arguments

Type IntentOptional Attributes Name
logical :: correct
integer :: rc

Source Code

subroutine  check_mesh_sph_1DeC_EM_file(correct, rc)
  type(ESMF_Mesh) :: mesh
  logical :: correct
  integer :: rc
  integer, pointer :: nodeIds(:),nodeOwners(:),nodeMask(:)
  real(ESMF_KIND_R8), pointer :: nodeCoords(:)
  real(ESMF_KIND_R8), pointer :: ownedNodeCoords(:)
  integer :: numNodes, numOwnedNodes, numOwnedNodesTst
  integer :: numElems,numOwnedElemsTst
  integer :: numElemConns, numTriElems, numQuadElems
  real(ESMF_KIND_R8), pointer :: elemCoords(:)
  integer, pointer :: elemIds(:),elemTypes(:),elemConn(:)
  integer, allocatable :: elemMask(:)
  real(ESMF_KIND_R8),allocatable :: elemArea(:)
  integer :: petCount, localPet
  type(ESMF_VM) :: vm
  type(ESMF_DistGrid) :: elemDistgrid
  type(ESMF_Field) :: nodeField, elemField
  integer :: i,j,k
  integer :: numNodesTst, numElemsTst,numElemConnsTst
  integer,allocatable :: elemIdsTst(:)
  integer,allocatable :: elemTypesTst(:)
  integer,allocatable :: elemConnTst(:)
  real(ESMF_KIND_R8),allocatable :: elemAreaTst(:)
  integer,allocatable :: nodeIdsTst(:)
  real(ESMF_KIND_R8),allocatable :: nodeCoordsTst(:)
  integer,allocatable :: nodeOwnersTst(:)
  integer,allocatable :: nodeMaskTst(:)
  integer,allocatable :: elemMaskTst(:)
  real(ESMF_KIND_R8), allocatable :: elemCoordsTst(:)
  logical :: nodeMaskIsPresentTst
  logical :: elemMaskIsPresentTst
  logical :: elemAreaIsPresentTst
  logical :: elemCoordsIsPresentTst             
  integer :: pdim, sdim
  type(ESMF_CoordSys_Flag) :: coordSys

  ! get global VM
  call ESMF_VMGetGlobal(vm, rc=rc)
  if (rc /= ESMF_SUCCESS) return
  call ESMF_VMGet(vm, localPet=localPet, petCount=petCount, rc=rc)
  if (rc /= ESMF_SUCCESS) return


  ! return with an error if not 1 or 4 PETs
  if ((petCount /= 1) .and. (petCount /=4)) then
     rc=ESMF_FAILURE
     return
  endif


  ! Setup mesh info depending on the
  ! number of PETs
  if (petCount .eq. 1) then

     ! Fill in node data
     numNodes=16

     !! node ids
     allocate(nodeIds(numNodes))
     nodeIds=(/1,2,3,4,5,6,7,8, &
               9,10,11,12,13,14,&
               15,16/)

     !! node Mask
     allocate(nodeMask(numNodes))
     nodeMask=(/1,0,1,0,1,0,1,0, &
                1,0,1,0,1,0,&
                1,0/)


     !! node Coords
     allocate(nodeCoords(numNodes*2))
     nodeCoords=(/0.0,0.0, & ! 1
                 1.0,0.0, &  ! 2
                 2.0,0.0, &  ! 3
                 3.0,0.0, &  ! 4
                 0.0,1.0, &  ! 5
                 1.0,1.0, &  ! 6
                 2.0,1.0, &  ! 7
                 3.0,1.0, &  ! 8
                 0.0,2.0, &  ! 9
                 1.0,2.0, &  ! 10
                 2.0,2.0, &  ! 11
                 3.0,2.0, &  ! 12
                 0.0,3.0, &  ! 13
                 1.0,3.0, &  ! 14
                 2.0,3.0, &  ! 15
                 3.0,3.0 /)  ! 16


      !! node owners
      allocate(nodeOwners(numNodes))
      nodeOwners=0 ! everything on proc 0

      ! Fill in elem data
      numTriElems=6
      numQuadElems=6
      numElems=numTriElems+numQuadElems
      numElemConns=3*numTriElems+4*numQuadElems

      !! elem ids
      allocate(elemIds(numElems))
      elemIds=(/1,2,3,4,5,6,7,8,9,10,11,12/)

      !! elem types
      allocate(elemTypes(numElems))
      elemTypes=(/ESMF_MESHELEMTYPE_TRI,  & ! 1
                  ESMF_MESHELEMTYPE_TRI,  & ! 2
                  ESMF_MESHELEMTYPE_TRI,  & ! 3
                  ESMF_MESHELEMTYPE_TRI,  & ! 4
                  ESMF_MESHELEMTYPE_QUAD, & ! 5
                  ESMF_MESHELEMTYPE_QUAD, & ! 6
                  ESMF_MESHELEMTYPE_QUAD, & ! 7
                  ESMF_MESHELEMTYPE_QUAD, & ! 8
                  ESMF_MESHELEMTYPE_QUAD, & ! 9
                  ESMF_MESHELEMTYPE_QUAD, & ! 10
                  ESMF_MESHELEMTYPE_TRI,  & ! 11
                  ESMF_MESHELEMTYPE_TRI/)   ! 12

      !! elem coords
      allocate(elemCoords(2*numElems))
      elemCoords=(/0.25,0.75, & ! 1
                   0.75,0.25, & ! 2
                   1.25,0.75, & ! 3
                   1.75,0.25, & ! 4
                   2.5,0.5,   & ! 5
                   0.5,1.5,   & ! 6
                   1.5,1.5,   & ! 7
                   2.5,1.5,   & ! 8
                   0.5,2.5,   & ! 9
                   1.5,2.5,   & ! 10
                   2.75,2.25, & ! 11
                   2.25,2.75/)  ! 12

      !! elem conn
      allocate(elemConn(numElemConns))
      elemConn=(/1,6,5,     & ! 1
                 1,2,6,     & ! 2
                 2,7,6,     & ! 3
                 2,3,7,     & ! 4
                 3,4,8,7,   & ! 5
                 5,6,10,9,  & ! 6
                 6,7,11,10, & ! 7
                 7,8,12,11, & ! 8
                 9,10,14,13, & ! 9
                 10,11,15,14, & ! 10
                 11,12,16, & ! 11
                 11,16,15/) ! 12


      !! elem mask
      allocate(elemMask(numElems))
      elemMask=(/1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1/)

      !! elem area
      !! (need kinds because not R4 doesn't rep. close enough)
      allocate(elemArea(numElems))
      elemArea=(/0.51_ESMF_KIND_R8, 0.52_ESMF_KIND_R8, 0.53_ESMF_KIND_R8, 0.54_ESMF_KIND_R8, &
           1.5_ESMF_KIND_R8, 1.6_ESMF_KIND_R8, 1.7_ESMF_KIND_R8, 1.8_ESMF_KIND_R8, &
           1.9_ESMF_KIND_R8, 1.1_ESMF_KIND_R8, 0.511_ESMF_KIND_R8,  0.512_ESMF_KIND_R8/)


   else if (petCount .eq. 4) then
     ! Setup mesh data depending on PET
     if (localPet .eq. 0) then

     ! Fill in node data
     numNodes=4

     !! node ids
     allocate(nodeIds(numNodes))
     nodeIds=(/1,2,5,6/)

     !! node Mask
     allocate(nodeMask(numNodes))
     nodeMask=(/1,0,1,0/)


     !! node Coords
     allocate(nodeCoords(numNodes*2))
     nodeCoords=(/0.0,0.0, & ! 1
                 1.0,0.0, &  ! 2
                 0.0,1.0, &  ! 5
                 1.0,1.0 /)  ! 6

      !! node owners
      allocate(nodeOwners(numNodes))
      nodeOwners=0 ! everything on proc 0


      ! Fill in elem data
      numTriElems=2
      numQuadElems=0
      numElems=numTriElems+numQuadElems
      numElemConns=3*numTriElems+4*numQuadElems

      !! elem ids
      allocate(elemIds(numElems))
      elemIds=(/1,2/)

      !! elem types
      allocate(elemTypes(numElems))
      elemTypes=(/ESMF_MESHELEMTYPE_TRI, & ! 1
                  ESMF_MESHELEMTYPE_TRI/)  ! 2

      !! elem coords
      allocate(elemCoords(2*numElems))
      elemCoords=(/0.25,0.75, &  ! 1
                   0.75,0.25/) ! 2

      !! elem conn
      allocate(elemConn(numElemConns))
      elemConn=(/1,4,3, &  ! 1
                 1,2,4 /)  ! 2

      !! elem mask
      allocate(elemMask(numElems))
      elemMask=(/1,0/)

      !! elem area
      !! (need kinds because not R4 doesn't rep. close enough)
      allocate(elemArea(numElems))
      elemArea=(/0.51_ESMF_KIND_R8, 0.52_ESMF_KIND_R8/)

     else if (localPet .eq. 1) then

     ! Fill in node data
     numNodes=6

     !! node ids
     allocate(nodeIds(numNodes))
     nodeIds=(/2,3,4,6,7,8/)

     !! node Mask
     allocate(nodeMask(numNodes))
     nodeMask=(/0,1,0,0,1,0/)

     !! node Coords
     allocate(nodeCoords(numNodes*2))
     nodeCoords=(/1.0,0.0, &  ! 2
                  2.0,0.0, &  ! 3
                  3.0,0.0, &  ! 4
                  1.0,1.0, &  ! 6
                  2.0,1.0, &  ! 7
                  3.0,1.0 /)  ! 8



      !! node owners
      allocate(nodeOwners(numNodes))
      nodeOwners=(/0, & ! 2
                   1, & ! 3
                   1, & ! 4
                   0, & ! 6
                   1, & ! 7
                   1/)  ! 8

      ! Fill in elem data
      numTriElems=2
      numQuadElems=1
      numElems=numTriElems+numQuadElems
      numElemConns=3*numTriElems+4*numQuadElems

      !! elem ids
      allocate(elemIds(numElems))
      elemIds=(/3,4,5/)

      !! elem types
      allocate(elemTypes(numElems))
       elemTypes=(/ESMF_MESHELEMTYPE_TRI, &  ! 3
                   ESMF_MESHELEMTYPE_TRI, &  ! 4
                   ESMF_MESHELEMTYPE_QUAD/)  ! 5

      !! elem coords
      allocate(elemCoords(2*numElems))
      elemCoords=(/1.25,0.75, & ! 3
                   1.75,0.25, & ! 4
                   2.50,0.50/)  ! 5

      !! elem conn
      allocate(elemConn(numElemConns))
      elemConn=(/1,5,4,   &   ! 3
                 1,2,5,   &   ! 4
                 2,3,6,5/)    ! 5

      !! elem mask
      allocate(elemMask(numElems))
      elemMask=(/1,1,1/)

      !! elem area
      !! (need kinds because not R4 doesn't rep. close enough)
      allocate(elemArea(numElems))
      elemArea=(/0.53_ESMF_KIND_R8,0.54_ESMF_KIND_R8,1.5_ESMF_KIND_R8/)

     else if (localPet .eq. 2) then

     ! Fill in node data
     numNodes=9

     !! node ids
     allocate(nodeIds(numNodes))
     nodeIds=(/5,6,7,   &
               9,10,11, &
               13,14,15/)

     !! node Mask
     allocate(nodeMask(numNodes))
     nodeMask=(/1,0,1, &
                1,0,1, &
                1,0,1/)

     !! node Coords
     allocate(nodeCoords(numNodes*2))
     nodeCoords=(/0.0,1.0, &  ! 5
                   1.0,1.0, &  ! 6
                  2.0,1.0, &  ! 7
                  0.0,2.0, &  ! 9
                  1.0,2.0, &  ! 10
                  2.0,2.0, &  ! 11
                  0.0,3.0, &  ! 13
                  1.0,3.0, &  ! 14
                  2.0,3.0/)  ! 15


      !! node owners
      allocate(nodeOwners(numNodes))
      nodeOwners=(/0, & ! 5
                   0, & ! 6
                   1, & ! 7
                   2, & ! 9
                   2, & ! 10
                   2, & ! 11
                   2, & ! 13
                   2, & ! 14
                   2/)  ! 15

      ! Fill in elem data
      numTriElems=0
      numQuadElems=4
      numElems=numTriElems+numQuadElems
      numElemConns=3*numTriElems+4*numQuadElems

       !! elem ids
      allocate(elemIds(numElems))
      elemIds=(/6,7,9,10/)

      !! elem types
      allocate(elemTypes(numElems))
      elemTypes=(/ESMF_MESHELEMTYPE_QUAD, & ! 6
                  ESMF_MESHELEMTYPE_QUAD, & ! 7
                  ESMF_MESHELEMTYPE_QUAD, & ! 9
                  ESMF_MESHELEMTYPE_QUAD/)  ! 10


      !! elem coords
      allocate(elemCoords(2*numElems))
      elemCoords=(/0.5,1.5, & ! 6
                   1.5,1.5, & ! 7
                   0.5,2.5, & ! 9
                   1.5,2.5/)  ! 10

      !! elem conn
      allocate(elemConn(numElemConns))
      elemConn=(/1,2,5,4,  & ! 6
                 2,3,6,5,  & ! 7
                 4,5,8,7,  & ! 9
                 5,6,9,8/)   ! 10

      !! elem mask
      allocate(elemMask(numElems))
      elemMask=(/1,1,1,1/)

      !! elem area
      !! (need kinds because not R4 doesn't rep. close enough)
      allocate(elemArea(numElems))
      elemArea=(/1.6_ESMF_KIND_R8,1.7_ESMF_KIND_R8,1.9_ESMF_KIND_R8,1.1_ESMF_KIND_R8/)

     else if (localPet .eq. 3) then

     ! Fill in node data
     numNodes=6

     !! node ids
     allocate(nodeIds(numNodes))
     nodeIds=(/7,8,11,12,15,16/)

     !! node Mask
     allocate(nodeMask(numNodes))
     nodeMask=(/1,0,1,0,1,0/)

     !! node Coords
     allocate(nodeCoords(numNodes*2))
     nodeCoords=(/2.0,1.0, &  ! 7
                  3.0,1.0, &  ! 8
                  2.0,2.0, &  ! 11
                  3.0,2.0, &  ! 12
                  2.0,3.0, &  ! 15
                  3.0,3.0 /)  ! 16


      !! node owners
      allocate(nodeOwners(numNodes))
      nodeOwners=(/1, & ! 7
                   1, & ! 8
                   2, & ! 11
                   3, & ! 12
                   2, & ! 15
                   3/)  ! 16

      ! Fill in elem data
      numTriElems=2
      numQuadElems=1
      numElems=numTriElems+numQuadElems
      numElemConns=3*numTriElems+4*numQuadElems

      !! elem ids
      allocate(elemIds(numElems))
      elemIds=(/8,11,12/)

      !! elem types
      allocate(elemTypes(numElems))
      elemTypes=(/ESMF_MESHELEMTYPE_QUAD, & ! 8
                  ESMF_MESHELEMTYPE_TRI,  & ! 11
                  ESMF_MESHELEMTYPE_TRI/)   ! 12

      !! elem coords
      allocate(elemCoords(2*numElems))
      elemCoords=(/2.5,1.5, & ! 8
                   2.75,2.25,& ! 11
                   2.25,2.75/)  ! 12

      !! elem conn
      allocate(elemConn(numElemConns))
      elemConn=(/1,2,4,3, & ! 8
                 3,4,6, & ! 11
                 3,6,5/) ! 12

      !! elem mask
      allocate(elemMask(numElems))
      elemMask=(/1,0,1/)

      !! elem area
      !! (need kinds because not R4 doesn't rep. close enough)
      allocate(elemArea(numElems))
      elemArea=(/1.8_ESMF_KIND_R8,0.511_ESMF_KIND_R8,0.512_ESMF_KIND_R8/)
     endif
   endif


   ! Create distgrid to ensure that elements are created on the
   ! same PET as described above
   elemdistgrid=ESMF_DistGridCreate(elemIds, rc=rc)
   if (rc /= ESMF_SUCCESS) return
  

   ! Read mesh from file that's the same as the one described by info set up above
   mesh=ESMF_MeshCreate("data/test_sph_1D_elemConn_esmf.nc", &
        fileformat=ESMF_FILEFORMAT_ESMFMESH, &
        addUserArea=.true., &
        elementDistgrid=elemDistgrid, &
        rc=rc)
   if (rc /= ESMF_SUCCESS) return

   ! DEBUG OUTPUT
   !call ESMF_MeshWrite(mesh,"mesh_1Dec",rc=rc)
   !if (rc .ne. ESMF_SUCCESS) return

   ! Init correct to true before looking for problems
   correct=.true.

   ! Get dim and coord info
   call ESMF_MeshGet(mesh, &
        parametricDim=pdim, &
        spatialDim=sdim, &
        coordSys=coordSys, &
        rc=rc)
   if (rc /= ESMF_SUCCESS) return

   ! Check Counts
   if (pdim .ne. 2) correct=.false.
   if (sdim .ne. 2) correct=.false.
   if (coordSys .ne. ESMF_COORDSYS_SPH_DEG) correct=.false.


   ! Get counts 
   call ESMF_MeshGet(mesh, &
        nodeCount=numNodesTst, &
        elementCount=numElemsTst, &
        elementConnCount=numElemConnsTst, &
        rc=rc)
   if (rc /= ESMF_SUCCESS) return

   ! Check Counts
   if (numNodes .ne. numNodesTst) correct=.false.
   if (numElems .ne. numElemsTst) correct=.false.
   if (numElemConns .ne. numElemConnsTst) correct=.false.

   ! Debugging
   !write(*,*) "numNodes=",numNodes,numNodesTst
   !write(*,*) "numElems=",numElems,numElemsTst
   !write(*,*) "numElemConns=",numElemConns,numElemConnsTst


   ! Get is present information
   call ESMF_MeshGet(mesh, &
        nodeMaskIsPresent=nodeMaskIsPresentTst, &
        elementMaskIsPresent=elemMaskIsPresentTst, &
        elementAreaIsPresent=elemAreaIsPresentTst, &
        elementCoordsIsPresent=elemCoordsIsPresentTst, &
        rc=rc)
   if (rc /= ESMF_SUCCESS) return

   ! Debug
   ! write(*,*) "nodeMaskIsPresent=",nodeMaskIsPresentTst
   ! write(*,*) "elemMaskIsPresent=",elemMaskIsPresentTst
   ! write(*,*) "elemAreaIsPresent=",elemAreaIsPresentTst
   ! write(*,*) "elemCoordsIsPresent=",elemCoordsIsPresentTst

   ! Check is present info
   if (.not. nodeMaskIsPresentTst) correct=.false.
   if (.not. elemMaskIsPresentTst) correct=.false.
   if (.not. elemAreaIsPresentTst) correct=.false.
   if (.not. elemCoordsIsPresentTst) correct=.false.


   ! Allocate space for tst arrays
   allocate(nodeIdsTst(numNodesTst))
   allocate(nodeCoordsTst(2*numNodesTst))
   allocate(nodeOwnersTst(numNodesTst))
   allocate(nodeMaskTst(numNodesTst))
   allocate(elemIdsTst(numElemsTst))
   allocate(elemTypesTst(numElemsTst))
   allocate(elemConnTst(numElemConnsTst))
   allocate(elemMaskTst(numElemsTst))
   allocate(elemAreaTst(numElemsTst))
   allocate(elemCoordsTst(2*numElemsTst))

 ! XMRKX
   ! Get Information
   call ESMF_MeshGet(mesh, &
        nodeIds=nodeIdsTst, &
        nodeCoords=nodeCoordsTst, &
        nodeOwners=nodeOwnersTst, &
        nodeMask=nodeMaskTst, &
        elementIds=elemIdsTst, &
        elementTypes=elemTypesTst, &
        elementConn=elemConnTst, &
        elementMask=elemMaskTst, & 
        elementArea=elemAreaTst, & 
        elementCoords=elemCoordsTst, &
        rc=rc)
   if (rc /= ESMF_SUCCESS) return

!   do i=1,numElemsTst
!      write(*,*) localPet,"# ",elemIdsTst(i),"ea=",elemAreaTst(i),"ec=",elemCoordsTst(2*i-1),elemCoordsTst(2*i),"em=",elemMaskTst(i)
!   enddo


   ! Check node ids
   do i=1,numNodesTst
      if (nodeIds(i) .ne. nodeIdsTst(i)) correct=.false.
   enddo

   ! Check node Coords
   k=1
   do i=1,numNodesTst ! Loop over nodes
      do j=1,2 ! Loop over coord spatial dim
         if (nodeCoords(k) .ne. nodeCoordsTst(k)) correct=.false.
         k=k+1
      enddo
   enddo

   ! Check node Owners
   do i=1,numNodesTst
      if (nodeOwners(i) .ne. nodeOwnersTst(i)) correct=.false.
   enddo

   ! Check node Mask
   do i=1,numNodesTst
      if (nodeMask(i) .ne. nodeMaskTst(i)) correct=.false.
   enddo


   ! Debug output
   ! write(*,*) "nodeMask   =",nodeMask
   ! write(*,*) "nodeMaskTst=",nodeMaskTst

   ! Check elem ids
   do i=1,numElemsTst
      if (elemIds(i) .ne. elemIdsTst(i)) correct=.false.
   enddo

   ! Check elem Types
   do i=1,numElemsTst
      if (elemTypes(i) .ne. elemTypesTst(i)) correct=.false.
   enddo

   ! Check elem Connections
   do i=1,numElemConnsTst
      if (elemConn(i) .ne. elemConnTst(i)) correct=.false.
   enddo

   ! Check elem mask
   do i=1,numElems
      if (elemMask(i) .ne. elemMaskTst(i)) correct=.false.
   enddo

   ! Check elem area
   do i=1,numElems
      if (elemArea(i) .ne. elemAreaTst(i)) correct=.false.
   enddo

   ! Check elem Coords
   k=1
   do i=1,numElemsTst ! Loop over nodes
      do j=1,2 ! Loop over coord spatial dim
         if (elemCoords(k) .ne. elemCoordsTst(k)) correct=.false.
         k=k+1
      enddo
   enddo

   ! Debug output
   ! write(*,*) "elemCoords   =",elemCoords
   ! write(*,*) "elemCoordsTst=",elemCoordsTst

   ! deallocate node data
   deallocate(nodeIds)
   deallocate(nodeCoords)
   deallocate(nodeOwners)
   deallocate(nodeMask)

   ! deallocate elem data
   deallocate(elemIds)
   deallocate(elemTypes)
   deallocate(elemCoords)
   deallocate(elemConn)
   deallocate(elemMask)
   deallocate(elemArea)

   ! Deallocate tst Arrays
   deallocate(nodeIdsTst)
   deallocate(nodeCoordsTst)
   deallocate(nodeOwnersTst)
   deallocate(nodeMaskTst)
   deallocate(elemIdsTst)
   deallocate(elemTypesTst)
   deallocate(elemConnTst)
   deallocate(elemMaskTst)
   deallocate(elemAreaTst)
   deallocate(elemCoordsTst)

   ! Get rid of Mesh
   call ESMF_MeshDestroy(mesh, rc=rc)
   if (rc /= ESMF_SUCCESS) return

  ! Get rid of elemDistgrid
  call ESMF_DistGridDestroy(elemDistgrid, rc=rc)
  if (rc /= ESMF_SUCCESS) return

   ! Return success
   rc=ESMF_SUCCESS

end subroutine check_mesh_sph_1DeC_EM_file