Simple Migration

src/pmpo_MPMesh.cpp

@@ -74,6 +74,7 @@ void MPMesh::CVTTrackingEdgeCenterBased(Vec2dView dx){
 
 
 void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){
+    Kokkos::Timer timer;
     int numVtxs = p_mesh->getNumVertices();
     int numElms = p_mesh->getNumElements();
     auto numMPs = p_MPs->getCount();
@@ -84,7 +85,9 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){
 
     auto mpPositions = p_MPs->getData<MPF_Cur_Pos_XYZ>();
     auto mpTgtPos = p_MPs->getData<MPF_Tgt_Pos_XYZ>();
-    auto MPs2Elm = p_MPs->getData<MPF_Tgt_Elm_ID>();;
+    auto MPs2Elm = p_MPs->getData<MPF_Tgt_Elm_ID>();
+    auto MPs2Proc = p_MPs->getData<MPF_Tgt_Proc_ID>();
+    auto elm2Process = p_mesh->getElm2Process();
 
     if(printVTPIndex>=0) {
       printVTP_mesh(printVTPIndex);
@@ -132,6 +135,8 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){
                 }
                 if(closestElm<0){
                     MPs2Elm(mp) = iElm;
+                    if (elm2Process.size() > 0) 
+                        MPs2Proc(mp) = elm2Process(iElm);
                     break;
                 }else{
                     iElm = closestElm;
@@ -193,6 +198,7 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){
         fprintf(pFile,"        </DataArray>\n      </Lines>\n    </Piece>\n  </PolyData>\n</VTKFile>\n");
         fclose(pFile);
     }
+    pumipic::RecordTime("PolyMPO_CVTTrackingElmCenterBased", timer.seconds());
 }
 
 void MPMesh::T2LTracking(Vec2dView dx){
@@ -258,17 +264,37 @@ void MPMesh::T2LTracking(Vec2dView dx){
     p_MPs->parallel_for(T2LCalc,"T2lTrackingCalc");
 }
 
+bool getAnyIsMigrating(bool isMigrating) {
+  Kokkos::Timer timer;
+  int comm_rank;
+  MPI_Comm_rank(MPI_COMM_WORLD, &comm_rank);
+  int comm_size;
+  MPI_Comm_size(MPI_COMM_WORLD, &comm_size);
+
+  bool anyIsMigrating = false;
+  MPI_Allreduce(&isMigrating, &anyIsMigrating, 1, MPI_C_BOOL, MPI_LOR, MPI_COMM_WORLD);
+  pumipic::RecordTime("PolyMPO_getAnyIsMigrating", timer.seconds());
+  return anyIsMigrating;
+}
+
 void MPMesh::push(){
+  Kokkos::Timer timer;
   p_mesh->computeRotLatLonIncr();
   sphericalInterpolation<MeshF_RotLatLonIncr, MPF_Rot_Lat_Lon_Incr>(*this);
   p_MPs->updateRotLatLonAndXYZ2Tgt(p_mesh->getSphereRadius()); // set Tgt_XYZ
 
-  CVTTrackingElmCenterBased(); // move to Tgt_XYZ
+  bool anyIsMigrating = false;
+  do {
+    CVTTrackingElmCenterBased(); // move to Tgt_XYZ
+    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
+    p_MPs->updateMPElmID(); //update mpElm IDs slices
+    bool isMigrating = p_MPs->migrate();
+    anyIsMigrating = getAnyIsMigrating(isMigrating);
+  } 
+  while (anyIsMigrating);
 
-  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
-  p_MPs->rebuild(); //rebuild pumi-pic
-  p_MPs->updateMPElmID(); //update mpElm IDs slices
+  pumipic::RecordTime("PolyMPO_push", timer.seconds());
 }
 
 void MPMesh::printVTP_mesh(int printVTPIndex){

src/pmpo_c.cpp

@@ -595,7 +595,7 @@ void polympo_getMeshOnSurfVeloIncr_f(MPMesh_ptr p_mpmesh, const int nComps, cons
   Kokkos::deep_copy(arrayHost, vtxField);
 }
 
-void polympo_setMeshOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, const double* array) {
+void polympo_setMeshVtxOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, const double* array) {
   //check mpMesh is valid
   checkMPMeshValid(p_mpmesh);
   auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
@@ -642,3 +642,29 @@ void polympo_push_f(MPMesh_ptr p_mpmesh){
   checkMPMeshValid(p_mpmesh);
   ((polyMPO::MPMesh*)p_mpmesh) ->push();
 }
+
+void polympo_setOwningProc_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());
+  kkViewHostU<const int*> arrayHost(array,nCells); 
+
+  //check the size
+  PMT_ALWAYS_ASSERT(nCells == p_mesh->getNumElements());
+
+  Kokkos::View<int*> owningProc("owningProc",nCells);
+  Kokkos::deep_copy(owningProc, arrayHost);
+  p_mesh->setOwningProc(owningProc);
+}
+
+void polympo_enableTiming_f(){
+  pumipic::EnableTiming();
+}
+
+void polympo_summarizeTime_f(){
+  pumipic::SummarizeTime();
+}
+
+void polympo_setTimingVerbosity_f(int v){
+  pumipic::SetTimingVerbosity(v);
+}

src/pmpo_c.h

@@ -50,6 +50,7 @@ void polympo_getMeshNumElms_f(MPMesh_ptr p_mpmesh, int & numElms);
 void polympo_setMeshNumEdgesPerElm_f(MPMesh_ptr p_mpmesh, const int nCells, const int* array);
 void polympo_setMeshElm2VtxConn_f(MPMesh_ptr p_mpmesh, const int maxEdges, const int nCells, const int* array);
 void polympo_setMeshElm2ElmConn_f(MPMesh_ptr p_mpmesh, const int maxEdges, const int nCells, const int* array);
+void polympo_setOwningProc_f(MPMesh_ptr p_mpmesh, const int nCells, const int* array);
 
 //Mesh fields
 void polympo_setMeshVtxCoords_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* xArray, const double* yArray, const double* zArray);
@@ -58,10 +59,15 @@ void polympo_setMeshVtxRotLat_f(MPMesh_ptr p_mpmesh, const int nVertices, const
 void polympo_getMeshVtxRotLat_f(MPMesh_ptr p_mpmesh, const int nVertices, double* latitude);
 void polympo_setMeshOnSurfVeloIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, const double* array);//vec2d
 void polympo_getMeshOnSurfVeloIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, double* array);//vec2d
-void polympo_setMeshOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, const double* array);//vec2d
+void polympo_setMeshVtxOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, const double* array);//vec2d
 void polympo_getMeshOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, double* array);//vec2d
 
 // calculations
 void polympo_push_f(MPMesh_ptr p_mpmesh);
+
+// Timing
+void polympo_enableTiming_f();
+void polympo_summarizeTime_f();
+void polympo_setTimingVerbosity_f(int v);
 }
 #endif

src/pmpo_fortran.f90

@@ -483,8 +483,8 @@ subroutine polympo_getMeshOnSurfVeloIncr(mpMesh, nComps, nVertices, array) &
   !> @param nVertices(in) numVertices
   !> @param array(in) input mesh velocity 2D array (2,numVtx)
   !---------------------------------------------------------------------------
-  subroutine polympo_setMeshOnSurfDispIncr(mpMesh, nComps, nVertices, array) &
-             bind(C, NAME='polympo_setMeshOnSurfDispIncr_f')
+  subroutine polympo_setMeshVtxOnSurfDispIncr(mpMesh, nComps, nVertices, array) &
+             bind(C, NAME='polympo_setMeshVtxOnSurfDispIncr_f')
     use :: iso_c_binding
     type(c_ptr), value :: mpMesh
     integer(c_int), value :: nComps, nVertices
@@ -507,6 +507,19 @@ subroutine polympo_getMeshOnSurfDispIncr(mpMesh, nComps, nVertices, array) &
     type(c_ptr), value :: array
   end subroutine
   !---------------------------------------------------------------------------
+  !> @brief set the owning process array
+  !> @param mpmesh(in/out) MPMesh object
+  !> @param nCells(in) number of cells
+  !> @param array(in) input mesh cell to process array
+  !---------------------------------------------------------------------------
+  subroutine polympo_setOwningProc(mpMesh, nCells, array) &
+    bind(C, NAME='polympo_setOwningProc_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: nCells
+    type(c_ptr), intent(in), value :: array
+  end subroutine
+  !---------------------------------------------------------------------------
   !> @brief calculate the MPs from given mesh vertices rotational latitude
   !>        longitude, update the MP slices
   !>        MPs MUST have rotated flag set to True(>0)
@@ -517,6 +530,25 @@ subroutine polympo_push(mpMesh) &
     use :: iso_c_binding
     type(c_ptr), value :: mpMesh
   end subroutine
+  !---------------------------------------------------------------------------
+  !> @brief enable timing functions
+  !---------------------------------------------------------------------------
+  subroutine polympo_enableTiming() bind(C, NAME='polympo_enableTiming_f')
+    use :: iso_c_binding
+  end subroutine
+  !---------------------------------------------------------------------------
+  !> @brief print summary of timings collected
+  !---------------------------------------------------------------------------
+  subroutine polympo_summarizeTime() bind(C, NAME='polympo_summarizeTime_f')
+    use :: iso_c_binding
+  end subroutine
+  !---------------------------------------------------------------------------
+  !> @brief set verbosity of timings collected
+  !---------------------------------------------------------------------------
+  subroutine polympo_setTimingVerbosity(v) bind(C, NAME='polympo_setTimingVerbosity_f')
+    use :: iso_c_binding
+    integer(c_int), value :: v
+  end subroutine
   end interface
   contains
   !---------------------------------------------------------------------------

src/pmpo_materialPoints.cpp

@@ -101,6 +101,33 @@ void MaterialPoints::finishRebuild() {
   rebuildFields.ongoing = false;
 }
 
+bool MaterialPoints::migrate() {
+  Kokkos::Timer timer;
+  auto MPs2Elm = getData<MPF_Tgt_Elm_ID>();
+  auto MPs2Proc = getData<MPF_Tgt_Proc_ID>();
+
+  IntView new_elem("new_elem", MPs->capacity());
+  IntView new_process("new_process", MPs->capacity());
+  IntView isMigrating("isMigrating", 1);
+
+  int rank;
+  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+  auto setMigrationFields = PS_LAMBDA(const int& e, const int& mp, const bool& mask) {
+    if (mask) {
+      new_elem(mp) = MPs2Elm(mp);
+      new_process(mp) = MPs2Proc(mp);
+      if (new_process(mp) != rank) isMigrating(0) = 1;
+    }
+  };
+  parallel_for(setMigrationFields, "setMigrationFields");
+  MPs->migrate(new_elem, new_process);
+
+  if (getOpMode() == polyMPO::MP_DEBUG)
+    printf("Material point migration: %f\n", timer.seconds());
+  pumipic::RecordTime("PolyMPO_migrate", timer.seconds());
+  return pumipic::getLastValue(isMigrating) > 0;
+}
+
 bool MaterialPoints::rebuildOngoing() { return rebuildFields.ongoing; }
 
 void MaterialPoints::setAppIDFunc(IntFunc getAppIDIn) { getAppID = getAppIDIn; }

src/pmpo_materialPoints.hpp

@@ -21,6 +21,7 @@ using defaultSpace = Kokkos::DefaultExecutionSpace::memory_space;
 typedef int mp_flag_t;
 typedef int mp_id_t;
 typedef int  mp_elm_id_t;
+typedef int mp_proc_id_t;
 typedef double mp_sclr_t[1];//TODO
 typedef vec2d_t mp_vec2d_t;
 typedef double mp_vec3d_t[3];//TODO
@@ -49,7 +50,8 @@ enum MaterialPointSlice {
   MPF_Stress_Div,
   MPF_Shear_Traction,
   MPF_Constv_Mdl_Param,
-  MPF_MP_APP_ID
+  MPF_MP_APP_ID,
+  MPF_Tgt_Proc_ID
 };
 
 enum Operating_Mode{
@@ -76,7 +78,8 @@ const static std::map<MaterialPointSlice, std::pair<int,MeshFieldIndex>>
                            {MPF_Stress_Div,      {3,MeshF_Unsupported}},
                            {MPF_Shear_Traction,  {3,MeshF_Unsupported}},
                            {MPF_Constv_Mdl_Param,{12,MeshF_Unsupported}},
-                           {MPF_MP_APP_ID,       {1,MeshF_Invalid}}};
+                           {MPF_MP_APP_ID,       {1,MeshF_Invalid}},
+                           {MPF_Tgt_Proc_ID,     {0,MeshF_Invalid}}}; //TODO: What does this integer mean?
 
 const static std::vector<std::pair<MaterialPointSlice, MaterialPointSlice>>
         mpSliceSwap = {{MPF_Cur_Elm_ID, MPF_Tgt_Elm_ID},
@@ -101,7 +104,8 @@ typedef MemberTypes<mp_flag_t,              //MP_Status
                     mp_vec3d_t,             //MP_Stress_Div
                     mp_vec3d_t,             //MP_Shear_Traction
                     mp_constv_mdl_param_t,  //MP_Constv_Mdl_Param
-                    mp_id_t                 //MP_APP_ID
+                    mp_id_t,                //MP_APP_ID
+                    mp_proc_id_t            //MPF_Tgt_Proc_ID
                     >MaterialPointTypes;
 typedef ps::ParticleStructure<MaterialPointTypes> PS;
 
@@ -134,7 +138,9 @@ class MaterialPoints {
     void startRebuild(IntView tgtElm, int addedNumMPs, IntView addedMP2elm, IntView addedMPAppID, Kokkos::View<const int*> addedMPMask);
     void finishRebuild();
     bool rebuildOngoing();
-
+
+    bool migrate();
+
     template<int mpSliceIndex, typename mpSliceData>
     typename std::enable_if<mpSliceData::rank==1>::type
     setRebuildMPSlice(mpSliceData mpSliceIn);
@@ -202,6 +208,7 @@ class MaterialPoints {
         updateMPSlice<MPF_Cur_Pos_XYZ,MPF_Tgt_Pos_XYZ>();
     }
     void updateRotLatLonAndXYZ2Tgt(const double radius){
+        Kokkos::Timer timer;
         auto curPosRotLatLon = MPs->get<MPF_Cur_Pos_Rot_Lat_Lon>();
         auto tgtPosRotLatLon = MPs->get<MPF_Tgt_Pos_Rot_Lat_Lon>();
         auto tgtPosXYZ = MPs->get<MPF_Tgt_Pos_XYZ>();
@@ -227,6 +234,7 @@ class MaterialPoints {
             PMT_ALWAYS_ASSERT(false);
         } 
         ps::parallel_for(MPs, updateRotLatLon,"updateRotationalLatitudeLongitude"); 
+        pumipic::RecordTime("PolyMPO_updateRotLatLonAndXYZ2Tgt", timer.seconds()); 
     } 
 
     template <int index>

src/pmpo_mesh.cpp

@@ -40,6 +40,7 @@ namespace polyMPO{
     }
 
     void Mesh::computeRotLatLonIncr(){
+        Kokkos::Timer timer;
         PMT_ALWAYS_ASSERT(geomType_ == geom_spherical_surf);
 
         auto dispIncr = getMeshField<MeshF_OnSurfDispIncr>();
@@ -53,6 +54,11 @@ namespace polyMPO{
             rotLatLonIncr(iVtx, 0) = dispIncr(iVtx, 1)/sphereRadius;
             rotLatLonIncr(iVtx, 1) = dispIncr(iVtx, 0)/(sphereRadius * std::cos(lat(iVtx)));
         });
+        pumipic::RecordTime("PolyMPO_computeRotLatLonIncr", timer.seconds());
+    }
+
+    IntView Mesh::getElm2Process() {
+        return owningProc_;
     }
 
 } // namespace polyMPO

src/pmpo_mesh.hpp

@@ -2,6 +2,7 @@
 #define POLYMPO_MESH_H
 
 #include "pmpo_utils.hpp"
+#include <pumipic_kktypes.hpp>
 
 namespace polyMPO{
 
@@ -66,7 +67,9 @@ class Mesh {
     //IntView nEdgesPerElm_;
     IntVtx2ElmView elm2VtxConn_;
     IntElm2ElmView elm2ElmConn_;
-
+
+    IntView owningProc_;
+
     //start of meshFields
     DoubleVec3dView vtxCoords_;
     DoubleSclrView vtxRotLat_;
@@ -103,6 +106,8 @@ class Mesh {
     bool checkMeshType(int meshType);
     bool checkGeomType(int geomType);
 
+    IntView getElm2Process();
+
     mesh_type getMeshType() { return meshType_; }
     geom_type getGeomType() { return geomType_; }
     double getSphereRadius() { return sphereRadius_; }
@@ -129,6 +134,8 @@ class Mesh {
                                                      elm2VtxConn_ = elm2VtxConn; }
     void setElm2ElmConn(IntElm2ElmView elm2ElmConn) {PMT_ALWAYS_ASSERT(meshEdit_);
                                                      elm2ElmConn_ = elm2ElmConn; }
+    void setOwningProc(IntView owningProc) {PMT_ALWAYS_ASSERT(meshEdit_);
+                                            owningProc_ = owningProc; }
 
     void computeRotLatLonIncr();
 };

src/pmpo_wachspressBasis.hpp

@@ -306,6 +306,7 @@ void getBasisByAreaGblForm_1(Vec2d MP, int numVtxs, Vec2d* vtxCoords, double* ba
 // spherical interpolation of values from mesh vertices to MPsi
 template <MeshFieldIndex mfIndex, MaterialPointSlice mpfIndex>
 void sphericalInterpolation(MPMesh& mpMesh){
+    Kokkos::Timer timer;
     auto p_mesh = mpMesh.p_mesh;
     auto vtxCoords = p_mesh->getMeshField<polyMPO::MeshF_VtxCoords>();
     int numVtxs = p_mesh->getNumVertices();
@@ -355,6 +356,7 @@ void sphericalInterpolation(MPMesh& mpMesh){
         }
     };
     p_MPs->parallel_for(interpolation, "interpolation");
+    pumipic::RecordTime("PolyMPO_sphericalInterpolation", timer.seconds());
 }
 
 } //namespace polyMPO end