From a591f1136df62afc658d208b366e3e2861204a77 Mon Sep 17 00:00:00 2001 From: dhyan1272 Date: Tue, 25 Mar 2025 10:36:15 -0700 Subject: [PATCH 01/17] Multi process debugging --- src/pmpo_MPMesh.cpp | 32 +++++++++++++++++++++++++++++++- src/pmpo_MPMesh_assembly.hpp | 6 ++++-- src/pmpo_c.cpp | 33 ++++++++++++++++++++++++--------- 3 files changed, 59 insertions(+), 12 deletions(-) diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index a5f709a..3c5fec3 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -127,7 +127,12 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){ auto MPs2Elm = p_MPs->getData(); auto MPs2Proc = p_MPs->getData(); auto elm2Process = p_mesh->getElm2Process(); - + + Kokkos::parallel_for("countProcess", numElms, KOKKOS_LAMBDA(const int iElm){ + int pp_id=elm2Process(iElm); + printf("Mesh elm %d owning element %d \n", iElm, pp_id); + }); + if(printVTPIndex>=0) { printVTP_mesh(printVTPIndex); } @@ -324,15 +329,37 @@ bool getAnyIsMigrating(MaterialPoints* p_MPs, bool isMigrating) { } void MPMesh::push(){ + static int count=0; + std::cout<<__FUNCTION__<<" "<1) exit(1); Kokkos::Timer timer; p_mesh->computeRotLatLonIncr(); + + assert(cudaDeviceSynchronize() == cudaSuccess); + MPI_Barrier(MPI_COMM_WORLD); + printf("FooPush\n"); + sphericalInterpolation(*this); + assert(cudaDeviceSynchronize() == cudaSuccess); + MPI_Barrier(MPI_COMM_WORLD); + printf("FooPush \n"); + p_MPs->updateRotLatLonAndXYZ2Tgt(p_mesh->getSphereRadius()); // set Tgt_XYZ + assert(cudaDeviceSynchronize() == cudaSuccess); + MPI_Barrier(MPI_COMM_WORLD); + printf("FooPush\n"); + auto elm2Process = p_mesh->getElm2Process(); + assert(cudaDeviceSynchronize() == cudaSuccess); + MPI_Barrier(MPI_COMM_WORLD); + printf("FooPush \n"); bool anyIsMigrating = false; do { CVTTrackingElmCenterBased(); // move to Tgt_XYZ + assert(cudaDeviceSynchronize() == cudaSuccess); + printf("FooPushInside\n"); p_MPs->updateMPSlice(); // Tgt_XYZ becomes Cur_XYZ p_MPs->updateMPSlice(); // Tgt becomes Cur if (elm2Process.size() > 0) @@ -343,6 +370,9 @@ void MPMesh::push(){ reconstructSlices(); } while (anyIsMigrating); + assert(cudaDeviceSynchronize() == cudaSuccess); + MPI_Barrier(MPI_COMM_WORLD); + printf("FooPush\n"); pumipic::RecordTime("PolyMPO_push", timer.seconds()); } diff --git a/src/pmpo_MPMesh_assembly.hpp b/src/pmpo_MPMesh_assembly.hpp index 25899b9..69e02e6 100644 --- a/src/pmpo_MPMesh_assembly.hpp +++ b/src/pmpo_MPMesh_assembly.hpp @@ -104,7 +104,7 @@ void MPMesh::resetPreComputeFlag(){ } void MPMesh::computeMatricesAndSolve(){ - + Kokkos::Timer timer; //Mesh Information auto elm2VtxConn = p_mesh->getElm2VtxConn(); int numVtx = p_mesh->getNumVertices(); @@ -205,11 +205,12 @@ void MPMesh::computeMatricesAndSolve(){ VtxCoeffs(vtx,i)=coeff[i]; }); this->precomputedVtxCoeffs = VtxCoeffs; + pumipic::RecordTime("PolyMPO_Calculate_MLS_Coeff", timer.seconds()); } template void MPMesh::assemblyVtx1() { - + Kokkos::Timer timer; //If no reconstruction till now calculate the coeffs if (!isPreComputed) { computeMatricesAndSolve(); @@ -270,6 +271,7 @@ void MPMesh::assemblyVtx1() { for(int k=0; k diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index 67653b5..57a07ad 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -61,17 +61,19 @@ void polympo_setMPICommunicator_f(MPMesh_ptr p_mpmesh, MPI_Fint fcomm){ void polympo_createMPs_f(MPMesh_ptr p_mpmesh, const int numElms, - const int numMPs, // total number of MPs which is GREATER than or equal to number of active MPs + const int numMPs, // total nof of MPs which is >= no of active MPs int* mpsPerElm, const int* mp2Elm, const int* isMPActive) { checkMPMeshValid(p_mpmesh); - + std::cout<<__FUNCTION__<p_mesh; PMT_ALWAYS_ASSERT(!p_mesh->meshEditable()); PMT_ALWAYS_ASSERT(p_mesh->getNumElements() == numElms); + //Find the total no of MPs across all ranks + //And loop over all MPs and find the smallest element id associated across a MP int numActiveMPs = 0; int minElmID = numElms+1; for(int i = 0; i < numMPs; i++) { @@ -82,20 +84,27 @@ void polympo_createMPs_f(MPMesh_ptr p_mpmesh, } } } - //TODO do we care about empty ranks? check just in case... - PMT_ALWAYS_ASSERT(numActiveMPs>0); - - int firstElmWithMPs=-1; + long long globalNumActiveMPs = 0; + int globalMinElmID; + MPI_Allreduce(&numActiveMPs, &globalNumActiveMPs, 1, MPI_LONG_LONG_INT, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(&minElmID, &globalMinElmID, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD); + PMT_ALWAYS_ASSERT(globalNumActiveMPs>0); + + //Loop over all mesh elements 0,1,... and find the first element that has an associated MP + int firstElmWithMPs=numElms+1; for (int i=0; ip_MPs; ((polyMPO::MPMesh*)p_mpmesh)->p_MPs = new polyMPO::MaterialPoints(numElms, numActiveMPs, mpsPerElm_d, active_mp2Elm_d, active_mpIDs_d); + auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; p_MPs->setElmIDoffset(offset); + + assert(cudaDeviceSynchronize() == cudaSuccess); + MPI_Barrier(MPI_COMM_WORLD); + printf("Foo1\n"); + } void polympo_startRebuildMPs_f(MPMesh_ptr p_mpmesh, From 7be04f94d08b42554986d4ce7738eff801049647 Mon Sep 17 00:00:00 2001 From: dhyan1272 Date: Thu, 27 Mar 2025 12:41:35 -0700 Subject: [PATCH 02/17] More debugging --- src/pmpo_MPMesh.cpp | 8 ++++++-- src/pmpo_c.cpp | 4 +++- 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index 3c5fec3..d213649 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -127,7 +127,11 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){ auto MPs2Elm = p_MPs->getData(); auto MPs2Proc = p_MPs->getData(); auto elm2Process = p_mesh->getElm2Process(); - + + assert(cudaDeviceSynchronize() == cudaSuccess); + MPI_Barrier(MPI_COMM_WORLD); + printf("FooTracking \n"); + Kokkos::parallel_for("countProcess", numElms, KOKKOS_LAMBDA(const int iElm){ int pp_id=elm2Process(iElm); printf("Mesh elm %d owning element %d \n", iElm, pp_id); @@ -332,7 +336,7 @@ void MPMesh::push(){ static int count=0; std::cout<<__FUNCTION__<<" "<1) exit(1); + //if(count>1) exit(1); Kokkos::Timer timer; p_mesh->computeRotLatLonIncr(); diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index 57a07ad..406e88e 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -70,7 +70,7 @@ void polympo_createMPs_f(MPMesh_ptr p_mpmesh, //the mesh must be fixed/set before adding MPs auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; PMT_ALWAYS_ASSERT(!p_mesh->meshEditable()); - PMT_ALWAYS_ASSERT(p_mesh->getNumElements() == numElms); + //PMT_ALWAYS_ASSERT(p_mesh->getNumElements() == numElms); //Find the total no of MPs across all ranks //And loop over all MPs and find the smallest element id associated across a MP @@ -84,6 +84,7 @@ void polympo_createMPs_f(MPMesh_ptr p_mpmesh, } } } + printf("Num Active MPs and minElmId %d %d\n", numActiveMPs, minElmID); long long globalNumActiveMPs = 0; int globalMinElmID; MPI_Allreduce(&numActiveMPs, &globalNumActiveMPs, 1, MPI_LONG_LONG_INT, MPI_SUM, MPI_COMM_WORLD); @@ -98,6 +99,7 @@ void polympo_createMPs_f(MPMesh_ptr p_mpmesh, break; } } + printf("First elem with MP %d\n", firstElmWithMPs); int globalFirstElmWithMPs; MPI_Allreduce(&firstElmWithMPs, &globalFirstElmWithMPs, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD); From fefd5e58295fa96e10633f4b77ef87e8cd36f48f Mon Sep 17 00:00:00 2001 From: nathd2 Date: Mon, 31 Mar 2025 11:05:13 -0400 Subject: [PATCH 03/17] Added global ids to polyMPO from MPAS --- src/pmpo_MPMesh.cpp | 40 +++++++++++++-------------------- src/pmpo_c.cpp | 45 ++++++++++++++++++++++++------------- src/pmpo_c.h | 2 ++ src/pmpo_fortran.f90 | 15 +++++++++++++ src/pmpo_materialPoints.cpp | 10 ++++++--- src/pmpo_materialPoints.hpp | 2 +- src/pmpo_mesh.cpp | 6 +++++ src/pmpo_mesh.hpp | 6 ++++- 8 files changed, 81 insertions(+), 45 deletions(-) diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index d213649..b282231 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -127,15 +127,21 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){ auto MPs2Elm = p_MPs->getData(); auto MPs2Proc = p_MPs->getData(); auto elm2Process = p_mesh->getElm2Process(); - + auto elm2global = p_mesh->getElmGlobal(); + + + MPI_Comm comm = p_MPs->getMPIComm(); + int comm_rank; + MPI_Comm_rank(comm, &comm_rank); assert(cudaDeviceSynchronize() == cudaSuccess); MPI_Barrier(MPI_COMM_WORLD); printf("FooTracking \n"); - Kokkos::parallel_for("countProcess", numElms, KOKKOS_LAMBDA(const int iElm){ int pp_id=elm2Process(iElm); - printf("Mesh elm %d owning element %d \n", iElm, pp_id); + if(iElm<3 || iElm==1470 || iElm==1471 || iElm==1472 ) + printf("Mesh elm %d Process %d owning proc %d global %d \n", iElm, comm_rank, pp_id, elm2global(iElm)); }); + if(printVTPIndex>=0) { printVTP_mesh(printVTPIndex); @@ -333,51 +339,37 @@ bool getAnyIsMigrating(MaterialPoints* p_MPs, bool isMigrating) { } void MPMesh::push(){ + static int count=0; std::cout<<__FUNCTION__<<" "<1) exit(1); Kokkos::Timer timer; + p_mesh->computeRotLatLonIncr(); - assert(cudaDeviceSynchronize() == cudaSuccess); - MPI_Barrier(MPI_COMM_WORLD); - printf("FooPush\n"); - sphericalInterpolation(*this); - assert(cudaDeviceSynchronize() == cudaSuccess); - MPI_Barrier(MPI_COMM_WORLD); - printf("FooPush \n"); - + p_MPs->updateRotLatLonAndXYZ2Tgt(p_mesh->getSphereRadius()); // set Tgt_XYZ - assert(cudaDeviceSynchronize() == cudaSuccess); - MPI_Barrier(MPI_COMM_WORLD); - printf("FooPush\n"); - + auto elm2Process = p_mesh->getElm2Process(); - assert(cudaDeviceSynchronize() == cudaSuccess); - MPI_Barrier(MPI_COMM_WORLD); - printf("FooPush \n"); bool anyIsMigrating = false; do { CVTTrackingElmCenterBased(); // move to Tgt_XYZ assert(cudaDeviceSynchronize() == cudaSuccess); - printf("FooPushInside\n"); + printf("FooPushDoWhile\n"); p_MPs->updateMPSlice(); // Tgt_XYZ becomes Cur_XYZ p_MPs->updateMPSlice(); // Tgt becomes Cur if (elm2Process.size() > 0) anyIsMigrating = getAnyIsMigrating(p_MPs, p_MPs->migrate()); else p_MPs->rebuild(); //rebuild pumi-pic + printf("Is migrating %d \n", anyIsMigrating); p_MPs->updateMPElmID(); //update mpElm IDs slices reconstructSlices(); } while (anyIsMigrating); - assert(cudaDeviceSynchronize() == cudaSuccess); - MPI_Barrier(MPI_COMM_WORLD); - printf("FooPush\n"); - + pumipic::RecordTime("PolyMPO_push", timer.seconds()); } diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index 406e88e..b60a63a 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -75,34 +75,33 @@ void polympo_createMPs_f(MPMesh_ptr p_mpmesh, //Find the total no of MPs across all ranks //And loop over all MPs and find the smallest element id associated across a MP int numActiveMPs = 0; - int minElmID = numElms+1; + int minElmID = INT_MAX; for(int i = 0; i < numMPs; i++) { if(isMPActive[i] == MP_ACTIVE) { - if(mp2Elm[i] < minElmID) { + numActiveMPs++; + if(mp2Elm[i] < minElmID) minElmID = mp2Elm[i]; - numActiveMPs++; - } } } - printf("Num Active MPs and minElmId %d %d\n", numActiveMPs, minElmID); - long long globalNumActiveMPs = 0; + int globalNumActiveMPs = 0; int globalMinElmID; - MPI_Allreduce(&numActiveMPs, &globalNumActiveMPs, 1, MPI_LONG_LONG_INT, MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(&minElmID, &globalMinElmID, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD); + MPI_Allreduce(&numActiveMPs, &globalNumActiveMPs, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(&minElmID, &globalMinElmID, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD); PMT_ALWAYS_ASSERT(globalNumActiveMPs>0); //Loop over all mesh elements 0,1,... and find the first element that has an associated MP - int firstElmWithMPs=numElms+1; + int firstElmWithMPs=INT_MAX; for (int i=0; igetElmGlobal(); auto mpsPerElm_d = create_mirror_view_and_copy(mpsPerElm, numElms); auto active_mp2Elm_d = create_mirror_view_and_copy(active_mp2Elm.data(), numActiveMPs); @@ -130,19 +130,15 @@ void polympo_createMPs_f(MPMesh_ptr p_mpmesh, delete ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; ((polyMPO::MPMesh*)p_mpmesh)->p_MPs = - new polyMPO::MaterialPoints(numElms, numActiveMPs, mpsPerElm_d, active_mp2Elm_d, active_mpIDs_d); + new polyMPO::MaterialPoints(numElms, numActiveMPs, mpsPerElm_d, active_mp2Elm_d, active_mpIDs_d, elm2global); auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; p_MPs->setElmIDoffset(offset); - assert(cudaDeviceSynchronize() == cudaSuccess); - MPI_Barrier(MPI_COMM_WORLD); - printf("Foo1\n"); - } void polympo_startRebuildMPs_f(MPMesh_ptr p_mpmesh, - const int numMPs, // total number of MPs which is GREATER than or equal to number of active MPs + const int numMPs, // Total # MPs which is GREATER than or equal to number of active MPs const int* allMP2Elm, const int* addedMPMask) { checkMPMeshValid(p_mpmesh); @@ -1046,6 +1042,23 @@ void polympo_setOwningProc_f(MPMesh_ptr p_mpmesh, const int nCells, const int* a p_mesh->setOwningProc(owningProc); } +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 arrayHost("arrayHost", nCells); + for (int i = 0; i < nCells; i++) { + arrayHost(i) = array[i] - 1; // Decrease each value by 1 + } + //check the size + PMT_ALWAYS_ASSERT(nCells == p_mesh->getNumElements()); + + Kokkos::View elmGlobal("elmGlobal",nCells); + Kokkos::deep_copy(elmGlobal, arrayHost); + p_mesh->setElmGlobal(elmGlobal); +} + + void polympo_enableTiming_f(){ pumipic::EnableTiming(); } diff --git a/src/pmpo_c.h b/src/pmpo_c.h index 174f862..f62ffdc 100644 --- a/src/pmpo_c.h +++ b/src/pmpo_c.h @@ -57,6 +57,8 @@ void polympo_setMeshNumEdgesPerElm_f(MPMesh_ptr p_mpmesh, const int nCells, cons 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); +void polympo_setElmGlobal_f(MPMesh_ptr p_mpmesh, const int nCells, const int* array); + //Mesh fields int polympo_getMeshFVtxType_f(); diff --git a/src/pmpo_fortran.f90 b/src/pmpo_fortran.f90 index cf013ca..f386058 100644 --- a/src/pmpo_fortran.f90 +++ b/src/pmpo_fortran.f90 @@ -698,6 +698,21 @@ subroutine polympo_setOwningProc(mpMesh, nCells, array) & integer(c_int), value :: nCells type(c_ptr), intent(in), 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_setElmGlobal(mpMesh, nCells, array) & + bind(C, NAME='polympo_setElmGlobal_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 diff --git a/src/pmpo_materialPoints.cpp b/src/pmpo_materialPoints.cpp index c20fe00..5e4cdae 100644 --- a/src/pmpo_materialPoints.cpp +++ b/src/pmpo_materialPoints.cpp @@ -40,8 +40,12 @@ PS* createDPS(int numElms, int numMPs, MPSView positions, IntVi return dps; } -PS* createDPS(int numElms, int numMPs, IntView mpsPerElm, IntView mp2elm, IntView mpAppID) { +PS* createDPS(int numElms, int numMPs, IntView mpsPerElm, IntView mp2elm, IntView mpAppID, IntView elm2global) { PS::kkGidView elmGids("elementGlobalIds", numElms); //TODO - initialize this to [0..numElms) + Kokkos::parallel_for("setGids", numElms, KOKKOS_LAMBDA(const int elm){ + elmGids(elm) = elm2global(elm); + }); + auto mpInfo = createInternalMemberViews(numMPs, mp2elm, mpAppID); Kokkos::TeamPolicy policy(numElms,Kokkos::AUTO); auto dps = new DPS(policy, numElms, numMPs, mpsPerElm, elmGids, mp2elm, mpInfo); @@ -57,8 +61,8 @@ MaterialPoints::MaterialPoints(int numElms, int numMPs, MPSView operating_mode = MP_RELEASE; }; -MaterialPoints::MaterialPoints(int numElms, int numMPs, IntView mpsPerElm, IntView mp2elm, IntView mpAppID) { - MPs = createDPS(numElms, numMPs, mpsPerElm, mp2elm, mpAppID); +MaterialPoints::MaterialPoints(int numElms, int numMPs, IntView mpsPerElm, IntView mp2elm, IntView mpAppID, IntView elm2global){ + MPs = createDPS(numElms, numMPs, mpsPerElm, mp2elm, mpAppID, elm2global); updateMaxAppID(); operating_mode = MP_RELEASE; }; diff --git a/src/pmpo_materialPoints.hpp b/src/pmpo_materialPoints.hpp index 99d2964..1701985 100644 --- a/src/pmpo_materialPoints.hpp +++ b/src/pmpo_materialPoints.hpp @@ -130,7 +130,7 @@ class MaterialPoints { public: MaterialPoints() : MPs(nullptr) {}; MaterialPoints(int numElms, int numMPs, MPSView positions, IntView mpsPerElm, IntView mp2elm); - MaterialPoints(int numElms, int numMPs, IntView mpsPerElm, IntView mp2elm, IntView mpAppID); + MaterialPoints(int numElms, int numMPs, IntView mpsPerElm, IntView mp2elm, IntView mpAppID, IntView elm2global); ~MaterialPoints(); void rebuild(IntView addedMP2elm, IntView addedMPAppID); diff --git a/src/pmpo_mesh.cpp b/src/pmpo_mesh.cpp index 9232ad1..57b838f 100644 --- a/src/pmpo_mesh.cpp +++ b/src/pmpo_mesh.cpp @@ -77,5 +77,11 @@ namespace polyMPO{ IntView Mesh::getElm2Process() { return owningProc_; } + + IntView Mesh::getElmGlobal() { + return globalElm_; + } + + } // namespace polyMPO diff --git a/src/pmpo_mesh.hpp b/src/pmpo_mesh.hpp index 2058e13..0819da2 100644 --- a/src/pmpo_mesh.hpp +++ b/src/pmpo_mesh.hpp @@ -88,7 +88,7 @@ class Mesh { IntVtx2ElmView elm2VtxConn_; IntElm2ElmView elm2ElmConn_; IntView owningProc_; - + IntView globalElm_; //start of meshFields MeshFView vtxCoords_; MeshFView vtxRotLat_; @@ -163,6 +163,10 @@ class Mesh { void setOwningProc(IntView owningProc) {PMT_ALWAYS_ASSERT(meshEdit_); owningProc_ = owningProc; } + void setElmGlobal(IntView globalElm) {PMT_ALWAYS_ASSERT(meshEdit_); + globalElm_ = globalElm; } + IntView getElmGlobal(); + void computeRotLatLonIncr(); }; From 8bc07803aebf619050c46c59086d727cb7b15e9c Mon Sep 17 00:00:00 2001 From: nathd2 Date: Tue, 1 Apr 2025 15:20:35 -0400 Subject: [PATCH 04/17] Running on multiple processors --- src/pmpo_MPMesh.cpp | 22 +++++++++------------- src/pmpo_c.cpp | 5 ++--- src/pmpo_materialPoints.cpp | 26 ++++++++++++++++++++++---- src/pmpo_materialPoints.hpp | 3 ++- 4 files changed, 35 insertions(+), 21 deletions(-) diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index b282231..f2f0929 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -129,20 +129,13 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){ auto elm2Process = p_mesh->getElm2Process(); auto elm2global = p_mesh->getElmGlobal(); - MPI_Comm comm = p_MPs->getMPIComm(); int comm_rank; MPI_Comm_rank(comm, &comm_rank); - assert(cudaDeviceSynchronize() == cudaSuccess); - MPI_Barrier(MPI_COMM_WORLD); - printf("FooTracking \n"); Kokkos::parallel_for("countProcess", numElms, KOKKOS_LAMBDA(const int iElm){ int pp_id=elm2Process(iElm); - if(iElm<3 || iElm==1470 || iElm==1471 || iElm==1472 ) - printf("Mesh elm %d Process %d owning proc %d global %d \n", iElm, comm_rank, pp_id, elm2global(iElm)); }); - if(printVTPIndex>=0) { printVTP_mesh(printVTPIndex); } @@ -341,8 +334,9 @@ bool getAnyIsMigrating(MaterialPoints* p_MPs, bool isMigrating) { void MPMesh::push(){ static int count=0; - std::cout<<__FUNCTION__<<" "<computeRotLatLonIncr(); @@ -357,14 +351,16 @@ void MPMesh::push(){ do { CVTTrackingElmCenterBased(); // move to Tgt_XYZ assert(cudaDeviceSynchronize() == cudaSuccess); - printf("FooPushDoWhile\n"); p_MPs->updateMPSlice(); // Tgt_XYZ becomes Cur_XYZ p_MPs->updateMPSlice(); // Tgt becomes Cur - if (elm2Process.size() > 0) - anyIsMigrating = getAnyIsMigrating(p_MPs, p_MPs->migrate()); + + bool anyIsMigrating = getAnyIsMigrating(p_MPs, p_MPs->check_migrate()); + + if(anyIsMigrating) + p_MPs->migrate(); else - p_MPs->rebuild(); //rebuild pumi-pic - printf("Is migrating %d \n", anyIsMigrating); + p_MPs->rebuild(); + p_MPs->updateMPElmID(); //update mpElm IDs slices reconstructSlices(); } diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index b60a63a..19f1307 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -66,11 +66,10 @@ void polympo_createMPs_f(MPMesh_ptr p_mpmesh, const int* mp2Elm, const int* isMPActive) { checkMPMeshValid(p_mpmesh); - std::cout<<__FUNCTION__<p_mesh; PMT_ALWAYS_ASSERT(!p_mesh->meshEditable()); - //PMT_ALWAYS_ASSERT(p_mesh->getNumElements() == numElms); + PMT_ALWAYS_ASSERT(p_mesh->getNumElements() == numElms); //Find the total no of MPs across all ranks //And loop over all MPs and find the smallest element id associated across a MP @@ -100,7 +99,7 @@ void polympo_createMPs_f(MPMesh_ptr p_mpmesh, int globalFirstElmWithMPs; MPI_Allreduce(&firstElmWithMPs, &globalFirstElmWithMPs, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD); - printf("With a MP, globally smallest mesh elm %d and first elm %d \n", globalMinElmID, globalFirstElmWithMPs); + //printf("With a MP, globally smallest mesh elm %d and first elm %d \n", globalMinElmID, globalFirstElmWithMPs); int offset = -1; if(globalMinElmID-globalFirstElmWithMPs==1) { diff --git a/src/pmpo_materialPoints.cpp b/src/pmpo_materialPoints.cpp index 5e4cdae..1d99c95 100644 --- a/src/pmpo_materialPoints.cpp +++ b/src/pmpo_materialPoints.cpp @@ -113,14 +113,34 @@ void MaterialPoints::setMPIComm(MPI_Comm comm) { mpi_comm = comm; } -bool MaterialPoints::migrate() { +bool MaterialPoints::check_migrate(){ + Kokkos::Timer timer; + auto MPs2Elm = getData(); + auto MPs2Proc = getData(); + + IntView isMigrating("isMigrating", 1); + + int rank; + MPI_Comm_rank(mpi_comm, &rank); + auto setMigrationFields = PS_LAMBDA(const int& e, const int& mp, const bool& mask) { + if (mask) { + if (MPs2Proc(mp) != rank) isMigrating(0) = 1; + } + }; + parallel_for(setMigrationFields, "setMigrationFields"); + if (getOpMode() == polyMPO::MP_DEBUG) + printf("Material point check migration: %f\n", timer.seconds()); + pumipic::RecordTime("PolyMPO_check_migrate", timer.seconds()); + return pumipic::getLastValue(isMigrating) > 0; +} + +void MaterialPoints::migrate() { Kokkos::Timer timer; auto MPs2Elm = getData(); auto MPs2Proc = getData(); 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, &rank); @@ -128,7 +148,6 @@ bool MaterialPoints::migrate() { if (mask) { new_elem(mp) = MPs2Elm(mp); new_process(mp) = MPs2Proc(mp); - if (new_process(mp) != rank) isMigrating(0) = 1; } }; parallel_for(setMigrationFields, "setMigrationFields"); @@ -137,7 +156,6 @@ bool MaterialPoints::migrate() { 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; } diff --git a/src/pmpo_materialPoints.hpp b/src/pmpo_materialPoints.hpp index 1701985..7511a1d 100644 --- a/src/pmpo_materialPoints.hpp +++ b/src/pmpo_materialPoints.hpp @@ -138,7 +138,8 @@ class MaterialPoints { void finishRebuild(); bool rebuildOngoing(); - bool migrate(); + bool check_migrate(); + void migrate(); MPI_Comm getMPIComm(); void setMPIComm(MPI_Comm comm); From e004bc53c13e9cc617cfb504ca069ecd906b03f7 Mon Sep 17 00:00:00 2001 From: nathd2 Date: Mon, 7 Apr 2025 11:39:50 -0400 Subject: [PATCH 05/17] Vtp files issue for multi-process --- src/pmpo_MPMesh.cpp | 50 ++++++++++++++++++++++++++++++--------------- 1 file changed, 34 insertions(+), 16 deletions(-) diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index f2f0929..db27d06 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -135,16 +135,21 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){ Kokkos::parallel_for("countProcess", numElms, KOKKOS_LAMBDA(const int iElm){ int pp_id=elm2Process(iElm); }); - - if(printVTPIndex>=0) { - printVTP_mesh(printVTPIndex); + + //Since Mesh is static print pnly for 1 time step + if(printVTPIndex==0) { + printVTP_mesh(comm_rank); } + + assert(cudaDeviceSynchronize()==cudaSuccess); + MPI_Barrier(MPI_COMM_WORLD); Vec3dView history("positionHistory",numMPs); Vec3dView resultLeft("positionResult",numMPs); Vec3dView resultRight("positionResult",numMPs); Vec3dView mpTgtPosArray("positionTarget",numMPs); - + Kokkos::View 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)); if(mask){ @@ -163,7 +168,7 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){ for(int i=1; i<=numConnElms; i++){ int elmID = elm2ElmConn(iElm,i)-1; - //New delta + //New delta Vec3d center(elmCenter(elmID, 0), elmCenter(elmID, 1), elmCenter(elmID, 2)); delta = MPnew - center; @@ -183,7 +188,9 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){ iElm = closestElm; } } - if(printVTPIndex>=0){ + + if(printVTPIndex>=0 && numMPs>0){ + //printf("Rank %d mp %d counter %d \n", comm_rank, mp, counter); double d1 = dx[0]; double d2 = dx[2]; double d3 = dx[3]; @@ -195,29 +202,36 @@ 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; - history(mp) = MP; - resultLeft(mp) = MPLeft; - resultRight(mp) = MPRight; - mpTgtPosArray(mp) = MPnew; + auto xx=Kokkos::atomic_fetch_add(&counter(0), 1); + history(xx) = MP; + resultLeft(xx) = MPLeft; + resultRight(xx) = MPRight; + mpTgtPosArray(xx) = MPnew; } } }; p_MPs->parallel_for(CVTElmCalc,"CVTTrackingElmCenterBasedCalc"); - if(printVTPIndex>=0){ + assert(cudaDeviceSynchronize()==cudaSuccess); + MPI_Barrier(MPI_COMM_WORLD); + printf("After Tracking \n"); + + if(printVTPIndex>=0 && numMPs>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::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); + printf("Host counter value: %d\n", h_counter(0)); // printVTP file char* fileOutput = (char *)malloc(sizeof(char) * 256); - sprintf(fileOutput, "polyMPOCVTTrackingElmCenter_MPtracks_%d.vtp", printVTPIndex); + sprintf(fileOutput, "polyMPOCVTTrackingElmCenter_MPtracks_%d_%d.vtp", comm_rank, printVTPIndex); FILE * pFile = fopen(fileOutput,"w"); free(fileOutput); fprintf(pFile, "\n \n \n \n \n",numMPs*4,numMPs*2); @@ -239,6 +253,11 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){ fprintf(pFile," \n \n \n \n\n"); fclose(pFile); } + assert(cudaDeviceSynchronize()==cudaSuccess); + MPI_Barrier(MPI_COMM_WORLD); + printf("After printing particle paths \n"); + + pumipic::RecordTime("PolyMPO_CVTTrackingElmCenterBased", timer.seconds()); } @@ -335,7 +354,6 @@ void MPMesh::push(){ static int count=0; std::cout<<"Push"<<" "<updateMPSlice(); // Tgt_XYZ becomes Cur_XYZ p_MPs->updateMPSlice(); // Tgt becomes Cur @@ -365,7 +383,7 @@ void MPMesh::push(){ reconstructSlices(); } while (anyIsMigrating); - + count ++; pumipic::RecordTime("PolyMPO_push", timer.seconds()); } From 4274cd68a8245c71607f5b03a420c27bbc780d05 Mon Sep 17 00:00:00 2001 From: nathd2 Date: Wed, 9 Apr 2025 22:19:34 -0400 Subject: [PATCH 06/17] Towards multi Process --- src/pmpo_MPMesh.cpp | 9 --------- src/pmpo_c.cpp | 9 +++++++++ src/pmpo_c.h | 4 ++++ src/pmpo_defines.h | 1 + src/pmpo_fortran.f90 | 22 ++++++++++++++++++++++ src/pmpo_materialPoints.cpp | 30 ++++++++++++++++++++++++++++++ 6 files changed, 66 insertions(+), 9 deletions(-) diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index db27d06..0509043 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -212,10 +212,6 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){ }; p_MPs->parallel_for(CVTElmCalc,"CVTTrackingElmCenterBasedCalc"); - assert(cudaDeviceSynchronize()==cudaSuccess); - MPI_Barrier(MPI_COMM_WORLD); - printf("After Tracking \n"); - if(printVTPIndex>=0 && numMPs>0){ Vec3dView::HostMirror h_history = Kokkos::create_mirror_view(history); Vec3dView::HostMirror h_resultLeft = Kokkos::create_mirror_view(resultLeft); @@ -228,7 +224,6 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){ Kokkos::deep_copy(h_resultRight, resultRight); Kokkos::deep_copy(h_mpTgtPos, mpTgtPosArray); Kokkos::deep_copy(h_counter, counter); - printf("Host counter value: %d\n", h_counter(0)); // printVTP file char* fileOutput = (char *)malloc(sizeof(char) * 256); sprintf(fileOutput, "polyMPOCVTTrackingElmCenter_MPtracks_%d_%d.vtp", comm_rank, printVTPIndex); @@ -253,10 +248,6 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){ fprintf(pFile," \n \n \n \n\n"); fclose(pFile); } - assert(cudaDeviceSynchronize()==cudaSuccess); - MPI_Barrier(MPI_COMM_WORLD); - printf("After printing particle paths \n"); - pumipic::RecordTime("PolyMPO_CVTTrackingElmCenterBased", timer.seconds()); } diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index 19f1307..1f9ddea 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -212,6 +212,15 @@ void polympo_setAppIDFunc_f(MPMesh_ptr p_mpmesh, IntVoidFunc getNext, void* appI p_MPs->setAppIDFunc(getNextAppID); } +void polympo_setMPASAppIDFunc_f(MPMesh_ptr p_mpmesh, VoidVoidFunc getMPASAppID, void* arg1, void*arg2, void*arg3, + void* arg4, void* arg5, void*arg6) { + checkMPMeshValid(p_mpmesh); + auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; + std::function polyMPO_getMPASAppID = [getMPASAppID, arg1, arg2, arg3, arg4, arg5, arg6]() + {getMPASAppID(arg1, arg2, arg3, arg4, arg5, arg6); return *(static_cast(arg5));}; + p_MPs->setAppIDFunc(polyMPO_getMPASAppID); +} + void polympo_getMPCurElmID_f(MPMesh_ptr p_mpmesh, const int numMPs, int* elmIDs){ diff --git a/src/pmpo_c.h b/src/pmpo_c.h index f62ffdc..d749745 100644 --- a/src/pmpo_c.h +++ b/src/pmpo_c.h @@ -22,7 +22,11 @@ void polympo_setMPICommunicator_f(MPMesh_ptr p_mpmesh, MPI_Fint fcomm); void polympo_createMPs_f(MPMesh_ptr p_mpmesh, const int numElms, const int numMPs, int* mpsPerElm, const int* mp2Elm, const int* isMPActive); void polympo_startRebuildMPs_f(MPMesh_ptr p_mpmesh, const int numMPs, const int* allTgtMpElmIn, const int* addedMPMask); void polympo_finishRebuildMPs_f(MPMesh_ptr p_mpmesh); + void polympo_setAppIDFunc_f(MPMesh_ptr p_mpmesh, IntVoidFunc getNext, void* appIDs); +void polympo_setMPASAppIDFunc_f(MPMesh_ptr p_mpmesh, VoidVoidFunc getMPASAppID, void* arg1, void* arg2, void*arg3, + void* arg4, void*arg5, void* arg6); + void polympo_getMPCurElmID_f(MPMesh_ptr p_mpmesh, const int numMPs, int* elmIDs); void polympo_setMPLatLonRotatedFlag_f(MPMesh_ptr p_mpmesh, const int isRotateFlag); diff --git a/src/pmpo_defines.h b/src/pmpo_defines.h index 34857b9..3bb8898 100644 --- a/src/pmpo_defines.h +++ b/src/pmpo_defines.h @@ -8,6 +8,7 @@ typedef void* MPMesh_ptr; //Function that receives void* and returns an int typedef int (*IntVoidFunc)(void*); +typedef void (*VoidVoidFunc)(void*, void*, void*, void*, void*, void*); using space_t = Kokkos::DefaultExecutionSpace::memory_space; diff --git a/src/pmpo_fortran.f90 b/src/pmpo_fortran.f90 index f386058..1c9f473 100644 --- a/src/pmpo_fortran.f90 +++ b/src/pmpo_fortran.f90 @@ -110,6 +110,28 @@ subroutine polympo_setAppIDFunc(mpMesh, getNext, appIDs) & type(c_funptr), value :: getNext type(c_ptr), value :: appIDs end subroutine + + !--------------------------------------------------------------------------- + !> @brief Stores pointer to appID data structure and a function to retrieve them used in migration + !> @param mpmesh(in/out) MPMesh object + !> @param getNext(in) Pointer to function that returns next App IDs + !> @param appIDs(in) Pointer to opaque data application data structure (that may contain all available app IDs) + !--------------------------------------------------------------------------- + subroutine polympo_setMPASAppIDFunc(mpMesh, getMPASAppID, & + arg1, arg2, arg3, arg4, arg5, arg6) & + bind(C, NAME='polympo_setMPASAppIDFunc_f') + use :: iso_c_binding + type(c_ptr), value :: mpMesh + type(c_funptr), value :: getMPASAppID + type(c_ptr), value :: arg1 + type(c_ptr), value :: arg2 + type(c_ptr), value :: arg3 + type(c_ptr), value :: arg4 + type(c_ptr), value :: arg5 + type(c_ptr), value :: arg6 + end subroutine + + !--------------------------------------------------------------------------- !> @brief get the current element ID MP array from a polympo array !> @param mpmesh(in/out) MPMesh object diff --git a/src/pmpo_materialPoints.cpp b/src/pmpo_materialPoints.cpp index 1d99c95..8b7f7a6 100644 --- a/src/pmpo_materialPoints.cpp +++ b/src/pmpo_materialPoints.cpp @@ -138,6 +138,7 @@ void MaterialPoints::migrate() { Kokkos::Timer timer; auto MPs2Elm = getData(); auto MPs2Proc = getData(); + auto mpAppID = getData(); IntView new_elem("new_elem", MPs->capacity()); IntView new_process("new_process", MPs->capacity()); @@ -148,11 +149,40 @@ void MaterialPoints::migrate() { if (mask) { new_elem(mp) = MPs2Elm(mp); new_process(mp) = MPs2Proc(mp); + if(rank!=new_process(mp)){ + mpAppID(mp)=-1; + printf("Particle migrated and so its AppID is -1\n"); + } } }; parallel_for(setMigrationFields, "setMigrationFields"); MPs->migrate(new_elem, new_process); + //Since rebuilt + mpAppID = getData(); + //Count MPs that have -1 appID, so that we can count no of MPs received + Kokkos::View numReceivedMPs("numReceivedMPs", 1); + Kokkos::deep_copy(numReceivedMPs, 0); + auto countnewMPs = PS_LAMBDA(const int& e, const int& mp, const bool& mask) { + if(mask){ + if (mpAppID(mp) == -1) + Kokkos::atomic_add(&numReceivedMPs(0), 1); + } + }; + parallel_for(countnewMPs, "countReceivedPtcls"); + Kokkos::fence(); + + auto numReceivedMPs_host = Kokkos::create_mirror_view(numReceivedMPs); + Kokkos::deep_copy(numReceivedMPs_host, numReceivedMPs); + if(numReceivedMPs_host(0)) + std::cout <<"Rank "< added_mpIDs(numReceivedMPs_host(0)); + for(int i=0; i Date: Thu, 10 Apr 2025 12:17:05 -0400 Subject: [PATCH 07/17] One more way of trying MPAS AppIDs --- src/pmpo_c.cpp | 6 +++--- src/pmpo_c.h | 2 +- src/pmpo_defines.h | 2 +- src/pmpo_fortran.f90 | 6 ++---- src/pmpo_materialPoints.cpp | 2 +- 5 files changed, 8 insertions(+), 10 deletions(-) diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index 1f9ddea..c705493 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -213,11 +213,11 @@ void polympo_setAppIDFunc_f(MPMesh_ptr p_mpmesh, IntVoidFunc getNext, void* appI } void polympo_setMPASAppIDFunc_f(MPMesh_ptr p_mpmesh, VoidVoidFunc getMPASAppID, void* arg1, void*arg2, void*arg3, - void* arg4, void* arg5, void*arg6) { + const int arg4) { checkMPMeshValid(p_mpmesh); auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; - std::function polyMPO_getMPASAppID = [getMPASAppID, arg1, arg2, arg3, arg4, arg5, arg6]() - {getMPASAppID(arg1, arg2, arg3, arg4, arg5, arg6); return *(static_cast(arg5));}; + std::function polyMPO_getMPASAppID = [getMPASAppID, arg1, arg2, arg3, arg4]() + { int iParticleNew; getMPASAppID(arg1, arg2, arg3, arg4, iParticleNew); return iParticleNew;}; p_MPs->setAppIDFunc(polyMPO_getMPASAppID); } diff --git a/src/pmpo_c.h b/src/pmpo_c.h index d749745..be83a83 100644 --- a/src/pmpo_c.h +++ b/src/pmpo_c.h @@ -25,7 +25,7 @@ void polympo_finishRebuildMPs_f(MPMesh_ptr p_mpmesh); void polympo_setAppIDFunc_f(MPMesh_ptr p_mpmesh, IntVoidFunc getNext, void* appIDs); void polympo_setMPASAppIDFunc_f(MPMesh_ptr p_mpmesh, VoidVoidFunc getMPASAppID, void* arg1, void* arg2, void*arg3, - void* arg4, void*arg5, void* arg6); + const int arg4); void polympo_getMPCurElmID_f(MPMesh_ptr p_mpmesh, const int numMPs, int* elmIDs); void polympo_setMPLatLonRotatedFlag_f(MPMesh_ptr p_mpmesh, const int isRotateFlag); diff --git a/src/pmpo_defines.h b/src/pmpo_defines.h index 3bb8898..f4abcb0 100644 --- a/src/pmpo_defines.h +++ b/src/pmpo_defines.h @@ -8,7 +8,7 @@ typedef void* MPMesh_ptr; //Function that receives void* and returns an int typedef int (*IntVoidFunc)(void*); -typedef void (*VoidVoidFunc)(void*, void*, void*, void*, void*, void*); +typedef void (*VoidVoidFunc)(void*, void*, void*, const int, int); using space_t = Kokkos::DefaultExecutionSpace::memory_space; diff --git a/src/pmpo_fortran.f90 b/src/pmpo_fortran.f90 index 1c9f473..646a07d 100644 --- a/src/pmpo_fortran.f90 +++ b/src/pmpo_fortran.f90 @@ -118,7 +118,7 @@ subroutine polympo_setAppIDFunc(mpMesh, getNext, appIDs) & !> @param appIDs(in) Pointer to opaque data application data structure (that may contain all available app IDs) !--------------------------------------------------------------------------- subroutine polympo_setMPASAppIDFunc(mpMesh, getMPASAppID, & - arg1, arg2, arg3, arg4, arg5, arg6) & + arg1, arg2, arg3, arg4) & bind(C, NAME='polympo_setMPASAppIDFunc_f') use :: iso_c_binding type(c_ptr), value :: mpMesh @@ -126,9 +126,7 @@ subroutine polympo_setMPASAppIDFunc(mpMesh, getMPASAppID, & type(c_ptr), value :: arg1 type(c_ptr), value :: arg2 type(c_ptr), value :: arg3 - type(c_ptr), value :: arg4 - type(c_ptr), value :: arg5 - type(c_ptr), value :: arg6 + integer(c_int), intent(in), value :: arg4 end subroutine diff --git a/src/pmpo_materialPoints.cpp b/src/pmpo_materialPoints.cpp index 8b7f7a6..394163e 100644 --- a/src/pmpo_materialPoints.cpp +++ b/src/pmpo_materialPoints.cpp @@ -151,7 +151,7 @@ void MaterialPoints::migrate() { new_process(mp) = MPs2Proc(mp); if(rank!=new_process(mp)){ mpAppID(mp)=-1; - printf("Particle migrated and so its AppID is -1\n"); + printf("Particle migrated and so its AppID is -1\n"); } } }; From c52bfa73167516379b7a10ac949916c1f55de011 Mon Sep 17 00:00:00 2001 From: nathd2 Date: Fri, 11 Apr 2025 10:54:44 -0400 Subject: [PATCH 08/17] MPASAppID now takes cell no as input --- src/pmpo_c.cpp | 6 +++--- src/pmpo_defines.h | 2 +- src/pmpo_materialPoints.cpp | 27 +++++++++++++++++++++------ src/pmpo_materialPoints.hpp | 12 +++++++++--- 4 files changed, 34 insertions(+), 13 deletions(-) diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index c705493..6616f9e 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -209,15 +209,15 @@ void polympo_setAppIDFunc_f(MPMesh_ptr p_mpmesh, IntVoidFunc getNext, void* appI checkMPMeshValid(p_mpmesh); auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; polyMPO::IntFunc getNextAppID = [getNext, appIDs]() { return getNext(appIDs); }; - p_MPs->setAppIDFunc(getNextAppID); + //p_MPs->setAppIDFunc(getNextAppID); } void polympo_setMPASAppIDFunc_f(MPMesh_ptr p_mpmesh, VoidVoidFunc getMPASAppID, void* arg1, void*arg2, void*arg3, const int arg4) { checkMPMeshValid(p_mpmesh); auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; - std::function polyMPO_getMPASAppID = [getMPASAppID, arg1, arg2, arg3, arg4]() - { int iParticleNew; getMPASAppID(arg1, arg2, arg3, arg4, iParticleNew); return iParticleNew;}; + std::function polyMPO_getMPASAppID = [getMPASAppID, arg1, arg2, arg3, arg4](int iCell) + { int iParticleNew; getMPASAppID(arg1, arg2, arg3, arg4, iCell, iParticleNew); return iParticleNew;}; p_MPs->setAppIDFunc(polyMPO_getMPASAppID); } diff --git a/src/pmpo_defines.h b/src/pmpo_defines.h index f4abcb0..21a7d49 100644 --- a/src/pmpo_defines.h +++ b/src/pmpo_defines.h @@ -8,7 +8,7 @@ typedef void* MPMesh_ptr; //Function that receives void* and returns an int typedef int (*IntVoidFunc)(void*); -typedef void (*VoidVoidFunc)(void*, void*, void*, const int, int); +typedef void (*VoidVoidFunc)(void*, void*, void*, const int, int, int); using space_t = Kokkos::DefaultExecutionSpace::memory_space; diff --git a/src/pmpo_materialPoints.cpp b/src/pmpo_materialPoints.cpp index 394163e..9bea799 100644 --- a/src/pmpo_materialPoints.cpp +++ b/src/pmpo_materialPoints.cpp @@ -158,7 +158,7 @@ void MaterialPoints::migrate() { parallel_for(setMigrationFields, "setMigrationFields"); MPs->migrate(new_elem, new_process); - //Since rebuilt + //SINCE REBUILT, mpAppID needs to be be recalled mpAppID = getData(); //Count MPs that have -1 appID, so that we can count no of MPs received Kokkos::View numReceivedMPs("numReceivedMPs", 1); @@ -171,16 +171,31 @@ void MaterialPoints::migrate() { }; parallel_for(countnewMPs, "countReceivedPtcls"); Kokkos::fence(); - + //Bring them to CPU and print particles received auto numReceivedMPs_host = Kokkos::create_mirror_view(numReceivedMPs); Kokkos::deep_copy(numReceivedMPs_host, numReceivedMPs); if(numReceivedMPs_host(0)) std::cout <<"Rank "< added_mpIDs(numReceivedMPs_host(0)); + //Another array that contains element id of the the MPs that have migrated + Kokkos::View receivedMPs2Elm("ReceivedMPs2Elm", numReceivedMPs_host(0)); + Kokkos::View counter("counter", 1); + auto set_new_elem = PS_LAMBDA(const int& e, const int& mp, const bool& mask) { + if(mask){ + if (mpAppID(mp) == -1){ + auto xx=Kokkos::atomic_fetch_add(&counter(0), 1); + receivedMPs2Elm(xx) = e; + } + } + }; + parallel_for(set_new_elem, "countReceivedPtcls"); + //NEED AN ASSERT ELEMNT DOING COUNTER=#RECEIVED MPs + auto receivedMPs2Elm_host = Kokkos::create_mirror_view(receivedMPs2Elm); + Kokkos::deep_copy(receivedMPs2Elm_host, receivedMPs2Elm); + for(int i=0; i IntFunc; +typedef std::function IntIntFunc; + + enum MaterialPointSlice { MPF_Status = 0, @@ -124,7 +127,8 @@ class MaterialPoints { bool isRotatedFlag = false; Operating_Mode operating_mode; RebuildHelper rebuildFields; - IntFunc getAppID; + //IntFunc getAppID; + IntIntFunc getAppID; MPI_Comm mpi_comm; public: @@ -151,8 +155,10 @@ class MaterialPoints { typename std::enable_if::type setRebuildMPSlice(mpSliceData mpSliceIn); - void setAppIDFunc(IntFunc getAppIDIn); - int getNextAppID(); + //void setAppIDFunc(IntFunc getAppIDIn); + void setAppIDFunc(IntIntFunc getAppIDIn); + + int getNextAppID(int iElm); void rebuild() { IntView tgtElm("tgtElm", MPs->capacity()); From 8e310cfd67a772ce07fcad5e3c1f3523e21de62f Mon Sep 17 00:00:00 2001 From: nathd2 Date: Sat, 12 Apr 2025 07:59:05 -0400 Subject: [PATCH 09/17] iParticleNew Issue --- src/pmpo_MPMesh.cpp | 2 +- src/pmpo_c.cpp | 6 +++--- src/pmpo_c.h | 2 +- src/pmpo_defines.h | 2 +- src/pmpo_fortran.f90 | 3 ++- src/pmpo_materialPoints.cpp | 2 +- 6 files changed, 9 insertions(+), 8 deletions(-) diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index 0509043..cd3fd33 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -345,7 +345,7 @@ void MPMesh::push(){ static int count=0; std::cout<<"Push"<<" "<473) exit(-1); Kokkos::Timer timer; p_mesh->computeRotLatLonIncr(); diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index 6616f9e..81df542 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -213,11 +213,11 @@ void polympo_setAppIDFunc_f(MPMesh_ptr p_mpmesh, IntVoidFunc getNext, void* appI } void polympo_setMPASAppIDFunc_f(MPMesh_ptr p_mpmesh, VoidVoidFunc getMPASAppID, void* arg1, void*arg2, void*arg3, - const int arg4) { + const int arg4, const int arg5) { checkMPMeshValid(p_mpmesh); auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; - std::function polyMPO_getMPASAppID = [getMPASAppID, arg1, arg2, arg3, arg4](int iCell) - { int iParticleNew; getMPASAppID(arg1, arg2, arg3, arg4, iCell, iParticleNew); return iParticleNew;}; + std::function polyMPO_getMPASAppID = [getMPASAppID, arg1, arg2, arg3, arg4, arg5](int iCell) + { int iParticleNew; getMPASAppID(arg1, arg2, arg3, arg4, iCell, iParticleNew, arg5); return iParticleNew;}; p_MPs->setAppIDFunc(polyMPO_getMPASAppID); } diff --git a/src/pmpo_c.h b/src/pmpo_c.h index be83a83..13f36a8 100644 --- a/src/pmpo_c.h +++ b/src/pmpo_c.h @@ -25,7 +25,7 @@ void polympo_finishRebuildMPs_f(MPMesh_ptr p_mpmesh); void polympo_setAppIDFunc_f(MPMesh_ptr p_mpmesh, IntVoidFunc getNext, void* appIDs); void polympo_setMPASAppIDFunc_f(MPMesh_ptr p_mpmesh, VoidVoidFunc getMPASAppID, void* arg1, void* arg2, void*arg3, - const int arg4); + const int arg4, const int arg5); void polympo_getMPCurElmID_f(MPMesh_ptr p_mpmesh, const int numMPs, int* elmIDs); void polympo_setMPLatLonRotatedFlag_f(MPMesh_ptr p_mpmesh, const int isRotateFlag); diff --git a/src/pmpo_defines.h b/src/pmpo_defines.h index 21a7d49..05b301b 100644 --- a/src/pmpo_defines.h +++ b/src/pmpo_defines.h @@ -8,7 +8,7 @@ typedef void* MPMesh_ptr; //Function that receives void* and returns an int typedef int (*IntVoidFunc)(void*); -typedef void (*VoidVoidFunc)(void*, void*, void*, const int, int, int); +typedef void (*VoidVoidFunc)(void*, void*, void*, const int, int, int, const int); using space_t = Kokkos::DefaultExecutionSpace::memory_space; diff --git a/src/pmpo_fortran.f90 b/src/pmpo_fortran.f90 index 646a07d..e5670f3 100644 --- a/src/pmpo_fortran.f90 +++ b/src/pmpo_fortran.f90 @@ -118,7 +118,7 @@ subroutine polympo_setAppIDFunc(mpMesh, getNext, appIDs) & !> @param appIDs(in) Pointer to opaque data application data structure (that may contain all available app IDs) !--------------------------------------------------------------------------- subroutine polympo_setMPASAppIDFunc(mpMesh, getMPASAppID, & - arg1, arg2, arg3, arg4) & + arg1, arg2, arg3, arg4, arg5) & bind(C, NAME='polympo_setMPASAppIDFunc_f') use :: iso_c_binding type(c_ptr), value :: mpMesh @@ -127,6 +127,7 @@ subroutine polympo_setMPASAppIDFunc(mpMesh, getMPASAppID, & type(c_ptr), value :: arg2 type(c_ptr), value :: arg3 integer(c_int), intent(in), value :: arg4 + integer(c_int), intent(in), value :: arg5 end subroutine diff --git a/src/pmpo_materialPoints.cpp b/src/pmpo_materialPoints.cpp index 9bea799..2a90ba1 100644 --- a/src/pmpo_materialPoints.cpp +++ b/src/pmpo_materialPoints.cpp @@ -195,7 +195,7 @@ void MaterialPoints::migrate() { for(int i=0; i Date: Wed, 16 Apr 2025 12:56:15 -0400 Subject: [PATCH 10/17] New particleID can be found from MPAS in polyMPO --- src/pmpo_MPMesh.cpp | 22 +++++++++++++--- src/pmpo_c.cpp | 9 ++++--- src/pmpo_c.h | 3 +-- src/pmpo_defines.h | 2 +- src/pmpo_fortran.f90 | 8 +++--- src/pmpo_materialPoints.cpp | 50 ++++++++++++++++++++++++++----------- 6 files changed, 64 insertions(+), 30 deletions(-) diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index cd3fd33..57afa6b 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -2,7 +2,7 @@ #include "pmpo_utils.hpp" #include "pmpo_MPMesh.hpp" #include "pmpo_wachspressBasis.hpp" - +#include namespace polyMPO{ void printVTP_mesh(MPMesh& mpMesh, int printVTPIndex=-1); @@ -345,7 +345,7 @@ void MPMesh::push(){ static int count=0; std::cout<<"Push"<<" "<473) exit(-1); + Kokkos::Timer timer; p_mesh->computeRotLatLonIncr(); @@ -358,7 +358,7 @@ void MPMesh::push(){ bool anyIsMigrating = false; do { - CVTTrackingElmCenterBased(count); // move to Tgt_XYZ + CVTTrackingElmCenterBased(); // move to Tgt_XYZ assert(cudaDeviceSynchronize() == cudaSuccess); p_MPs->updateMPSlice(); // Tgt_XYZ becomes Cur_XYZ p_MPs->updateMPSlice(); // Tgt becomes Cur @@ -369,12 +369,26 @@ void MPMesh::push(){ p_MPs->migrate(); else p_MPs->rebuild(); + printf("Done till here 0\n"); p_MPs->updateMPElmID(); //update mpElm IDs slices reconstructSlices(); } while (anyIsMigrating); - count ++; + + count ++; + + volatile int i = 0; + char hostname[256]; + gethostname(hostname, sizeof(hostname)); + MPI_Comm comm = p_MPs->getMPIComm(); + int comm_rank; + MPI_Comm_rank(comm, &comm_rank); + printf("Rank %d PID %d on %s ready for attach\n", comm_rank, getpid(), hostname); + fflush(stdout); + if(count==480) sleep(100); + + printf("Done till here\n"); pumipic::RecordTime("PolyMPO_push", timer.seconds()); } diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index 81df542..94e1589 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -212,12 +212,13 @@ void polympo_setAppIDFunc_f(MPMesh_ptr p_mpmesh, IntVoidFunc getNext, void* appI //p_MPs->setAppIDFunc(getNextAppID); } -void polympo_setMPASAppIDFunc_f(MPMesh_ptr p_mpmesh, VoidVoidFunc getMPASAppID, void* arg1, void*arg2, void*arg3, - const int arg4, const int arg5) { +//arg1 is blockPtr, arg2 is blockSize and arg3 is nCells +void polympo_setMPASAppIDFunc_f(MPMesh_ptr p_mpmesh, VoidVoidFunc getMPASAppID, void*arg1, + const int arg2, const int arg3) { checkMPMeshValid(p_mpmesh); auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; - std::function polyMPO_getMPASAppID = [getMPASAppID, arg1, arg2, arg3, arg4, arg5](int iCell) - { int iParticleNew; getMPASAppID(arg1, arg2, arg3, arg4, iCell, iParticleNew, arg5); return iParticleNew;}; + std::function polyMPO_getMPASAppID = [getMPASAppID, arg1, arg2, arg3](int iCell) + { int iParticleNew=-1; getMPASAppID(arg1, arg2, arg3, iCell, iParticleNew); return iParticleNew;}; p_MPs->setAppIDFunc(polyMPO_getMPASAppID); } diff --git a/src/pmpo_c.h b/src/pmpo_c.h index 13f36a8..6dcbf4e 100644 --- a/src/pmpo_c.h +++ b/src/pmpo_c.h @@ -24,8 +24,7 @@ void polympo_startRebuildMPs_f(MPMesh_ptr p_mpmesh, const int numMPs, const int* void polympo_finishRebuildMPs_f(MPMesh_ptr p_mpmesh); void polympo_setAppIDFunc_f(MPMesh_ptr p_mpmesh, IntVoidFunc getNext, void* appIDs); -void polympo_setMPASAppIDFunc_f(MPMesh_ptr p_mpmesh, VoidVoidFunc getMPASAppID, void* arg1, void* arg2, void*arg3, - const int arg4, const int arg5); +void polympo_setMPASAppIDFunc_f(MPMesh_ptr p_mpmesh, VoidVoidFunc getMPASAppID, void* arg1, const int arg2, const int arg3); void polympo_getMPCurElmID_f(MPMesh_ptr p_mpmesh, const int numMPs, int* elmIDs); void polympo_setMPLatLonRotatedFlag_f(MPMesh_ptr p_mpmesh, const int isRotateFlag); diff --git a/src/pmpo_defines.h b/src/pmpo_defines.h index 05b301b..14435eb 100644 --- a/src/pmpo_defines.h +++ b/src/pmpo_defines.h @@ -8,7 +8,7 @@ typedef void* MPMesh_ptr; //Function that receives void* and returns an int typedef int (*IntVoidFunc)(void*); -typedef void (*VoidVoidFunc)(void*, void*, void*, const int, int, int, const int); +typedef void (*VoidVoidFunc)(void*, const int, const int, const int, int&); using space_t = Kokkos::DefaultExecutionSpace::memory_space; diff --git a/src/pmpo_fortran.f90 b/src/pmpo_fortran.f90 index e5670f3..c0cd8ca 100644 --- a/src/pmpo_fortran.f90 +++ b/src/pmpo_fortran.f90 @@ -118,16 +118,14 @@ subroutine polympo_setAppIDFunc(mpMesh, getNext, appIDs) & !> @param appIDs(in) Pointer to opaque data application data structure (that may contain all available app IDs) !--------------------------------------------------------------------------- subroutine polympo_setMPASAppIDFunc(mpMesh, getMPASAppID, & - arg1, arg2, arg3, arg4, arg5) & + arg1, arg2, arg3) & bind(C, NAME='polympo_setMPASAppIDFunc_f') use :: iso_c_binding type(c_ptr), value :: mpMesh type(c_funptr), value :: getMPASAppID type(c_ptr), value :: arg1 - type(c_ptr), value :: arg2 - type(c_ptr), value :: arg3 - integer(c_int), intent(in), value :: arg4 - integer(c_int), intent(in), value :: arg5 + integer(c_int), intent(in), value :: arg2 + integer(c_int), intent(in), value :: arg3 end subroutine diff --git a/src/pmpo_materialPoints.cpp b/src/pmpo_materialPoints.cpp index 2a90ba1..0aa7f40 100644 --- a/src/pmpo_materialPoints.cpp +++ b/src/pmpo_materialPoints.cpp @@ -158,9 +158,9 @@ void MaterialPoints::migrate() { parallel_for(setMigrationFields, "setMigrationFields"); MPs->migrate(new_elem, new_process); - //SINCE REBUILT, mpAppID needs to be be recalled + //AS REBUILT, mpAppID needs to be be recalled mpAppID = getData(); - //Count MPs that have -1 appID, so that we can count no of MPs received + //Count MPs that have -1 appID, so that we can count no of MPs received by a rank Kokkos::View numReceivedMPs("numReceivedMPs", 1); Kokkos::deep_copy(numReceivedMPs, 0); auto countnewMPs = PS_LAMBDA(const int& e, const int& mp, const bool& mask) { @@ -170,34 +170,56 @@ void MaterialPoints::migrate() { } }; parallel_for(countnewMPs, "countReceivedPtcls"); - Kokkos::fence(); - //Bring them to CPU and print particles received auto numReceivedMPs_host = Kokkos::create_mirror_view(numReceivedMPs); Kokkos::deep_copy(numReceivedMPs_host, numReceivedMPs); - if(numReceivedMPs_host(0)) - std::cout <<"Rank "< receivedMPs2Elm("ReceivedMPs2Elm", numReceivedMPs_host(0)); Kokkos::View counter("counter", 1); + Kokkos::deep_copy(counter, 0); auto set_new_elem = PS_LAMBDA(const int& e, const int& mp, const bool& mask) { if(mask){ if (mpAppID(mp) == -1){ - auto xx=Kokkos::atomic_fetch_add(&counter(0), 1); - receivedMPs2Elm(xx) = e; + auto count_temp=Kokkos::atomic_fetch_add(&counter(0), 1); + receivedMPs2Elm(count_temp) = e; } } }; - parallel_for(set_new_elem, "countReceivedPtcls"); - //NEED AN ASSERT ELEMNT DOING COUNTER=#RECEIVED MPs + parallel_for(set_new_elem, "countReceivedPtcls"); + //Bring them to CPU so that elm_id can be passed and a new mpAppID can be found from CPU auto receivedMPs2Elm_host = Kokkos::create_mirror_view(receivedMPs2Elm); Kokkos::deep_copy(receivedMPs2Elm_host, receivedMPs2Elm); - + auto counter_host = Kokkos::create_mirror_view(counter); + Kokkos::deep_copy(counter_host, counter); + assert(numReceivedMPs_host(0)==counter_host(0)); + if(counter_host(0)) + std::cout <<"Rank "< appIDs; for(int i=0; i appIDs_host(appIDs.data(), appIDs.size()); + Kokkos::View appIDs_d("appIDsDevice", appIDs.size()); + Kokkos::deep_copy(appIDs_d, appIDs_host); + //If the mpAppID is -1 assign a new mpAppID + Kokkos::deep_copy(counter, 0); + auto set_appID = PS_LAMBDA(const int& e, const int& mp, const bool& mask) { + if(mask){ + if (mpAppID(mp) == -1){ + auto count_temp=Kokkos::atomic_fetch_add(&counter(0), 1); + mpAppID(mp)=appIDs_d(count_temp)-1; + printf("Count_temp %d and new appID %d \n", count_temp, mpAppID(mp)); + } + } + }; + parallel_for(set_appID, "setApplicationIDs"); + if (getOpMode() == polyMPO::MP_DEBUG) printf("Material point migration: %f\n", timer.seconds()); pumipic::RecordTime("PolyMPO_migrate", timer.seconds()); From a8582e86e14493c0fe052918b64025e8cd169d86 Mon Sep 17 00:00:00 2001 From: nathd2 Date: Thu, 24 Apr 2025 09:37:22 -0400 Subject: [PATCH 11/17] Before Phase 1 --- src/pmpo_MPMesh.cpp | 11 +++++------ src/pmpo_materialPoints.cpp | 4 +--- 2 files changed, 6 insertions(+), 9 deletions(-) diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index 57afa6b..06b13d8 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -369,7 +369,6 @@ void MPMesh::push(){ p_MPs->migrate(); else p_MPs->rebuild(); - printf("Done till here 0\n"); p_MPs->updateMPElmID(); //update mpElm IDs slices reconstructSlices(); @@ -384,11 +383,11 @@ void MPMesh::push(){ MPI_Comm comm = p_MPs->getMPIComm(); int comm_rank; MPI_Comm_rank(comm, &comm_rank); - printf("Rank %d PID %d on %s ready for attach\n", comm_rank, getpid(), hostname); - fflush(stdout); - if(count==480) sleep(100); - - printf("Done till here\n"); + if(count==480){ + printf("Rank %d PID %d on %s ready for attach\n", comm_rank, getpid(), hostname); + fflush(stdout); + sleep(100); + } pumipic::RecordTime("PolyMPO_push", timer.seconds()); } diff --git a/src/pmpo_materialPoints.cpp b/src/pmpo_materialPoints.cpp index 0aa7f40..76218aa 100644 --- a/src/pmpo_materialPoints.cpp +++ b/src/pmpo_materialPoints.cpp @@ -151,7 +151,6 @@ void MaterialPoints::migrate() { new_process(mp) = MPs2Proc(mp); if(rank!=new_process(mp)){ mpAppID(mp)=-1; - printf("Particle migrated and so its AppID is -1\n"); } } }; @@ -200,7 +199,7 @@ void MaterialPoints::migrate() { std::vector appIDs; for(int i=0; i appIDs_host(appIDs.data(), appIDs.size()); @@ -214,7 +213,6 @@ void MaterialPoints::migrate() { if (mpAppID(mp) == -1){ auto count_temp=Kokkos::atomic_fetch_add(&counter(0), 1); mpAppID(mp)=appIDs_d(count_temp)-1; - printf("Count_temp %d and new appID %d \n", count_temp, mpAppID(mp)); } } }; From 3988a6aee621216a30668deb225854443ca556b1 Mon Sep 17 00:00:00 2001 From: nathd2 Date: Wed, 30 Apr 2025 08:58:42 -0400 Subject: [PATCH 12/17] Push operation from MPAS in while loop --- src/pmpo_MPMesh.cpp | 35 ++++++++++++++++++++-- src/pmpo_MPMesh.hpp | 1 + src/pmpo_c.cpp | 59 +++++++++++++++++++++++++++++++++++-- src/pmpo_c.h | 7 ++++- src/pmpo_fortran.f90 | 35 +++++++++++++++++++++- src/pmpo_materialPoints.cpp | 19 +++++++++--- src/pmpo_materialPoints.hpp | 7 ++++- 7 files changed, 151 insertions(+), 12 deletions(-) diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index 06b13d8..f15f2ff 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -141,9 +141,9 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){ printVTP_mesh(comm_rank); } - assert(cudaDeviceSynchronize()==cudaSuccess); - MPI_Barrier(MPI_COMM_WORLD); - + //assert(cudaDeviceSynchronize()==cudaSuccess); + //MPI_Barrier(MPI_COMM_WORLD); + printf("NumMPs %d \n", numMPs); Vec3dView history("positionHistory",numMPs); Vec3dView resultLeft("positionResult",numMPs); Vec3dView resultRight("positionResult",numMPs); @@ -341,6 +341,35 @@ bool getAnyIsMigrating(MaterialPoints* p_MPs, bool isMigrating) { return anyIsMigrating; } +bool MPMesh::push1P(){ + + static int count=0; + std::cout<<"Push1P"<<" "<computeRotLatLonIncr(); + sphericalInterpolation(*this); + + //Push the MPs + p_MPs->updateRotLatLonAndXYZ2Tgt(p_mesh->getSphereRadius()); + + //Given target location find the new MP element and the process it belongs to + CVTTrackingElmCenterBased(); + //From the above two inputs find if any particle needs to be migrated + bool anyIsMigrating = getAnyIsMigrating(p_MPs, p_MPs->check_migrate()); + + //New element (maynot be final) so that MPAS can do the migration + p_MPs->updateMPElmID(); + + + p_MPs->updateMPSlice(); + p_MPs->updateMPSlice(); + + count++; + return anyIsMigrating; + +} + void MPMesh::push(){ static int count=0; diff --git a/src/pmpo_MPMesh.hpp b/src/pmpo_MPMesh.hpp index 33f0599..b30bcf3 100644 --- a/src/pmpo_MPMesh.hpp +++ b/src/pmpo_MPMesh.hpp @@ -43,6 +43,7 @@ class MPMesh{ void CVTTrackingEdgeCenterBased(Vec2dView dx); void CVTTrackingElmCenterBased(const int printVTPIndex = -1); void T2LTracking(Vec2dView dx); + bool push1P(); void push(); void calcBasis(); diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index 94e1589..85046a8 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -110,6 +110,7 @@ void polympo_createMPs_f(MPMesh_ptr p_mpmesh, fprintf(stderr,"The minElmID is incorrect! Offset is wrong!\n"); exit(1); } + printf("Offset %d \n", offset); std::vector active_mpIDs(numMPs); std::vector active_mp2Elm(numMPs); @@ -198,6 +199,43 @@ void polympo_startRebuildMPs_f(MPMesh_ptr p_mpmesh, } } +void polympo_startRebuildMPs_f2(MPMesh_ptr p_mpmesh, + const int sizeMP2elm, + const int* elem_ids, + const int nMPs_delete, + const int nMPs_add, + const int* recvMPs_elm, + const int* recvMPs_ids) { + checkMPMeshValid(p_mpmesh); + auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; + int offset = p_MPs->getElmIDoffset(); + printf("To delete %d MPs add %d MPs\n", nMPs_delete, nMPs_add); + auto elem_ids_d = create_mirror_view_and_copy(elem_ids, sizeMP2elm); + auto recvMPs_elm_d = create_mirror_view_and_copy(recvMPs_elm, nMPs_add); + auto recvMPs_ids_d = create_mirror_view_and_copy(recvMPs_ids, nMPs_add); + + Kokkos::View mp2Elm("mp2Elm", p_MPs->getCapacity()); + Kokkos::View numDeletedMPs_d("numDeletedMPs", 1); + auto mpAppID = p_MPs->getData(); + + auto setMP2Elm = PS_LAMBDA(const int& e, const int& mp, const int& mask) { + if(mask) { + mp2Elm(mp) = elem_ids_d(mpAppID(mp)); + if (mp2Elm(mp) == MP_DELETE){ + Kokkos::atomic_increment(&numDeletedMPs_d(0)); + printf("Deleted particle is %d it was in %d \n", mpAppID(mp), e); + } + } + }; + p_MPs->parallel_for(setMP2Elm, "setMP2Elm"); + int numDeletedMPs = pumipic::getLastValue(numDeletedMPs_d); + if(numDeletedMPs) printf("Dleted %d particles \n", numDeletedMPs); + + p_MPs->startRebuild(mp2Elm, nMPs_add, recvMPs_elm_d, recvMPs_ids_d); +} + + + void polympo_finishRebuildMPs_f(MPMesh_ptr p_mpmesh) { checkMPMeshValid(p_mpmesh); @@ -255,13 +293,14 @@ void polympo_setMPLatLonRotatedFlag_f(MPMesh_ptr p_mpmesh, const int isRotateFla void polympo_setMPPositions_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, - const double* mpPositionsIn){ + const double* mpPositionsIn, + const int setTgtPos){ Kokkos::Timer timer; checkMPMeshValid(p_mpmesh); auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; PMT_ALWAYS_ASSERT(nComps == vec3d_nEntries); PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount()); - PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID()); + //PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID()); kkViewHostU mpPositionsIn_h(mpPositionsIn,nComps,numMPs); if (p_MPs->rebuildOngoing()) { @@ -270,6 +309,10 @@ void polympo_setMPPositions_f(MPMesh_ptr p_mpmesh, } auto mpPositions = p_MPs->getData(); + if(setTgtPos){ + printf("Setting the TGT position instead\n"); + mpPositions = p_MPs->getData(); + } auto mpAppID = p_MPs->getData(); Kokkos::View mpPositionsIn_d("mpPositionsDevice",vec3d_nEntries,numMPs); Kokkos::deep_copy(mpPositionsIn_d, mpPositionsIn_h); @@ -281,6 +324,7 @@ void polympo_setMPPositions_f(MPMesh_ptr p_mpmesh, } }; p_MPs->parallel_for(setPos, "setMPPositions"); + printf("Done setting MP pos\n"); pumipic::RecordTime("PolyMPO_setMPPositions", timer.seconds()); } @@ -987,6 +1031,12 @@ void polympo_getMeshVtxOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, c Kokkos::deep_copy(arrayHost, array_d); } +bool polympo_push1P_f(MPMesh_ptr p_mpmesh){ + checkMPMeshValid(p_mpmesh); + bool is_migrating=((polyMPO::MPMesh*)p_mpmesh) ->push1P(); + return is_migrating; +} + void polympo_push_f(MPMesh_ptr p_mpmesh){ checkMPMeshValid(p_mpmesh); ((polyMPO::MPMesh*)p_mpmesh) ->push(); @@ -1068,6 +1118,11 @@ void polympo_setElmGlobal_f(MPMesh_ptr p_mpmesh, const int nCells, const int* ar } +int polympo_getMPCount_f(MPMesh_ptr p_mpmesh) { + auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; + return p_MPs->getCount(); // This returns the int you're after +} + void polympo_enableTiming_f(){ pumipic::EnableTiming(); } diff --git a/src/pmpo_c.h b/src/pmpo_c.h index 6dcbf4e..3ab2331 100644 --- a/src/pmpo_c.h +++ b/src/pmpo_c.h @@ -21,6 +21,10 @@ void polympo_setMPICommunicator_f(MPMesh_ptr p_mpmesh, MPI_Fint fcomm); //MP info void polympo_createMPs_f(MPMesh_ptr p_mpmesh, const int numElms, const int numMPs, int* mpsPerElm, const int* mp2Elm, const int* isMPActive); void polympo_startRebuildMPs_f(MPMesh_ptr p_mpmesh, const int numMPs, const int* allTgtMpElmIn, const int* addedMPMask); +void polympo_startRebuildMPs_f2(MPMesh_ptr p_mpmesh, const int size1, const int* arg1, const int size2, const int size3, const int* arg2, const int* arg3); + +int polympo_getMPCount_f(MPMesh_ptr p_mpmesh); + void polympo_finishRebuildMPs_f(MPMesh_ptr p_mpmesh); void polympo_setAppIDFunc_f(MPMesh_ptr p_mpmesh, IntVoidFunc getNext, void* appIDs); @@ -30,7 +34,7 @@ void polympo_getMPCurElmID_f(MPMesh_ptr p_mpmesh, const int numMPs, int* elmIDs) void polympo_setMPLatLonRotatedFlag_f(MPMesh_ptr p_mpmesh, const int isRotateFlag); //MP slices -void polympo_setMPPositions_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpPositionsIn); +void polympo_setMPPositions_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpPositionsIn, const int setTgtPos); void polympo_getMPPositions_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpPositionsIn); void polympo_setMPRotLatLon_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpRotLatLonIn); void polympo_getMPRotLatLon_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpRotLatLonHost); @@ -90,6 +94,7 @@ void polympo_getMeshDualTriangleArea_f(MPMesh_ptr p_mpmesh, const int nVertices, // calculations void polympo_push_f(MPMesh_ptr p_mpmesh); +bool polympo_push1P_f(MPMesh_ptr p_mpmesh); // Reconstruction of variables from MPs to mesh vertices void polympo_setReconstructionOfMass_f(MPMesh_ptr p_mpmesh, const int order, const int meshEntType); diff --git a/src/pmpo_fortran.f90 b/src/pmpo_fortran.f90 index c0cd8ca..4efe641 100644 --- a/src/pmpo_fortran.f90 +++ b/src/pmpo_fortran.f90 @@ -87,6 +87,20 @@ subroutine polympo_startRebuildMPs(mpMesh, numMPs, allMP2Elm, addedMPMask) & type(c_ptr), intent(in), value :: allMP2Elm type(c_ptr), intent(in), value :: addedMPMask end subroutine + + + subroutine polympo_startRebuildMPs2(mpMesh, size1, arg1, size2, size3, arg2, arg3) & + bind(C, NAME='polympo_startRebuildMPs_f2') + use :: iso_c_binding + type(c_ptr), value :: mpMesh + integer(c_int), value :: size1 + type(c_ptr), intent(in), value :: arg1 + integer(c_int), value :: size2 + integer(c_int), value :: size3 + type(c_ptr), intent(in), value :: arg2 + type(c_ptr), intent(in), value :: arg3 + end subroutine + ! !--------------------------------------------------------------------------- !> @brief called after startRebuild() !> @brief called after initializing MP fields @@ -160,12 +174,13 @@ subroutine polympo_setMPLatLonRotatedFlag(mpMesh, isRotateFlag) & !> @param numMPs(in) number of the MPs !> @param array(in) MP current position 2D array (3,numMPs), allocated by user on host !--------------------------------------------------------------------------- - subroutine polympo_setMPPositions(mpMesh, nComps, numMPs, array) & + subroutine polympo_setMPPositions(mpMesh, nComps, numMPs, array, setTgtPos) & bind(C, NAME='polympo_setMPPositions_f') use :: iso_c_binding type(c_ptr), value :: mpMesh integer(c_int), value :: nComps, numMPs type(c_ptr), value :: array + integer(c_int), value:: setTgtPos end subroutine !--------------------------------------------------------------------------- !> @brief get the MP positions array from a polympo array @@ -743,6 +758,16 @@ subroutine polympo_push(mpMesh) & use :: iso_c_binding type(c_ptr), value :: mpMesh end subroutine + + + !--------------------------------------------------------------------------- + !> @brief calculate the MPs from given mesh vertices rotational latitude + logical function polympo_push1P(mpMesh) & + bind(C, NAME='polympo_push1P_f') + use :: iso_c_binding + type(c_ptr), value :: mpMesh + end function + !--------------------------------------------------------------------------- !> @brief start the reconstruction of MP Mass to Mesh Vertices !> @param mpmesh(in/out) MPMesh object @@ -823,6 +848,13 @@ subroutine polympo_setTimingVerbosity(v) bind(C, NAME='polympo_setTimingVerbosit use :: iso_c_binding integer(c_int), value :: v end subroutine + + integer function polympo_getMPCount(mpMesh) bind(C, name="polympo_getMPCount_f") + use :: iso_c_binding + implicit none + type(c_ptr), value :: mpMesh + end function + end interface contains !--------------------------------------------------------------------------- @@ -838,4 +870,5 @@ subroutine polympo_checkPrecisionForRealKind(APP_RKIND) call exit(1) end if end subroutine + end module diff --git a/src/pmpo_materialPoints.cpp b/src/pmpo_materialPoints.cpp index 76218aa..aca7f2e 100644 --- a/src/pmpo_materialPoints.cpp +++ b/src/pmpo_materialPoints.cpp @@ -95,6 +95,16 @@ void MaterialPoints::startRebuild(IntView tgtElm, int addedNumMPs, IntView added rebuildFields.addedSlices_h = createInternalMemberViews(addedNumMPs, addedMP2elm_h, addedMPAppID_h); } +void MaterialPoints::startRebuild(IntView tgtElm, int addedNumMPs, IntView addedMP2elm, IntView addedMPAppID) { + rebuildFields.ongoing = true; + rebuildFields.addedNumMPs = addedNumMPs; + rebuildFields.addedMP2elm = addedMP2elm; + rebuildFields.allTgtElm = tgtElm; + auto addedMP2elm_h = Kokkos::create_mirror_view_and_copy(hostSpace(), addedMP2elm); + auto addedMPAppID_h = Kokkos::create_mirror_view_and_copy(hostSpace(), addedMPAppID); + rebuildFields.addedSlices_h = createInternalMemberViews(addedNumMPs, addedMP2elm_h, addedMPAppID_h); +} + void MaterialPoints::finishRebuild() { auto addedSlices_d = ps::createMemberViews(rebuildFields.addedNumMPs); ps::CopyMemSpaceToMemSpace(addedSlices_d, rebuildFields.addedSlices_h); @@ -114,20 +124,20 @@ void MaterialPoints::setMPIComm(MPI_Comm comm) { } bool MaterialPoints::check_migrate(){ - Kokkos::Timer timer; - auto MPs2Elm = getData(); + + Kokkos::Timer timer; auto MPs2Proc = getData(); - IntView isMigrating("isMigrating", 1); - int rank; MPI_Comm_rank(mpi_comm, &rank); + auto setMigrationFields = PS_LAMBDA(const int& e, const int& mp, const bool& mask) { if (mask) { if (MPs2Proc(mp) != rank) isMigrating(0) = 1; } }; parallel_for(setMigrationFields, "setMigrationFields"); + if (getOpMode() == polyMPO::MP_DEBUG) printf("Material point check migration: %f\n", timer.seconds()); pumipic::RecordTime("PolyMPO_check_migrate", timer.seconds()); @@ -150,6 +160,7 @@ void MaterialPoints::migrate() { new_elem(mp) = MPs2Elm(mp); new_process(mp) = MPs2Proc(mp); if(rank!=new_process(mp)){ + printf("Particle %d in rank %d to be moved from %d %d \n", mpAppID(mp), rank, e, new_elem(mp) ); mpAppID(mp)=-1; } } diff --git a/src/pmpo_materialPoints.hpp b/src/pmpo_materialPoints.hpp index 13a9a86..c2a037f 100644 --- a/src/pmpo_materialPoints.hpp +++ b/src/pmpo_materialPoints.hpp @@ -139,6 +139,8 @@ class MaterialPoints { void rebuild(IntView addedMP2elm, IntView addedMPAppID); void startRebuild(IntView tgtElm, int addedNumMPs, IntView addedMP2elm, IntView addedMPAppID, Kokkos::View addedMPMask); + void startRebuild(IntView tgtElm, int addedNumMPs, IntView addedMP2elm, IntView addedMPAppID); + void finishRebuild(); bool rebuildOngoing(); @@ -163,9 +165,12 @@ class MaterialPoints { void rebuild() { IntView tgtElm("tgtElm", MPs->capacity()); auto tgtMpElm = MPs->get(); - auto setTgtElm = PS_LAMBDA(const int&, const int& mp, const int& mask) { + auto mpAppID = MPs->get(); + auto setTgtElm = PS_LAMBDA(const int& e, const int& mp, const int& mask) { if(mask) { + auto app_id=mpAppID(mp); tgtElm(mp) = tgtMpElm(mp); + if(app_id==191) printf("Going from elm %d to %d \n", e, tgtMpElm(mp)); } }; ps::parallel_for(MPs, setTgtElm, "setTargetElement"); From 77cde5b1ae1d00be4e1198149ee4dc1e272f3d85 Mon Sep 17 00:00:00 2001 From: nathd2 Date: Mon, 5 May 2025 22:14:27 -0400 Subject: [PATCH 13/17] Woking 2 proc case --- src/pmpo_MPMesh.cpp | 54 ++++++---- src/pmpo_MPMesh.hpp | 3 + src/pmpo_c.cpp | 204 ++++++++++++++++++++++++++++++++---- src/pmpo_c.h | 15 ++- src/pmpo_fortran.f90 | 105 ++++++++++++++++++- src/pmpo_materialPoints.cpp | 18 ++++ 6 files changed, 355 insertions(+), 44 deletions(-) diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index f15f2ff..f2fa0c6 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -128,7 +128,8 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){ auto MPs2Proc = p_MPs->getData(); auto elm2Process = p_mesh->getElm2Process(); auto elm2global = p_mesh->getElmGlobal(); - + auto mpAppID = p_MPs->getData(); + MPI_Comm comm = p_MPs->getMPIComm(); int comm_rank; MPI_Comm_rank(comm, &comm_rank); @@ -143,12 +144,14 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){ //assert(cudaDeviceSynchronize()==cudaSuccess); //MPI_Barrier(MPI_COMM_WORLD); - printf("NumMPs %d \n", numMPs); Vec3dView history("positionHistory",numMPs); Vec3dView resultLeft("positionResult",numMPs); Vec3dView resultRight("positionResult",numMPs); Vec3dView mpTgtPosArray("positionTarget",numMPs); Kokkos::View counter("counter",1); + + assert(cudaDeviceSynchronize() == cudaSuccess); + //printf("Rank %d Foo4 Begin\n", comm_rank); auto CVTElmCalc = PS_LAMBDA(const int& elm, const int& mp, const int&mask){ Vec3d MP(mpPositions(mp,0),mpPositions(mp,1),mpPositions(mp,2)); @@ -181,13 +184,15 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){ if(closestElm<0){ MPs2Elm(mp) = iElm; - if (elm2Process.size() > 0) - MPs2Proc(mp) = elm2Process(iElm); + MPs2Proc(mp) = elm2Process(iElm); break; }else{ iElm = closestElm; } } + if(mpAppID(mp)==0 || mpAppID(mp)==191) + printf("Pos %.15e %.15e %.15e => %.15e %.15e %.15e\n", mpPositions(mp,0), mpPositions(mp,1), + mpPositions(mp,2), mpTgtPos(mp,0), mpTgtPos(mp,1), mpTgtPos(mp, 2)); if(printVTPIndex>=0 && numMPs>0){ //printf("Rank %d mp %d counter %d \n", comm_rank, mp, counter); @@ -211,6 +216,10 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){ } }; p_MPs->parallel_for(CVTElmCalc,"CVTTrackingElmCenterBasedCalc"); + + assert(cudaDeviceSynchronize() == cudaSuccess); + //printf("Rank %d Foo4 End\n", comm_rank); + if(printVTPIndex>=0 && numMPs>0){ Vec3dView::HostMirror h_history = Kokkos::create_mirror_view(history); @@ -341,35 +350,44 @@ bool getAnyIsMigrating(MaterialPoints* p_MPs, bool isMigrating) { return anyIsMigrating; } -bool MPMesh::push1P(){ - - static int count=0; - std::cout<<"Push1P"<<" "<computeRotLatLonIncr(); sphericalInterpolation(*this); - //Push the MPs p_MPs->updateRotLatLonAndXYZ2Tgt(p_mesh->getSphereRadius()); - - //Given target location find the new MP element and the process it belongs to + count0 ++; +} + +bool MPMesh::push1P(){ + //Given target location find the new element or the last element in a partioned mesh + //and the process it belongs to so that migration can be checked + static int count_p=0; CVTTrackingElmCenterBased(); + //From the above two inputs find if any particle needs to be migrated bool anyIsMigrating = getAnyIsMigrating(p_MPs, p_MPs->check_migrate()); + count_p=count_p+1; + return anyIsMigrating; - //New element (maynot be final) so that MPAS can do the migration +} + +void MPMesh::push_swap(){ + //current becomes target, target becomes -1 p_MPs->updateMPElmID(); - +} + +void MPMesh::push_swap_pos(){ + //current becomes target, target becomes -1 + //Making read for next push_ahead p_MPs->updateMPSlice(); p_MPs->updateMPSlice(); - - count++; - return anyIsMigrating; - } + void MPMesh::push(){ static int count=0; diff --git a/src/pmpo_MPMesh.hpp b/src/pmpo_MPMesh.hpp index b30bcf3..7b54afb 100644 --- a/src/pmpo_MPMesh.hpp +++ b/src/pmpo_MPMesh.hpp @@ -44,6 +44,9 @@ class MPMesh{ void CVTTrackingElmCenterBased(const int printVTPIndex = -1); void T2LTracking(Vec2dView dx); bool push1P(); + void push_ahead(); + void push_swap(); + void push_swap_pos(); void push(); void calcBasis(); diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index 85046a8..ce7b7ba 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -204,12 +204,24 @@ void polympo_startRebuildMPs_f2(MPMesh_ptr p_mpmesh, const int* elem_ids, const int nMPs_delete, const int nMPs_add, - const int* recvMPs_elm, - const int* recvMPs_ids) { + int* recvMPs_elm, + int* recvMPs_ids) { checkMPMeshValid(p_mpmesh); auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; int offset = p_MPs->getElmIDoffset(); - printf("To delete %d MPs add %d MPs\n", nMPs_delete, nMPs_add); + + int rank; + auto mpi_comm=p_MPs->getMPIComm(); + MPI_Comm_rank(mpi_comm, &rank); + + printf("Rank %d to add %d and delete %d ptcls \n", rank, nMPs_add, nMPs_delete); + for (int k=0; k mp2Elm("mp2Elm", p_MPs->getCapacity()); Kokkos::View numDeletedMPs_d("numDeletedMPs", 1); auto mpAppID = p_MPs->getData(); - auto setMP2Elm = PS_LAMBDA(const int& e, const int& mp, const int& mask) { if(mask) { - mp2Elm(mp) = elem_ids_d(mpAppID(mp)); - if (mp2Elm(mp) == MP_DELETE){ + mp2Elm(mp) = elem_ids_d(mpAppID(mp))-offset; + if (mp2Elm(mp) == MP_DELETE){ + //printf("Rank %d to delete %d particle in element %d \n", rank, mpAppID(mp), e); Kokkos::atomic_increment(&numDeletedMPs_d(0)); - printf("Deleted particle is %d it was in %d \n", mpAppID(mp), e); - } + } } }; p_MPs->parallel_for(setMP2Elm, "setMP2Elm"); int numDeletedMPs = pumipic::getLastValue(numDeletedMPs_d); - if(numDeletedMPs) printf("Dleted %d particles \n", numDeletedMPs); - + assert(nMPs_delete==numDeletedMPs); p_MPs->startRebuild(mp2Elm, nMPs_add, recvMPs_elm_d, recvMPs_ids_d); } @@ -260,13 +270,42 @@ void polympo_setMPASAppIDFunc_f(MPMesh_ptr p_mpmesh, VoidVoidFunc getMPASAppID, p_MPs->setAppIDFunc(polyMPO_getMPASAppID); } +void polympo_getMPTgtElmID_f(MPMesh_ptr p_mpmesh, + const int numMPs, + int* elmIDs){ + checkMPMeshValid(p_mpmesh); + auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; + PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount()); + //PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID()); + auto mpTgtElmID = p_MPs->getData(); + auto mpAppID = p_MPs->getData(); + auto elmIDoffset = p_MPs->getElmIDoffset(); + + kkIntViewHostU arrayHost(elmIDs,numMPs); + polyMPO::IntView mpTgtElmIDCopy("mpTgtElmIDNewValue",numMPs); + + int rank; + auto mpi_comm=p_MPs->getMPIComm(); + MPI_Comm_rank(mpi_comm, &rank); + + auto setTgtElmId = PS_LAMBDA(const int&, const int& mp, const int& mask){ + if(mask){ + mpTgtElmIDCopy(mpAppID(mp)) = mpTgtElmID(mp)+elmIDoffset; + } + }; + p_MPs->parallel_for(setTgtElmId, "set mpTgtElmID"); + Kokkos::deep_copy( arrayHost, mpTgtElmIDCopy); + +} + void polympo_getMPCurElmID_f(MPMesh_ptr p_mpmesh, const int numMPs, int* elmIDs){ + checkMPMeshValid(p_mpmesh); auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount()); - PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID()); + //PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID()); auto mpCurElmID = p_MPs->getData(); auto mpAppID = p_MPs->getData(); auto elmIDoffset = p_MPs->getElmIDoffset(); @@ -274,13 +313,18 @@ void polympo_getMPCurElmID_f(MPMesh_ptr p_mpmesh, kkIntViewHostU arrayHost(elmIDs,numMPs); polyMPO::IntView mpCurElmIDCopy("mpCurElmIDNewValue",numMPs); + int rank; + auto mpi_comm=p_MPs->getMPIComm(); + MPI_Comm_rank(mpi_comm, &rank); + auto getElmId = PS_LAMBDA(const int&, const int& mp, const int& mask){ if(mask){ - mpCurElmIDCopy(mpAppID(mp)) = mpCurElmID(mp)+elmIDoffset; + mpCurElmIDCopy(mpAppID(mp)) = mpCurElmID(mp)+elmIDoffset; } }; p_MPs->parallel_for(getElmId, "get mpCurElmID"); Kokkos::deep_copy( arrayHost, mpCurElmIDCopy); + } void polympo_setMPLatLonRotatedFlag_f(MPMesh_ptr p_mpmesh, const int isRotateFlag){ @@ -293,8 +337,7 @@ void polympo_setMPLatLonRotatedFlag_f(MPMesh_ptr p_mpmesh, const int isRotateFla void polympo_setMPPositions_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, - const double* mpPositionsIn, - const int setTgtPos){ + const double* mpPositionsIn){ Kokkos::Timer timer; checkMPMeshValid(p_mpmesh); auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; @@ -309,11 +352,8 @@ void polympo_setMPPositions_f(MPMesh_ptr p_mpmesh, } auto mpPositions = p_MPs->getData(); - if(setTgtPos){ - printf("Setting the TGT position instead\n"); - mpPositions = p_MPs->getData(); - } auto mpAppID = p_MPs->getData(); + Kokkos::View mpPositionsIn_d("mpPositionsDevice",vec3d_nEntries,numMPs); Kokkos::deep_copy(mpPositionsIn_d, mpPositionsIn_h); auto setPos = PS_LAMBDA(const int&, const int& mp, const int& mask){ @@ -324,7 +364,6 @@ void polympo_setMPPositions_f(MPMesh_ptr p_mpmesh, } }; p_MPs->parallel_for(setPos, "setMPPositions"); - printf("Done setting MP pos\n"); pumipic::RecordTime("PolyMPO_setMPPositions", timer.seconds()); } @@ -336,7 +375,7 @@ void polympo_getMPPositions_f(MPMesh_ptr p_mpmesh, auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; PMT_ALWAYS_ASSERT(nComps == vec3d_nEntries); PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount()); - PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID()); + //PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID()); auto mpPositions = p_MPs->getData(); auto mpAppID = p_MPs->getData(); @@ -353,6 +392,66 @@ void polympo_getMPPositions_f(MPMesh_ptr p_mpmesh, Kokkos::deep_copy(arrayHost, mpPositionsCopy); } + +void polympo_setMPTgtPositions_f(MPMesh_ptr p_mpmesh, + const int nComps, + const int numMPs, + const double* mpPositionsIn){ + Kokkos::Timer timer; + checkMPMeshValid(p_mpmesh); + auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; + PMT_ALWAYS_ASSERT(nComps == vec3d_nEntries); + PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount()); + //printf("NumMPs %d %d \n", numMPs, p_MPs->getCount()); + //PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID()); + kkViewHostU mpPositionsIn_h(mpPositionsIn,nComps,numMPs); + + if (p_MPs->rebuildOngoing()) { + p_MPs->setRebuildMPSlice(mpPositionsIn_h); + return; + } + + auto mpPositions = p_MPs->getData(); + auto mpAppID = p_MPs->getData(); + Kokkos::View mpPositionsIn_d("mpPositionsDevice",vec3d_nEntries,numMPs); + Kokkos::deep_copy(mpPositionsIn_d, mpPositionsIn_h); + auto setPos = PS_LAMBDA(const int&, const int& mp, const int& mask){ + if(mask){ + mpPositions(mp,0) = mpPositionsIn_d(0, mpAppID(mp)); + mpPositions(mp,1) = mpPositionsIn_d(1, mpAppID(mp)); + mpPositions(mp,2) = mpPositionsIn_d(2, mpAppID(mp)); + } + }; + p_MPs->parallel_for(setPos, "setMPPositions"); + pumipic::RecordTime("PolyMPO_setMPPositions", timer.seconds()); +} + +void polympo_getMPTgtPositions_f(MPMesh_ptr p_mpmesh, + const int nComps, + const int numMPs, + double* mpPositionsHost){ + checkMPMeshValid(p_mpmesh); + auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; + PMT_ALWAYS_ASSERT(nComps == vec3d_nEntries); + PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount()); + //PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID()); + + auto mpPositions = p_MPs->getData(); + auto mpAppID = p_MPs->getData(); + Kokkos::View mpPositionsCopy("mpPositionsCopy",vec3d_nEntries,numMPs); + auto getPos = PS_LAMBDA(const int&, const int& mp, const int& mask){ + if(mask){ + mpPositionsCopy(0,mpAppID(mp)) = mpPositions(mp,0); + mpPositionsCopy(1,mpAppID(mp)) = mpPositions(mp,1); + mpPositionsCopy(2,mpAppID(mp)) = mpPositions(mp,2); + } + }; + p_MPs->parallel_for(getPos, "getMPPositions"); + kkDbl2dViewHostU arrayHost(mpPositionsHost,nComps,numMPs); + Kokkos::deep_copy(arrayHost, mpPositionsCopy); +} + + void polympo_setMPRotLatLon_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, @@ -404,6 +503,56 @@ void polympo_getMPRotLatLon_f(MPMesh_ptr p_mpmesh, Kokkos::deep_copy(arrayHost, mpRotLatLonCopy); } + +void polympo_setMPTgtRotLatLon_f(MPMesh_ptr p_mpmesh, + const int nComps, + const int numMPs, + const double* mpRotLatLonIn){ + checkMPMeshValid(p_mpmesh); + auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; + PMT_ALWAYS_ASSERT(nComps == vec2d_nEntries); + PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount()); + //PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID()); + + auto mpRotLatLon = p_MPs->getData(); + auto mpAppID = p_MPs->getData(); + kkViewHostU mpRotLatLonIn_h(mpRotLatLonIn,nComps,numMPs); + Kokkos::View mpRotLatLonIn_d("mpRotLatLonDevice",vec2d_nEntries,numMPs); + Kokkos::deep_copy(mpRotLatLonIn_d, mpRotLatLonIn_h); + auto setPos = PS_LAMBDA(const int& elm, const int& mp, const int& mask){ + if(mask){ + mpRotLatLon(mp,0) = mpRotLatLonIn_d(0, mpAppID(mp)); + mpRotLatLon(mp,1) = mpRotLatLonIn_d(1, mpAppID(mp)); + } + }; + p_MPs->parallel_for(setPos, "setMPRotLatLon"); +} + +void polympo_getMPTgtRotLatLon_f(MPMesh_ptr p_mpmesh, + const int nComps, + const int numMPs, + double* mpRotLatLonHost){ + checkMPMeshValid(p_mpmesh); + auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; + PMT_ALWAYS_ASSERT(nComps == vec2d_nEntries); + PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount()); + //PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID()); + + auto mpRotLatLon = p_MPs->getData(); + auto mpAppID = p_MPs->getData(); + Kokkos::View mpRotLatLonCopy("mpRotLatLonCopy",vec2d_nEntries,numMPs); + auto getPos = PS_LAMBDA(const int& elm, const int& mp, const int& mask){ + if(mask){ + mpRotLatLonCopy(0,mpAppID(mp)) = mpRotLatLon(mp,0); + mpRotLatLonCopy(1,mpAppID(mp)) = mpRotLatLon(mp,1); + } + }; + p_MPs->parallel_for(getPos, "getMPRotLatLon"); + kkDbl2dViewHostU arrayHost(mpRotLatLonHost,nComps,numMPs); + Kokkos::deep_copy(arrayHost, mpRotLatLonCopy); +} + + void polympo_setMPMass_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpMassIn) { Kokkos::Timer timer; checkMPMeshValid(p_mpmesh); @@ -1037,6 +1186,21 @@ bool polympo_push1P_f(MPMesh_ptr p_mpmesh){ return is_migrating; } +void polympo_push_ahead_f(MPMesh_ptr p_mpmesh){ + checkMPMeshValid(p_mpmesh); + ((polyMPO::MPMesh*)p_mpmesh) ->push_ahead(); +} + +void polympo_push_swap_f(MPMesh_ptr p_mpmesh){ + checkMPMeshValid(p_mpmesh); + ((polyMPO::MPMesh*)p_mpmesh) ->push_swap(); +} + +void polympo_push_swap_pos_f(MPMesh_ptr p_mpmesh){ + checkMPMeshValid(p_mpmesh); + ((polyMPO::MPMesh*)p_mpmesh) ->push_swap_pos(); +} + void polympo_push_f(MPMesh_ptr p_mpmesh){ checkMPMeshValid(p_mpmesh); ((polyMPO::MPMesh*)p_mpmesh) ->push(); diff --git a/src/pmpo_c.h b/src/pmpo_c.h index 3ab2331..cfbe628 100644 --- a/src/pmpo_c.h +++ b/src/pmpo_c.h @@ -21,7 +21,7 @@ void polympo_setMPICommunicator_f(MPMesh_ptr p_mpmesh, MPI_Fint fcomm); //MP info void polympo_createMPs_f(MPMesh_ptr p_mpmesh, const int numElms, const int numMPs, int* mpsPerElm, const int* mp2Elm, const int* isMPActive); void polympo_startRebuildMPs_f(MPMesh_ptr p_mpmesh, const int numMPs, const int* allTgtMpElmIn, const int* addedMPMask); -void polympo_startRebuildMPs_f2(MPMesh_ptr p_mpmesh, const int size1, const int* arg1, const int size2, const int size3, const int* arg2, const int* arg3); +void polympo_startRebuildMPs_f2(MPMesh_ptr p_mpmesh, const int size1, const int* arg1, const int size2, const int size3, int* arg2, int* arg3); int polympo_getMPCount_f(MPMesh_ptr p_mpmesh); @@ -30,14 +30,22 @@ void polympo_finishRebuildMPs_f(MPMesh_ptr p_mpmesh); void polympo_setAppIDFunc_f(MPMesh_ptr p_mpmesh, IntVoidFunc getNext, void* appIDs); void polympo_setMPASAppIDFunc_f(MPMesh_ptr p_mpmesh, VoidVoidFunc getMPASAppID, void* arg1, const int arg2, const int arg3); +void polympo_getMPTgtElmID_f(MPMesh_ptr p_mpmesh, const int numMPs, int* elmIDs); void polympo_getMPCurElmID_f(MPMesh_ptr p_mpmesh, const int numMPs, int* elmIDs); void polympo_setMPLatLonRotatedFlag_f(MPMesh_ptr p_mpmesh, const int isRotateFlag); //MP slices -void polympo_setMPPositions_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpPositionsIn, const int setTgtPos); +//Positions +void polympo_setMPPositions_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpPositionsIn); void polympo_getMPPositions_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpPositionsIn); +void polympo_setMPTgtPositions_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpPositionsIn); +void polympo_getMPTgtPositions_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpPositionsIn); +//LatLon void polympo_setMPRotLatLon_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpRotLatLonIn); void polympo_getMPRotLatLon_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpRotLatLonHost); +void polympo_setMPTgtRotLatLon_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpRotLatLonIn); +void polympo_getMPTgtRotLatLon_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpRotLatLonHost); + void polympo_setMPMass_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpMassIn); void polympo_getMPMass_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpMassHost); void polympo_setMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpVelIn); @@ -94,7 +102,10 @@ void polympo_getMeshDualTriangleArea_f(MPMesh_ptr p_mpmesh, const int nVertices, // calculations void polympo_push_f(MPMesh_ptr p_mpmesh); +void polympo_push_ahead_f(MPMesh_ptr p_mpmesh); bool polympo_push1P_f(MPMesh_ptr p_mpmesh); +void polympo_push_swap_f(MPMesh_ptr p_mpmesh); +void polympo_push_swap_pos_f(MPMesh_ptr p_mpmesh); // Reconstruction of variables from MPs to mesh vertices void polympo_setReconstructionOfMass_f(MPMesh_ptr p_mpmesh, const int order, const int meshEntType); diff --git a/src/pmpo_fortran.f90 b/src/pmpo_fortran.f90 index 4efe641..b05ce49 100644 --- a/src/pmpo_fortran.f90 +++ b/src/pmpo_fortran.f90 @@ -97,8 +97,8 @@ subroutine polympo_startRebuildMPs2(mpMesh, size1, arg1, size2, size3, arg2, arg type(c_ptr), intent(in), value :: arg1 integer(c_int), value :: size2 integer(c_int), value :: size3 - type(c_ptr), intent(in), value :: arg2 - type(c_ptr), intent(in), value :: arg3 + type(c_ptr), value :: arg2 + type(c_ptr), value :: arg3 end subroutine ! !--------------------------------------------------------------------------- @@ -156,6 +156,16 @@ subroutine polympo_getMPCurElmID(mpMesh, numMPs, array) & integer(c_int), value :: numMPs type(c_ptr), value :: array end subroutine + + + subroutine polympo_getMPTgtElmID(mpMesh, numMPs, array) & + bind(C, NAME='polympo_getMPTgtElmID_f') + use :: iso_c_binding + type(c_ptr), value :: mpMesh + integer(c_int), value :: numMPs + type(c_ptr), value :: array + end subroutine + ! !--------------------------------------------------------------------------- !> @brief set the mp lat lon is rotational or normal !> @param mpmesh(in/out) MPMesh object @@ -167,6 +177,8 @@ subroutine polympo_setMPLatLonRotatedFlag(mpMesh, isRotateFlag) & type(c_ptr), value :: mpMesh integer(c_int), value :: isRotateFlag end subroutine + + !--------------------------------------------------------------------------- !> @brief set the MP positions array from a host array !> @param mpmesh(in/out) MPMesh object @@ -174,13 +186,12 @@ subroutine polympo_setMPLatLonRotatedFlag(mpMesh, isRotateFlag) & !> @param numMPs(in) number of the MPs !> @param array(in) MP current position 2D array (3,numMPs), allocated by user on host !--------------------------------------------------------------------------- - subroutine polympo_setMPPositions(mpMesh, nComps, numMPs, array, setTgtPos) & + subroutine polympo_setMPPositions(mpMesh, nComps, numMPs, array) & bind(C, NAME='polympo_setMPPositions_f') use :: iso_c_binding type(c_ptr), value :: mpMesh integer(c_int), value :: nComps, numMPs type(c_ptr), value :: array - integer(c_int), value:: setTgtPos end subroutine !--------------------------------------------------------------------------- !> @brief get the MP positions array from a polympo array @@ -197,6 +208,37 @@ subroutine polympo_getMPPositions(mpMesh, nComps, numMPs, array) & integer(c_int), value :: nComps, numMPs type(c_ptr), value :: array end subroutine + !--------------------------------------------------------------------------- + !> @brief set the MP positions array from a host array + !> @param mpmesh(in/out) MPMesh object + !> @param nComps(in) number of components, should always be 3 + !> @param numMPs(in) number of the MPs + !> @param array(in) MP current position 2D array (3,numMPs), allocated by user on host + !--------------------------------------------------------------------------- + subroutine polympo_setMPTgtPositions(mpMesh, nComps, numMPs, array) & + bind(C, NAME='polympo_setMPTgtPositions_f') + use :: iso_c_binding + type(c_ptr), value :: mpMesh + integer(c_int), value :: nComps, numMPs + type(c_ptr), value :: array + end subroutine + !--------------------------------------------------------------------------- + !> @brief get the MP positions array from a polympo array + !> @param mpmesh(in/out) MPMesh object + !> @param nComps(in) number of components, should always be 3 + !> @param numMPs(in) number of the MPs + !> @param array(in/out) output MP current position 2D array (3,numMPs), + !> allocated by user + !--------------------------------------------------------------------------- + subroutine polympo_getMPTgtPositions(mpMesh, nComps, numMPs, array) & + bind(C, NAME='polympo_getMPTgtPositions_f') + use :: iso_c_binding + type(c_ptr), value :: mpMesh + integer(c_int), value :: nComps, numMPs + type(c_ptr), value :: array + end subroutine + + !--------------------------------------------------------------------------- !> @brief set the MP latitude and longtitude array from a host array !> @param mpmesh(in/out) MPMesh object @@ -227,6 +269,40 @@ subroutine polympo_getMPRotLatLon(mpMesh, nComps, numMPs, array) & integer(c_int), value :: nComps, numMPs type(c_ptr), value :: array end subroutine + !--------------------------------------------------------------------------- + !> @brief set the MP latitude and longtitude array from a host array + !> @param mpmesh(in/out) MPMesh object + !> @param nComps(in) number of components, should always be 2 + !> @param numMPs(in) number of the MPs + !> @param array(in) input MP current lat and lon 2D array (2,numMPs), + !> allocated by user + !--------------------------------------------------------------------------- + subroutine polympo_setMPTgtRotLatLon(mpMesh, nComps, numMPs, array) & + bind(C, NAME='polympo_setMPTgtRotLatLon_f') + use :: iso_c_binding + type(c_ptr), value :: mpMesh + integer(c_int), value :: nComps, numMPs + type(c_ptr), intent(in), value :: array + end subroutine + !--------------------------------------------------------------------------- + !> @brief get the MP latitude and longtitude array from a polympo array + !> @param mpmesh(in/out) MPMesh object + !> @param nComps(in) number of components, should always be 2 + !> @param numMPs(in) number of the MPs + !> @param array(in/out) output MP current lat and lon 2D array (2,numMPs), + !> allocated by user + !--------------------------------------------------------------------------- + subroutine polympo_getMPTgtRotLatLon(mpMesh, nComps, numMPs, array) & + bind(C, NAME='polympo_getMPTgtRotLatLon_f') + use :: iso_c_binding + type(c_ptr), value :: mpMesh + integer(c_int), value :: nComps, numMPs + type(c_ptr), value :: array + end subroutine + + + + !--------------------------------------------------------------------------- !> @brief set the Mass MP array from a host array !> @param mpmesh(in/out) MPMesh object @@ -768,6 +844,27 @@ logical function polympo_push1P(mpMesh) & type(c_ptr), value :: mpMesh end function + subroutine polympo_push_ahead(mpMesh) & + bind(C, NAME='polympo_push_ahead_f') + use :: iso_c_binding + type(c_ptr), value :: mpMesh + end subroutine + + subroutine polympo_push_swap(mpMesh) & + bind(C, NAME='polympo_push_swap_f') + use :: iso_c_binding + type(c_ptr), value :: mpMesh + end subroutine + + subroutine polympo_push_swap_pos(mpMesh) & + bind(C, NAME='polympo_push_swap_pos_f') + use :: iso_c_binding + type(c_ptr), value :: mpMesh + end subroutine + + + + !--------------------------------------------------------------------------- !> @brief start the reconstruction of MP Mass to Mesh Vertices !> @param mpmesh(in/out) MPMesh object diff --git a/src/pmpo_materialPoints.cpp b/src/pmpo_materialPoints.cpp index aca7f2e..2b05f92 100644 --- a/src/pmpo_materialPoints.cpp +++ b/src/pmpo_materialPoints.cpp @@ -8,11 +8,13 @@ template pumipic::MemberTypeViews createInternalMemberViews(int numMPs, View mp2elm, View mpAppID){ auto mpInfo = ps::createMemberViews(numMPs); auto mpCurElmPos_m = ps::getMemberView(mpInfo); + auto mpTgtElmPos_m = ps::getMemberView(mpInfo); auto mpAppID_m = ps::getMemberView(mpInfo); auto mpStatus_m = ps::getMemberView(mpInfo); auto policy = Kokkos::RangePolicy(typename MemSpace::execution_space(), 0, numMPs); Kokkos::parallel_for("setMPinfo", policy, KOKKOS_LAMBDA(int i) { mpCurElmPos_m(i) = mp2elm(i); + mpTgtElmPos_m(i) = -1; mpStatus_m(i) = MP_ACTIVE; mpAppID_m(i) = mpAppID(i); }); @@ -113,6 +115,22 @@ void MaterialPoints::finishRebuild() { ps::destroyViews(rebuildFields.addedSlices_h); ps::destroyViews(addedSlices_d); rebuildFields.ongoing = false; + + //Debug + /* + int rank; + MPI_Comm_rank(mpi_comm, &rank); + if (rank==0) return; + auto curr_elm=getData(); + auto tgt_elm =getData(); + auto mpAppID = getData(); + auto testElm = PS_LAMBDA(const int& e, const int& mp, const bool& mask) { + if (mask) { + printf("R1: finishRebuild AppID %d Curr %d Tgt %d e %d \n", mpAppID(mp), curr_elm(mp), tgt_elm(mp), e); + } + }; + parallel_for(testElm, "curr_elm"); + */ } MPI_Comm MaterialPoints::getMPIComm() { From a9e74ab5300e295914e8e8ab12fc965a3a87ef82 Mon Sep 17 00:00:00 2001 From: Nath Date: Tue, 20 May 2025 14:26:36 -0400 Subject: [PATCH 14/17] Removed print annd some cleaning --- src/pmpo_MPMesh.cpp | 45 ++------------------------------------------- src/pmpo_c.cpp | 8 -------- 2 files changed, 2 insertions(+), 51 deletions(-) diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index f2fa0c6..df9f3bb 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -133,26 +133,18 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){ MPI_Comm comm = p_MPs->getMPIComm(); int comm_rank; MPI_Comm_rank(comm, &comm_rank); - Kokkos::parallel_for("countProcess", numElms, KOKKOS_LAMBDA(const int iElm){ - int pp_id=elm2Process(iElm); - }); - + //Since Mesh is static print pnly for 1 time step if(printVTPIndex==0) { printVTP_mesh(comm_rank); } - //assert(cudaDeviceSynchronize()==cudaSuccess); - //MPI_Barrier(MPI_COMM_WORLD); Vec3dView history("positionHistory",numMPs); Vec3dView resultLeft("positionResult",numMPs); Vec3dView resultRight("positionResult",numMPs); Vec3dView mpTgtPosArray("positionTarget",numMPs); Kokkos::View counter("counter",1); - assert(cudaDeviceSynchronize() == cudaSuccess); - //printf("Rank %d Foo4 Begin\n", comm_rank); - auto CVTElmCalc = PS_LAMBDA(const int& elm, const int& mp, const int&mask){ Vec3d MP(mpPositions(mp,0),mpPositions(mp,1),mpPositions(mp,2)); if(mask){ @@ -190,12 +182,7 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){ iElm = closestElm; } } - if(mpAppID(mp)==0 || mpAppID(mp)==191) - printf("Pos %.15e %.15e %.15e => %.15e %.15e %.15e\n", mpPositions(mp,0), mpPositions(mp,1), - mpPositions(mp,2), mpTgtPos(mp,0), mpTgtPos(mp,1), mpTgtPos(mp, 2)); - if(printVTPIndex>=0 && numMPs>0){ - //printf("Rank %d mp %d counter %d \n", comm_rank, mp, counter); double d1 = dx[0]; double d2 = dx[2]; double d3 = dx[3]; @@ -217,10 +204,6 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){ }; p_MPs->parallel_for(CVTElmCalc,"CVTTrackingElmCenterBasedCalc"); - assert(cudaDeviceSynchronize() == cudaSuccess); - //printf("Rank %d Foo4 End\n", comm_rank); - - if(printVTPIndex>=0 && numMPs>0){ Vec3dView::HostMirror h_history = Kokkos::create_mirror_view(history); Vec3dView::HostMirror h_resultLeft = Kokkos::create_mirror_view(resultLeft); @@ -351,27 +334,20 @@ bool getAnyIsMigrating(MaterialPoints* p_MPs, bool isMigrating) { } void MPMesh::push_ahead(){ - static int count0=0; - std::cout<<"Push_ahead"<<" "<computeRotLatLonIncr(); sphericalInterpolation(*this); //Push the MPs p_MPs->updateRotLatLonAndXYZ2Tgt(p_mesh->getSphereRadius()); - count0 ++; } bool MPMesh::push1P(){ //Given target location find the new element or the last element in a partioned mesh //and the process it belongs to so that migration can be checked - static int count_p=0; - CVTTrackingElmCenterBased(); - + CVTTrackingElmCenterBased(); //From the above two inputs find if any particle needs to be migrated bool anyIsMigrating = getAnyIsMigrating(p_MPs, p_MPs->check_migrate()); - count_p=count_p+1; return anyIsMigrating; - } void MPMesh::push_swap(){ @@ -379,7 +355,6 @@ void MPMesh::push_swap(){ p_MPs->updateMPElmID(); } - void MPMesh::push_swap_pos(){ //current becomes target, target becomes -1 //Making read for next push_ahead @@ -390,9 +365,6 @@ void MPMesh::push_swap_pos(){ void MPMesh::push(){ - static int count=0; - std::cout<<"Push"<<" "<computeRotLatLonIncr(); @@ -422,19 +394,6 @@ void MPMesh::push(){ } while (anyIsMigrating); - count ++; - - volatile int i = 0; - char hostname[256]; - gethostname(hostname, sizeof(hostname)); - MPI_Comm comm = p_MPs->getMPIComm(); - int comm_rank; - MPI_Comm_rank(comm, &comm_rank); - if(count==480){ - printf("Rank %d PID %d on %s ready for attach\n", comm_rank, getpid(), hostname); - fflush(stdout); - sleep(100); - } pumipic::RecordTime("PolyMPO_push", timer.seconds()); } diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index ce7b7ba..5b0cba8 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -210,16 +210,9 @@ void polympo_startRebuildMPs_f2(MPMesh_ptr p_mpmesh, auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; int offset = p_MPs->getElmIDoffset(); - int rank; - auto mpi_comm=p_MPs->getMPIComm(); - MPI_Comm_rank(mpi_comm, &rank); - - printf("Rank %d to add %d and delete %d ptcls \n", rank, nMPs_add, nMPs_delete); for (int k=0; k Date: Tue, 20 May 2025 17:35:12 -0400 Subject: [PATCH 15/17] Some more redudancy removed like trying to get new AppID from MPAS --- src/pmpo_c.cpp | 14 +----- src/pmpo_c.h | 1 - src/pmpo_fortran.f90 | 18 -------- src/pmpo_materialPoints.cpp | 88 ++----------------------------------- src/pmpo_materialPoints.hpp | 12 ++--- 5 files changed, 8 insertions(+), 125 deletions(-) diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index 5b0cba8..71473fd 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -249,17 +249,7 @@ void polympo_setAppIDFunc_f(MPMesh_ptr p_mpmesh, IntVoidFunc getNext, void* appI checkMPMeshValid(p_mpmesh); auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; polyMPO::IntFunc getNextAppID = [getNext, appIDs]() { return getNext(appIDs); }; - //p_MPs->setAppIDFunc(getNextAppID); -} - -//arg1 is blockPtr, arg2 is blockSize and arg3 is nCells -void polympo_setMPASAppIDFunc_f(MPMesh_ptr p_mpmesh, VoidVoidFunc getMPASAppID, void*arg1, - const int arg2, const int arg3) { - checkMPMeshValid(p_mpmesh); - auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; - std::function polyMPO_getMPASAppID = [getMPASAppID, arg1, arg2, arg3](int iCell) - { int iParticleNew=-1; getMPASAppID(arg1, arg2, arg3, iCell, iParticleNew); return iParticleNew;}; - p_MPs->setAppIDFunc(polyMPO_getMPASAppID); + p_MPs->setAppIDFunc(getNextAppID); } void polympo_getMPTgtElmID_f(MPMesh_ptr p_mpmesh, @@ -1263,7 +1253,7 @@ void polympo_setElmGlobal_f(MPMesh_ptr p_mpmesh, const int nCells, const int* ar PMT_ALWAYS_ASSERT(p_mesh->meshEditable()); Kokkos::View arrayHost("arrayHost", nCells); for (int i = 0; i < nCells; i++) { - arrayHost(i) = array[i] - 1; // Decrease each value by 1 + arrayHost(i) = array[i] - 1; // TODO right now elmID offset is set after MPs initialized } //check the size PMT_ALWAYS_ASSERT(nCells == p_mesh->getNumElements()); diff --git a/src/pmpo_c.h b/src/pmpo_c.h index cfbe628..ead18d6 100644 --- a/src/pmpo_c.h +++ b/src/pmpo_c.h @@ -28,7 +28,6 @@ int polympo_getMPCount_f(MPMesh_ptr p_mpmesh); void polympo_finishRebuildMPs_f(MPMesh_ptr p_mpmesh); void polympo_setAppIDFunc_f(MPMesh_ptr p_mpmesh, IntVoidFunc getNext, void* appIDs); -void polympo_setMPASAppIDFunc_f(MPMesh_ptr p_mpmesh, VoidVoidFunc getMPASAppID, void* arg1, const int arg2, const int arg3); void polympo_getMPTgtElmID_f(MPMesh_ptr p_mpmesh, const int numMPs, int* elmIDs); void polympo_getMPCurElmID_f(MPMesh_ptr p_mpmesh, const int numMPs, int* elmIDs); diff --git a/src/pmpo_fortran.f90 b/src/pmpo_fortran.f90 index b05ce49..694b6e0 100644 --- a/src/pmpo_fortran.f90 +++ b/src/pmpo_fortran.f90 @@ -125,24 +125,6 @@ subroutine polympo_setAppIDFunc(mpMesh, getNext, appIDs) & type(c_ptr), value :: appIDs end subroutine - !--------------------------------------------------------------------------- - !> @brief Stores pointer to appID data structure and a function to retrieve them used in migration - !> @param mpmesh(in/out) MPMesh object - !> @param getNext(in) Pointer to function that returns next App IDs - !> @param appIDs(in) Pointer to opaque data application data structure (that may contain all available app IDs) - !--------------------------------------------------------------------------- - subroutine polympo_setMPASAppIDFunc(mpMesh, getMPASAppID, & - arg1, arg2, arg3) & - bind(C, NAME='polympo_setMPASAppIDFunc_f') - use :: iso_c_binding - type(c_ptr), value :: mpMesh - type(c_funptr), value :: getMPASAppID - type(c_ptr), value :: arg1 - integer(c_int), intent(in), value :: arg2 - integer(c_int), intent(in), value :: arg3 - end subroutine - - !--------------------------------------------------------------------------- !> @brief get the current element ID MP array from a polympo array !> @param mpmesh(in/out) MPMesh object diff --git a/src/pmpo_materialPoints.cpp b/src/pmpo_materialPoints.cpp index 2b05f92..46fb41a 100644 --- a/src/pmpo_materialPoints.cpp +++ b/src/pmpo_materialPoints.cpp @@ -114,23 +114,7 @@ void MaterialPoints::finishRebuild() { updateMaxAppID(); ps::destroyViews(rebuildFields.addedSlices_h); ps::destroyViews(addedSlices_d); - rebuildFields.ongoing = false; - - //Debug - /* - int rank; - MPI_Comm_rank(mpi_comm, &rank); - if (rank==0) return; - auto curr_elm=getData(); - auto tgt_elm =getData(); - auto mpAppID = getData(); - auto testElm = PS_LAMBDA(const int& e, const int& mp, const bool& mask) { - if (mask) { - printf("R1: finishRebuild AppID %d Curr %d Tgt %d e %d \n", mpAppID(mp), curr_elm(mp), tgt_elm(mp), e); - } - }; - parallel_for(testElm, "curr_elm"); - */ + rebuildFields.ongoing = false; } MPI_Comm MaterialPoints::getMPIComm() { @@ -166,7 +150,6 @@ void MaterialPoints::migrate() { Kokkos::Timer timer; auto MPs2Elm = getData(); auto MPs2Proc = getData(); - auto mpAppID = getData(); IntView new_elem("new_elem", MPs->capacity()); IntView new_process("new_process", MPs->capacity()); @@ -177,76 +160,11 @@ void MaterialPoints::migrate() { if (mask) { new_elem(mp) = MPs2Elm(mp); new_process(mp) = MPs2Proc(mp); - if(rank!=new_process(mp)){ - printf("Particle %d in rank %d to be moved from %d %d \n", mpAppID(mp), rank, e, new_elem(mp) ); - mpAppID(mp)=-1; - } } }; parallel_for(setMigrationFields, "setMigrationFields"); MPs->migrate(new_elem, new_process); - //AS REBUILT, mpAppID needs to be be recalled - mpAppID = getData(); - //Count MPs that have -1 appID, so that we can count no of MPs received by a rank - Kokkos::View numReceivedMPs("numReceivedMPs", 1); - Kokkos::deep_copy(numReceivedMPs, 0); - auto countnewMPs = PS_LAMBDA(const int& e, const int& mp, const bool& mask) { - if(mask){ - if (mpAppID(mp) == -1) - Kokkos::atomic_add(&numReceivedMPs(0), 1); - } - }; - parallel_for(countnewMPs, "countReceivedPtcls"); - auto numReceivedMPs_host = Kokkos::create_mirror_view(numReceivedMPs); - Kokkos::deep_copy(numReceivedMPs_host, numReceivedMPs); - - //Another array that contains new element id of the the MPs that have migrated - //Array size is #elements received - Kokkos::View receivedMPs2Elm("ReceivedMPs2Elm", numReceivedMPs_host(0)); - Kokkos::View counter("counter", 1); - Kokkos::deep_copy(counter, 0); - auto set_new_elem = PS_LAMBDA(const int& e, const int& mp, const bool& mask) { - if(mask){ - if (mpAppID(mp) == -1){ - auto count_temp=Kokkos::atomic_fetch_add(&counter(0), 1); - receivedMPs2Elm(count_temp) = e; - } - } - }; - parallel_for(set_new_elem, "countReceivedPtcls"); - //Bring them to CPU so that elm_id can be passed and a new mpAppID can be found from CPU - auto receivedMPs2Elm_host = Kokkos::create_mirror_view(receivedMPs2Elm); - Kokkos::deep_copy(receivedMPs2Elm_host, receivedMPs2Elm); - auto counter_host = Kokkos::create_mirror_view(counter); - Kokkos::deep_copy(counter_host, counter); - assert(numReceivedMPs_host(0)==counter_host(0)); - if(counter_host(0)) - std::cout <<"Rank "< appIDs; - for(int i=0; i appIDs_host(appIDs.data(), appIDs.size()); - Kokkos::View appIDs_d("appIDsDevice", appIDs.size()); - Kokkos::deep_copy(appIDs_d, appIDs_host); - - //If the mpAppID is -1 assign a new mpAppID - Kokkos::deep_copy(counter, 0); - auto set_appID = PS_LAMBDA(const int& e, const int& mp, const bool& mask) { - if(mask){ - if (mpAppID(mp) == -1){ - auto count_temp=Kokkos::atomic_fetch_add(&counter(0), 1); - mpAppID(mp)=appIDs_d(count_temp)-1; - } - } - }; - parallel_for(set_appID, "setApplicationIDs"); - if (getOpMode() == polyMPO::MP_DEBUG) printf("Material point migration: %f\n", timer.seconds()); pumipic::RecordTime("PolyMPO_migrate", timer.seconds()); @@ -254,8 +172,8 @@ void MaterialPoints::migrate() { bool MaterialPoints::rebuildOngoing() { return rebuildFields.ongoing; } -void MaterialPoints::setAppIDFunc(IntIntFunc getAppIDIn) { getAppID = getAppIDIn; } +void MaterialPoints::setAppIDFunc(IntFunc getAppIDIn) { getAppID = getAppIDIn; } -int MaterialPoints::getNextAppID(int iElm) { return getAppID(iElm); } +int MaterialPoints::getNextAppID() { return getAppID(); } } diff --git a/src/pmpo_materialPoints.hpp b/src/pmpo_materialPoints.hpp index c2a037f..70c71dc 100644 --- a/src/pmpo_materialPoints.hpp +++ b/src/pmpo_materialPoints.hpp @@ -127,8 +127,7 @@ class MaterialPoints { bool isRotatedFlag = false; Operating_Mode operating_mode; RebuildHelper rebuildFields; - //IntFunc getAppID; - IntIntFunc getAppID; + IntFunc getAppID; MPI_Comm mpi_comm; public: @@ -157,20 +156,15 @@ class MaterialPoints { typename std::enable_if::type setRebuildMPSlice(mpSliceData mpSliceIn); - //void setAppIDFunc(IntFunc getAppIDIn); - void setAppIDFunc(IntIntFunc getAppIDIn); - - int getNextAppID(int iElm); + void setAppIDFunc(IntFunc getAppIDIn); + int getNextAppID(); void rebuild() { IntView tgtElm("tgtElm", MPs->capacity()); auto tgtMpElm = MPs->get(); - auto mpAppID = MPs->get(); auto setTgtElm = PS_LAMBDA(const int& e, const int& mp, const int& mask) { if(mask) { - auto app_id=mpAppID(mp); tgtElm(mp) = tgtMpElm(mp); - if(app_id==191) printf("Going from elm %d to %d \n", e, tgtMpElm(mp)); } }; ps::parallel_for(MPs, setTgtElm, "setTargetElement"); From 4b6930f709eb022e79052883b282dcaffa969cc8 Mon Sep 17 00:00:00 2001 From: Nath Date: Wed, 21 May 2025 12:12:30 -0400 Subject: [PATCH 16/17] All tests on remus pass --- src/pmpo_MPMesh.cpp | 29 ++++++++++------------------ src/pmpo_c.cpp | 2 -- src/pmpo_mesh.hpp | 3 +-- test/testFortranCreateRebuildMPs.f90 | 5 ++++- test/testFortranMPAdvection.f90 | 23 ++++++---------------- test/testFortranMPAppIDs.f90 | 15 +++++++++++--- test/testFortranMPReconstruction.f90 | 8 +++++--- test/testMPAppIDs.cpp | 2 +- 8 files changed, 39 insertions(+), 48 deletions(-) diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index df9f3bb..d525b04 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -128,22 +128,16 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){ auto MPs2Proc = p_MPs->getData(); auto elm2Process = p_mesh->getElm2Process(); auto elm2global = p_mesh->getElmGlobal(); - auto mpAppID = p_MPs->getData(); - 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 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::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, "\n \n \n \n \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(); // Tgt_XYZ becomes Cur_XYZ p_MPs->updateMPSlice(); // Tgt becomes Cur diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index 71473fd..4ac7275 100644 --- a/src/pmpo_c.cpp +++ b/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 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 diff --git a/src/pmpo_mesh.hpp b/src/pmpo_mesh.hpp index 0819da2..22e6ddd 100644 --- a/src/pmpo_mesh.hpp +++ b/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(); diff --git a/test/testFortranCreateRebuildMPs.f90 b/test/testFortranCreateRebuildMPs.f90 index 61ea24a..478708a 100644 --- a/test/testFortranCreateRebuildMPs.f90 +++ b/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) diff --git a/test/testFortranMPAdvection.f90 b/test/testFortranMPAdvection.f90 index abdedf7..292d806 100644 --- a/test/testFortranMPAdvection.f90 +++ b/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 diff --git a/test/testFortranMPAppIDs.f90 b/test/testFortranMPAppIDs.f90 index bcae9d7..bdcbc8e 100644 --- a/test/testFortranMPAppIDs.f90 +++ b/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 \ No newline at end of file +end program diff --git a/test/testFortranMPReconstruction.f90 b/test/testFortranMPReconstruction.f90 index 03a4165..4a9e910 100644 --- a/test/testFortranMPReconstruction.f90 +++ b/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 diff --git a/test/testMPAppIDs.cpp b/test/testMPAppIDs.cpp index 4d33658..c3064cc 100644 --- a/test/testMPAppIDs.cpp +++ b/test/testMPAppIDs.cpp @@ -67,4 +67,4 @@ void testAppIDPointer(MPMesh_ptr p_mpmesh) { PMT_ALWAYS_ASSERT(numAddedMPsAfter_h == numAddedMPs); } -#endif \ No newline at end of file +#endif From 00643ff10962e3f642072bdffa646ad73c6e3831 Mon Sep 17 00:00:00 2001 From: Nath Date: Wed, 21 May 2025 12:21:57 -0400 Subject: [PATCH 17/17] Removed print statement --- test/testFortranMPAppIDs.f90 | 1 - 1 file changed, 1 deletion(-) diff --git a/test/testFortranMPAppIDs.f90 b/test/testFortranMPAppIDs.f90 index bdcbc8e..a7cf7f2 100644 --- a/test/testFortranMPAppIDs.f90 +++ b/test/testFortranMPAppIDs.f90 @@ -68,7 +68,6 @@ subroutine testAppIDPointer(mpMesh) & !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));