All tests on remus pass

Author: dhyan1272

src/pmpo_MPMesh.cpp

@@ -128,22 +128,16 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){
     auto MPs2Proc = p_MPs->getData<MPF_Tgt_Proc_ID>();
     auto elm2Process = p_mesh->getElm2Process();
     auto elm2global = p_mesh->getElmGlobal();
-    auto mpAppID = p_MPs->getData<polyMPO::MPF_MP_APP_ID>();
 
-    MPI_Comm comm = p_MPs->getMPIComm(); 
-    int comm_rank;
-    MPI_Comm_rank(comm, &comm_rank);
-    
     //Since Mesh is static print pnly for 1 time step 
-    if(printVTPIndex==0) {
-      printVTP_mesh(comm_rank);
+    if(printVTPIndex>=0) {
+      printVTP_mesh(printVTPIndex);
     }
 
     Vec3dView history("positionHistory",numMPs);
     Vec3dView resultLeft("positionResult",numMPs);
     Vec3dView resultRight("positionResult",numMPs);
     Vec3dView mpTgtPosArray("positionTarget",numMPs);
-    Kokkos::View<int*> counter("counter",1);
 
     auto CVTElmCalc = PS_LAMBDA(const int& elm, const int& mp, const int&mask){
         Vec3d MP(mpPositions(mp,0),mpPositions(mp,1),mpPositions(mp,2));
@@ -176,7 +170,8 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){
 
                 if(closestElm<0){
                     MPs2Elm(mp) = iElm;
-                    MPs2Proc(mp) = elm2Process(iElm);
+                    if (elm2Process.size() > 0)
+                        MPs2Proc(mp) = elm2Process(iElm);
                     break;
                 }else{
                     iElm = closestElm;
@@ -194,31 +189,28 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){
                 Vec3d shift = dx.cross(r) * ((1.0-0.7)*dx.magnitude()/(dx.cross(r)).magnitude());
                 Vec3d MPLeft = MParrow + shift;
                 Vec3d MPRight = MParrow - shift;
-                auto xx=Kokkos::atomic_fetch_add(&counter(0), 1);
-                history(xx) = MP;
-                resultLeft(xx) = MPLeft;
-                resultRight(xx) = MPRight;
-                mpTgtPosArray(xx) = MPnew;
+                history(mp) = MP;
+                resultLeft(mp) = MPLeft;
+                resultRight(mp) = MPRight;
+                mpTgtPosArray(mp) = MPnew;
             }
         }
     };
     p_MPs->parallel_for(CVTElmCalc,"CVTTrackingElmCenterBasedCalc");
 
-    if(printVTPIndex>=0 && numMPs>0){
+    if(printVTPIndex>=0){
         Vec3dView::HostMirror h_history = Kokkos::create_mirror_view(history);
         Vec3dView::HostMirror h_resultLeft = Kokkos::create_mirror_view(resultLeft);
         Vec3dView::HostMirror h_resultRight = Kokkos::create_mirror_view(resultRight);
         Vec3dView::HostMirror h_mpTgtPos = Kokkos::create_mirror_view(mpTgtPosArray);
-	Kokkos::View<int*>::HostMirror h_counter = Kokkos::create_mirror_view(counter);
 
         Kokkos::deep_copy(h_history, history);
         Kokkos::deep_copy(h_resultLeft, resultLeft);
         Kokkos::deep_copy(h_resultRight, resultRight);
         Kokkos::deep_copy(h_mpTgtPos, mpTgtPosArray);
-	Kokkos::deep_copy(h_counter, counter);
         // printVTP file
         char* fileOutput = (char *)malloc(sizeof(char) * 256); 
-        sprintf(fileOutput, "polyMPOCVTTrackingElmCenter_MPtracks_%d_%d.vtp", comm_rank, printVTPIndex);
+        sprintf(fileOutput, "polyMPOCVTTrackingElmCenter_MPtracks_%d.vtp", printVTPIndex);
         FILE * pFile = fopen(fileOutput,"w");
         free(fileOutput);   
         fprintf(pFile, "<VTKFile type=\"PolyData\" version=\"1.0\" byte_order=\"LittleEndian\" header_type=\"UInt64\">\n  <PolyData>\n    <Piece NumberOfPoints=\"%d\" NumberOfVerts=\"0\" NumberOfLines=\"%d\" NumberOfStrips=\"0\" NumberOfPolys=\"0\">\n      <Points>\n        <DataArray type=\"Float32\" Name=\"Points\" NumberOfComponents=\"3\" format=\"ascii\">\n",numMPs*4,numMPs*2); 
@@ -378,7 +370,6 @@ void MPMesh::push(){
   bool anyIsMigrating = false;
   do {
     CVTTrackingElmCenterBased(); // move to Tgt_XYZ
-    assert(cudaDeviceSynchronize() == cudaSuccess); 
     p_MPs->updateMPSlice<MPF_Cur_Pos_XYZ, MPF_Tgt_Pos_XYZ>(); // Tgt_XYZ becomes Cur_XYZ
     p_MPs->updateMPSlice<MPF_Cur_Pos_Rot_Lat_Lon, MPF_Tgt_Pos_Rot_Lat_Lon>(); // Tgt becomes Cur
 

src/pmpo_c.cpp

@@ -123,7 +123,6 @@ void polympo_createMPs_f(MPMesh_ptr p_mpmesh,
     }
   }
   auto elm2global = p_mesh->getElmGlobal();
-
   auto mpsPerElm_d = create_mirror_view_and_copy(mpsPerElm, numElms);
   auto active_mp2Elm_d = create_mirror_view_and_copy(active_mp2Elm.data(), numActiveMPs);
   auto active_mpIDs_d = create_mirror_view_and_copy(active_mpIDs.data(), numActiveMPs);
@@ -1250,7 +1249,6 @@ void polympo_setOwningProc_f(MPMesh_ptr p_mpmesh, const int nCells, const int* a
 void polympo_setElmGlobal_f(MPMesh_ptr p_mpmesh, const int nCells, const int* array){
   checkMPMeshValid(p_mpmesh);
   auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; 
-  PMT_ALWAYS_ASSERT(p_mesh->meshEditable());
   Kokkos::View<int*, Kokkos::HostSpace> arrayHost("arrayHost", nCells);
   for (int i = 0; i < nCells; i++) {
     arrayHost(i) = array[i] - 1;  // TODO right now elmID offset is set after MPs initialized

src/pmpo_mesh.hpp

@@ -163,8 +163,7 @@ class Mesh {
     void setOwningProc(IntView owningProc) {PMT_ALWAYS_ASSERT(meshEdit_);
                                             owningProc_ = owningProc; }
 
-    void setElmGlobal(IntView globalElm) {PMT_ALWAYS_ASSERT(meshEdit_);
-                                            globalElm_ = globalElm; }
+    void setElmGlobal(IntView globalElm) {globalElm_ = globalElm;}
     IntView getElmGlobal();
 
     void computeRotLatLonIncr();

test/testFortranCreateRebuildMPs.f90

@@ -25,7 +25,7 @@ subroutine createMPsTest(mpMesh, nCells, numMPs, mp2Elm, isMPActive, mpPosition)
     real(kind=MPAS_RKIND) :: ptOne = 0.1_MPAS_RKIND
     integer, parameter :: MP_ACTIVE = 1
     integer, parameter :: MP_INACTIVE = 0
-    integer, dimension(:), pointer :: mpsPerElm, mp2Elm, isMPActive
+    integer, dimension(:), pointer :: mpsPerElm, mp2Elm, isMPActive, globalElms
     real(kind=MPAS_RKIND), dimension(:,:), pointer :: mpPosition
 
     isMPActive = MP_ACTIVE !no inactive MPs and some changed below
@@ -42,11 +42,13 @@ subroutine createMPsTest(mpMesh, nCells, numMPs, mp2Elm, isMPActive, mpPosition)
     end do
 
     allocate(mpsPerElm(nCells))
+    allocate(globalElms(nCells))
     mpsPerElm = 1 !all elements have 1 MP and some changed below
     mpsPerElm(1) = 0 !1st element has 0 MPs
     mpsPerElm(2) = 2 !2nd element has 2 MPs
     mpsPerElm(3) = 2 !3rd element has 2 MPs 
 
+    call polympo_setElmGlobal(mpMesh, nCells, c_loc(globalElms))
     call polympo_createMPs(mpMesh,nCells,numMPs,c_loc(mpsPerElm),c_loc(mp2Elm),c_loc(isMPActive))
 
     !set mp positions
@@ -80,6 +82,7 @@ subroutine createMPsTest(mpMesh, nCells, numMPs, mp2Elm, isMPActive, mpPosition)
 
     !deallocate MP variables
     deallocate(mpsPerElm)
+    deallocate(globalElms)
 end subroutine
 
 subroutine rebuildMPsTests(mpMesh, numMPs, mp2Elm, isMPActive, mpPosition)

test/testFortranMPAdvection.f90

@@ -218,7 +218,7 @@ program main
   real(kind=MPAS_RKIND), dimension(:), pointer :: xCell, yCell, zCell
   integer, dimension(:,:), pointer :: verticesOnCell, cellsOnCell
   integer :: numMPs, numMPsCount, numPush
-  integer, dimension(:), pointer :: mpsPerElm, mp2Elm, isMPActive, mp2Elm_new
+  integer, dimension(:), pointer :: mpsPerElm, mp2Elm, isMPActive, mp2Elm_new, globalElms
   real(kind=MPAS_RKIND), dimension(:,:), pointer :: mpPosition, mpLatLon, mpPositions_new, mpLatLon_new
   integer, parameter :: MP_ACTIVE = 1
   integer, parameter :: MP_INACTIVE = 0
@@ -284,6 +284,7 @@ program main
 
   allocate(lonCell(nCells))
   allocate(mpsPerElm(nCells))
+  allocate(globalElms(nCells))
   allocate(mp2Elm(numMPs))
   allocate(mp2Elm_new(numMPs))
 
@@ -301,6 +302,7 @@ program main
     mp2Elm(numMPsCount+1:numMPsCount+localNumMPs) = i
     mpsPerElm(i) = localNumMPs
     numMPsCount = numMPsCount + localNumMPs
+    globalElms(i)=i
   end do
 
   call assert(numMPsCount == numMPs, "num mps miscounted")
@@ -356,24 +358,11 @@ program main
   end do
 
   call assert(numMPsCount == numMPs, "num mps miscounted")
-
+  call polympo_setElmGlobal(mpMesh, nCells, c_loc(globalElms))
   call polympo_createMPs(mpMesh,nCells,numMPs,c_loc(mpsPerElm),c_loc(mp2Elm),c_loc(isMPActive))
   call polympo_setMPRotLatLon(mpMesh,2,numMPs,c_loc(mpLatLon))
   call polympo_setMPPositions(mpMesh,3,numMPs,c_loc(mpPosition))
-
-  !Another advection test to test if material poins come back to the same position
-  call runAdvectionTest2(mpMesh, numPush, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius)
-  call polympo_getMPPositions(mpMesh, 3, numMPs, c_loc(mpPositions_new))
-  call polympo_getMPRotLatLon(mpMesh, 2, numMPs, c_loc(mpLatLon_new))
-    call polympo_getMPCurElmID(mpMesh, numMPS, c_loc(mp2Elm_new))
-
-  do i = 1, numMPs
-    if ( abs(mpLatLon_new(2,i)-mpLatLon(2,i)) > max_push_diff ) then
-      max_push_diff = abs(mpLatLon_new(2,i)-mpLatLon(2,i))
-    end if
-  end do
-  call assert(max_push_diff.le.TOLERANCE_PUSH , "MPs donot come back check push!")
-
+ 
   if (testType == "API") then
     call runApiTest(mpMesh, numMPs, nVertices, nCells, numPush, mpLatLon, mpPosition, xVertex, yVertex, zVertex, latVertex)
   else if (testType == "MIGRATION") then
@@ -412,7 +401,7 @@ program main
   deallocate(isMPActive)
   deallocate(mpPosition)
   deallocate(mpLatLon)
-
+  deallocate(globalElms)
   stop
 
 end program

test/testFortranMPAppIDs.f90

@@ -40,7 +40,7 @@ subroutine testAppIDPointer(mpMesh) &
     integer :: setMeshOption, setMPOption
     integer :: ierr, self
     integer :: mpi_comm_handle = MPI_COMM_WORLD
-    integer, dimension(:), pointer :: mpsPerElm, mp2Elm, isMPActive
+    integer, dimension(:), pointer :: mpsPerElm, mp2Elm, isMPActive, globalElms
     type(c_ptr) :: mpMesh
     integer :: nCells, numMPs, appID
     integer, parameter :: MP_ACTIVE = 1
@@ -56,16 +56,19 @@ subroutine testAppIDPointer(mpMesh) &
     setMeshOption = 1 !create a hard coded planar test mesh
     setMPOption = 0   !create an empty set of MPs
     mpMesh = polympo_createMPMesh(setMeshOption, setMPOption)
-
     numMPs = 1
     allocate(mpsPerElm(nCells))
+    allocate(globalElms(nCells))
     allocate(mp2Elm(numMPs))
     allocate(isMPActive(numMPs))
     call polympo_getMeshNumElms(mpMesh, nCells)
     mpsPerElm = 1
     mp2Elm = 1
     isMPActive = MP_ACTIVE
+    !This is a dummy array of global Cell IDs for each mesh element
+    call polympo_setElmGlobal(mpMesh, nCells, c_loc(globalElms))
     call polympo_createMPs(mpMesh, nCells, numMPs, c_loc(mpsPerElm), c_loc(mp2Elm), c_loc(isMPActive))
+    print *, "Done here"
     call polympo_setMPICommunicator(mpMesh, mpi_comm_handle)
     ! Set function and opaque data structure(list/queue) used to retrieve appIDS
     call polympo_setAppIDFunc(mpMesh, c_funloc(GetAppID), c_loc(queue));
@@ -74,9 +77,15 @@ subroutine testAppIDPointer(mpMesh) &
     call testAppIDPointer(mpMesh)
 
     ! Clean Up
+    deallocate(mpsPerElm)
+    deallocate(globalElms)
+    deallocate(mp2Elm)
+    deallocate(isMPActive)
+    
+    
     call polympo_deleteMPMesh(mpMesh)
     call queue_destroy(queue)
     call polympo_finalize()
     call mpi_finalize(ierr)
 
-end program
+end program

test/testFortranMPReconstruction.f90

@@ -32,7 +32,7 @@ program main
 
   integer, dimension(:,:), pointer :: verticesOnCell, cellsOnCell
   integer :: numMPs, vID
-  integer, dimension(:), pointer :: mpsPerElm, mp2Elm, isMPActive
+  integer, dimension(:), pointer :: mpsPerElm, mp2Elm, isMPActive, globalElms
   real(kind=MPAS_RKIND), dimension(:,:), pointer :: mpPosition, mpLatLon
   real(kind=MPAS_RKIND), dimension(:,:), pointer :: mpMass, mpVel
   real(kind=MPAS_RKIND), dimension(:), pointer :: meshVtxMass, meshElmMass, meshVtxMass1
@@ -84,6 +84,7 @@ program main
   !createMPs
   numMPs = nCells
   allocate(mpsPerElm(nCells))
+  allocate(globalElms(nCells))
   allocate(mp2Elm(numMPs))
   allocate(isMPActive(numMPs))
   allocate(mpPosition(3,numMPs))
@@ -109,7 +110,8 @@ program main
     mpPosition(2,i) = yCell(i)
     mpPosition(3,i) = zCell(i)
   end do
-
+  
+  call polympo_setElmGlobal(mpMesh, nCells, c_loc(globalElms))
   call polympo_createMPs(mpMesh,nCells,numMPs,c_loc(mpsPerElm),c_loc(mp2Elm),c_loc(isMPActive))
   call polympo_setMPICommunicator(mpMesh, mpi_comm_handle)
   call polympo_setMPRotLatLon(mpMesh,2,numMPs,c_loc(mpLatLon))
@@ -188,7 +190,7 @@ program main
   deallocate(meshElmMass)
   deallocate(meshVtxVelu)
   deallocate(meshVtxVelv)
-
+  deallocate(globalElms)
   stop
 
   contains

test/testMPAppIDs.cpp

@@ -67,4 +67,4 @@ void testAppIDPointer(MPMesh_ptr p_mpmesh) {
   PMT_ALWAYS_ASSERT(numAddedMPsAfter_h == numAddedMPs);
 }
 
-#endif
+#endif