Ctests passing in remus

Author: dhyan1272

src/pmpo_MPMesh_assembly.hpp

@@ -28,6 +28,7 @@ DoubleView MPMesh::assemblyV0(){
 template <MeshFieldIndex meshFieldIndex>
 void MPMesh::assemblyVtx0(){
   Kokkos::Timer timer;
+
   constexpr MaterialPointSlice mpfIndex = meshFieldIndexToMPSlice<meshFieldIndex>;
   auto elm2VtxConn = p_mesh->getElm2VtxConn();  
   auto mpData = p_MPs->getData<mpfIndex>();
@@ -94,7 +95,8 @@ void MPMesh::assemblyElm0() {
   pumipic::RecordTime("PolyMPO_Reconstruct_Elm0", timer.seconds());
 }
 
-void MPMesh::reconstruct_coeff_full(){  
+void MPMesh::reconstruct_coeff_full(){
+  std::cout<<__FUNCTION__<<std::endl;
   Kokkos::Timer timer;
   int self, numProcsTot;
   MPI_Comm comm = p_MPs->getMPIComm();
@@ -142,7 +144,7 @@ void MPMesh::reconstruct_coeff_full(){
         Kokkos::atomic_add(&vtxMatrices(vID,6), w_vtx*(-vtxCoords(vID,0)+mpPos(mp,0))*(-vtxCoords(vID,2)+mpPos(mp,2))/(radius*radius));
         Kokkos::atomic_add(&vtxMatrices(vID,7), w_vtx*(-vtxCoords(vID,1)+mpPos(mp,1))*(-vtxCoords(vID,1)+mpPos(mp,1))/(radius*radius));
         Kokkos::atomic_add(&vtxMatrices(vID,8), w_vtx*(-vtxCoords(vID,1)+mpPos(mp,1))*(-vtxCoords(vID,2)+mpPos(mp,2))/(radius*radius));
-        Kokkos::atomic_add(&vtxMatrices(vID,9), w_vtx*(-vtxCoords(vID,2)+mpPos(mp,2))*(-vtxCoords(vID,2)+mpPos(mp,2))/(radius*radius));  
+        Kokkos::atomic_add(&vtxMatrices(vID,9), w_vtx*(-vtxCoords(vID,2)+mpPos(mp,2))*(-vtxCoords(vID,2)+mpPos(mp,2))/(radius*radius));
       }
     }
   };
@@ -278,9 +280,9 @@ void MPMesh::invertMatrix(const Kokkos::View<double**>& vtxMatrices, const doubl
     VtxCoeffs(vtx, 0, 2) = temp[1];
     VtxCoeffs(vtx, 0, 3) = temp[2];
 
-    //Debugging
+    //Debugging 
     /*
-    if(vtx == 2630){
+    if(vtx == 10){
       printf("Matrices %d vtx \n", vtx);
       printf("[ %.15e  %.15e  %.15e  %.15e ]\n", vtxMatrices(vtx,0), vtxMatrices(vtx,1), vtxMatrices(vtx,2), vtxMatrices(vtx,3));
       printf("[ %.15e  %.15e  %.15e  %.15e ]\n", vtxMatrices(vtx,1), vtxMatrices(vtx,4), vtxMatrices(vtx,5), vtxMatrices(vtx,6));
@@ -306,6 +308,7 @@ void MPMesh::invertMatrix(const Kokkos::View<double**>& vtxMatrices, const doubl
 
 template <MeshFieldIndex meshFieldIndex>
 void MPMesh::assemblyVtx1() {
+  std::cout<<__FUNCTION__<<std::endl;
   Kokkos::Timer timer;
 
   int self, numProcsTot;
@@ -368,22 +371,23 @@ void MPMesh::assemblyVtx1() {
     communicate_and_take_halo_contributions(meshField, numVertices, numEntries, 0, 0);
   pumipic::RecordTime("Communicate Field Values" + std::to_string(self), timer.seconds());
 
-  //Debug
-  //Kokkos::fence();
-  //assert(cudaDeviceSynchronize() == cudaSuccess);
+  //Debugging
+  /*
+  assert(cudaDeviceSynchronize() == cudaSuccess);
   Kokkos::parallel_for("printSymmetricBlock", numVertices, KOKKOS_LAMBDA(const int vtx){
     if (vtx >= 10 && vtx <= 10) {
       printf("Field in %d:  ", vtx);
       for (int k=0; k<numEntries; k++)
         printf("  %.15e ", meshField(vtx, k));
       printf("\n");
     }
-  });  
+  });
+  */  
 }
 
 //Start Communication routine
 void MPMesh::startCommunication(){
-
+  std::cout<<__FUNCTION__<<std::endl;
   Kokkos::Timer timer;
   int self, numProcsTot;
   MPI_Comm comm = p_MPs->getMPIComm(); 
@@ -431,7 +435,7 @@ void MPMesh::startCommunication(){
                          ent2global);
 
   //Do Map of Global To Local ID
-  //TODO make ordered map; which faster?
+  //Check unordered vs ordered map; which faster?
   std::map<int, int> global2local;
   //std::unordered_map<int, int> global2local;
   for (int iEnt = 0; iEnt < numEntities; iEnt++) {
@@ -527,12 +531,7 @@ void MPMesh::startCommunication(){
 
   pumipic::RecordTime("Start Communication" + std::to_string(self), timer.seconds());
 
-  bool isRotated = p_mesh->getRotatedFlag();
-  p_mesh->setGnomonicProjection(isRotated);
-
   /*
-  if (p_MPs->getOpMode() != polyMPO::MP_DEBUG) 
-   return;
   printf("Rank %d Owners %d Halos %d Total %d \n", self, numOwnersTot, numHalosTot, numEntities);
   for (int i=0; i<numProcsTot; i++){
     printf("Rank %d has %d halos which are owners in other rank %d \n", self, numOwnersOnOtherProcs[i], i);
@@ -757,8 +756,6 @@ void MPMesh::communicateFields(const std::vector<std::vector<double>>& fieldData
   MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
 
   /*
-  if (p_MPs->getOpMode() != polyMPO::MP_DEBUG)
-    return;
   static int count_deb=0;
   if(self==0) std::cout<<"====================="<<count_deb<<"========================"<<std::endl;
   count_deb++;

src/pmpo_c.cpp

@@ -1199,6 +1199,13 @@ void polympo_getMeshDualTriangleArea_f(MPMesh_ptr p_mpmesh, const int nVertices,
     areaTriangle[i] = h_dualArea(i,0);
 }
 
+void polympo_setGnomonicProjection_f(MPMesh_ptr p_mpmesh){
+  checkMPMeshValid(p_mpmesh);
+  auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
+  bool isRotated = p_mesh->getRotatedFlag();
+  p_mesh->setGnomonicProjection(isRotated);
+}
+
 void polyMPO_setTanLatVertexRotatedOverRadius_f(MPMesh_ptr p_mpmesh, const int nVertices, double* array){
   //chech validity
   checkMPMeshValid(p_mpmesh);

src/pmpo_c.h

@@ -93,6 +93,7 @@ void polympo_setMeshElmCenter_f(MPMesh_ptr p_mpmesh, const int nCells, const dou
 void polympo_getMeshElmCenter_f(MPMesh_ptr p_mpmesh, const int nCells, double* xArray, double* yArray, double* zArray);
 void polympo_setMeshDualTriangleArea_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* areaTriangle);
 void polympo_getMeshDualTriangleArea_f(MPMesh_ptr p_mpmesh, const int nVertices, double* areaTriangle);
+void polympo_setGnomonicProjection_f(MPMesh_ptr p_mpmesh);
 void polyMPO_setTanLatVertexRotatedOverRadius_f(MPMesh_ptr p_mpmesh, const int nVertices, double* array); 
 
 // Advection calculations

src/pmpo_fortran.f90

@@ -885,16 +885,21 @@ subroutine polympo_getMeshDualTriangleArea(mpMesh, nVertices, areaTriangle) &
     integer(c_int), value :: nVertices
     type(c_ptr), value :: areaTriangle
   end subroutine
-  
+
+  subroutine polympo_setGnomonicProjection(mpMesh) &
+             bind(C, NAME='polympo_setGnomonicProjection_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+  end subroutine
+
   subroutine  polyMPO_setTanLatVertexRotatedOverRadius(mpMesh, nVertices, array) &
-             bind(C, NAME=' polyMPO_setTanLatVertexRotatedOverRadius_f')
+              bind(C, NAME=' polyMPO_setTanLatVertexRotatedOverRadius_f')
     use :: iso_c_binding
     type(c_ptr), value :: mpMesh
     integer(c_int), value :: nVertices
     type(c_ptr), value :: array
   end subroutine
 
-
   !---------------------------------------------------------------------------
   !> @brief calculate the MPs from given mesh vertices rotational latitude
   !>        longitude, update the MP slices

test/CMakeLists.txt

@@ -112,7 +112,7 @@ else()
 endif()
 pmpo_add_test(testFortranReadMPAS 1 ./testFortranReadMPAS ${TEST_NC_FILE})
 pmpo_add_test(testFortranCreateRebuildMPs 1 ./testFortranCreateRebuildMPs ${TEST_NC_FILE})
-#pmpo_add_test(testFortranMPReconstruction 1 ./testFortranMPReconstruction ${TEST_NC_FILE})
+pmpo_add_test(testFortranMPReconstruction 1 ./testFortranMPReconstruction ${TEST_NC_FILE})
 
 pmpo_add_test(test_timing 1 ./timeAssmblyWachspress 1 ${TEST_NC_FILE})
 pmpo_add_test(test_wachspress 1 ./testWachspress ${TEST_NC_FILE})

test/testFortranMPAdvection.f90

@@ -109,13 +109,16 @@ subroutine runReconstructionTest(mpMesh, numMPs, numPush, nCells, nVertices, mp2
     call polympo_setMPMass(mpMesh,1,numMPs,c_loc(mpMass))
     call polympo_setMPVel(mpMesh,2,numMPs,c_loc(mpVel))
 
+    ! Although this test just does 0th order reconstruction testing, and just needs the BasisSlice,
+    ! calculating the coefficeints too as that will involve calculating the Basis Slice
+    call polympo_reconstruct_coeff_with_MPI(mpmesh)
     ! Test push reconstruction
-
     do j = 1, numPush
       call calcSurfDispIncr(mpMesh, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius)
       call polympo_setReconstructionOfMass(mpMesh,0,polympo_getMeshFElmType())
       call polympo_setReconstructionOfMass(mpMesh,0,polympo_getMeshFVtxType())
       call polympo_setReconstructionOfVel(mpMesh,0,polympo_getMeshFVtxType())
+      call polympo_applyReconstruction(mpMesh)
       call polympo_push(mpMesh)
       call polympo_getMeshElmMass(mpMesh,nCells,c_loc(meshElmMass))
       call polympo_getMeshVtxMass(mpMesh,nVertices,c_loc(meshVtxMass))
@@ -270,7 +273,8 @@ program main
                         latVertex, &
                         xCell, yCell, zCell, &
                         verticesOnCell, cellsOnCell)
-
+  
+  call polympo_setGnomonicProjection(mpMesh)
   call polympo_setMPICommunicator(mpMesh, mpi_comm_handle);
 
   !createMPs

test/testFortranMPReconstruction.f90

@@ -78,7 +78,8 @@ program main
                         latVertex, &
                         xCell, yCell, zCell, &
                         verticesOnCell, cellsOnCell, areaTriangle)
- 
+  
+  call polympo_setGnomonicProjection(mpMesh)
   nCompsDisp = 2
   allocate(dispIncr(nCompsDisp,nVertices))
   !createMPs
@@ -120,29 +121,30 @@ program main
   call polympo_setMPMass(mpMesh,1,numMPs,c_loc(mpMass))
   call polympo_setMPVel(mpMesh,2,numMPs,c_loc(mpVel))
 
-  ! Test vtx reconstruction
-  call polympo_setReconstructionOfMass(mpMesh,0,polympo_getMeshFVtxType())
-  call polympo_applyReconstruction(mpMesh)
-  call polympo_getMeshVtxMass(mpMesh,nVertices,c_loc(meshVtxMass))
-
-  do i = 1, nVertices
-    call assert(meshVtxMass(i) < TEST_VAL+TOLERANCE .and. meshVtxMass(i) > TEST_VAL-TOLERANCE, "Error: wrong vtx mass")
-    !write(*, *) 'The value of LRV is:', meshVtxMass(i)
-  end do
-  
+  !First Order reconstruction done before 0th order reconstruction as calculation of coefficients will
+  !fill the MpBasis slice
   !Test vtx order 1 reconstruction
+  call polympo_reconstruct_coeff_with_MPI(mpmesh)
   call polympo_setReconstructionOfMass(mpMesh,1,polympo_getMeshFVtxType())
   call polympo_setReconstructionOfVel(mpMesh, 1, polympo_getMeshFVtxType())
   call polympo_applyReconstruction(mpMesh)
   call polympo_getMeshVtxMass(mpMesh,nVertices,c_loc(meshVtxMass1))
   call polympo_getMeshVtxVel(mpMesh, nVertices, c_loc(meshVtxVelu), c_loc(meshVtxVelv))
   do i = 1, nVertices
     call assert(meshVtxMass1(i) < TEST_VAL+TOLERANCE1 .and. meshVtxMass1(i) > TEST_VAL-TOLERANCE1, "Error: wrong vtx mass order 1")
-    call assert(meshVtxVelu(i) < TEST_VAL+TOLERANCE1 .and. meshVtxVelu(i) > TEST_VAL-TOLERANCE1, "Error: wrong vtx velU order 1")
-    call assert(meshVtxVelv(i) < TEST_VAL+TOLERANCE1 .and. meshVtxVelv(i) > TEST_VAL-TOLERANCE1, "Error: wrong vtx velV order 1")
-    write(*, *) 'The value of LRV is:', meshVtxVelu(i)-TEST_VAL, meshVtxVelv(i)-TEST_VAL
+    call assert(meshVtxVelu(i)  < TEST_VAL+TOLERANCE1 .and. meshVtxVelu(i)  > TEST_VAL-TOLERANCE1, "Error: wrong vtx velU order 1")
+    call assert(meshVtxVelv(i)  < TEST_VAL+TOLERANCE1 .and. meshVtxVelv(i)  > TEST_VAL-TOLERANCE1, "Error: wrong vtx velV order 1")
   end do
 
+  ! Test vtx order 0  reconstruction
+  call polympo_setReconstructionOfMass(mpMesh,0,polympo_getMeshFVtxType())
+  call polympo_applyReconstruction(mpMesh)
+  call polympo_getMeshVtxMass(mpMesh,nVertices,c_loc(meshVtxMass))
+
+  do i = 1, nVertices
+    call assert(meshVtxMass(i) < TEST_VAL+TOLERANCE .and. meshVtxMass(i) > TEST_VAL-TOLERANCE, "Error: wrong vtx mass")
+  end do
+  
 
   ! Test push reconstruction
   do j = 1, 5