Simple Migration

src/pmpo_MPMesh.cpp

@@ -84,7 +84,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 +134,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;
@@ -258,17 +262,36 @@ void MPMesh::T2LTracking(Vec2dView dx){
     p_MPs->parallel_for(T2LCalc,"T2lTrackingCalc");
 }
 
+bool getAnyIsMigrating(bool isMigrating) {
+  int comm_rank;
+  MPI_Comm_rank(MPI_COMM_WORLD, &comm_rank);
+  int comm_size;
+  MPI_Comm_size(MPI_COMM_WORLD, &comm_size);
+
+  bool anyIsMigrating = isMigrating;
+  for (int i=0; i < comm_size; i++) {
+    if (i == comm_rank) MPI_Bcast(&anyIsMigrating, 1, MPI_C_BOOL, i, MPI_COMM_WORLD);
+    else MPI_Bcast(&isMigrating, 1, MPI_C_BOOL, i, MPI_COMM_WORLD);
+    anyIsMigrating |= isMigrating;
+  }
+  return anyIsMigrating;
+}
+
 void MPMesh::push(){
   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
+  while(true) {
+    CVTTrackingElmCenterBased(); // move to Tgt_XYZ
+    bool isMigrating = p_MPs->migrate();
+    bool anyIsMigrating = getAnyIsMigrating(isMigrating);
+    p_MPs->updateMPElmID(); //update mpElm IDs slices
+    if (!anyIsMigrating) break;
+  }
 
   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
 }
 
 void MPMesh::printVTP_mesh(int printVTPIndex){

src/pmpo_c.cpp

@@ -642,3 +642,17 @@ 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);
+}

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);

src/pmpo_fortran.f90

@@ -507,6 +507,19 @@ subroutine polympo_getMeshOnSurfDispIncr(mpMesh, nComps, nVertices, array) &
     type(c_ptr), value :: array
   end subroutine
   !---------------------------------------------------------------------------
+  !> @brief set the index to cell id array from a host 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)

src/pmpo_materialPoints.cpp

@@ -101,6 +101,28 @@ void MaterialPoints::finishRebuild() {
   rebuildFields.ongoing = false;
 }
 
+bool MaterialPoints::migrate() {
+  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);
+  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);

src/pmpo_mesh.cpp

@@ -55,4 +55,8 @@ namespace polyMPO{
         });
     }
 
+    IntView Mesh::getElm2Process() {
+        return owningProc_;
+    }
+
 } // namespace polyMPO

src/pmpo_mesh.hpp

@@ -66,7 +66,9 @@ class Mesh {
     //IntView nEdgesPerElm_;
     IntVtx2ElmView elm2VtxConn_;
     IntElm2ElmView elm2ElmConn_;
-
+
+    IntView owningProc_;
+
     //start of meshFields
     DoubleVec3dView vtxCoords_;
     DoubleSclrView vtxRotLat_;
@@ -103,6 +105,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 +133,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();
 };

test/CMakeLists.txt

@@ -48,6 +48,7 @@ pmpo_add_exe(testReconstruction testReconstruction.cpp)
 pmpo_add_exe(testTracking testTracking.cpp)
 pmpo_add_exe(testRebuild testRebuild.cpp)
 pmpo_add_exe(timeAssmblyWachspress testTiming.cpp)
+pmpo_add_exe(testMigration testMigration.cpp)
 pmpo_add_fortran_exe(testFortranInit testFortranInit.f90)
 pmpo_add_fortran_exe(testFortran testFortran.f90)
 pmpo_add_fortran_exe(testFortranReadMPAS testFortranReadMPAS.f90)
@@ -59,24 +60,24 @@ pmpo_add_fortran_exe(testFortranCreateRebuildMPs testFortranCreateRebuildMPs.f90
 #add test
 find_program(MPIRUN_EXECUTABLE NAMES mpirun) 
 
-function(pmpo_add_test testname command)
+function(pmpo_add_test testname procs command)
   if(MPIRUN_EXECUTABLE)
     add_test(NAME ${testname} 
-             COMMAND mpirun --bind-to core -np 1 ${command} ${ARGN})
+             COMMAND mpirun --bind-to core -np ${procs} ${command} ${ARGN})
   else()
     add_test(NAME ${testname} COMMAND ${command} ${ARGN})
   endif()
 endfunction()
 
-pmpo_add_test(unit_test ./unitTest)
-pmpo_add_test(mp_unit_test ./mpUnitTest)
-pmpo_add_test(test_reconstruction ./testReconstruction)
-pmpo_add_test(test_tracking ./testTracking)
-pmpo_add_test(test_rebuild ./testRebuild)
-pmpo_add_test(testFortranInit ./testFortranInit)
-pmpo_add_test(testFortran ./testFortran)
-pmpo_add_test(testFortranMPMeshModule ./testFortranMPMeshModule)
-pmpo_add_test(testFortranMPAppIDs ./testFortranMPAppIDs)
+pmpo_add_test(unit_test 1 ./unitTest)
+pmpo_add_test(mp_unit_test 1 ./mpUnitTest)
+pmpo_add_test(test_reconstruction 1 ./testReconstruction)
+pmpo_add_test(test_tracking 1 ./testTracking)
+pmpo_add_test(test_rebuild 1 ./testRebuild)
+pmpo_add_test(testFortranInit 1 ./testFortranInit)
+pmpo_add_test(testFortran 1 ./testFortran)
+pmpo_add_test(testFortranMPMeshModule 1 ./testFortranMPMeshModule)
+pmpo_add_test(testFortranMPAppIDs 1 ./testFortranMPAppIDs)
 
 #set NC file for test
 set(TEST_NC_FILE_PLANAR "${CMAKE_SOURCE_DIR}/test/sample_mpas_meshes/planar_nonuniform_cvt_for_square_673elms.nc")
@@ -86,12 +87,13 @@ if(EXISTS ${TEST_NC_FILE_SPHERICAL})
     message(STATUS "Found Spherical NC File: ${TEST_NC_FILE_SPHERICAL}") 
     set(TEST_NC_FILE ${TEST_NC_FILE_SPHERICAL})
     # this test only works for spherical NC files!
-    pmpo_add_test(testFortranMPAdvection ./testFortranMPAdvection ${TEST_NC_FILE}) 
+    pmpo_add_test(testFortranMPAdvection 1 ./testFortranMPAdvection ${TEST_NC_FILE}) 
 else()
     set(TEST_NC_FILE ${TEST_NC_FILE_PLANAR})
 endif()
-pmpo_add_test(testFortranReadMPAS ./testFortranReadMPAS ${TEST_NC_FILE})
-pmpo_add_test(testFortranCreateRebuildMPs ./testFortranCreateRebuildMPs ${TEST_NC_FILE})
+pmpo_add_test(testFortranReadMPAS 1 ./testFortranReadMPAS ${TEST_NC_FILE})
+pmpo_add_test(testFortranCreateRebuildMPs 1 ./testFortranCreateRebuildMPs ${TEST_NC_FILE})
 
-pmpo_add_test(test_timing ./timeAssmblyWachspress 1 ${TEST_NC_FILE})
-pmpo_add_test(test_wachspress ./testWachspress ${TEST_NC_FILE})
+pmpo_add_test(test_timing 1 ./timeAssmblyWachspress 1 ${TEST_NC_FILE})
+pmpo_add_test(test_wachspress 1 ./testWachspress ${TEST_NC_FILE})
+pmpo_add_test(test_migration 8 ./testMigration ${TEST_NC_FILE}) 

test/testMigration.cpp

@@ -0,0 +1,83 @@
+#include "pmpo_MPMesh.hpp"
+#include "testUtils.hpp"
+
+using namespace polyMPO;
+
+IntView procWedges(Mesh* mesh, int numElms, int comm_size) {
+    auto elm2VtxConn = mesh->getElm2VtxConn();
+    auto vtxLat = mesh->getMeshField<polyMPO::MeshF_VtxRotLat>(); 
+
+    DoubleView elmLat("elmLat",numElms);
+    DoubleView max("max", 1);
+    DoubleView min("min", 1);
+    Kokkos::parallel_for("calcElementLat", numElms, KOKKOS_LAMBDA(const int elm){  
+        int numVtx = elm2VtxConn(elm,0);
+        double sum = 0.0;
+        for(int i=1; i<= numVtx; i++){
+            sum += vtxLat(elm2VtxConn(elm,i)-1);
+        }
+        elmLat(elm) = sum/numVtx;
+        Kokkos::atomic_max(&max(0), elmLat(elm));
+        Kokkos::atomic_min(&min(0), elmLat(elm));
+    });
+
+    IntView owningProc("owningProc", numElms);
+    Kokkos::parallel_for("setOwningProc", numElms, KOKKOS_LAMBDA(const int elm){
+        double normalizedLat = (elmLat(elm) - min(0)) / (max(0) - min(0));
+        owningProc(elm) = normalizedLat * comm_size;
+    });
+    return owningProc;
+}
+
+int main(int argc, char* argv[] ) {
+    MPI_Init(&argc, &argv);
+    Kokkos::initialize(argc,argv);
+
+    int comm_rank;
+    MPI_Comm_rank(MPI_COMM_WORLD, &comm_rank);
+    int comm_size;
+    MPI_Comm_size(MPI_COMM_WORLD, &comm_size);
+
+    {
+        MPMesh* mpmesh = NULL;
+        void* p;
+        setWithMPASMeshByFortran(&p, argv[1], (int)strlen(argv[1]));
+        mpmesh = (MPMesh*)p;
+        Mesh* mesh = mpmesh->p_mesh;
+        int numElms = mesh->getNumElements();
+        int numMPs = numElms;
+
+        IntView mpsPerElm("mpsPerElm", numElms);
+        IntView activeMP2Elm("activeMP2Elm", numMPs);
+        IntView activeMPIDs("activeMPIDs", numMPs);
+
+        Kokkos::parallel_for("setMPsPerElm", numElms, KOKKOS_LAMBDA(const int elm){
+            mpsPerElm(elm) = 1;
+        });
+        Kokkos::parallel_for("setMPs", numMPs, KOKKOS_LAMBDA(const int mp){
+            activeMP2Elm(mp) = mp;
+            activeMPIDs(mp) = mp;
+        });
+
+        MaterialPoints p_MPs(numElms, numMPs, mpsPerElm, activeMP2Elm, activeMPIDs);
+        mpmesh->p_MPs = &p_MPs;
+
+        IntView owningProc = procWedges(mesh, numElms, comm_size);
+        mesh->setMeshEdit(true);
+        mesh->setOwningProc(owningProc);
+        mpmesh->push();
+
+        auto MPs2Proc = p_MPs.getData<MPF_Tgt_Proc_ID>();
+        auto checkPostBackMigrate = PS_LAMBDA(const int& e, const int& mp, const bool& mask) {
+            if (mask) {
+                assert(MPs2Proc(mp) == comm_rank);
+                assert(owningProc(e) == comm_rank);
+            }
+        };
+        p_MPs.parallel_for(checkPostBackMigrate, "checkPostBackMigrate");
+    } 
+
+    Kokkos::finalize();
+    MPI_Finalize();
+    return 0;
+}