From 7ba0b24475b55d1ba46d99dd71005796b05b8aff Mon Sep 17 00:00:00 2001 From: Shahrear Jahan Santho Date: Thu, 24 Sep 2026 18:31:29 -0700 Subject: [PATCH 01/12] Add CUDA-aware MPI halo exchange --- CMakeLists.txt | 2 +- src/CMakeLists.txt | 1 + src/pmpo_MPMesh.cpp | 2 +- src/pmpo_MPMesh.hpp | 820 +++++++++++++++++++++++++++++++---- src/pmpo_MPMesh_assembly.hpp | 6 +- src/pmpo_c.cpp | 1 + 6 files changed, 736 insertions(+), 96 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 080a4186..3f627daa 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -50,4 +50,4 @@ if(IS_TESTING) add_subdirectory (test) endif() -bob_end_package() \ No newline at end of file +bob_end_package() diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 76059965..d0fb8e74 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -21,6 +21,7 @@ set(SOURCES ) add_library(polyMPO-core ${SOURCES}) +target_compile_definitions(polyMPO-core PUBLIC CUDA_AWARE_MPI) set_property(TARGET polyMPO-core PROPERTY CXX_STANDARD "17") set_property(TARGET polyMPO-core PROPERTY CXX_STANDARD_REQUIRED ON) set_property(TARGET polyMPO-core PROPERTY CXX_EXTENSIONS OFF) diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index 3801dcd9..837b896e 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -171,7 +171,7 @@ void MPMesh::calculateStressDivergence(){ timer.reset(); if(numProcsTot>1){ //Takes contribution of halo vertices and adds it in owner procs - communicate_and_take_halo_contributions1(stress_divUV, numVertices, 2, 0, 0); + communicate_and_take_halo_contributions1_improved(stress_divUV, numVertices, 2, 0, 0); //Transfer the correct values at owned vertices to halo vertices //communicate_and_take_halo_contributions(stress_divUV, numVertices, 2, 1, 1); } diff --git a/src/pmpo_MPMesh.hpp b/src/pmpo_MPMesh.hpp index d3aac1b7..dfd922e5 100644 --- a/src/pmpo_MPMesh.hpp +++ b/src/pmpo_MPMesh.hpp @@ -4,6 +4,9 @@ #include "pmpo_utils.hpp" #include "pmpo_mesh.hpp" #include "pmpo_materialPoints.hpp" +#include +#include +#include namespace polyMPO{ @@ -28,22 +31,27 @@ class MPMesh{ int numOwnersTot, numHalosTot; std::vector numOwnersOnOtherProcs; std::vector numHalosOnOtherProcs; - std::vectorhaloOwnerProcs; + std::vector haloOwnerProcs; std::vector> haloOwnerLocalIDs; std::vector> ownerOwnerLocalIDs; std::vector> ownerHaloLocalIDs; void startCommunication(); - void communicate_and_take_halo_contributions(const Kokkos::View& meshField, int nEntities, int numEntries, int mode, int op); + void communicate_and_take_halo_contributions( + const Kokkos::View& meshField, + int nEntities, + int numEntries, + int mode, + int op); - //Now Kokkos views are made 1D + // Original CPU-staging function template void communicate_and_take_halo_contributions1( const ViewType& meshField, int nEntities, int numEntries, - int mode , + int mode, int op){ int self; @@ -51,95 +59,155 @@ class MPMesh{ MPI_Comm_rank(comm, &self); Kokkos::Timer timer; - auto reconVals_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), meshField); + auto reconVals_host = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), meshField); + pumipic::RecordTime("SD: GPU-CPU copy-" + std::to_string(self), timer.seconds()); timer.reset(); - std::vector> recvIDVec; + std::vector> recvIDVec; std::vector> recvDataVec; + pumipic::RecordTime("SD: Recv Vec Allocation-" + std::to_string(self), timer.seconds()); timer.reset(); - //communicateFields1(fieldData1, nEntities, numEntries, mode, recvIDVec, recvDataVec); - communicateFields1(reconVals_host, nEntities, numEntries, mode, recvIDVec, recvDataVec); + + communicateFields1( + reconVals_host, + nEntities, + numEntries, + mode, + recvIDVec, + recvDataVec); + pumipic::RecordTime("SD: IP Comm-" + std::to_string(self), timer.seconds()); timer.reset(); + int numProcsTot = recvIDVec.size(); - //Flatten IDs + int totalSize = 0; - std::vector offsets(numProcsTot, 0); - for(int i=0; i offsets(numProcsTot, 0); + + for(int i = 0; i < numProcsTot; i++){ offsets[i] = totalSize; totalSize += recvIDVec[i].size(); } - std::vector flatIDVec(totalSize, 0); - for(int i=0; i recvIDGPU("recvIDGPU", totalSize); + auto hostView = + Kokkos::View("recvIDCPU", totalSize); + + for(int i = 0; i < numProcsTot; i++){ + std::copy( + recvIDVec[i].begin(), + recvIDVec[i].end(), + hostView.data() + offsets[i]); } + pumipic::RecordTime("SD: Flatten IDs-" + std::to_string(self), timer.seconds()); timer.reset(); - Kokkos::View recvIDGPU("recvIDGPU", totalSize); - auto hostView = Kokkos::View("recvIDCPU", totalSize); - std::copy(flatIDVec.begin(), flatIDVec.end(), hostView.data()); + Kokkos::deep_copy(recvIDGPU, hostView); Kokkos::fence(); + pumipic::RecordTime("SD: Copy CPU-GPU-" + std::to_string(self), timer.seconds()); - //Flatten Data timer.reset(); - int totalSize_data=0; + + int totalSize_data = 0; std::vector offsets_data(numProcsTot, 0); - for(int i=0; i flatDataVec(totalSize_data, 0); - for(int i=0; i recvDataGPU("recvDataGPU", totalSize_data); + auto hostView_data = + Kokkos::View("recvDataCPU", totalSize_data); + + for(int i = 0; i < numProcsTot; i++){ + std::copy( + recvDataVec[i].begin(), + recvDataVec[i].end(), + hostView_data.data() + offsets_data[i]); } + pumipic::RecordTime("SD: Flatten Data-" + std::to_string(self), timer.seconds()); timer.reset(); - Kokkos::View recvDataGPU("recvDataGPU", totalSize_data); - auto hostView_data= Kokkos::View("recvDataCPU", totalSize_data); - std::copy(flatDataVec.begin(), flatDataVec.end(), hostView_data.data()); + Kokkos::deep_copy(recvDataGPU, hostView_data); Kokkos::fence(); - assert(totalSize_data == totalSize*numEntries); - for (int i=0; i>& fieldData, const int numEntities, const int numEntries, int mode, - std::vector>& recvIDVec, std::vector>& recvDataVec); + void communicateFields( + const std::vector>& fieldData, + const int numEntities, + const int numEntries, + int mode, + std::vector>& recvIDVec, + std::vector>& recvDataVec); template void communicateFields1( - const ViewType& fieldData, - const int numEntities, const int numEntries, int mode, + const ViewType& fieldData, + const int numEntities, + const int numEntries, + int mode, std::vector>& recvIDVec, std::vector>& recvDataVec){ int self, numProcsTot; + MPI_Comm comm = p_MPs->getMPIComm(); + MPI_Comm_rank(comm, &self); MPI_Comm_size(comm, &numProcsTot); @@ -151,98 +219,191 @@ class MPMesh{ recvDataVec.resize(numProcsTot); for(int i = 0; i < numProcsTot; i++){ - if(i==self) continue; + if(i == self) continue; + + int numToSend = 0; + int numToRecv = 0; - int numToSend = 0, numToRecv = 0; - if(mode == 0) { - //gather (halos send to owners) + if(mode == 0){ numToSend = numOwnersOnOtherProcs[i]; numToRecv = numHalosOnOtherProcs[i]; } - else{ - //scatter (owners send to halos) + else{ numToSend = numHalosOnOtherProcs[i]; numToRecv = numOwnersOnOtherProcs[i]; } if(numToSend > 0){ - sendDataVec[i].reserve(numToSend*numEntries); + sendDataVec[i].reserve(numToSend * numEntries); } + if(numToRecv > 0){ - recvDataVec[i].resize(numToRecv*numEntries); + recvDataVec[i].resize(numToRecv * numEntries); recvIDVec[i].resize(numToRecv); } } if(mode == 0){ - // Halos sends to owners - for (int iEnt = 0; iEnt < numHalosTot; iEnt++){ + for(int iEnt = 0; iEnt < numHalosTot; iEnt++){ auto ownerProc = haloOwnerProcs[iEnt]; - for (int iDouble = 0; iDouble < numEntries; iDouble++) - sendDataVec[ownerProc].push_back(fieldData(numOwnersTot+iEnt, iDouble)); + + for(int iDouble = 0; iDouble < numEntries; iDouble++){ + sendDataVec[ownerProc].push_back( + fieldData(numOwnersTot + iEnt, iDouble)); + } } } else if(mode == 1){ - // Owner sends to halos - for (size_t iProc=0; iProc requests; - requests.reserve(4*numProcsTot); + requests.reserve(4 * numProcsTot); + for(int proc = 0; proc < numProcsTot; proc++){ - if(proc == self) continue; + if(proc == self) continue; + if(mode == 0 && numHalosOnOtherProcs[proc]){ - assert(recvIDVec[proc].size() == (size_t)numHalosOnOtherProcs[proc]); - assert(recvDataVec[proc].size() == recvIDVec[proc].size() * (size_t)numEntries); - MPI_Request req3, req4; - MPI_Irecv(recvIDVec[proc].data(), recvIDVec[proc].size(), MPI_INT, proc, 1, comm, &req3); - MPI_Irecv(recvDataVec[proc].data(), recvDataVec[proc].size(), MPI_DOUBLE, proc, 2, comm, &req4); + assert(recvIDVec[proc].size() == + static_cast(numHalosOnOtherProcs[proc])); + + assert(recvDataVec[proc].size() == + recvIDVec[proc].size() * static_cast(numEntries)); + + MPI_Request req3; + MPI_Request req4; + + MPI_Irecv( + recvIDVec[proc].data(), + recvIDVec[proc].size(), + MPI_INT, + proc, + 1, + comm, + &req3); + + MPI_Irecv( + recvDataVec[proc].data(), + recvDataVec[proc].size(), + MPI_DOUBLE, + proc, + 2, + comm, + &req4); + requests.push_back(req3); requests.push_back(req4); } - if(mode == 0 && numOwnersOnOtherProcs[proc]) { - assert(haloOwnerLocalIDs[proc].size() == (size_t)numOwnersOnOtherProcs[proc]); - assert(sendDataVec[proc].size() == haloOwnerLocalIDs[proc].size() * (size_t)numEntries); - MPI_Request req1, req2; - MPI_Isend(haloOwnerLocalIDs[proc].data(), haloOwnerLocalIDs[proc].size(), MPI_INT, proc, 1, comm, &req1); - MPI_Isend(sendDataVec[proc].data(), sendDataVec[proc].size(), MPI_DOUBLE, proc, 2, comm, &req2); + + if(mode == 0 && numOwnersOnOtherProcs[proc]){ + assert(haloOwnerLocalIDs[proc].size() == + static_cast(numOwnersOnOtherProcs[proc])); + + assert(sendDataVec[proc].size() == + haloOwnerLocalIDs[proc].size() * static_cast(numEntries)); + + MPI_Request req1; + MPI_Request req2; + + MPI_Isend( + haloOwnerLocalIDs[proc].data(), + haloOwnerLocalIDs[proc].size(), + MPI_INT, + proc, + 1, + comm, + &req1); + + MPI_Isend( + sendDataVec[proc].data(), + sendDataVec[proc].size(), + MPI_DOUBLE, + proc, + 2, + comm, + &req2); + requests.push_back(req1); requests.push_back(req2); } if(mode == 1 && numOwnersOnOtherProcs[proc]){ - MPI_Request req3, req4; - MPI_Irecv(recvIDVec[proc].data(), recvIDVec[proc].size(), MPI_INT, proc, 1, comm, &req3); - MPI_Irecv(recvDataVec[proc].data(), recvDataVec[proc].size(), MPI_DOUBLE, proc, 2, comm, &req4); + MPI_Request req3; + MPI_Request req4; + + MPI_Irecv( + recvIDVec[proc].data(), + recvIDVec[proc].size(), + MPI_INT, + proc, + 1, + comm, + &req3); + + MPI_Irecv( + recvDataVec[proc].data(), + recvDataVec[proc].size(), + MPI_DOUBLE, + proc, + 2, + comm, + &req4); + requests.push_back(req3); requests.push_back(req4); } - if(mode == 1 && numHalosOnOtherProcs[proc]) { - MPI_Request req1, req2; - MPI_Isend(ownerHaloLocalIDs[proc].data(), ownerHaloLocalIDs[proc].size(), MPI_INT, proc, 1, comm, &req1); - MPI_Isend(sendDataVec[proc].data(), sendDataVec[proc].size(), MPI_DOUBLE, proc, 2, comm, &req2); + + if(mode == 1 && numHalosOnOtherProcs[proc]){ + MPI_Request req1; + MPI_Request req2; + + MPI_Isend( + ownerHaloLocalIDs[proc].data(), + ownerHaloLocalIDs[proc].size(), + MPI_INT, + proc, + 1, + comm, + &req1); + + MPI_Isend( + sendDataVec[proc].data(), + sendDataVec[proc].size(), + MPI_DOUBLE, + proc, + 2, + comm, + &req2); + requests.push_back(req1); requests.push_back(req2); } } + MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE); } + MPMesh(Mesh* inMesh, MaterialPoints* inMPs): - p_mesh(inMesh), p_MPs(inMPs) { + p_mesh(inMesh), + p_MPs(inMPs) { }; - ~MPMesh() { + + ~MPMesh(){ delete p_mesh; delete p_MPs; } - //MP advection and tracking + + // MP advection and tracking void CVTTrackingEdgeCenterBased(Vec2dView dx); void CVTTrackingElmCenterBased(const int printVTPIndex = -1); void T2LTracking(Vec2dView dx); @@ -252,44 +413,521 @@ class MPMesh{ void push_swap_pos(); void push(); - //Used before advection to interpolate fields from mesh to MPs - //And also before reconstruction + + // Used before advection to interpolate fields from mesh to MPs + // And also before reconstruction void calcBasis(); - //Reconstruction + + // Reconstruction DoubleView assemblyV0(); + template void assemblyVtx0(); + template void assemblyElm0(); + template void assemblyVtx1(); + void reconstruct_coeff_full(); - void invertMatrix(const Kokkos::View& vtxMatrices, const double& radius); + + void invertMatrix( + const Kokkos::View& vtxMatrices, + const double& radius); + Kokkos::View precomputedVtxCoeffs_new; - Kokkos::View nearAnEdge; + Kokkos::View nearAnEdge; Kokkos::View vtxMatrixMass; - //Not used currently - std::map> reconstructSlice = std::map>(); + + // Not used currently + std::map> reconstructSlice = + std::map>(); + template DoubleView wtScaAssembly(); + template Vec2dView wtVec2Assembly(); + template - void assembly(int order, MeshFieldType type, bool basisWeightFlag, bool massWeightFlag); + void assembly( + int order, + MeshFieldType type, + bool basisWeightFlag, + bool massWeightFlag); + template - void setReconstructSlice(int order, MeshFieldType type); + void setReconstructSlice( + int order, + MeshFieldType type); + void reconstructSlices(); void printVTP_mesh(int printVTPIndex); - void writeMPTrackingVTP(int printVTPIndex, int numMPs, const Vec3dView& history, const Vec3dView& resultLeft, - const Vec3dView& resultRight, const Vec3dView& mpTgtPosArray); + void writeMPTrackingVTP( + int printVTPIndex, + int numMPs, + const Vec3dView& history, + const Vec3dView& resultLeft, + const Vec3dView& resultRight, + const Vec3dView& mpTgtPosArray); void calculateStrain(); void calculateStress(const int constitutive_relation); void calculateStressDivergence(); + + + +#ifdef CUDA_AWARE_MPI + + // Cached CUDA-aware MPI communication metadata and per-neighbor GPU buffers. + // Important change from the previous version: + // MPI is always given the base pointer of a Kokkos allocation, not + // "base pointer + offset". This avoids Cray MPICH/GTL CUDA IPC problems. + bool cudaAwareMPICacheValid = true; + bool cudaAwareMPIDisabled = false; + bool cudaAwareMPIEnvChecked = false; + bool cudaAwareMPIForceCPU = false; + bool cudaAwareMPILogged = false; + + struct CudaAwareMPIFieldCache{ + bool valid = false; + int cachedNumProcs = -1; + + std::vector sendCounts; + std::vector recvCounts; + + std::vector> sendEntityGPUPerProc; + std::vector> recvIDGPUPerProc; + + std::vector> sendDataGPUPerProc; + std::vector> recvDataGPUPerProc; + }; + + std::map, CudaAwareMPIFieldCache> cudaAwareMPICaches; + + bool cudaAwareMPIForceDisabled(){ + if(!cudaAwareMPIEnvChecked){ + const char* value = std::getenv("POLYMPO_DISABLE_CUDA_AWARE_MPI"); + cudaAwareMPIForceCPU = + value != nullptr && value[0] != '\0' && value[0] != '0'; + cudaAwareMPIEnvChecked = true; + } + + return cudaAwareMPIForceCPU; + } + + // Fully CUDA-aware MPI version: + // Field data is sent/received using GPU pointers. Receive IDs are cached + // once from the fixed halo/owner mapping and are not sent every call. + // + // Important: + // This function caches communication metadata and GPU buffers per + // (mode, numEntries). If the communication pattern changes, clear + // cudaAwareMPICaches before the next call. + template + void communicate_and_take_halo_contributions1_improved( + const ViewType& meshField, + int nEntities, + int numEntries, + int mode, + int op){ + + int self, numProcsTot; + + MPI_Comm comm = p_MPs->getMPIComm(); + + MPI_Comm_rank(comm, &self); + MPI_Comm_size(comm, &numProcsTot); + + assert(mode == 0 || mode == 1); + assert(op == 0 || op == 1); + assert(nEntities == numOwnersTot + numHalosTot); + + if(cudaAwareMPIDisabled || cudaAwareMPIForceDisabled()){ + communicate_and_take_halo_contributions1( + meshField, + nEntities, + numEntries, + mode, + op); + return; + } + +#ifdef POLYMPO_VERBOSE_MPI + if(self == 0 && !cudaAwareMPILogged){ + std::cout + << "[CUDA_AWARE_MPI] Using per-proc cached full GPU-aware MPI path in communicate_and_take_halo_contributions1_improved()" + << "\n"; + cudaAwareMPILogged = true; + } +#endif + + Kokkos::Timer timer; + + if(!cudaAwareMPICacheValid){ + cudaAwareMPICaches.clear(); + cudaAwareMPICacheValid = true; + } + + auto& cudaAwareCache = + cudaAwareMPICaches[std::make_pair(mode, numEntries)]; + + const bool needRebuild = + (!cudaAwareCache.valid) || + (cudaAwareCache.cachedNumProcs != numProcsTot); + + if(needRebuild){ + + cudaAwareCache.cachedNumProcs = numProcsTot; + + cudaAwareCache.sendCounts.assign(numProcsTot, 0); + cudaAwareCache.recvCounts.assign(numProcsTot, 0); + + cudaAwareCache.sendEntityGPUPerProc.clear(); + cudaAwareCache.recvIDGPUPerProc.clear(); + cudaAwareCache.sendDataGPUPerProc.clear(); + cudaAwareCache.recvDataGPUPerProc.clear(); + + cudaAwareCache.sendEntityGPUPerProc.resize(numProcsTot); + cudaAwareCache.recvIDGPUPerProc.resize(numProcsTot); + cudaAwareCache.sendDataGPUPerProc.resize(numProcsTot); + cudaAwareCache.recvDataGPUPerProc.resize(numProcsTot); + + for(int proc = 0; proc < numProcsTot; proc++){ + if(proc == self) continue; + + if(mode == 0){ + cudaAwareCache.sendCounts[proc] = numOwnersOnOtherProcs[proc]; + cudaAwareCache.recvCounts[proc] = numHalosOnOtherProcs[proc]; + } + else{ + cudaAwareCache.sendCounts[proc] = numHalosOnOtherProcs[proc]; + cudaAwareCache.recvCounts[proc] = numOwnersOnOtherProcs[proc]; + } + } + + for(int proc = 0; proc < numProcsTot; proc++){ + if(proc == self) continue; + + const int sendCount = cudaAwareCache.sendCounts[proc]; + const int recvCount = cudaAwareCache.recvCounts[proc]; + + if(sendCount > 0){ + cudaAwareCache.sendEntityGPUPerProc[proc] = + Kokkos::View( + "cudaAwareMPISendEntityGPUPerProc", + sendCount); + + cudaAwareCache.sendDataGPUPerProc[proc] = + Kokkos::View( + "cudaAwareMPISendDataGPUPerProc", + sendCount * numEntries); + + auto sendEntityCPU = + Kokkos::View( + "sendEntityCPU", + sendCount); + + if(mode == 0){ + assert(haloOwnerLocalIDs[proc].size() == + static_cast(sendCount)); + + int localIndex = 0; + + for(int iEnt = 0; iEnt < numHalosTot; iEnt++){ + int ownerProc = haloOwnerProcs[iEnt]; + + if(ownerProc != proc) continue; + + assert(localIndex < sendCount); + + sendEntityCPU(localIndex) = numOwnersTot + iEnt; + + localIndex++; + } + + assert(localIndex == sendCount); + } + else{ + assert(ownerOwnerLocalIDs[proc].size() == + static_cast(sendCount)); + + for(int i = 0; i < sendCount; i++){ + sendEntityCPU(i) = ownerOwnerLocalIDs[proc][i]; + } + } + + Kokkos::deep_copy( + cudaAwareCache.sendEntityGPUPerProc[proc], + sendEntityCPU); + } + + if(recvCount > 0){ + cudaAwareCache.recvIDGPUPerProc[proc] = + Kokkos::View( + "cudaAwareMPIRecvIDGPUPerProc", + recvCount); + + cudaAwareCache.recvDataGPUPerProc[proc] = + Kokkos::View( + "cudaAwareMPIRecvDataGPUPerProc", + recvCount * numEntries); + + auto recvIDCPU = + Kokkos::View( + "recvIDCPU", + recvCount); + + if(mode == 0){ + assert(ownerOwnerLocalIDs[proc].size() == + static_cast(recvCount)); + + for(int i = 0; i < recvCount; i++){ + recvIDCPU(i) = ownerOwnerLocalIDs[proc][i]; + } + } + else{ + int localIndex = 0; + + for(int iEnt = 0; iEnt < numHalosTot; iEnt++){ + if(haloOwnerProcs[iEnt] != proc) continue; + + assert(localIndex < recvCount); + + recvIDCPU(localIndex) = numOwnersTot + iEnt; + + localIndex++; + } + + assert(localIndex == recvCount); + } + + Kokkos::deep_copy( + cudaAwareCache.recvIDGPUPerProc[proc], + recvIDCPU); + } + } + + Kokkos::fence(); + + cudaAwareCache.valid = true; + + pumipic::RecordTime( + "SD: CUDA-aware MPI Cache Build m" + std::to_string(mode) + + " e" + std::to_string(numEntries) + "-" + std::to_string(self), + timer.seconds()); + + timer.reset(); + } + + for(int proc = 0; proc < numProcsTot; proc++){ + if(proc == self) continue; + if(cudaAwareCache.sendCounts[proc] <= 0) continue; + + auto sendEntityGPU = cudaAwareCache.sendEntityGPUPerProc[proc]; + auto sendDataGPU = cudaAwareCache.sendDataGPUPerProc[proc]; + int sendCount = cudaAwareCache.sendCounts[proc]; + + Kokkos::parallel_for( + "pack cached cuda-aware mpi send buffer per proc", + sendCount, + KOKKOS_LAMBDA(const int i){ + int entity = sendEntityGPU(i); + + for(int k = 0; k < numEntries; k++){ + sendDataGPU(i * numEntries + k) = + meshField(entity, k); + } + }); + } + + Kokkos::fence(); + + pumipic::RecordTime( + "SD: CUDA-aware MPI Pack m" + std::to_string(mode) + + " e" + std::to_string(numEntries) + "-" + std::to_string(self), + timer.seconds()); + + timer.reset(); + + Kokkos::Timer mpiTotalTimer; + std::vector requests; + requests.reserve(2 * numProcsTot); + int mpiError = MPI_SUCCESS; + + for(int proc = 0; proc < numProcsTot; proc++){ + if(proc == self) continue; + + if(cudaAwareCache.recvCounts[proc] > 0){ + MPI_Request reqData; + + mpiError = MPI_Irecv( + cudaAwareCache.recvDataGPUPerProc[proc].data(), + cudaAwareCache.recvCounts[proc] * numEntries, + MPI_DOUBLE, + proc, + 2, + comm, + &reqData); + if(mpiError != MPI_SUCCESS) break; + requests.push_back(reqData); + } + + if(cudaAwareCache.sendCounts[proc] > 0){ + MPI_Request reqData; + + mpiError = MPI_Isend( + cudaAwareCache.sendDataGPUPerProc[proc].data(), + cudaAwareCache.sendCounts[proc] * numEntries, + MPI_DOUBLE, + proc, + 2, + comm, + &reqData); + if(mpiError != MPI_SUCCESS) break; + requests.push_back(reqData); + } + } + + pumipic::RecordTime( + "SD: CUDA-aware MPI Post m" + std::to_string(mode) + + " e" + std::to_string(numEntries) + "-" + std::to_string(self), + timer.seconds()); + + timer.reset(); + + if(mpiError == MPI_SUCCESS && !requests.empty()){ + mpiError = MPI_Waitall( + static_cast(requests.size()), + requests.data(), + MPI_STATUSES_IGNORE); + } + + pumipic::RecordTime( + "SD: CUDA-aware MPI Wait m" + std::to_string(mode) + + " e" + std::to_string(numEntries) + "-" + std::to_string(self), + timer.seconds()); + + if(mpiError != MPI_SUCCESS){ + cudaAwareMPIDisabled = true; + + if(self == 0){ + std::cout + << "[CUDA_AWARE_MPI] Device-pointer MPI failed." + << std::endl; + } + + if(requests.empty()){ + if(self == 0){ + std::cout + << "[CUDA_AWARE_MPI] Falling back to CPU-staged communication." + << std::endl; + } + + communicate_and_take_halo_contributions1( + meshField, + nEntities, + numEntries, + mode, + op); + return; + } + + if(self == 0){ + std::cout + << "[CUDA_AWARE_MPI] Failure happened after MPI requests were posted. " + << "Set POLYMPO_DISABLE_CUDA_AWARE_MPI=1 before running to force the CPU-staged path." + << std::endl; + } + + MPI_Abort(comm, mpiError); + return; + } + + pumipic::RecordTime( + "SD: CUDA-aware MPI Comm m" + std::to_string(mode) + + " e" + std::to_string(numEntries) + "-" + std::to_string(self), + mpiTotalTimer.seconds()); + + timer.reset(); + + for(int proc = 0; proc < numProcsTot; proc++){ + if(proc == self) continue; + if(cudaAwareCache.recvCounts[proc] <= 0) continue; + + auto recvIDGPU = cudaAwareCache.recvIDGPUPerProc[proc]; + auto recvDataGPU = cudaAwareCache.recvDataGPUPerProc[proc]; + int recvCount = cudaAwareCache.recvCounts[proc]; + + if(op == 0){ + Kokkos::parallel_for( + "halo add cached cuda-aware mpi per proc", + recvCount, + KOKKOS_LAMBDA(const int i){ + const int vertex = recvIDGPU(i); + + for(int k = 0; k < numEntries; k++){ +#ifdef POLYMPO_ASSUME_UNIQUE_HALO_CONTRIBS + meshField(vertex, k) += + recvDataGPU(i * numEntries + k); +#else + Kokkos::atomic_add( + &meshField(vertex, k), + recvDataGPU(i * numEntries + k)); +#endif + } + }); + } + else{ + Kokkos::parallel_for( + "halo assign cached cuda-aware mpi per proc", + recvCount, + KOKKOS_LAMBDA(const int i){ + const int vertex = recvIDGPU(i); + + for(int k = 0; k < numEntries; k++){ + meshField(vertex, k) = + recvDataGPU(i * numEntries + k); + } + }); + } + } + + Kokkos::fence(); + + pumipic::RecordTime( + "SD: CUDA-aware MPI Contribution m" + std::to_string(mode) + + " e" + std::to_string(numEntries) + "-" + std::to_string(self), + timer.seconds()); + } + +#else + + // Fallback path: + // if CUDA_AWARE_MPI is not defined, use the original GPU-CPU staging function. + template + void communicate_and_take_halo_contributions1_improved( + const ViewType& meshField, + int nEntities, + int numEntries, + int mode, + int op){ + + communicate_and_take_halo_contributions1( + meshField, + nEntities, + numEntries, + mode, + op); + } + +#endif + }; }//namespace polyMPO end diff --git a/src/pmpo_MPMesh_assembly.hpp b/src/pmpo_MPMesh_assembly.hpp index 9f184e05..7705c8f7 100644 --- a/src/pmpo_MPMesh_assembly.hpp +++ b/src/pmpo_MPMesh_assembly.hpp @@ -161,10 +161,10 @@ void MPMesh::reconstruct_coeff_full(){ int mode = 0; int op = 0; if (numProcsTot >1){ - communicate_and_take_halo_contributions1(vtxMatrices, numVertices, numEntriesMatrix, mode, op); + communicate_and_take_halo_contributions1_improved(vtxMatrices, numVertices, numEntriesMatrix, mode, op); mode=1; op=1; - communicate_and_take_halo_contributions1(vtxMatrices, numVertices, numEntriesMatrix, mode, op); + communicate_and_take_halo_contributions1_improved(vtxMatrices, numVertices, numEntriesMatrix, mode, op); } pumipic::RecordTime("Communicate Matrix Values" + std::to_string(self), timer.seconds()); @@ -380,7 +380,7 @@ void MPMesh::assemblyVtx1(){ timer.reset(); if(numProcsTot>1){ - communicate_and_take_halo_contributions1(meshField, numVertices, numEntries, 0, 0); + communicate_and_take_halo_contributions1_improved(meshField, numVertices, numEntries, 0, 0); } pumipic::RecordTime("Communicate Field Values" + std::to_string(self), timer.seconds()); } diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index 53c330e1..9c88f672 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -1609,6 +1609,7 @@ void polympo_set_halo_vel_from_owner_f(MPMesh_ptr p_mpmesh){ auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; int numVertices = p_mesh->getNumVertices(); auto vtxFieldVel = p_mesh->getMeshField(); + mpMesh->communicate_and_take_halo_contributions1(vtxFieldVel, numVertices, 2, 1, 1); } From ce6cfc966a540c6d67ca1d1ebde8130a534671c0 Mon Sep 17 00:00:00 2001 From: Shahrear Jahan Santho Date: Thu, 24 Sep 2026 18:03:12 -0700 Subject: [PATCH 02/12] Add GPU-aware halo exchange, unit test and cleanup --- src/CMakeLists.txt | 20 +- src/pmpo_MPMesh.cpp | 9 +- src/pmpo_MPMesh.hpp | 916 +++++++++++------------------- src/pmpo_MPMesh_assembly.hpp | 46 +- src/pmpo_c.cpp | 123 ++-- test/CMakeLists.txt | 2 + test/testGPUAwareHaloExchange.cpp | 284 +++++++++ 7 files changed, 727 insertions(+), 673 deletions(-) create mode 100644 test/testGPUAwareHaloExchange.cpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index d0fb8e74..32026a3e 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -21,7 +21,25 @@ set(SOURCES ) add_library(polyMPO-core ${SOURCES}) -target_compile_definitions(polyMPO-core PUBLIC CUDA_AWARE_MPI) + +if(DEFINED ENV{PE_MPICH_GTL_DIR_nvidia80} + AND DEFINED ENV{PE_MPICH_GTL_LIBS_nvidia80}) + + target_link_options(polyMPO-core PUBLIC + "$ENV{PE_MPICH_GTL_DIR_nvidia80}" + "$ENV{PE_MPICH_GTL_LIBS_nvidia80}" + ) + + message(STATUS + "polyMPO: enabling Cray MPICH CUDA GTL linkage") +else() + message(FATAL_ERROR + "Cray MPICH CUDA GTL environment variables are unavailable. " + "Load cudatoolkit and craype-accel-nvidia80, or set " + "CRAY_ACCEL_TARGET=nvidia80 before configuring.") +endif() + +target_compile_definitions(polyMPO-core PUBLIC GPU_AWARE_MPI) set_property(TARGET polyMPO-core PROPERTY CXX_STANDARD "17") set_property(TARGET polyMPO-core PROPERTY CXX_STANDARD_REQUIRED ON) set_property(TARGET polyMPO-core PROPERTY CXX_EXTENSIONS OFF) diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index 837b896e..fbb16485 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -89,7 +89,6 @@ void MPMesh::calculateStress(const int constitutive_relation){ void MPMesh::calculateStressDivergence(){ - Kokkos::Timer timer; int self, numProcsTot; MPI_Comm comm = p_MPs->getMPIComm(); MPI_Comm_rank(comm, &self); @@ -166,17 +165,13 @@ void MPMesh::calculateStressDivergence(){ }; p_MPs->parallel_for(stress_div, " stress_div_assembly"); Kokkos::fence(); - pumipic::RecordTime("Stress_Divergence_Reconstruction" + std::to_string(self), timer.seconds()); - - timer.reset(); - if(numProcsTot>1){ + if(numProcsTot>1){ //Takes contribution of halo vertices and adds it in owner procs - communicate_and_take_halo_contributions1_improved(stress_divUV, numVertices, 2, 0, 0); + communicate_and_take_halo_contributions_gpu_aware(stress_divUV, numVertices, 2, 0, 0); //Transfer the correct values at owned vertices to halo vertices //communicate_and_take_halo_contributions(stress_divUV, numVertices, 2, 1, 1); } Kokkos::fence(); - pumipic::RecordTime("Stress_Divergence Communication" + std::to_string(self), timer.seconds()); } void MPMesh::calcBasis() { diff --git a/src/pmpo_MPMesh.hpp b/src/pmpo_MPMesh.hpp index dfd922e5..ba75b56d 100644 --- a/src/pmpo_MPMesh.hpp +++ b/src/pmpo_MPMesh.hpp @@ -7,6 +7,11 @@ #include #include #include +#include +#include +#ifdef KOKKOS_ENABLE_CUDA +#include +#endif namespace polyMPO{ @@ -37,21 +42,15 @@ class MPMesh{ std::vector> ownerHaloLocalIDs; void startCommunication(); - - void communicate_and_take_halo_contributions( - const Kokkos::View& meshField, - int nEntities, - int numEntries, - int mode, - int op); - - // Original CPU-staging function + void communicate_and_take_halo_contributions(const Kokkos::View& meshField, int nEntities, + int numEntries, int mode, int op); + template - void communicate_and_take_halo_contributions1( + void communicate_and_take_halo_contributions_staged( const ViewType& meshField, int nEntities, int numEntries, - int mode, + int mode , int op){ int self; @@ -59,155 +58,94 @@ class MPMesh{ MPI_Comm_rank(comm, &self); Kokkos::Timer timer; - auto reconVals_host = - Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), meshField); - + auto reconVals_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), meshField); pumipic::RecordTime("SD: GPU-CPU copy-" + std::to_string(self), timer.seconds()); timer.reset(); - std::vector> recvIDVec; + std::vector> recvIDVec; std::vector> recvDataVec; - pumipic::RecordTime("SD: Recv Vec Allocation-" + std::to_string(self), timer.seconds()); - - timer.reset(); - - communicateFields1( - reconVals_host, - nEntities, - numEntries, - mode, - recvIDVec, - recvDataVec); + timer.reset(); + //communicateFieldsFromHostView(fieldData1, nEntities, numEntries, mode, recvIDVec, recvDataVec); + communicateFieldsFromHostView(reconVals_host, nEntities, numEntries, mode, recvIDVec, recvDataVec); pumipic::RecordTime("SD: IP Comm-" + std::to_string(self), timer.seconds()); timer.reset(); - int numProcsTot = recvIDVec.size(); - + //Flatten IDs int totalSize = 0; std::vector offsets(numProcsTot, 0); - - for(int i = 0; i < numProcsTot; i++){ + for(int i=0; i recvIDGPU("recvIDGPU", totalSize); - auto hostView = - Kokkos::View("recvIDCPU", totalSize); - - for(int i = 0; i < numProcsTot; i++){ - std::copy( - recvIDVec[i].begin(), - recvIDVec[i].end(), - hostView.data() + offsets[i]); + std::vector flatIDVec(totalSize, 0); + for(int i=0; i recvIDGPU("recvIDGPU", totalSize); + auto hostView = Kokkos::View("recvIDCPU", totalSize); + std::copy(flatIDVec.begin(), flatIDVec.end(), hostView.data()); Kokkos::deep_copy(recvIDGPU, hostView); Kokkos::fence(); - pumipic::RecordTime("SD: Copy CPU-GPU-" + std::to_string(self), timer.seconds()); + //Flatten Data timer.reset(); - - int totalSize_data = 0; + int totalSize_data=0; std::vector offsets_data(numProcsTot, 0); - - for(int i = 0; i < numProcsTot; i++){ + for(int i=0; i recvDataGPU("recvDataGPU", totalSize_data); - auto hostView_data = - Kokkos::View("recvDataCPU", totalSize_data); - - for(int i = 0; i < numProcsTot; i++){ - std::copy( - recvDataVec[i].begin(), - recvDataVec[i].end(), - hostView_data.data() + offsets_data[i]); + std::vector flatDataVec(totalSize_data, 0); + for(int i=0; i recvDataGPU("recvDataGPU", totalSize_data); + auto hostView_data= Kokkos::View("recvDataCPU", totalSize_data); + std::copy(flatDataVec.begin(), flatDataVec.end(), hostView_data.data()); Kokkos::deep_copy(recvDataGPU, hostView_data); Kokkos::fence(); - - assert(totalSize_data == totalSize * numEntries); - - for(int i = 0; i < numProcsTot; i++){ + assert(totalSize_data == totalSize*numEntries); + for (int i=0; i>& fieldData, const int numEntities, + const int numEntries, int mode, std::vector>& recvIDVec, + std::vector>& recvDataVec); - void communicateFields( - const std::vector>& fieldData, - const int numEntities, - const int numEntries, - int mode, - std::vector>& recvIDVec, - std::vector>& recvDataVec); - - - template - void communicateFields1( + template + void communicateFieldsFromHostView( const ViewType& fieldData, - const int numEntities, - const int numEntries, - int mode, + const int numEntities, const int numEntries, int mode, std::vector>& recvIDVec, std::vector>& recvDataVec){ - - int self, numProcsTot; + int self, numProcsTot; MPI_Comm comm = p_MPs->getMPIComm(); - MPI_Comm_rank(comm, &self); MPI_Comm_size(comm, &numProcsTot); @@ -219,191 +157,98 @@ class MPMesh{ recvDataVec.resize(numProcsTot); for(int i = 0; i < numProcsTot; i++){ - if(i == self) continue; - - int numToSend = 0; - int numToRecv = 0; + if(i==self) continue; - if(mode == 0){ + int numToSend = 0, numToRecv = 0; + if(mode == 0) { + //gather (halos send to owners) numToSend = numOwnersOnOtherProcs[i]; numToRecv = numHalosOnOtherProcs[i]; } else{ + //scatter (owners send to halos) numToSend = numHalosOnOtherProcs[i]; numToRecv = numOwnersOnOtherProcs[i]; } if(numToSend > 0){ - sendDataVec[i].reserve(numToSend * numEntries); + sendDataVec[i].reserve(numToSend*numEntries); } - if(numToRecv > 0){ - recvDataVec[i].resize(numToRecv * numEntries); + recvDataVec[i].resize(numToRecv*numEntries); recvIDVec[i].resize(numToRecv); } } if(mode == 0){ - for(int iEnt = 0; iEnt < numHalosTot; iEnt++){ + // Halos sends to owners + for (int iEnt = 0; iEnt < numHalosTot; iEnt++){ auto ownerProc = haloOwnerProcs[iEnt]; - - for(int iDouble = 0; iDouble < numEntries; iDouble++){ - sendDataVec[ownerProc].push_back( - fieldData(numOwnersTot + iEnt, iDouble)); - } + for (int iDouble = 0; iDouble < numEntries; iDouble++) + sendDataVec[ownerProc].push_back(fieldData(numOwnersTot+iEnt, iDouble)); } } else if(mode == 1){ - for(size_t iProc = 0; iProc < ownerOwnerLocalIDs.size(); iProc++){ - for(auto& ownerID : ownerOwnerLocalIDs[iProc]){ - for(int iDouble = 0; iDouble < numEntries; iDouble++){ - sendDataVec[iProc].push_back( - fieldData(ownerID, iDouble)); - } + // Owner sends to halos + for (size_t iProc=0; iProc requests; - requests.reserve(4 * numProcsTot); - + requests.reserve(4*numProcsTot); for(int proc = 0; proc < numProcsTot; proc++){ if(proc == self) continue; - if(mode == 0 && numHalosOnOtherProcs[proc]){ - assert(recvIDVec[proc].size() == - static_cast(numHalosOnOtherProcs[proc])); - - assert(recvDataVec[proc].size() == - recvIDVec[proc].size() * static_cast(numEntries)); - - MPI_Request req3; - MPI_Request req4; - - MPI_Irecv( - recvIDVec[proc].data(), - recvIDVec[proc].size(), - MPI_INT, - proc, - 1, - comm, - &req3); - - MPI_Irecv( - recvDataVec[proc].data(), - recvDataVec[proc].size(), - MPI_DOUBLE, - proc, - 2, - comm, - &req4); - + assert(recvIDVec[proc].size() == (size_t)numHalosOnOtherProcs[proc]); + assert(recvDataVec[proc].size() == recvIDVec[proc].size() * (size_t)numEntries); + MPI_Request req3, req4; + MPI_Irecv(recvIDVec[proc].data(), recvIDVec[proc].size(), MPI_INT, proc, 1, comm, &req3); + MPI_Irecv(recvDataVec[proc].data(), recvDataVec[proc].size(), MPI_DOUBLE, proc, 2, comm, &req4); requests.push_back(req3); requests.push_back(req4); } - - if(mode == 0 && numOwnersOnOtherProcs[proc]){ - assert(haloOwnerLocalIDs[proc].size() == - static_cast(numOwnersOnOtherProcs[proc])); - - assert(sendDataVec[proc].size() == - haloOwnerLocalIDs[proc].size() * static_cast(numEntries)); - - MPI_Request req1; - MPI_Request req2; - - MPI_Isend( - haloOwnerLocalIDs[proc].data(), - haloOwnerLocalIDs[proc].size(), - MPI_INT, - proc, - 1, - comm, - &req1); - - MPI_Isend( - sendDataVec[proc].data(), - sendDataVec[proc].size(), - MPI_DOUBLE, - proc, - 2, - comm, - &req2); - + if(mode == 0 && numOwnersOnOtherProcs[proc]) { + assert(haloOwnerLocalIDs[proc].size() == (size_t)numOwnersOnOtherProcs[proc]); + assert(sendDataVec[proc].size() == haloOwnerLocalIDs[proc].size() * (size_t)numEntries); + MPI_Request req1, req2; + MPI_Isend(haloOwnerLocalIDs[proc].data(), haloOwnerLocalIDs[proc].size(), MPI_INT, proc, 1, comm, &req1); + MPI_Isend(sendDataVec[proc].data(), sendDataVec[proc].size(), MPI_DOUBLE, proc, 2, comm, &req2); requests.push_back(req1); requests.push_back(req2); } if(mode == 1 && numOwnersOnOtherProcs[proc]){ - MPI_Request req3; - MPI_Request req4; - - MPI_Irecv( - recvIDVec[proc].data(), - recvIDVec[proc].size(), - MPI_INT, - proc, - 1, - comm, - &req3); - - MPI_Irecv( - recvDataVec[proc].data(), - recvDataVec[proc].size(), - MPI_DOUBLE, - proc, - 2, - comm, - &req4); - + MPI_Request req3, req4; + MPI_Irecv(recvIDVec[proc].data(), recvIDVec[proc].size(), MPI_INT, proc, 1, comm, &req3); + MPI_Irecv(recvDataVec[proc].data(), recvDataVec[proc].size(), MPI_DOUBLE, proc, 2, comm, &req4); requests.push_back(req3); requests.push_back(req4); } - - if(mode == 1 && numHalosOnOtherProcs[proc]){ - MPI_Request req1; - MPI_Request req2; - - MPI_Isend( - ownerHaloLocalIDs[proc].data(), - ownerHaloLocalIDs[proc].size(), - MPI_INT, - proc, - 1, - comm, - &req1); - - MPI_Isend( - sendDataVec[proc].data(), - sendDataVec[proc].size(), - MPI_DOUBLE, - proc, - 2, - comm, - &req2); - + if(mode == 1 && numHalosOnOtherProcs[proc]) { + MPI_Request req1, req2; + MPI_Isend(ownerHaloLocalIDs[proc].data(), ownerHaloLocalIDs[proc].size(), MPI_INT, proc, 1, comm, &req1); + MPI_Isend(sendDataVec[proc].data(), sendDataVec[proc].size(), MPI_DOUBLE, proc, 2, comm, &req2); requests.push_back(req1); requests.push_back(req2); } } - MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE); } - MPMesh(Mesh* inMesh, MaterialPoints* inMPs): - p_mesh(inMesh), - p_MPs(inMPs) { + p_mesh(inMesh), p_MPs(inMPs) { }; - - ~MPMesh(){ + ~MPMesh() { delete p_mesh; delete p_MPs; } - - // MP advection and tracking + //MP advection and tracking void CVTTrackingEdgeCenterBased(Vec2dView dx); void CVTTrackingElmCenterBased(const int printVTPIndex = -1); void T2LTracking(Vec2dView dx); @@ -413,86 +258,139 @@ class MPMesh{ void push_swap_pos(); void push(); - - // Used before advection to interpolate fields from mesh to MPs - // And also before reconstruction + //Used before advection to interpolate fields from mesh to MPs + //And also before reconstruction void calcBasis(); - - // Reconstruction + //Reconstruction DoubleView assemblyV0(); - template void assemblyVtx0(); - template void assemblyElm0(); - template void assemblyVtx1(); - void reconstruct_coeff_full(); - - void invertMatrix( - const Kokkos::View& vtxMatrices, - const double& radius); - + void invertMatrix(const Kokkos::View& vtxMatrices, const double& radius); Kokkos::View precomputedVtxCoeffs_new; Kokkos::View nearAnEdge; Kokkos::View vtxMatrixMass; - - // Not used currently - std::map> reconstructSlice = - std::map>(); - + //Not used currently + std::map> reconstructSlice = std::map>(); template DoubleView wtScaAssembly(); - template Vec2dView wtVec2Assembly(); - template - void assembly( - int order, - MeshFieldType type, - bool basisWeightFlag, - bool massWeightFlag); - + void assembly(int order, MeshFieldType type, bool basisWeightFlag, bool massWeightFlag); template - void setReconstructSlice( - int order, - MeshFieldType type); - + void setReconstructSlice(int order, MeshFieldType type); void reconstructSlices(); void printVTP_mesh(int printVTPIndex); - - void writeMPTrackingVTP( - int printVTPIndex, - int numMPs, - const Vec3dView& history, - const Vec3dView& resultLeft, - const Vec3dView& resultRight, - const Vec3dView& mpTgtPosArray); + void writeMPTrackingVTP(int printVTPIndex, int numMPs, const Vec3dView& history, const Vec3dView& resultLeft, + const Vec3dView& resultRight, const Vec3dView& mpTgtPosArray); void calculateStrain(); void calculateStress(const int constitutive_relation); void calculateStressDivergence(); + //Prints the selected halo-exchange path once, on rank 0 + bool haloExchangePathReported = false; + void reportHaloExchangePath(const char* path){ + if(haloExchangePathReported) return; + haloExchangePathReported = true; + int self; + MPI_Comm_rank(p_MPs->getMPIComm(), &self); + if(self == 0){ + std::cout << "polyMPO halo exchange: " << path << std::endl; + } + } + +#ifdef GPU_AWARE_MPI + +#ifdef KOKKOS_ENABLE_CUDA + //Device buffer allocated with plain cudaMalloc and exposed as an + //unmanaged Kokkos::View. Used for buffers passed directly to MPI. + // + //Kokkos (>= 4.2, CUDA_MALLOC_ASYNC=ON by default) allocates Views from + //cudaMallocAsync memory pools, which cuIpcGetMemHandle rejects + //(CUDA_ERROR_INVALID_VALUE). Cray MPICH/GTL relies on CUDA IPC for + //intra-node GPU-to-GPU transfers, so these buffers bypass the Kokkos + //allocator regardless of how Kokkos was built. + template + struct RawCudaMPIBuffer{ + T* ptr = nullptr; + size_t count = 0; + + void allocate(size_t n){ + free(); + count = n; + if(n > 0){ + cudaError_t err = cudaMalloc(&ptr, n * sizeof(T)); + if(err != cudaSuccess){ + throw std::runtime_error(std::string("RawCudaMPIBuffer: cudaMalloc failed: ") + cudaGetErrorString(err)); + } + } + } + + void free(){ + if(ptr != nullptr){ cudaFree(ptr); ptr = nullptr; } + count = 0; + } + + T* data() const{ return ptr; } + size_t size() const{ return count; } + Kokkos::View view() const{ + return Kokkos::View(ptr, count); + } + + RawCudaMPIBuffer() = default; + ~RawCudaMPIBuffer(){ free(); } + RawCudaMPIBuffer(const RawCudaMPIBuffer&) = delete; + RawCudaMPIBuffer& operator=(const RawCudaMPIBuffer&) = delete; + RawCudaMPIBuffer(RawCudaMPIBuffer&& other) noexcept{ + ptr = other.ptr; count = other.count; + other.ptr = nullptr; other.count = 0; + } -#ifdef CUDA_AWARE_MPI + RawCudaMPIBuffer& operator=(RawCudaMPIBuffer&& other) noexcept{ + if(this != &other){ + free(); + ptr = other.ptr; count = other.count; + other.ptr = nullptr; other.count = 0; + } + return *this; + } + }; + + using CudaAwareMPIIntBuffer = RawCudaMPIBuffer; + using CudaAwareMPIDoubleBuffer = RawCudaMPIBuffer; +#else + //Non-CUDA backend: plain Kokkos::View with the same interface as RawCudaMPIBuffer + template + struct KokkosMPIBuffer{ + Kokkos::View v; + + void allocate(size_t n){ + v = Kokkos::View("cudaAwareMPIBuffer_batched", n); + } + + T* data() const{ return v.data(); } + size_t size() const{ return v.extent(0); } + Kokkos::View view() const{ return v; } + }; + + using CudaAwareMPIIntBuffer = KokkosMPIBuffer; + using CudaAwareMPIDoubleBuffer = KokkosMPIBuffer; +#endif - // Cached CUDA-aware MPI communication metadata and per-neighbor GPU buffers. - // Important change from the previous version: - // MPI is always given the base pointer of a Kokkos allocation, not - // "base pointer + offset". This avoids Cray MPICH/GTL CUDA IPC problems. bool cudaAwareMPICacheValid = true; bool cudaAwareMPIDisabled = false; - bool cudaAwareMPIEnvChecked = false; - bool cudaAwareMPIForceCPU = false; - bool cudaAwareMPILogged = false; + bool mpiGpuSupportChecked = false; + bool mpiGpuSupportEnabled = false; struct CudaAwareMPIFieldCache{ bool valid = false; @@ -500,37 +398,46 @@ class MPMesh{ std::vector sendCounts; std::vector recvCounts; - - std::vector> sendEntityGPUPerProc; - std::vector> recvIDGPUPerProc; - - std::vector> sendDataGPUPerProc; - std::vector> recvDataGPUPerProc; + std::vector sendOffsets; //prefix sum of sendCounts, in entities + std::vector recvOffsets; //prefix sum of recvCounts, in entities + + int totalSendCount = 0; + int totalRecvCount = 0; + + //Batched buffers for all neighbors. Per-proc slices are + //[offset, offset + count) for IDs and + //[offset * numEntries, (offset + count) * numEntries) for data. + CudaAwareMPIIntBuffer sendEntityGPU; + CudaAwareMPIIntBuffer recvIDGPU; + CudaAwareMPIDoubleBuffer sendDataGPU; + CudaAwareMPIDoubleBuffer recvDataGPU; }; std::map, CudaAwareMPIFieldCache> cudaAwareMPICaches; - bool cudaAwareMPIForceDisabled(){ - if(!cudaAwareMPIEnvChecked){ - const char* value = std::getenv("POLYMPO_DISABLE_CUDA_AWARE_MPI"); - cudaAwareMPIForceCPU = - value != nullptr && value[0] != '\0' && value[0] != '0'; - cudaAwareMPIEnvChecked = true; + //The halo-exchange path follows MPICH_GPU_SUPPORT_ENABLED: + // 1 -> GPU-aware MPI + // anything else (0, unset, yes, ...) -> CPU-staged path + //MPI must not be given device pointers unless GPU support is enabled; + //doing so crashes instead of returning an MPI error. + bool mpiGpuSupport(){ + if(!mpiGpuSupportChecked){ + const char* value = std::getenv("MPICH_GPU_SUPPORT_ENABLED"); + mpiGpuSupportEnabled = value != nullptr && std::string(value) == "1"; + mpiGpuSupportChecked = true; } - - return cudaAwareMPIForceCPU; + return mpiGpuSupportEnabled; } - // Fully CUDA-aware MPI version: - // Field data is sent/received using GPU pointers. Receive IDs are cached - // once from the fixed halo/owner mapping and are not sent every call. + //GPU-aware MPI path: field data is sent/received directly from GPU + //buffers. Receive IDs come from the fixed halo/owner mapping and are + //not exchanged. All neighbors share one batched buffer per direction, + //packed/unpacked with a single kernel each. // - // Important: - // This function caches communication metadata and GPU buffers per - // (mode, numEntries). If the communication pattern changes, clear - // cudaAwareMPICaches before the next call. + //Metadata and buffers are cached per (mode, numEntries). If the + //communication pattern changes, clear cudaAwareMPICaches first. template - void communicate_and_take_halo_contributions1_improved( + void communicate_and_take_halo_contributions_gpu_aware( const ViewType& meshField, int nEntities, int numEntries, @@ -538,9 +445,7 @@ class MPMesh{ int op){ int self, numProcsTot; - MPI_Comm comm = p_MPs->getMPIComm(); - MPI_Comm_rank(comm, &self); MPI_Comm_size(comm, &numProcsTot); @@ -548,59 +453,31 @@ class MPMesh{ assert(op == 0 || op == 1); assert(nEntities == numOwnersTot + numHalosTot); - if(cudaAwareMPIDisabled || cudaAwareMPIForceDisabled()){ - communicate_and_take_halo_contributions1( - meshField, - nEntities, - numEntries, - mode, - op); + if(cudaAwareMPIDisabled || !mpiGpuSupport()){ + reportHaloExchangePath("CPU-staged (MPICH_GPU_SUPPORT_ENABLED is not 1)"); + communicate_and_take_halo_contributions_staged(meshField, nEntities, numEntries, mode, op); return; } -#ifdef POLYMPO_VERBOSE_MPI - if(self == 0 && !cudaAwareMPILogged){ - std::cout - << "[CUDA_AWARE_MPI] Using per-proc cached full GPU-aware MPI path in communicate_and_take_halo_contributions1_improved()" - << "\n"; - cudaAwareMPILogged = true; - } -#endif - - Kokkos::Timer timer; + reportHaloExchangePath("GPU-aware MPI"); if(!cudaAwareMPICacheValid){ cudaAwareMPICaches.clear(); cudaAwareMPICacheValid = true; } - auto& cudaAwareCache = - cudaAwareMPICaches[std::make_pair(mode, numEntries)]; - - const bool needRebuild = - (!cudaAwareCache.valid) || - (cudaAwareCache.cachedNumProcs != numProcsTot); + auto& cudaAwareCache = cudaAwareMPICaches[std::make_pair(mode, numEntries)]; + const bool needRebuild = (!cudaAwareCache.valid) || (cudaAwareCache.cachedNumProcs != numProcsTot); if(needRebuild){ - cudaAwareCache.cachedNumProcs = numProcsTot; - cudaAwareCache.sendCounts.assign(numProcsTot, 0); cudaAwareCache.recvCounts.assign(numProcsTot, 0); - - cudaAwareCache.sendEntityGPUPerProc.clear(); - cudaAwareCache.recvIDGPUPerProc.clear(); - cudaAwareCache.sendDataGPUPerProc.clear(); - cudaAwareCache.recvDataGPUPerProc.clear(); - - cudaAwareCache.sendEntityGPUPerProc.resize(numProcsTot); - cudaAwareCache.recvIDGPUPerProc.resize(numProcsTot); - cudaAwareCache.sendDataGPUPerProc.resize(numProcsTot); - cudaAwareCache.recvDataGPUPerProc.resize(numProcsTot); + cudaAwareCache.sendOffsets.assign(numProcsTot, 0); + cudaAwareCache.recvOffsets.assign(numProcsTot, 0); for(int proc = 0; proc < numProcsTot; proc++){ if(proc == self) continue; - if(mode == 0){ cudaAwareCache.sendCounts[proc] = numOwnersOnOtherProcs[proc]; cudaAwareCache.recvCounts[proc] = numHalosOnOtherProcs[proc]; @@ -611,319 +488,201 @@ class MPMesh{ } } + int totalSend = 0; + int totalRecv = 0; for(int proc = 0; proc < numProcsTot; proc++){ - if(proc == self) continue; - - const int sendCount = cudaAwareCache.sendCounts[proc]; - const int recvCount = cudaAwareCache.recvCounts[proc]; - - if(sendCount > 0){ - cudaAwareCache.sendEntityGPUPerProc[proc] = - Kokkos::View( - "cudaAwareMPISendEntityGPUPerProc", - sendCount); - - cudaAwareCache.sendDataGPUPerProc[proc] = - Kokkos::View( - "cudaAwareMPISendDataGPUPerProc", - sendCount * numEntries); - - auto sendEntityCPU = - Kokkos::View( - "sendEntityCPU", - sendCount); - - if(mode == 0){ - assert(haloOwnerLocalIDs[proc].size() == - static_cast(sendCount)); - - int localIndex = 0; - - for(int iEnt = 0; iEnt < numHalosTot; iEnt++){ - int ownerProc = haloOwnerProcs[iEnt]; - - if(ownerProc != proc) continue; - - assert(localIndex < sendCount); + cudaAwareCache.sendOffsets[proc] = totalSend; + totalSend += cudaAwareCache.sendCounts[proc]; + cudaAwareCache.recvOffsets[proc] = totalRecv; + totalRecv += cudaAwareCache.recvCounts[proc]; + } + cudaAwareCache.totalSendCount = totalSend; + cudaAwareCache.totalRecvCount = totalRecv; - sendEntityCPU(localIndex) = numOwnersTot + iEnt; + cudaAwareCache.sendEntityGPU.allocate(totalSend); + cudaAwareCache.sendDataGPU.allocate(totalSend * numEntries); + cudaAwareCache.recvIDGPU.allocate(totalRecv); + cudaAwareCache.recvDataGPU.allocate(totalRecv * numEntries); - localIndex++; - } + //Build flattened send-entity list on host, then copy to device + if(totalSend > 0){ + auto sendEntityCPU = Kokkos::View("sendEntityCPU_batched", totalSend); + if(mode == 0){ + for(int proc = 0; proc < numProcsTot; proc++){ + if(proc == self) continue; + if(cudaAwareCache.sendCounts[proc] <= 0) continue; + assert(haloOwnerLocalIDs[proc].size() == (size_t)cudaAwareCache.sendCounts[proc]); + } - assert(localIndex == sendCount); + std::vector cursor(cudaAwareCache.sendOffsets); + for(int iEnt = 0; iEnt < numHalosTot; iEnt++){ + int ownerProc = haloOwnerProcs[iEnt]; + if(ownerProc == self) continue; + sendEntityCPU(cursor[ownerProc]) = numOwnersTot + iEnt; + cursor[ownerProc]++; } - else{ - assert(ownerOwnerLocalIDs[proc].size() == - static_cast(sendCount)); + for(int proc = 0; proc < numProcsTot; proc++){ + if(proc == self) continue; + assert(cursor[proc] == cudaAwareCache.sendOffsets[proc] + cudaAwareCache.sendCounts[proc]); + } + } + else{ + for(int proc = 0; proc < numProcsTot; proc++){ + if(proc == self) continue; + int sendCount = cudaAwareCache.sendCounts[proc]; + if(sendCount <= 0) continue; + assert(ownerOwnerLocalIDs[proc].size() == (size_t)sendCount); + int base = cudaAwareCache.sendOffsets[proc]; for(int i = 0; i < sendCount; i++){ - sendEntityCPU(i) = ownerOwnerLocalIDs[proc][i]; + sendEntityCPU(base + i) = ownerOwnerLocalIDs[proc][i]; } } - - Kokkos::deep_copy( - cudaAwareCache.sendEntityGPUPerProc[proc], - sendEntityCPU); } + Kokkos::deep_copy(cudaAwareCache.sendEntityGPU.view(), sendEntityCPU); + } - if(recvCount > 0){ - cudaAwareCache.recvIDGPUPerProc[proc] = - Kokkos::View( - "cudaAwareMPIRecvIDGPUPerProc", - recvCount); - - cudaAwareCache.recvDataGPUPerProc[proc] = - Kokkos::View( - "cudaAwareMPIRecvDataGPUPerProc", - recvCount * numEntries); - - auto recvIDCPU = - Kokkos::View( - "recvIDCPU", - recvCount); - - if(mode == 0){ - assert(ownerOwnerLocalIDs[proc].size() == - static_cast(recvCount)); - + //Build flattened recv-ID list on host, then copy to device + if(totalRecv > 0){ + auto recvIDCPU = Kokkos::View("recvIDCPU_batched", totalRecv); + if(mode == 0){ + for(int proc = 0; proc < numProcsTot; proc++){ + if(proc == self) continue; + int recvCount = cudaAwareCache.recvCounts[proc]; + if(recvCount <= 0) continue; + assert(ownerOwnerLocalIDs[proc].size() == (size_t)recvCount); + int base = cudaAwareCache.recvOffsets[proc]; for(int i = 0; i < recvCount; i++){ - recvIDCPU(i) = ownerOwnerLocalIDs[proc][i]; + recvIDCPU(base + i) = ownerOwnerLocalIDs[proc][i]; } } - else{ - int localIndex = 0; - - for(int iEnt = 0; iEnt < numHalosTot; iEnt++){ - if(haloOwnerProcs[iEnt] != proc) continue; - - assert(localIndex < recvCount); - - recvIDCPU(localIndex) = numOwnersTot + iEnt; - - localIndex++; - } - - assert(localIndex == recvCount); + } + else{ + std::vector cursor(cudaAwareCache.recvOffsets); + for(int iEnt = 0; iEnt < numHalosTot; iEnt++){ + int ownerProc = haloOwnerProcs[iEnt]; + if(ownerProc == self) continue; + recvIDCPU(cursor[ownerProc]) = numOwnersTot + iEnt; + cursor[ownerProc]++; } - Kokkos::deep_copy( - cudaAwareCache.recvIDGPUPerProc[proc], - recvIDCPU); + for(int proc = 0; proc < numProcsTot; proc++){ + if(proc == self) continue; + assert(cursor[proc] == cudaAwareCache.recvOffsets[proc] + cudaAwareCache.recvCounts[proc]); + } } + Kokkos::deep_copy(cudaAwareCache.recvIDGPU.view(), recvIDCPU); } - Kokkos::fence(); - cudaAwareCache.valid = true; - - pumipic::RecordTime( - "SD: CUDA-aware MPI Cache Build m" + std::to_string(mode) + - " e" + std::to_string(numEntries) + "-" + std::to_string(self), - timer.seconds()); - - timer.reset(); } - for(int proc = 0; proc < numProcsTot; proc++){ - if(proc == self) continue; - if(cudaAwareCache.sendCounts[proc] <= 0) continue; - - auto sendEntityGPU = cudaAwareCache.sendEntityGPUPerProc[proc]; - auto sendDataGPU = cudaAwareCache.sendDataGPUPerProc[proc]; - int sendCount = cudaAwareCache.sendCounts[proc]; - - Kokkos::parallel_for( - "pack cached cuda-aware mpi send buffer per proc", - sendCount, - KOKKOS_LAMBDA(const int i){ - int entity = sendEntityGPU(i); - - for(int k = 0; k < numEntries; k++){ - sendDataGPU(i * numEntries + k) = - meshField(entity, k); - } - }); + //Pack send buffer + if(cudaAwareCache.totalSendCount > 0){ + auto sendEntityGPU = cudaAwareCache.sendEntityGPU.view(); + auto sendDataGPU = cudaAwareCache.sendDataGPU.view(); + Kokkos::parallel_for("pack cached gpu-aware mpi send buffer batched", cudaAwareCache.totalSendCount, KOKKOS_LAMBDA(const int i){ + int entity = sendEntityGPU(i); + for(int k = 0; k < numEntries; k++){ + sendDataGPU(i * numEntries + k) = meshField(entity, k); + } + }); } - + //MPI is not stream-aware: packing must finish before Isend Kokkos::fence(); - pumipic::RecordTime( - "SD: CUDA-aware MPI Pack m" + std::to_string(mode) + - " e" + std::to_string(numEntries) + "-" + std::to_string(self), - timer.seconds()); - - timer.reset(); - - Kokkos::Timer mpiTotalTimer; std::vector requests; requests.reserve(2 * numProcsTot); int mpiError = MPI_SUCCESS; + //Post Irecv/Isend per neighbor for(int proc = 0; proc < numProcsTot; proc++){ if(proc == self) continue; if(cudaAwareCache.recvCounts[proc] > 0){ MPI_Request reqData; - - mpiError = MPI_Irecv( - cudaAwareCache.recvDataGPUPerProc[proc].data(), - cudaAwareCache.recvCounts[proc] * numEntries, - MPI_DOUBLE, - proc, - 2, - comm, - &reqData); + double* recvPtr = cudaAwareCache.recvDataGPU.data() + (size_t)cudaAwareCache.recvOffsets[proc] * numEntries; + mpiError = MPI_Irecv(recvPtr, cudaAwareCache.recvCounts[proc] * numEntries, MPI_DOUBLE, proc, 2, comm, &reqData); if(mpiError != MPI_SUCCESS) break; requests.push_back(reqData); } if(cudaAwareCache.sendCounts[proc] > 0){ MPI_Request reqData; - - mpiError = MPI_Isend( - cudaAwareCache.sendDataGPUPerProc[proc].data(), - cudaAwareCache.sendCounts[proc] * numEntries, - MPI_DOUBLE, - proc, - 2, - comm, - &reqData); + double* sendPtr = cudaAwareCache.sendDataGPU.data() + (size_t)cudaAwareCache.sendOffsets[proc] * numEntries; + mpiError = MPI_Isend(sendPtr, cudaAwareCache.sendCounts[proc] * numEntries, MPI_DOUBLE, proc, 2, comm, &reqData); if(mpiError != MPI_SUCCESS) break; requests.push_back(reqData); } } - pumipic::RecordTime( - "SD: CUDA-aware MPI Post m" + std::to_string(mode) + - " e" + std::to_string(numEntries) + "-" + std::to_string(self), - timer.seconds()); - - timer.reset(); - if(mpiError == MPI_SUCCESS && !requests.empty()){ - mpiError = MPI_Waitall( - static_cast(requests.size()), - requests.data(), - MPI_STATUSES_IGNORE); + mpiError = MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE); } - pumipic::RecordTime( - "SD: CUDA-aware MPI Wait m" + std::to_string(mode) + - " e" + std::to_string(numEntries) + "-" + std::to_string(self), - timer.seconds()); - if(mpiError != MPI_SUCCESS){ cudaAwareMPIDisabled = true; - if(self == 0){ - std::cout - << "[CUDA_AWARE_MPI] Device-pointer MPI failed." - << std::endl; + std::cout << "[GPU_AWARE_MPI] Batched device-pointer MPI failed." << std::endl; } if(requests.empty()){ if(self == 0){ - std::cout - << "[CUDA_AWARE_MPI] Falling back to CPU-staged communication." - << std::endl; + std::cout << "[GPU_AWARE_MPI] Falling back to CPU-staged communication." << std::endl; } - - communicate_and_take_halo_contributions1( - meshField, - nEntities, - numEntries, - mode, - op); + communicate_and_take_halo_contributions_staged(meshField, nEntities, numEntries, mode, op); return; } if(self == 0){ - std::cout - << "[CUDA_AWARE_MPI] Failure happened after MPI requests were posted. " - << "Set POLYMPO_DISABLE_CUDA_AWARE_MPI=1 before running to force the CPU-staged path." - << std::endl; + std::cout << "[GPU_AWARE_MPI] Failure happened after MPI requests were posted. " + << "Set MPICH_GPU_SUPPORT_ENABLED=0 before running to use the CPU-staged path." << std::endl; } - MPI_Abort(comm, mpiError); return; } - pumipic::RecordTime( - "SD: CUDA-aware MPI Comm m" + std::to_string(mode) + - " e" + std::to_string(numEntries) + "-" + std::to_string(self), - mpiTotalTimer.seconds()); - - timer.reset(); - - for(int proc = 0; proc < numProcsTot; proc++){ - if(proc == self) continue; - if(cudaAwareCache.recvCounts[proc] <= 0) continue; - - auto recvIDGPU = cudaAwareCache.recvIDGPUPerProc[proc]; - auto recvDataGPU = cudaAwareCache.recvDataGPUPerProc[proc]; - int recvCount = cudaAwareCache.recvCounts[proc]; - + //Unpack received contributions + if(cudaAwareCache.totalRecvCount > 0){ + auto recvIDGPU = cudaAwareCache.recvIDGPU.view(); + auto recvDataGPU = cudaAwareCache.recvDataGPU.view(); if(op == 0){ - Kokkos::parallel_for( - "halo add cached cuda-aware mpi per proc", - recvCount, - KOKKOS_LAMBDA(const int i){ - const int vertex = recvIDGPU(i); - - for(int k = 0; k < numEntries; k++){ + Kokkos::parallel_for("halo add cached gpu-aware mpi batched", cudaAwareCache.totalRecvCount, KOKKOS_LAMBDA(const int i){ + const int vertex = recvIDGPU(i); + for(int k = 0; k < numEntries; k++){ #ifdef POLYMPO_ASSUME_UNIQUE_HALO_CONTRIBS - meshField(vertex, k) += - recvDataGPU(i * numEntries + k); + meshField(vertex, k) += recvDataGPU(i * numEntries + k); #else - Kokkos::atomic_add( - &meshField(vertex, k), - recvDataGPU(i * numEntries + k)); + Kokkos::atomic_add(&meshField(vertex, k), recvDataGPU(i * numEntries + k)); #endif - } - }); + } + }); } else{ - Kokkos::parallel_for( - "halo assign cached cuda-aware mpi per proc", - recvCount, - KOKKOS_LAMBDA(const int i){ - const int vertex = recvIDGPU(i); - - for(int k = 0; k < numEntries; k++){ - meshField(vertex, k) = - recvDataGPU(i * numEntries + k); - } - }); + Kokkos::parallel_for("halo assign cached gpu-aware mpi batched", cudaAwareCache.totalRecvCount, KOKKOS_LAMBDA(const int i){ + const int vertex = recvIDGPU(i); + for(int k = 0; k < numEntries; k++){ + meshField(vertex, k) = recvDataGPU(i * numEntries + k); + } + }); } } - Kokkos::fence(); - - pumipic::RecordTime( - "SD: CUDA-aware MPI Contribution m" + std::to_string(mode) + - " e" + std::to_string(numEntries) + "-" + std::to_string(self), - timer.seconds()); } #else - // Fallback path: - // if CUDA_AWARE_MPI is not defined, use the original GPU-CPU staging function. + //GPU_AWARE_MPI not defined: use the CPU-staged path template - void communicate_and_take_halo_contributions1_improved( + void communicate_and_take_halo_contributions_gpu_aware( const ViewType& meshField, int nEntities, int numEntries, int mode, int op){ - communicate_and_take_halo_contributions1( - meshField, - nEntities, - numEntries, - mode, - op); + reportHaloExchangePath("CPU-staged (built without GPU_AWARE_MPI)"); + communicate_and_take_halo_contributions_staged(meshField, nEntities, numEntries, mode, op); } #endif @@ -933,4 +692,3 @@ class MPMesh{ }//namespace polyMPO end #endif - diff --git a/src/pmpo_MPMesh_assembly.hpp b/src/pmpo_MPMesh_assembly.hpp index 7705c8f7..61609341 100644 --- a/src/pmpo_MPMesh_assembly.hpp +++ b/src/pmpo_MPMesh_assembly.hpp @@ -85,7 +85,7 @@ void MPMesh::assemblyElm0() { } }; p_MPs->parallel_for(assemble, "assembly"); - + Kokkos::MDRangePolicy> policy({0,0},{numElms, numEntries}); Kokkos::parallel_for("assembly average", policy, KOKKOS_LAMBDA(const int elm, const int entry){ if (mpsPerElm(elm) > 0){ @@ -96,12 +96,11 @@ void MPMesh::assemblyElm0() { } void MPMesh::reconstruct_coeff_full(){ - Kokkos::Timer timer; int self, numProcsTot; MPI_Comm comm = p_MPs->getMPIComm(); MPI_Comm_rank(comm, &self); MPI_Comm_size(comm, &numProcsTot); - + static int coeff_count=0; if(!self) std::cout<<"===="<<__FUNCTION__<<" "<parallel_for(assemble, "assembly"); Kokkos::fence(); - pumipic::RecordTime("Assemble Matrix Per Process" + std::to_string(self), timer.seconds()); //Mode 0 is Gather: Halos Send to Owners //Mode 1 is Scatter: Owners Send to Halos //Op 0 is addition //Op 1 is replacement - timer.reset(); int mode = 0; int op = 0; - if (numProcsTot >1){ - communicate_and_take_halo_contributions1_improved(vtxMatrices, numVertices, numEntriesMatrix, mode, op); - mode=1; + if (numProcsTot>1){ + communicate_and_take_halo_contributions_gpu_aware(vtxMatrices, numVertices, numEntriesMatrix, mode, op); + mode=1; op=1; - communicate_and_take_halo_contributions1_improved(vtxMatrices, numVertices, numEntriesMatrix, mode, op); + communicate_and_take_halo_contributions_gpu_aware(vtxMatrices, numVertices, numEntriesMatrix, mode, op); } - pumipic::RecordTime("Communicate Matrix Values" + std::to_string(self), timer.seconds()); - - //Stroe the 1st matrix element + //Store the 1st matrix element Kokkos::ViewvtxMatrixMass_l("vtxMass", numVertices); Kokkos::parallel_for("storeMatrixMass", numVertices, KOKKOS_LAMBDA(const int vtx){ vtxMatrixMass_l(vtx) = vtxMatrices(vtx, 0); - }); + }); + Kokkos::fence(); this->vtxMatrixMass = vtxMatrixMass_l; invertMatrix(vtxMatrices, radius); } void MPMesh::invertMatrix(const Kokkos::View& vtxMatrices, const double& radius){ - int nVertices = p_mesh->getNumVertices(); auto vtxCoords = p_mesh->getMeshField(); auto dual_triangle_area = p_mesh->getMeshField(); @@ -192,7 +187,7 @@ void MPMesh::invertMatrix(const Kokkos::View& vtxMatrices, const doubl Kokkos::View VtxCoeffs("VtxCoeffs", nVertices); Kokkos::deep_copy(VtxCoeffs, 0.0); Kokkos::View nearAnEdge_l("nearAnEdge_l", nVertices); - + Kokkos::parallel_for("invertMatrix", nVertices, KOKKOS_LAMBDA(const int vtx){ if(vtxMatrices(vtx, 0) < eps) return; @@ -209,7 +204,7 @@ void MPMesh::invertMatrix(const Kokkos::View& vtxMatrices, const doubl } auto cosLat = sqrt(pow(X, 2) + pow(Y, 2)); - auto invCosLat = 1.0/cosLat; + auto invCosLat = 1.0/cosLat; auto vtx_area_sqrt = sqrt(dual_triangle_area(vtx,0)/(radius*radius)); Vec3d v0 = { -Y * invCosLat, -Z * X * invCosLat, X / vtx_area_sqrt }; @@ -281,7 +276,7 @@ void MPMesh::invertMatrix(const Kokkos::View& vtxMatrices, const doubl iBlockC[0] = blockC[0] * invM2D[0] + blockC[1] * invM2D[1] + blockC[2] * invM2D[2]; iBlockC[1] = blockC[0] * invM2D[1] + blockC[1] * invM2D[3] + blockC[2] * invM2D[4]; iBlockC[2] = blockC[0] * invM2D[2] + blockC[1] * invM2D[4] + blockC[2] * invM2D[5]; - iBlockC = iBlockC*invM11; + iBlockC = iBlockC*invM11; VtxCoeffs(vtx, 0, 0) = invM11 + invM11*iBlockC.dot(blockC); auto temp = -rotateScaleM.rightMultiply(iBlockC); @@ -293,7 +288,7 @@ void MPMesh::invertMatrix(const Kokkos::View& vtxMatrices, const doubl VtxCoeffs(vtx, 1, 1) = invM2D[0] * rotateScaleM(0, 0) + invM2D[1] * rotateScaleM(0, 1) + invM2D[2] * rotateScaleM(0, 2); VtxCoeffs(vtx, 1, 2) = invM2D[0] * rotateScaleM(1, 0) + invM2D[1] * rotateScaleM(1, 1) + invM2D[2] * rotateScaleM(1, 2); VtxCoeffs(vtx, 1, 3) = invM2D[0] * rotateScaleM(2, 0) + invM2D[1] * rotateScaleM(2, 1) + invM2D[2] * rotateScaleM(2, 2); - + VtxCoeffs(vtx, 2, 0) = -iBlockC[1]; VtxCoeffs(vtx, 2, 1) = invM2D[1] * rotateScaleM(0, 0) + invM2D[3] * rotateScaleM(0, 1) + invM2D[4] * rotateScaleM(0, 2); VtxCoeffs(vtx, 2, 2) = invM2D[1] * rotateScaleM(1, 0) + invM2D[3] * rotateScaleM(1, 1) + invM2D[4] * rotateScaleM(1, 2); @@ -321,7 +316,6 @@ void MPMesh::invertMatrix(const Kokkos::View& vtxMatrices, const doubl template void MPMesh::assemblyVtx1(){ - Kokkos::Timer timer; int self, numProcsTot; MPI_Comm comm = p_MPs->getMPIComm(); @@ -358,7 +352,7 @@ void MPMesh::assemblyVtx1(){ int nVtxE = elm2VtxConn(elm,0); //number of vertices bounding the element for(int i=0; iparallel_for(reconstruct, "reconstruct"); Kokkos::fence(); - pumipic::RecordTime("Assemble Field per process" + std::to_string(self), timer.seconds()); - - timer.reset(); - if(numProcsTot>1){ - communicate_and_take_halo_contributions1_improved(meshField, numVertices, numEntries, 0, 0); + if(numProcsTot>1){ + communicate_and_take_halo_contributions_gpu_aware(meshField, numVertices, numEntries, 0, 0); } - pumipic::RecordTime("Communicate Field Values" + std::to_string(self), timer.seconds()); } template @@ -428,7 +418,7 @@ DoubleView MPMesh::wtScaAssembly(){ // last component of eVtxCoords stores the firs vertex (to avoid if-condition in the Wachspress computation) eVtxCoords[nElmVtxs][0] = vtxCoords(elm2VtxConn(elm,1)-1,0); eVtxCoords[nElmVtxs][1] = vtxCoords(elm2VtxConn(elm,1)-1,1); - + /* compute the values of basis functions at mp position */ double basisByArea[maxElmsPerVtx]; Vec2d mpCoord(mpPositions(mp,0), mpPositions(mp,1)); @@ -467,7 +457,7 @@ Vec2dView MPMesh::wtVec2Assembly(){ Vec2d eVtxCoords[maxVtxsPerElm + 1]; for (int i = 1; i <= nElmVtxs; i++) { // elm2VtxConn(elm,i) is the vertex ID (1-based index) of vertex #i of elm - eVtxCoords[i-1][0] = vtxCoords(elm2VtxConn(elm,i)-1,0); + eVtxCoords[i-1][0] = vtxCoords(elm2VtxConn(elm,i)-1,0); eVtxCoords[i-1][1] = vtxCoords(elm2VtxConn(elm,i)-1,1); } // last component of eVtxCoords stores the firs vertex (to avoid if-condition in the Wachspress computation) diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index 9c88f672..fd598cc1 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -6,7 +6,7 @@ namespace{ std::vector p_mpmeshes;////store the p_mpmeshes that is legal - + void checkMPMeshValid(MPMesh_ptr p_mpmesh){ auto p_mpmeshIter = std::find(p_mpmeshes.begin(),p_mpmeshes.end(),p_mpmesh); PMT_ALWAYS_ASSERT(p_mpmeshIter != p_mpmeshes.end()); @@ -39,7 +39,7 @@ MPMesh_ptr polympo_createMPMesh_f(const int testMeshOption, const int testMPOpti PMT_ALWAYS_ASSERT(testMeshOption >= 1); p_mps = polyMPO::initTestMPs(p_mesh, testMPOption); }else{ - p_mps = new polyMPO::MaterialPoints(); + p_mps = new polyMPO::MaterialPoints(); } MPMesh_ptr p_mpMeshReturn = (MPMesh_ptr) new polyMPO::MPMesh(p_mesh, p_mps); p_mpmeshes.push_back(p_mpMeshReturn); @@ -87,7 +87,7 @@ void polympo_createMPs_f(MPMesh_ptr p_mpmesh, for(int i = 0; i < numMPs; i++) { if(isMPActive[i] == MP_ACTIVE) { numActiveMPs++; - if(mp2Elm[i] < minElmID) + if(mp2Elm[i] < minElmID) minElmID = mp2Elm[i]; } } @@ -138,7 +138,7 @@ void polympo_createMPs_f(MPMesh_ptr p_mpmesh, 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); + p_MPs->setElmIDoffset(offset); } void polympo_startRebuildMPs_f(MPMesh_ptr p_mpmesh, @@ -151,7 +151,7 @@ void polympo_startRebuildMPs_f(MPMesh_ptr p_mpmesh, int self; MPI_Comm comm = p_MPs->getMPIComm(); MPI_Comm_rank(comm, &self); - + PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount()); //PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID()); @@ -225,7 +225,7 @@ void polympo_startRebuildMPs2_f(MPMesh_ptr p_mpmesh, recvMPs_elm[k]=recvMPs_elm[k]-offset; recvMPs_ids[k]=recvMPs_ids[k]-offset; } - + 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); @@ -309,7 +309,7 @@ void polympo_getMPCurElmID_f(MPMesh_ptr p_mpmesh, const int numMPs, int* elmIDs) } }; p_MPs->parallel_for(getElmId, "get mpCurElmID"); - Kokkos::deep_copy( arrayHost, mpCurElmIDCopy); + Kokkos::deep_copy( arrayHost, mpCurElmIDCopy); } void polympo_setMPLatLonRotatedFlag_f(MPMesh_ptr p_mpmesh, const int isRotateFlag){ @@ -335,7 +335,7 @@ void polympo_setMPPositions_f(MPMesh_ptr p_mpmesh, const int nComps, const int n 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){ @@ -526,7 +526,7 @@ void polympo_setMPMass_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs int self; MPI_Comm comm = p_MPs->getMPIComm(); MPI_Comm_rank(comm, &self); - + PMT_ALWAYS_ASSERT(nComps == 1); //TODO mp_sclr_t PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount()); //PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID()); @@ -616,7 +616,7 @@ void polympo_getMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, void polympo_calculateMPStrainRate_f(MPMesh_ptr p_mpmesh){ checkMPMeshValid(p_mpmesh); auto mpMesh = ((polyMPO::MPMesh*)p_mpmesh); - mpMesh->calculateStrain(); + mpMesh->calculateStrain(); } void polympo_setMPStrainRate_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpStrainRateIn){ @@ -674,7 +674,7 @@ void polympo_getMPStrainRate_f(MPMesh_ptr p_mpmesh, const int nComps, const int void polympo_calculateMPStress_f(MPMesh_ptr p_mpmesh, const int constitutive_model){ checkMPMeshValid(p_mpmesh); auto mpMesh = ((polyMPO::MPMesh*)p_mpmesh); - mpMesh->calculateStress(constitutive_model); + mpMesh->calculateStress(constitutive_model); } void polympo_setMPStress_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpStressIn) { @@ -786,7 +786,7 @@ void polympo_setIcePressureMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int void polympo_getReplacementPressureMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* replacementPressureMPHost){ Kokkos::Timer timer; checkMPMeshValid(p_mpmesh); - + auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; //Rank information int self; @@ -816,14 +816,14 @@ void polympo_getReplacementPressureMP_f(MPMesh_ptr p_mpmesh, const int nComps, c void polympo_startMeshFill_f(MPMesh_ptr p_mpmesh){ checkMPMeshValid(p_mpmesh); - ((polyMPO::MPMesh*)p_mpmesh)->p_mesh->setMeshEdit(true); + ((polyMPO::MPMesh*)p_mpmesh)->p_mesh->setMeshEdit(true); } void polympo_endMeshFill_f(MPMesh_ptr p_mpmesh){ checkMPMeshValid(p_mpmesh); - auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; + auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; PMT_ALWAYS_ASSERT(p_mesh->meshEditable()); - p_mesh->setMeshEdit(false); + p_mesh->setMeshEdit(false); } void polympo_checkMeshMaxSettings_f(MPMesh_ptr p_mpmesh, const int maxEdges, const int vertexDegree){ @@ -868,7 +868,7 @@ void polympo_setMeshNumVtxs_f(MPMesh_ptr p_mpmesh, const int numVtxs){ checkMPMeshValid(p_mpmesh); auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; p_mesh->setNumVtxs(numVtxs); - p_mesh->setMeshVtxBasedFieldSize(); + p_mesh->setMeshVtxBasedFieldSize(); } void polympo_setMeshNumVtxsOwned_f(MPMesh_ptr p_mpmesh, const int numVtxsOwned){ @@ -887,8 +887,8 @@ void polympo_setMeshNumElms_f(MPMesh_ptr p_mpmesh, const int numElms){ checkMPMeshValid(p_mpmesh); auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; - auto elm2Vtx = polyMPO::IntVtx2ElmView("MeshElementsToVertices",numElms); - auto elm2Elm = polyMPO::IntElm2ElmView("MeshElementsToElements",numElms); + auto elm2Vtx = polyMPO::IntVtx2ElmView("MeshElementsToVertices",numElms); + auto elm2Elm = polyMPO::IntElm2ElmView("MeshElementsToElements",numElms); p_mesh->setNumElms(numElms); p_mesh->setElm2VtxConn(elm2Vtx); @@ -922,8 +922,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){ //chech vailidity checkMPMeshValid(p_mpmesh); - kkViewHostU arrayHost(array,maxEdges,nCells); - auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; + kkViewHostU arrayHost(array,maxEdges,nCells); + auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; PMT_ALWAYS_ASSERT(p_mesh->meshEditable()); //check the size @@ -944,29 +944,29 @@ void polympo_setMeshElm2ElmConn_f(MPMesh_ptr p_mpmesh, const int maxEdges, const //chech vailidity checkMPMeshValid(p_mpmesh); kkViewHostU arrayHost(array,maxEdges,nCells); //Fortran is column-major - auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; + auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; PMT_ALWAYS_ASSERT(p_mesh->meshEditable()); //check the size PMT_ALWAYS_ASSERT(maxEdges <= maxVtxsPerElm); PMT_ALWAYS_ASSERT(nCells == p_mesh->getNumElements()); - + Kokkos::View elm2ElmArray("MeshElementsToVertices",maxEdges,nCells); Kokkos::deep_copy(elm2ElmArray, arrayHost); auto elm2ElmConn = p_mesh->getElm2ElmConn(); Kokkos::parallel_for("set elm2ElmConn", nCells, KOKKOS_LAMBDA(const int elm){ for(int i=0; ip_mesh; + auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; PMT_ALWAYS_ASSERT(p_mesh->meshEditable()); - kkViewHostU arrayHost(array,nCells); + kkViewHostU arrayHost(array,nCells); //check the size PMT_ALWAYS_ASSERT(nCells == p_mesh->getNumElements()); @@ -978,8 +978,8 @@ void polympo_setOwningProc_f(MPMesh_ptr p_mpmesh, const int nCells, const int* a void polympo_setOwningProcVertex_f(MPMesh_ptr p_mpmesh, const int nVertices, const int* array){ checkMPMeshValid(p_mpmesh); - auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; - kkViewHostU arrayHost(array,nVertices); + auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; + kkViewHostU arrayHost(array,nVertices); //check the size PMT_ALWAYS_ASSERT(nVertices == p_mesh->getNumVertices()); @@ -991,7 +991,7 @@ void polympo_setOwningProcVertex_f(MPMesh_ptr p_mpmesh, const int nVertices, con 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; + auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; 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 @@ -1006,7 +1006,7 @@ void polympo_setElmGlobal_f(MPMesh_ptr p_mpmesh, const int nCells, const int* ar void polympo_setVtxGlobal_f(MPMesh_ptr p_mpmesh, const int nVertices, const int* array){ checkMPMeshValid(p_mpmesh); - auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; + auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; PMT_ALWAYS_ASSERT(nVertices==p_mesh->getNumVertices()); Kokkos::View arrayHost("arrayHost", nVertices); for (int i = 0; i < nVertices; i++) { @@ -1020,7 +1020,7 @@ void polympo_setVtxGlobal_f(MPMesh_ptr p_mpmesh, const int nVertices, const int* void polympo_setInteriorVertex_f(MPMesh_ptr p_mpmesh, const int nVertices, const int* array){ checkMPMeshValid(p_mpmesh); - auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; + auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; PMT_ALWAYS_ASSERT(nVertices==p_mesh->getNumVertices()); Kokkos::View arrayHost("arrayHost", nVertices); for (int i = 0; i < nVertices; i++) { @@ -1048,7 +1048,7 @@ void polympo_setMeshVtxCoords_f(MPMesh_ptr p_mpmesh, const int nVertices, const auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; //check the size - PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); + PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); //copy the host array to the device auto coordsArray = p_mesh->getMeshField(); @@ -1068,9 +1068,9 @@ void polympo_getMeshVtxCoords_f(MPMesh_ptr p_mpmesh, const int nVertices, double auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; //check the size - PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); - - //copy the device to host + PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); + + //copy the device to host auto coordsArray = p_mesh->getMeshField(); auto h_coordsArray = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), coordsArray); @@ -1088,7 +1088,7 @@ void polympo_setMeshVtxRotLat_f(MPMesh_ptr p_mpmesh, const int nVertices, const auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; //check the size - PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); + PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); //copy the host array to the device auto coordsArray = p_mesh->getMeshField(); @@ -1106,9 +1106,9 @@ void polympo_getMeshVtxRotLat_f(MPMesh_ptr p_mpmesh, const int nVertices, double auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; //check the size - PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); - - //copy the device to host + PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); + + //copy the device to host auto coordsArray = p_mesh->getMeshField(); auto h_coordsArray = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), coordsArray); @@ -1123,7 +1123,7 @@ void polympo_setMeshVtxVel_f(MPMesh_ptr p_mpmesh, const int nVertices, const dou auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; //check the size - PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); + PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); //copy the host array to the device auto coordsArray = p_mesh->getMeshField(); @@ -1142,7 +1142,7 @@ void polympo_getMeshVtxVel_f(MPMesh_ptr p_mpmesh, const int nVertices, double* u auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; //check the size - PMT_ALWAYS_ASSERT(p_mesh->getNumVertices() == nVertices); + PMT_ALWAYS_ASSERT(p_mesh->getNumVertices() == nVertices); //copy the device array to the host auto coordsArray = p_mesh->getMeshField(); @@ -1160,7 +1160,7 @@ void polympo_setMeshVtxMass_f(MPMesh_ptr p_mpmesh, const int nVertices, const do auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; //check the size - PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); + PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); //copy the host array to the device auto coordsArray = p_mesh->getMeshField(); @@ -1180,9 +1180,9 @@ void polympo_getMeshVtxMass_f(MPMesh_ptr p_mpmesh, const int nVertices, double* int self; MPI_Comm comm = p_MPs->getMPIComm(); MPI_Comm_rank(comm, &self); - + //check the size - PMT_ALWAYS_ASSERT(p_mesh->getNumVertices() == nVertices); + PMT_ALWAYS_ASSERT(p_mesh->getNumVertices() == nVertices); //copy the device array to the host auto coordsArray = p_mesh->getMeshField(); @@ -1199,7 +1199,7 @@ void polympo_setMeshElmMass_f(MPMesh_ptr p_mpmesh, const int nCells, const doubl auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; //check the size - PMT_ALWAYS_ASSERT(p_mesh->getNumElements()==nCells); + PMT_ALWAYS_ASSERT(p_mesh->getNumElements()==nCells); //copy the host array to the device auto coordsArray = p_mesh->getMeshField(); @@ -1218,7 +1218,7 @@ void polympo_getMeshElmMass_f(MPMesh_ptr p_mpmesh, const int nCells, double* elm auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; //check the size - PMT_ALWAYS_ASSERT(p_mesh->getNumElements() == nCells); + PMT_ALWAYS_ASSERT(p_mesh->getNumElements() == nCells); //copy the device array to the host auto coordsArray = p_mesh->getMeshField(); @@ -1231,7 +1231,7 @@ void polympo_getMeshElmMass_f(MPMesh_ptr p_mpmesh, const int nCells, double* elm //Increments in vertex velcoity and displacement void polympo_setMeshVtxOnSurfVeloIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, const double* array) { - + Kokkos::Timer timer; //check mpMesh is valid checkMPMeshValid(p_mpmesh); @@ -1265,7 +1265,7 @@ void polympo_getMeshVtxOnSurfVeloIncr_f(MPMesh_ptr p_mpmesh, const int nComps, c //check the size PMT_ALWAYS_ASSERT(nComps == vec2d_nEntries); - PMT_ALWAYS_ASSERT(p_mesh->getNumVertices() == nVertices); + PMT_ALWAYS_ASSERT(p_mesh->getNumVertices() == nVertices); PMT_ALWAYS_ASSERT(static_cast(nVertices*vec2d_nEntries)==vtxField.size()); //copy the device array to the host @@ -1310,7 +1310,7 @@ void polympo_getMeshVtxOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, c //check the size PMT_ALWAYS_ASSERT(nComps == vec2d_nEntries); - PMT_ALWAYS_ASSERT(p_mesh->getNumVertices() == nVertices); + PMT_ALWAYS_ASSERT(p_mesh->getNumVertices() == nVertices); PMT_ALWAYS_ASSERT(static_cast(nVertices*vec2d_nEntries)==vtxField.size()); //copy the device array to the host @@ -1327,7 +1327,7 @@ void polympo_setMeshElmCenter_f(MPMesh_ptr p_mpmesh, const int nCells, const dou auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; //check the size - PMT_ALWAYS_ASSERT(p_mesh->getNumElements()==nCells); + PMT_ALWAYS_ASSERT(p_mesh->getNumElements()==nCells); //copy the host array to the device auto elmCenter = p_mesh->getMeshField(); @@ -1346,9 +1346,9 @@ void polympo_getMeshElmCenter_f(MPMesh_ptr p_mpmesh, const int nCells, double* x auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; //check the size - PMT_ALWAYS_ASSERT(p_mesh->getNumElements()==nCells); + PMT_ALWAYS_ASSERT(p_mesh->getNumElements()==nCells); - //copy the device to host + //copy the device to host auto elmCenter = p_mesh->getMeshField(); auto h_elmCenter = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), elmCenter); for(int i=0; ip_mesh; PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); - //copy the device to host + //copy the device to host auto dualArea = p_mesh->getMeshField(); auto h_dualArea = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), dualArea); for(int i=0; ip_mesh; //check the size - PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); + PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); - //copy the device to host + //copy the device to host auto stressDivergence = p_mesh->getMeshField(); auto h_stressDivergence = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), stressDivergence); for(int i=0; ip_MPs; + MPI_Comm comm = p_MPs->getMPIComm(); + MPI_Comm_size(comm, &numProcsTot); + if(numProcsTot == 1) return; + auto mpMesh = ((polyMPO::MPMesh*)p_mpmesh); auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; int numVertices = p_mesh->getNumVertices(); auto vtxFieldVel = p_mesh->getMeshField(); - mpMesh->communicate_and_take_halo_contributions1(vtxFieldVel, numVertices, 2, 1, 1); + mpMesh->communicate_and_take_halo_contributions_gpu_aware(vtxFieldVel, numVertices, 2, 1, 1); } //Advection Calcualtions @@ -1727,7 +1734,7 @@ void polympo_aggregate_deluDyn_f(MPMesh_ptr p_mpmesh){ } void polympo_finalize_deludelvDyn_f(MPMesh_ptr p_mpmesh){ - + checkMPMeshValid(p_mpmesh); auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; @@ -1738,7 +1745,7 @@ void polympo_finalize_deludelvDyn_f(MPMesh_ptr p_mpmesh){ auto vtxField = p_mesh->getMeshField(); auto vtxFieldVel = p_mesh->getMeshField(); auto vtxFieldVel_incr = p_mesh->getMeshField(); - + Kokkos::parallel_for("Finalize_increments", nVertices, KOKKOS_LAMBDA(const int vtx){ vtxField(vtx, 0) = vtxField(vtx, 0) * elasticTimeStep; vtxField(vtx, 1) = vtxField(vtx, 1) * elasticTimeStep; diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 3183cdc5..b4af325a 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -61,6 +61,7 @@ pmpo_add_exe(testReconstruction testReconstruction.cpp) pmpo_add_exe(testTracking testTracking.cpp) pmpo_add_exe(testRebuild testRebuild.cpp) pmpo_add_exe(timeAssmblyWachspress testTiming.cpp) +pmpo_add_exe(testGPUAwareHaloExchange testGPUAwareHaloExchange.cpp) pmpo_add_fortran_exe(testFortranInit testFortranInit.f90) pmpo_add_fortran_exe(testFortran testFortran.f90) pmpo_add_fortran_exe(testFortranReadMPAS testFortranReadMPAS.f90) @@ -91,6 +92,7 @@ pmpo_add_test(testFortranInit 1 ./testFortranInit) pmpo_add_test(testFortran 1 ./testFortran) pmpo_add_test(testFortranMPMeshModule 1 ./testFortranMPMeshModule) pmpo_add_test(testFortranMPAppIDs 1 ./testFortranMPAppIDs) +pmpo_add_test(test_gpu_aware_halo_exchange 4 ./testGPUAwareHaloExchange) #set NC file for test set(TEST_NC_FILE_PLANAR "${CMAKE_SOURCE_DIR}/test/sample_mpas_meshes/planar_nonuniform_cvt_for_square_673elms.nc") diff --git a/test/testGPUAwareHaloExchange.cpp b/test/testGPUAwareHaloExchange.cpp new file mode 100644 index 00000000..9c4902fa --- /dev/null +++ b/test/testGPUAwareHaloExchange.cpp @@ -0,0 +1,284 @@ +#include "pmpo_MPMesh.hpp" +#include "pmpo_createTestMPMesh.hpp" + +#include +#include + +#include + +int main(int argc, char** argv) +{ + MPI_Init(&argc, &argv); + Kokkos::initialize(argc, argv); + + int testResult = 0; + + { + int rank = -1; + int size = -1; + + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + MPI_Comm_size(MPI_COMM_WORLD, &size); + + if (rank == 0) { + std::cout + << "GPU-aware halo exchange test running with " + << size << " MPI ranks." + << std::endl; + } + + // This test uses the following communication pattern: + // + // Rank 0 ----\ + // Rank 1 ----- > Rank 3 + // Rank 2 ----/ + // | + // v + // Rank 0 + // + // Therefore: + // Rank 0 has one halo owned by Rank 3. + // Rank 1 has no halos. + // Rank 2 has no halos. + // Rank 3 has three halos owned by Ranks 0, 1, and 2. + + if (size != 4) { + if (rank == 0) { + std::cerr + << "This test requires exactly 4 MPI ranks." + << std::endl; + } + + MPI_Abort(MPI_COMM_WORLD, 1); + } + + // Create the existing polyMPO test mesh and MPMesh. + + polyMPO::Mesh* mesh = polyMPO::initTestMesh(1, 1); + + polyMPO::MPMesh mpMesh = polyMPO::initTestMPMesh(mesh, 1); + + // initTestMPMesh() creates MaterialPoints, but the test helper + // does not initialize its MPI communicator. + mpMesh.p_MPs->setMPIComm(MPI_COMM_WORLD); + + const int nVertices = mesh->getNumVertices(); + + // startCommunication() expects local vertices to be ordered: + + polyMPO::IntView owningProcVertex("testOwningProcVertex", nVertices); + + polyMPO::IntView globalVtx("testGlobalVtx",nVertices); + + auto owningProcHost = Kokkos::create_mirror_view(owningProcVertex); + + auto globalVtxHost = Kokkos::create_mirror_view(globalVtx); + + // Number of locally-owned vertices differs by rank. + // Rank 0: 18 owners + 1 halo + // Rank 1: 19 owners + 0 halos + // Rank 2: 19 owners + 0 halos + // Rank 3: 16 owners + 3 halos + + int numLocalOwners = nVertices; + + if (rank == 0) { + numLocalOwners = nVertices - 1; + } + else if (rank == 3) { + numLocalOwners = nVertices - 3; + } + + // Assign locally-owned vertices. + // Rank 0 global IDs: 0, 1, 2, ... + // Rank 1 global IDs: 100, 101, 102, ... + // Rank 2 global IDs: 200, 201, 202, ... + // Rank 3 global IDs: 300, 301, 302, ... + + for (int i = 0; i < numLocalOwners; ++i) { + owningProcHost(i) = rank; + globalVtxHost(i) = rank * 100 + i; + } + + // Define halo vertices. + + if (rank == 0) { + + // Rank 0 receives Rank 3 owner vertex 0. + // Rank 3 owner vertex 0 has global ID 300. + + owningProcHost(nVertices - 1) = 3; + globalVtxHost(nVertices - 1) = 300; + } + else if (rank == 3) { + + // Rank 3 receives Rank 0 owner vertex 0. + owningProcHost(nVertices - 3) = 0; + globalVtxHost(nVertices - 3) = 0; + + // Rank 3 receives Rank 1 owner vertex 0. + owningProcHost(nVertices - 2) = 1; + globalVtxHost(nVertices - 2) = 100; + + // Rank 3 receives Rank 2 owner vertex 0. + owningProcHost(nVertices - 1) = 2; + globalVtxHost(nVertices - 1) = 200; + } + + Kokkos::deep_copy(owningProcVertex, owningProcHost); + + Kokkos::deep_copy(globalVtx, globalVtxHost); + + mesh->setOwningProcVertex(owningProcVertex); + mesh->setVtxGlobal(globalVtx); + + // Build owner/halo communication metadata + mpMesh.startCommunication(); + + std::cout + << "Rank " << rank + << ": owners = " << mpMesh.numOwnersTot + << ", halos = " << mpMesh.numHalosTot + << std::endl; + + // Create one scalar field value per vertex. + + const int nEntities = mpMesh.numOwnersTot + mpMesh.numHalosTot; + + const int numEntries = 1; + + Kokkos::View field("gpuAwareHaloTestField", nEntities, numEntries); + + const int numOwners = mpMesh.numOwnersTot; + + // Owner values are deterministic: + // Rank 0 owner 0 = 0 + // Rank 1 owner 0 = 1000 + // Rank 2 owner 0 = 2000 + // Rank 3 owner 0 = 3000 + // All halos start at -1. + + Kokkos::parallel_for("initialize_gpu_aware_halo_test", nEntities,KOKKOS_LAMBDA(const int i) + { + if (i < numOwners) { + field(i, 0) = + 1000.0 * static_cast(rank) + + static_cast(i); + } + else { + field(i, 0) = -1.0; + } + }); + + Kokkos::fence(); + + // Perform owner -> halo assignment. + // mode = 1 : owner -> halo + // op = 1 : assignment + + const int mode = 1; + const int op = 1; + + mpMesh.communicate_and_take_halo_contributions_gpu_aware(field, nEntities, numEntries, mode, op); + + Kokkos::fence(); + + // Copy results to host and verify exact halo values. + + auto fieldHost = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), field); + + int localFailures = 0; + + if (rank == 0) { + + // Rank 0 receives Rank 3 owner vertex 0. + const int haloIndex = nVertices - 1; + const double expectedValue = 3000.0; + + if (fieldHost(haloIndex, 0) != expectedValue) { + + ++localFailures; + + std::cerr + << "Rank 0: halo value = " + << fieldHost(haloIndex, 0) + << ", expected = " + << expectedValue + << std::endl; + } + } + else if (rank == 3) { + + // Rank 3 receives owner vertex 0 from Ranks 0, 1, and 2. + + const int haloFromRank0 = nVertices - 3; + const int haloFromRank1 = nVertices - 2; + const int haloFromRank2 = nVertices - 1; + + if (fieldHost(haloFromRank0, 0) != 0.0) { + + ++localFailures; + + std::cerr + << "Rank 3: halo from Rank 0 = " + << fieldHost(haloFromRank0, 0) + << ", expected = 0" + << std::endl; + } + + if (fieldHost(haloFromRank1, 0) != 1000.0) { + + ++localFailures; + + std::cerr + << "Rank 3: halo from Rank 1 = " + << fieldHost(haloFromRank1, 0) + << ", expected = 1000" + << std::endl; + } + + if (fieldHost(haloFromRank2, 0) != 2000.0) { + + ++localFailures; + + std::cerr + << "Rank 3: halo from Rank 2 = " + << fieldHost(haloFromRank2, 0) + << ", expected = 2000" + << std::endl; + } + } + + int globalFailures = 0; + + MPI_Allreduce(&localFailures, &globalFailures, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); + + if (rank == 0) { + + if (globalFailures == 0) { + std::cout + << "GPU-aware halo exchange test PASSED." + << std::endl; + } + else { + std::cerr + << "GPU-aware halo exchange test FAILED with " + << globalFailures + << " halo errors." + << std::endl; + } + } + + if (globalFailures != 0) { + testResult = 1; + } + + // Do not delete mesh here. + // mpMesh owns the Mesh and MaterialPoints objects. + } + + Kokkos::finalize(); + MPI_Finalize(); + + return testResult; +} From 7760bae0294df26e31672d0e9991ea3dae364ff3 Mon Sep 17 00:00:00 2001 From: Shahrear Jahan Santho Date: Thu, 24 Sep 2026 22:27:19 -0700 Subject: [PATCH 03/12] Clean up unrelated formatting changes --- CMakeLists.txt | 2 +- src/pmpo_MPMesh.cpp | 5 ++ src/pmpo_MPMesh.hpp | 27 ++++----- src/pmpo_MPMesh_assembly.hpp | 34 ++++++++---- src/pmpo_c.cpp | 104 +++++++++++++++++------------------ 5 files changed, 91 insertions(+), 81 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 3f627daa..080a4186 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -50,4 +50,4 @@ if(IS_TESTING) add_subdirectory (test) endif() -bob_end_package() +bob_end_package() \ No newline at end of file diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index fbb16485..8eb41183 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -89,6 +89,7 @@ void MPMesh::calculateStress(const int constitutive_relation){ void MPMesh::calculateStressDivergence(){ + Kokkos::Timer timer; int self, numProcsTot; MPI_Comm comm = p_MPs->getMPIComm(); MPI_Comm_rank(comm, &self); @@ -165,6 +166,9 @@ void MPMesh::calculateStressDivergence(){ }; p_MPs->parallel_for(stress_div, " stress_div_assembly"); Kokkos::fence(); + pumipic::RecordTime("Stress_Divergence_Reconstruction" + std::to_string(self), timer.seconds()); + + timer.reset(); if(numProcsTot>1){ //Takes contribution of halo vertices and adds it in owner procs communicate_and_take_halo_contributions_gpu_aware(stress_divUV, numVertices, 2, 0, 0); @@ -172,6 +176,7 @@ void MPMesh::calculateStressDivergence(){ //communicate_and_take_halo_contributions(stress_divUV, numVertices, 2, 1, 1); } Kokkos::fence(); + pumipic::RecordTime("Stress_Divergence Communication" + std::to_string(self), timer.seconds()); } void MPMesh::calcBasis() { diff --git a/src/pmpo_MPMesh.hpp b/src/pmpo_MPMesh.hpp index ba75b56d..62411d0f 100644 --- a/src/pmpo_MPMesh.hpp +++ b/src/pmpo_MPMesh.hpp @@ -42,8 +42,7 @@ class MPMesh{ std::vector> ownerHaloLocalIDs; void startCommunication(); - void communicate_and_take_halo_contributions(const Kokkos::View& meshField, int nEntities, - int numEntries, int mode, int op); + void communicate_and_take_halo_contributions(const Kokkos::View& meshField, int nEntities, int numEntries, int mode, int op); template void communicate_and_take_halo_contributions_staged( @@ -65,7 +64,6 @@ class MPMesh{ std::vector> recvIDVec; std::vector> recvDataVec; pumipic::RecordTime("SD: Recv Vec Allocation-" + std::to_string(self), timer.seconds()); - timer.reset(); //communicateFieldsFromHostView(fieldData1, nEntities, numEntries, mode, recvIDVec, recvDataVec); communicateFieldsFromHostView(reconVals_host, nEntities, numEntries, mode, recvIDVec, recvDataVec); @@ -75,7 +73,7 @@ class MPMesh{ int numProcsTot = recvIDVec.size(); //Flatten IDs int totalSize = 0; - std::vector offsets(numProcsTot, 0); + std::vector offsets(numProcsTot, 0); for(int i=0; i recvDataGPU("recvDataGPU", totalSize_data); auto hostView_data= Kokkos::View("recvDataCPU", totalSize_data); - std::copy(flatDataVec.begin(), flatDataVec.end(), hostView_data.data()); + std::copy(flatDataVec.begin(), flatDataVec.end(), hostView_data.data()); Kokkos::deep_copy(recvDataGPU, hostView_data); Kokkos::fence(); assert(totalSize_data == totalSize*numEntries); @@ -119,7 +117,6 @@ class MPMesh{ assert(recvDataVec[i].size() == recvIDVec[i].size() * numEntries); } pumipic::RecordTime("SD: Copy CPU-GPU2-" + std::to_string(self), timer.seconds()); - //Take contributions from other procs timer.reset(); Kokkos::parallel_for("halo contribution", recvIDGPU.size(), KOKKOS_LAMBDA(const int i){ @@ -133,9 +130,8 @@ class MPMesh{ pumipic::RecordTime("SD: Contribution" + std::to_string(self), timer.seconds()); } - void communicateFields(const std::vector>& fieldData, const int numEntities, - const int numEntries, int mode, std::vector>& recvIDVec, - std::vector>& recvDataVec); + void communicateFields(const std::vector>& fieldData, const int numEntities, const int numEntries, int mode, + std::vector>& recvIDVec, std::vector>& recvDataVec); template void communicateFieldsFromHostView( @@ -143,7 +139,7 @@ class MPMesh{ const int numEntities, const int numEntries, int mode, std::vector>& recvIDVec, std::vector>& recvDataVec){ - + int self, numProcsTot; MPI_Comm comm = p_MPs->getMPIComm(); MPI_Comm_rank(comm, &self); @@ -159,13 +155,13 @@ class MPMesh{ for(int i = 0; i < numProcsTot; i++){ if(i==self) continue; - int numToSend = 0, numToRecv = 0; + int numToSend = 0, numToRecv = 0; if(mode == 0) { //gather (halos send to owners) numToSend = numOwnersOnOtherProcs[i]; numToRecv = numHalosOnOtherProcs[i]; } - else{ + else{ //scatter (owners send to halos) numToSend = numHalosOnOtherProcs[i]; numToRecv = numOwnersOnOtherProcs[i]; @@ -201,7 +197,7 @@ class MPMesh{ std::vector requests; requests.reserve(4*numProcsTot); for(int proc = 0; proc < numProcsTot; proc++){ - if(proc == self) continue; + if(proc == self) continue; if(mode == 0 && numHalosOnOtherProcs[proc]){ assert(recvIDVec[proc].size() == (size_t)numHalosOnOtherProcs[proc]); assert(recvDataVec[proc].size() == recvIDVec[proc].size() * (size_t)numEntries); @@ -238,7 +234,7 @@ class MPMesh{ } MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE); } - + MPMesh(Mesh* inMesh, MaterialPoints* inMPs): p_mesh(inMesh), p_MPs(inMPs) { }; @@ -273,7 +269,7 @@ class MPMesh{ void reconstruct_coeff_full(); void invertMatrix(const Kokkos::View& vtxMatrices, const double& radius); Kokkos::View precomputedVtxCoeffs_new; - Kokkos::View nearAnEdge; + Kokkos::View nearAnEdge; Kokkos::View vtxMatrixMass; //Not used currently @@ -291,7 +287,6 @@ class MPMesh{ void printVTP_mesh(int printVTPIndex); void writeMPTrackingVTP(int printVTPIndex, int numMPs, const Vec3dView& history, const Vec3dView& resultLeft, const Vec3dView& resultRight, const Vec3dView& mpTgtPosArray); - void calculateStrain(); void calculateStress(const int constitutive_relation); void calculateStressDivergence(); diff --git a/src/pmpo_MPMesh_assembly.hpp b/src/pmpo_MPMesh_assembly.hpp index 61609341..c4769aac 100644 --- a/src/pmpo_MPMesh_assembly.hpp +++ b/src/pmpo_MPMesh_assembly.hpp @@ -85,7 +85,7 @@ void MPMesh::assemblyElm0() { } }; p_MPs->parallel_for(assemble, "assembly"); - + Kokkos::MDRangePolicy> policy({0,0},{numElms, numEntries}); Kokkos::parallel_for("assembly average", policy, KOKKOS_LAMBDA(const int elm, const int entry){ if (mpsPerElm(elm) > 0){ @@ -96,11 +96,12 @@ void MPMesh::assemblyElm0() { } void MPMesh::reconstruct_coeff_full(){ + Kokkos::Timer timer; int self, numProcsTot; MPI_Comm comm = p_MPs->getMPIComm(); MPI_Comm_rank(comm, &self); MPI_Comm_size(comm, &numProcsTot); - + static int coeff_count=0; if(!self) std::cout<<"===="<<__FUNCTION__<<" "<parallel_for(assemble, "assembly"); Kokkos::fence(); + pumipic::RecordTime("Assemble Matrix Per Process" + std::to_string(self), timer.seconds()); //Mode 0 is Gather: Halos Send to Owners //Mode 1 is Scatter: Owners Send to Halos //Op 0 is addition //Op 1 is replacement + timer.reset(); int mode = 0; int op = 0; - if (numProcsTot>1){ + if (numProcsTot >1){ communicate_and_take_halo_contributions_gpu_aware(vtxMatrices, numVertices, numEntriesMatrix, mode, op); - mode=1; + mode=1; op=1; communicate_and_take_halo_contributions_gpu_aware(vtxMatrices, numVertices, numEntriesMatrix, mode, op); } + pumipic::RecordTime("Communicate Matrix Values" + std::to_string(self), timer.seconds()); + //Store the 1st matrix element Kokkos::ViewvtxMatrixMass_l("vtxMass", numVertices); Kokkos::parallel_for("storeMatrixMass", numVertices, KOKKOS_LAMBDA(const int vtx){ vtxMatrixMass_l(vtx) = vtxMatrices(vtx, 0); }); - Kokkos::fence(); this->vtxMatrixMass = vtxMatrixMass_l; invertMatrix(vtxMatrices, radius); } void MPMesh::invertMatrix(const Kokkos::View& vtxMatrices, const double& radius){ + int nVertices = p_mesh->getNumVertices(); auto vtxCoords = p_mesh->getMeshField(); auto dual_triangle_area = p_mesh->getMeshField(); @@ -187,7 +192,7 @@ void MPMesh::invertMatrix(const Kokkos::View& vtxMatrices, const doubl Kokkos::View VtxCoeffs("VtxCoeffs", nVertices); Kokkos::deep_copy(VtxCoeffs, 0.0); Kokkos::View nearAnEdge_l("nearAnEdge_l", nVertices); - + Kokkos::parallel_for("invertMatrix", nVertices, KOKKOS_LAMBDA(const int vtx){ if(vtxMatrices(vtx, 0) < eps) return; @@ -204,7 +209,7 @@ void MPMesh::invertMatrix(const Kokkos::View& vtxMatrices, const doubl } auto cosLat = sqrt(pow(X, 2) + pow(Y, 2)); - auto invCosLat = 1.0/cosLat; + auto invCosLat = 1.0/cosLat; auto vtx_area_sqrt = sqrt(dual_triangle_area(vtx,0)/(radius*radius)); Vec3d v0 = { -Y * invCosLat, -Z * X * invCosLat, X / vtx_area_sqrt }; @@ -276,7 +281,7 @@ void MPMesh::invertMatrix(const Kokkos::View& vtxMatrices, const doubl iBlockC[0] = blockC[0] * invM2D[0] + blockC[1] * invM2D[1] + blockC[2] * invM2D[2]; iBlockC[1] = blockC[0] * invM2D[1] + blockC[1] * invM2D[3] + blockC[2] * invM2D[4]; iBlockC[2] = blockC[0] * invM2D[2] + blockC[1] * invM2D[4] + blockC[2] * invM2D[5]; - iBlockC = iBlockC*invM11; + iBlockC = iBlockC*invM11; VtxCoeffs(vtx, 0, 0) = invM11 + invM11*iBlockC.dot(blockC); auto temp = -rotateScaleM.rightMultiply(iBlockC); @@ -288,7 +293,7 @@ void MPMesh::invertMatrix(const Kokkos::View& vtxMatrices, const doubl VtxCoeffs(vtx, 1, 1) = invM2D[0] * rotateScaleM(0, 0) + invM2D[1] * rotateScaleM(0, 1) + invM2D[2] * rotateScaleM(0, 2); VtxCoeffs(vtx, 1, 2) = invM2D[0] * rotateScaleM(1, 0) + invM2D[1] * rotateScaleM(1, 1) + invM2D[2] * rotateScaleM(1, 2); VtxCoeffs(vtx, 1, 3) = invM2D[0] * rotateScaleM(2, 0) + invM2D[1] * rotateScaleM(2, 1) + invM2D[2] * rotateScaleM(2, 2); - + VtxCoeffs(vtx, 2, 0) = -iBlockC[1]; VtxCoeffs(vtx, 2, 1) = invM2D[1] * rotateScaleM(0, 0) + invM2D[3] * rotateScaleM(0, 1) + invM2D[4] * rotateScaleM(0, 2); VtxCoeffs(vtx, 2, 2) = invM2D[1] * rotateScaleM(1, 0) + invM2D[3] * rotateScaleM(1, 1) + invM2D[4] * rotateScaleM(1, 2); @@ -316,6 +321,7 @@ void MPMesh::invertMatrix(const Kokkos::View& vtxMatrices, const doubl template void MPMesh::assemblyVtx1(){ + Kokkos::Timer timer; int self, numProcsTot; MPI_Comm comm = p_MPs->getMPIComm(); @@ -352,7 +358,7 @@ void MPMesh::assemblyVtx1(){ int nVtxE = elm2VtxConn(elm,0); //number of vertices bounding the element for(int i=0; iparallel_for(reconstruct, "reconstruct"); Kokkos::fence(); + pumipic::RecordTime("Assemble Field per process" + std::to_string(self), timer.seconds()); + + timer.reset(); if(numProcsTot>1){ communicate_and_take_halo_contributions_gpu_aware(meshField, numVertices, numEntries, 0, 0); } + pumipic::RecordTime("Communicate Field Values" + std::to_string(self), timer.seconds()); } template @@ -418,7 +428,7 @@ DoubleView MPMesh::wtScaAssembly(){ // last component of eVtxCoords stores the firs vertex (to avoid if-condition in the Wachspress computation) eVtxCoords[nElmVtxs][0] = vtxCoords(elm2VtxConn(elm,1)-1,0); eVtxCoords[nElmVtxs][1] = vtxCoords(elm2VtxConn(elm,1)-1,1); - + /* compute the values of basis functions at mp position */ double basisByArea[maxElmsPerVtx]; Vec2d mpCoord(mpPositions(mp,0), mpPositions(mp,1)); @@ -457,7 +467,7 @@ Vec2dView MPMesh::wtVec2Assembly(){ Vec2d eVtxCoords[maxVtxsPerElm + 1]; for (int i = 1; i <= nElmVtxs; i++) { // elm2VtxConn(elm,i) is the vertex ID (1-based index) of vertex #i of elm - eVtxCoords[i-1][0] = vtxCoords(elm2VtxConn(elm,i)-1,0); + eVtxCoords[i-1][0] = vtxCoords(elm2VtxConn(elm,i)-1,0); eVtxCoords[i-1][1] = vtxCoords(elm2VtxConn(elm,i)-1,1); } // last component of eVtxCoords stores the firs vertex (to avoid if-condition in the Wachspress computation) diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index fd598cc1..385111dd 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -6,7 +6,7 @@ namespace{ std::vector p_mpmeshes;////store the p_mpmeshes that is legal - + void checkMPMeshValid(MPMesh_ptr p_mpmesh){ auto p_mpmeshIter = std::find(p_mpmeshes.begin(),p_mpmeshes.end(),p_mpmesh); PMT_ALWAYS_ASSERT(p_mpmeshIter != p_mpmeshes.end()); @@ -39,7 +39,7 @@ MPMesh_ptr polympo_createMPMesh_f(const int testMeshOption, const int testMPOpti PMT_ALWAYS_ASSERT(testMeshOption >= 1); p_mps = polyMPO::initTestMPs(p_mesh, testMPOption); }else{ - p_mps = new polyMPO::MaterialPoints(); + p_mps = new polyMPO::MaterialPoints(); } MPMesh_ptr p_mpMeshReturn = (MPMesh_ptr) new polyMPO::MPMesh(p_mesh, p_mps); p_mpmeshes.push_back(p_mpMeshReturn); @@ -87,7 +87,7 @@ void polympo_createMPs_f(MPMesh_ptr p_mpmesh, for(int i = 0; i < numMPs; i++) { if(isMPActive[i] == MP_ACTIVE) { numActiveMPs++; - if(mp2Elm[i] < minElmID) + if(mp2Elm[i] < minElmID) minElmID = mp2Elm[i]; } } @@ -138,7 +138,7 @@ void polympo_createMPs_f(MPMesh_ptr p_mpmesh, 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); + p_MPs->setElmIDoffset(offset); } void polympo_startRebuildMPs_f(MPMesh_ptr p_mpmesh, @@ -151,7 +151,7 @@ void polympo_startRebuildMPs_f(MPMesh_ptr p_mpmesh, int self; MPI_Comm comm = p_MPs->getMPIComm(); MPI_Comm_rank(comm, &self); - + PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount()); //PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID()); @@ -225,7 +225,7 @@ void polympo_startRebuildMPs2_f(MPMesh_ptr p_mpmesh, recvMPs_elm[k]=recvMPs_elm[k]-offset; recvMPs_ids[k]=recvMPs_ids[k]-offset; } - + 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); @@ -309,7 +309,7 @@ void polympo_getMPCurElmID_f(MPMesh_ptr p_mpmesh, const int numMPs, int* elmIDs) } }; p_MPs->parallel_for(getElmId, "get mpCurElmID"); - Kokkos::deep_copy( arrayHost, mpCurElmIDCopy); + Kokkos::deep_copy( arrayHost, mpCurElmIDCopy); } void polympo_setMPLatLonRotatedFlag_f(MPMesh_ptr p_mpmesh, const int isRotateFlag){ @@ -335,7 +335,7 @@ void polympo_setMPPositions_f(MPMesh_ptr p_mpmesh, const int nComps, const int n 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){ @@ -526,7 +526,7 @@ void polympo_setMPMass_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs int self; MPI_Comm comm = p_MPs->getMPIComm(); MPI_Comm_rank(comm, &self); - + PMT_ALWAYS_ASSERT(nComps == 1); //TODO mp_sclr_t PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount()); //PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID()); @@ -616,7 +616,7 @@ void polympo_getMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, void polympo_calculateMPStrainRate_f(MPMesh_ptr p_mpmesh){ checkMPMeshValid(p_mpmesh); auto mpMesh = ((polyMPO::MPMesh*)p_mpmesh); - mpMesh->calculateStrain(); + mpMesh->calculateStrain(); } void polympo_setMPStrainRate_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpStrainRateIn){ @@ -674,7 +674,7 @@ void polympo_getMPStrainRate_f(MPMesh_ptr p_mpmesh, const int nComps, const int void polympo_calculateMPStress_f(MPMesh_ptr p_mpmesh, const int constitutive_model){ checkMPMeshValid(p_mpmesh); auto mpMesh = ((polyMPO::MPMesh*)p_mpmesh); - mpMesh->calculateStress(constitutive_model); + mpMesh->calculateStress(constitutive_model); } void polympo_setMPStress_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpStressIn) { @@ -786,7 +786,7 @@ void polympo_setIcePressureMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int void polympo_getReplacementPressureMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* replacementPressureMPHost){ Kokkos::Timer timer; checkMPMeshValid(p_mpmesh); - + auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; //Rank information int self; @@ -816,14 +816,14 @@ void polympo_getReplacementPressureMP_f(MPMesh_ptr p_mpmesh, const int nComps, c void polympo_startMeshFill_f(MPMesh_ptr p_mpmesh){ checkMPMeshValid(p_mpmesh); - ((polyMPO::MPMesh*)p_mpmesh)->p_mesh->setMeshEdit(true); + ((polyMPO::MPMesh*)p_mpmesh)->p_mesh->setMeshEdit(true); } void polympo_endMeshFill_f(MPMesh_ptr p_mpmesh){ checkMPMeshValid(p_mpmesh); - auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; + auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; PMT_ALWAYS_ASSERT(p_mesh->meshEditable()); - p_mesh->setMeshEdit(false); + p_mesh->setMeshEdit(false); } void polympo_checkMeshMaxSettings_f(MPMesh_ptr p_mpmesh, const int maxEdges, const int vertexDegree){ @@ -868,7 +868,7 @@ void polympo_setMeshNumVtxs_f(MPMesh_ptr p_mpmesh, const int numVtxs){ checkMPMeshValid(p_mpmesh); auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; p_mesh->setNumVtxs(numVtxs); - p_mesh->setMeshVtxBasedFieldSize(); + p_mesh->setMeshVtxBasedFieldSize(); } void polympo_setMeshNumVtxsOwned_f(MPMesh_ptr p_mpmesh, const int numVtxsOwned){ @@ -887,8 +887,8 @@ void polympo_setMeshNumElms_f(MPMesh_ptr p_mpmesh, const int numElms){ checkMPMeshValid(p_mpmesh); auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; - auto elm2Vtx = polyMPO::IntVtx2ElmView("MeshElementsToVertices",numElms); - auto elm2Elm = polyMPO::IntElm2ElmView("MeshElementsToElements",numElms); + auto elm2Vtx = polyMPO::IntVtx2ElmView("MeshElementsToVertices",numElms); + auto elm2Elm = polyMPO::IntElm2ElmView("MeshElementsToElements",numElms); p_mesh->setNumElms(numElms); p_mesh->setElm2VtxConn(elm2Vtx); @@ -950,23 +950,23 @@ void polympo_setMeshElm2ElmConn_f(MPMesh_ptr p_mpmesh, const int maxEdges, const //check the size PMT_ALWAYS_ASSERT(maxEdges <= maxVtxsPerElm); PMT_ALWAYS_ASSERT(nCells == p_mesh->getNumElements()); - + Kokkos::View elm2ElmArray("MeshElementsToVertices",maxEdges,nCells); Kokkos::deep_copy(elm2ElmArray, arrayHost); auto elm2ElmConn = p_mesh->getElm2ElmConn(); Kokkos::parallel_for("set elm2ElmConn", nCells, KOKKOS_LAMBDA(const int elm){ for(int i=0; ip_mesh; + auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; PMT_ALWAYS_ASSERT(p_mesh->meshEditable()); - kkViewHostU arrayHost(array,nCells); + kkViewHostU arrayHost(array,nCells); //check the size PMT_ALWAYS_ASSERT(nCells == p_mesh->getNumElements()); @@ -978,8 +978,8 @@ void polympo_setOwningProc_f(MPMesh_ptr p_mpmesh, const int nCells, const int* a void polympo_setOwningProcVertex_f(MPMesh_ptr p_mpmesh, const int nVertices, const int* array){ checkMPMeshValid(p_mpmesh); - auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; - kkViewHostU arrayHost(array,nVertices); + auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; + kkViewHostU arrayHost(array,nVertices); //check the size PMT_ALWAYS_ASSERT(nVertices == p_mesh->getNumVertices()); @@ -991,7 +991,7 @@ void polympo_setOwningProcVertex_f(MPMesh_ptr p_mpmesh, const int nVertices, con 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; + auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; 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 @@ -1006,7 +1006,7 @@ void polympo_setElmGlobal_f(MPMesh_ptr p_mpmesh, const int nCells, const int* ar void polympo_setVtxGlobal_f(MPMesh_ptr p_mpmesh, const int nVertices, const int* array){ checkMPMeshValid(p_mpmesh); - auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; + auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; PMT_ALWAYS_ASSERT(nVertices==p_mesh->getNumVertices()); Kokkos::View arrayHost("arrayHost", nVertices); for (int i = 0; i < nVertices; i++) { @@ -1020,7 +1020,7 @@ void polympo_setVtxGlobal_f(MPMesh_ptr p_mpmesh, const int nVertices, const int* void polympo_setInteriorVertex_f(MPMesh_ptr p_mpmesh, const int nVertices, const int* array){ checkMPMeshValid(p_mpmesh); - auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; + auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; PMT_ALWAYS_ASSERT(nVertices==p_mesh->getNumVertices()); Kokkos::View arrayHost("arrayHost", nVertices); for (int i = 0; i < nVertices; i++) { @@ -1048,7 +1048,7 @@ void polympo_setMeshVtxCoords_f(MPMesh_ptr p_mpmesh, const int nVertices, const auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; //check the size - PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); + PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); //copy the host array to the device auto coordsArray = p_mesh->getMeshField(); @@ -1069,8 +1069,8 @@ void polympo_getMeshVtxCoords_f(MPMesh_ptr p_mpmesh, const int nVertices, double //check the size PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); - - //copy the device to host + + //copy the device to host auto coordsArray = p_mesh->getMeshField(); auto h_coordsArray = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), coordsArray); @@ -1088,7 +1088,7 @@ void polympo_setMeshVtxRotLat_f(MPMesh_ptr p_mpmesh, const int nVertices, const auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; //check the size - PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); + PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); //copy the host array to the device auto coordsArray = p_mesh->getMeshField(); @@ -1106,9 +1106,9 @@ void polympo_getMeshVtxRotLat_f(MPMesh_ptr p_mpmesh, const int nVertices, double auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; //check the size - PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); - - //copy the device to host + PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); + + //copy the device to host auto coordsArray = p_mesh->getMeshField(); auto h_coordsArray = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), coordsArray); @@ -1123,7 +1123,7 @@ void polympo_setMeshVtxVel_f(MPMesh_ptr p_mpmesh, const int nVertices, const dou auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; //check the size - PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); + PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); //copy the host array to the device auto coordsArray = p_mesh->getMeshField(); @@ -1160,7 +1160,7 @@ void polympo_setMeshVtxMass_f(MPMesh_ptr p_mpmesh, const int nVertices, const do auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; //check the size - PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); + PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); //copy the host array to the device auto coordsArray = p_mesh->getMeshField(); @@ -1180,9 +1180,9 @@ void polympo_getMeshVtxMass_f(MPMesh_ptr p_mpmesh, const int nVertices, double* int self; MPI_Comm comm = p_MPs->getMPIComm(); MPI_Comm_rank(comm, &self); - + //check the size - PMT_ALWAYS_ASSERT(p_mesh->getNumVertices() == nVertices); + PMT_ALWAYS_ASSERT(p_mesh->getNumVertices() == nVertices); //copy the device array to the host auto coordsArray = p_mesh->getMeshField(); @@ -1199,7 +1199,7 @@ void polympo_setMeshElmMass_f(MPMesh_ptr p_mpmesh, const int nCells, const doubl auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; //check the size - PMT_ALWAYS_ASSERT(p_mesh->getNumElements()==nCells); + PMT_ALWAYS_ASSERT(p_mesh->getNumElements()==nCells); //copy the host array to the device auto coordsArray = p_mesh->getMeshField(); @@ -1218,7 +1218,7 @@ void polympo_getMeshElmMass_f(MPMesh_ptr p_mpmesh, const int nCells, double* elm auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; //check the size - PMT_ALWAYS_ASSERT(p_mesh->getNumElements() == nCells); + PMT_ALWAYS_ASSERT(p_mesh->getNumElements() == nCells); //copy the device array to the host auto coordsArray = p_mesh->getMeshField(); @@ -1231,7 +1231,7 @@ void polympo_getMeshElmMass_f(MPMesh_ptr p_mpmesh, const int nCells, double* elm //Increments in vertex velcoity and displacement void polympo_setMeshVtxOnSurfVeloIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, const double* array) { - + Kokkos::Timer timer; //check mpMesh is valid checkMPMeshValid(p_mpmesh); @@ -1265,7 +1265,7 @@ void polympo_getMeshVtxOnSurfVeloIncr_f(MPMesh_ptr p_mpmesh, const int nComps, c //check the size PMT_ALWAYS_ASSERT(nComps == vec2d_nEntries); - PMT_ALWAYS_ASSERT(p_mesh->getNumVertices() == nVertices); + PMT_ALWAYS_ASSERT(p_mesh->getNumVertices() == nVertices); PMT_ALWAYS_ASSERT(static_cast(nVertices*vec2d_nEntries)==vtxField.size()); //copy the device array to the host @@ -1310,7 +1310,7 @@ void polympo_getMeshVtxOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, c //check the size PMT_ALWAYS_ASSERT(nComps == vec2d_nEntries); - PMT_ALWAYS_ASSERT(p_mesh->getNumVertices() == nVertices); + PMT_ALWAYS_ASSERT(p_mesh->getNumVertices() == nVertices); PMT_ALWAYS_ASSERT(static_cast(nVertices*vec2d_nEntries)==vtxField.size()); //copy the device array to the host @@ -1327,7 +1327,7 @@ void polympo_setMeshElmCenter_f(MPMesh_ptr p_mpmesh, const int nCells, const dou auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; //check the size - PMT_ALWAYS_ASSERT(p_mesh->getNumElements()==nCells); + PMT_ALWAYS_ASSERT(p_mesh->getNumElements()==nCells); //copy the host array to the device auto elmCenter = p_mesh->getMeshField(); @@ -1346,9 +1346,9 @@ void polympo_getMeshElmCenter_f(MPMesh_ptr p_mpmesh, const int nCells, double* x auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; //check the size - PMT_ALWAYS_ASSERT(p_mesh->getNumElements()==nCells); + PMT_ALWAYS_ASSERT(p_mesh->getNumElements()==nCells); - //copy the device to host + //copy the device to host auto elmCenter = p_mesh->getMeshField(); auto h_elmCenter = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), elmCenter); for(int i=0; ip_mesh; PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); - //copy the device to host + //copy the device to host auto dualArea = p_mesh->getMeshField(); auto h_dualArea = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), dualArea); for(int i=0; ip_mesh; //check the size - PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); + PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); - //copy the device to host + //copy the device to host auto stressDivergence = p_mesh->getMeshField(); auto h_stressDivergence = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), stressDivergence); for(int i=0; ip_mesh; @@ -1745,7 +1745,7 @@ void polympo_finalize_deludelvDyn_f(MPMesh_ptr p_mpmesh){ auto vtxField = p_mesh->getMeshField(); auto vtxFieldVel = p_mesh->getMeshField(); auto vtxFieldVel_incr = p_mesh->getMeshField(); - + Kokkos::parallel_for("Finalize_increments", nVertices, KOKKOS_LAMBDA(const int vtx){ vtxField(vtx, 0) = vtxField(vtx, 0) * elasticTimeStep; vtxField(vtx, 1) = vtxField(vtx, 1) * elasticTimeStep; From b993971833994cd313d546d52430f5dc27cbffff Mon Sep 17 00:00:00 2001 From: Shahrear Jahan Santho Date: Thu, 24 Sep 2026 23:04:27 -0700 Subject: [PATCH 04/12] Clean up remaining unrelated formatting changes --- src/pmpo_MPMesh.cpp | 6 +++--- src/pmpo_MPMesh.hpp | 20 ++++++++++++-------- src/pmpo_MPMesh_assembly.hpp | 14 +++++++------- src/pmpo_c.cpp | 12 ++++++------ 4 files changed, 28 insertions(+), 24 deletions(-) diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index 8eb41183..008a42e3 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -166,17 +166,17 @@ void MPMesh::calculateStressDivergence(){ }; p_MPs->parallel_for(stress_div, " stress_div_assembly"); Kokkos::fence(); - pumipic::RecordTime("Stress_Divergence_Reconstruction" + std::to_string(self), timer.seconds()); + pumipic::RecordTime("Stress_Divergence_Reconstruction" + std::to_string(self), timer.seconds()); timer.reset(); - if(numProcsTot>1){ + if(numProcsTot>1){ //Takes contribution of halo vertices and adds it in owner procs communicate_and_take_halo_contributions_gpu_aware(stress_divUV, numVertices, 2, 0, 0); //Transfer the correct values at owned vertices to halo vertices //communicate_and_take_halo_contributions(stress_divUV, numVertices, 2, 1, 1); } Kokkos::fence(); - pumipic::RecordTime("Stress_Divergence Communication" + std::to_string(self), timer.seconds()); + pumipic::RecordTime("Stress_Divergence Communication" + std::to_string(self), timer.seconds()); } void MPMesh::calcBasis() { diff --git a/src/pmpo_MPMesh.hpp b/src/pmpo_MPMesh.hpp index 62411d0f..0369bc6a 100644 --- a/src/pmpo_MPMesh.hpp +++ b/src/pmpo_MPMesh.hpp @@ -42,8 +42,10 @@ class MPMesh{ std::vector> ownerHaloLocalIDs; void startCommunication(); - void communicate_and_take_halo_contributions(const Kokkos::View& meshField, int nEntities, int numEntries, int mode, int op); + void communicate_and_take_halo_contributions(const Kokkos::View& meshField, int nEntities, int numEntries, int mode, int op); + + //Now Kokkos views are made 1D template void communicate_and_take_halo_contributions_staged( const ViewType& meshField, @@ -131,15 +133,15 @@ class MPMesh{ } void communicateFields(const std::vector>& fieldData, const int numEntities, const int numEntries, int mode, - std::vector>& recvIDVec, std::vector>& recvDataVec); + std::vector>& recvIDVec, std::vector>& recvDataVec); - template + template void communicateFieldsFromHostView( - const ViewType& fieldData, + const ViewType& fieldData, const int numEntities, const int numEntries, int mode, std::vector>& recvIDVec, std::vector>& recvDataVec){ - + int self, numProcsTot; MPI_Comm comm = p_MPs->getMPIComm(); MPI_Comm_rank(comm, &self); @@ -197,7 +199,7 @@ class MPMesh{ std::vector requests; requests.reserve(4*numProcsTot); for(int proc = 0; proc < numProcsTot; proc++){ - if(proc == self) continue; + if(proc == self) continue; if(mode == 0 && numHalosOnOtherProcs[proc]){ assert(recvIDVec[proc].size() == (size_t)numHalosOnOtherProcs[proc]); assert(recvDataVec[proc].size() == recvIDVec[proc].size() * (size_t)numEntries); @@ -234,7 +236,7 @@ class MPMesh{ } MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE); } - + MPMesh(Mesh* inMesh, MaterialPoints* inMPs): p_mesh(inMesh), p_MPs(inMPs) { }; @@ -269,7 +271,7 @@ class MPMesh{ void reconstruct_coeff_full(); void invertMatrix(const Kokkos::View& vtxMatrices, const double& radius); Kokkos::View precomputedVtxCoeffs_new; - Kokkos::View nearAnEdge; + Kokkos::View nearAnEdge; Kokkos::View vtxMatrixMass; //Not used currently @@ -287,6 +289,8 @@ class MPMesh{ void printVTP_mesh(int printVTPIndex); void writeMPTrackingVTP(int printVTPIndex, int numMPs, const Vec3dView& history, const Vec3dView& resultLeft, const Vec3dView& resultRight, const Vec3dView& mpTgtPosArray); + + void calculateStrain(); void calculateStress(const int constitutive_relation); void calculateStressDivergence(); diff --git a/src/pmpo_MPMesh_assembly.hpp b/src/pmpo_MPMesh_assembly.hpp index c4769aac..c26d3899 100644 --- a/src/pmpo_MPMesh_assembly.hpp +++ b/src/pmpo_MPMesh_assembly.hpp @@ -167,19 +167,19 @@ void MPMesh::reconstruct_coeff_full(){ communicate_and_take_halo_contributions_gpu_aware(vtxMatrices, numVertices, numEntriesMatrix, mode, op); } pumipic::RecordTime("Communicate Matrix Values" + std::to_string(self), timer.seconds()); - + //Store the 1st matrix element Kokkos::ViewvtxMatrixMass_l("vtxMass", numVertices); Kokkos::parallel_for("storeMatrixMass", numVertices, KOKKOS_LAMBDA(const int vtx){ vtxMatrixMass_l(vtx) = vtxMatrices(vtx, 0); - }); + }); this->vtxMatrixMass = vtxMatrixMass_l; invertMatrix(vtxMatrices, radius); } void MPMesh::invertMatrix(const Kokkos::View& vtxMatrices, const double& radius){ - + int nVertices = p_mesh->getNumVertices(); auto vtxCoords = p_mesh->getMeshField(); auto dual_triangle_area = p_mesh->getMeshField(); @@ -209,7 +209,7 @@ void MPMesh::invertMatrix(const Kokkos::View& vtxMatrices, const doubl } auto cosLat = sqrt(pow(X, 2) + pow(Y, 2)); - auto invCosLat = 1.0/cosLat; + auto invCosLat = 1.0/cosLat; auto vtx_area_sqrt = sqrt(dual_triangle_area(vtx,0)/(radius*radius)); Vec3d v0 = { -Y * invCosLat, -Z * X * invCosLat, X / vtx_area_sqrt }; @@ -379,7 +379,7 @@ void MPMesh::assemblyVtx1(){ pumipic::RecordTime("Assemble Field per process" + std::to_string(self), timer.seconds()); timer.reset(); - if(numProcsTot>1){ + if(numProcsTot>1){ communicate_and_take_halo_contributions_gpu_aware(meshField, numVertices, numEntries, 0, 0); } pumipic::RecordTime("Communicate Field Values" + std::to_string(self), timer.seconds()); @@ -428,7 +428,7 @@ DoubleView MPMesh::wtScaAssembly(){ // last component of eVtxCoords stores the firs vertex (to avoid if-condition in the Wachspress computation) eVtxCoords[nElmVtxs][0] = vtxCoords(elm2VtxConn(elm,1)-1,0); eVtxCoords[nElmVtxs][1] = vtxCoords(elm2VtxConn(elm,1)-1,1); - + /* compute the values of basis functions at mp position */ double basisByArea[maxElmsPerVtx]; Vec2d mpCoord(mpPositions(mp,0), mpPositions(mp,1)); @@ -467,7 +467,7 @@ Vec2dView MPMesh::wtVec2Assembly(){ Vec2d eVtxCoords[maxVtxsPerElm + 1]; for (int i = 1; i <= nElmVtxs; i++) { // elm2VtxConn(elm,i) is the vertex ID (1-based index) of vertex #i of elm - eVtxCoords[i-1][0] = vtxCoords(elm2VtxConn(elm,i)-1,0); + eVtxCoords[i-1][0] = vtxCoords(elm2VtxConn(elm,i)-1,0); eVtxCoords[i-1][1] = vtxCoords(elm2VtxConn(elm,i)-1,1); } // last component of eVtxCoords stores the firs vertex (to avoid if-condition in the Wachspress computation) diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index 385111dd..2ee56258 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -6,7 +6,7 @@ namespace{ std::vector p_mpmeshes;////store the p_mpmeshes that is legal - + void checkMPMeshValid(MPMesh_ptr p_mpmesh){ auto p_mpmeshIter = std::find(p_mpmeshes.begin(),p_mpmeshes.end(),p_mpmesh); PMT_ALWAYS_ASSERT(p_mpmeshIter != p_mpmeshes.end()); @@ -922,8 +922,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){ //chech vailidity checkMPMeshValid(p_mpmesh); - kkViewHostU arrayHost(array,maxEdges,nCells); - auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; + kkViewHostU arrayHost(array,maxEdges,nCells); + auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; PMT_ALWAYS_ASSERT(p_mesh->meshEditable()); //check the size @@ -944,7 +944,7 @@ void polympo_setMeshElm2ElmConn_f(MPMesh_ptr p_mpmesh, const int maxEdges, const //chech vailidity checkMPMeshValid(p_mpmesh); kkViewHostU arrayHost(array,maxEdges,nCells); //Fortran is column-major - auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; + auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; PMT_ALWAYS_ASSERT(p_mesh->meshEditable()); //check the size @@ -1068,7 +1068,7 @@ void polympo_getMeshVtxCoords_f(MPMesh_ptr p_mpmesh, const int nVertices, double auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; //check the size - PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); + PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); //copy the device to host auto coordsArray = p_mesh->getMeshField(); @@ -1142,7 +1142,7 @@ void polympo_getMeshVtxVel_f(MPMesh_ptr p_mpmesh, const int nVertices, double* u auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; //check the size - PMT_ALWAYS_ASSERT(p_mesh->getNumVertices() == nVertices); + PMT_ALWAYS_ASSERT(p_mesh->getNumVertices() == nVertices); //copy the device array to the host auto coordsArray = p_mesh->getMeshField(); From e5096121343fdf1a3cb0ec1f1460a7dbbe945698 Mon Sep 17 00:00:00 2001 From: Shahrear Jahan Santho Date: Thu, 24 Sep 2026 23:13:05 -0700 Subject: [PATCH 05/12] Clean up remaining formatting changes --- src/pmpo_MPMesh.cpp | 2 +- src/pmpo_MPMesh.hpp | 7 +++++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index 008a42e3..f044d26f 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -166,7 +166,7 @@ void MPMesh::calculateStressDivergence(){ }; p_MPs->parallel_for(stress_div, " stress_div_assembly"); Kokkos::fence(); - pumipic::RecordTime("Stress_Divergence_Reconstruction" + std::to_string(self), timer.seconds()); + pumipic::RecordTime("Stress_Divergence_Reconstruction" + std::to_string(self), timer.seconds()); timer.reset(); if(numProcsTot>1){ diff --git a/src/pmpo_MPMesh.hpp b/src/pmpo_MPMesh.hpp index 0369bc6a..2cf71ede 100644 --- a/src/pmpo_MPMesh.hpp +++ b/src/pmpo_MPMesh.hpp @@ -42,7 +42,7 @@ class MPMesh{ std::vector> ownerHaloLocalIDs; void startCommunication(); - + void communicate_and_take_halo_contributions(const Kokkos::View& meshField, int nEntities, int numEntries, int mode, int op); //Now Kokkos views are made 1D @@ -66,6 +66,7 @@ class MPMesh{ std::vector> recvIDVec; std::vector> recvDataVec; pumipic::RecordTime("SD: Recv Vec Allocation-" + std::to_string(self), timer.seconds()); + timer.reset(); //communicateFieldsFromHostView(fieldData1, nEntities, numEntries, mode, recvIDVec, recvDataVec); communicateFieldsFromHostView(reconVals_host, nEntities, numEntries, mode, recvIDVec, recvDataVec); @@ -119,6 +120,7 @@ class MPMesh{ assert(recvDataVec[i].size() == recvIDVec[i].size() * numEntries); } pumipic::RecordTime("SD: Copy CPU-GPU2-" + std::to_string(self), timer.seconds()); + //Take contributions from other procs timer.reset(); Kokkos::parallel_for("halo contribution", recvIDGPU.size(), KOKKOS_LAMBDA(const int i){ @@ -132,8 +134,9 @@ class MPMesh{ pumipic::RecordTime("SD: Contribution" + std::to_string(self), timer.seconds()); } + void communicateFields(const std::vector>& fieldData, const int numEntities, const int numEntries, int mode, - std::vector>& recvIDVec, std::vector>& recvDataVec); + std::vector>& recvIDVec, std::vector>& recvDataVec); template void communicateFieldsFromHostView( From b096020666f48c06259efd3191e4068aa936f1e4 Mon Sep 17 00:00:00 2001 From: Shahrear Jahan Santho Date: Thu, 24 Sep 2026 23:18:39 -0700 Subject: [PATCH 06/12] Clean up --- src/pmpo_MPMesh.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index f044d26f..86f3cfd1 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -166,7 +166,7 @@ void MPMesh::calculateStressDivergence(){ }; p_MPs->parallel_for(stress_div, " stress_div_assembly"); Kokkos::fence(); - pumipic::RecordTime("Stress_Divergence_Reconstruction" + std::to_string(self), timer.seconds()); + pumipic::RecordTime("Stress_Divergence_Reconstruction" + std::to_string(self), timer.seconds()); timer.reset(); if(numProcsTot>1){ From 7a46fd520c32f27aef000deeb7c00c1e4f5d3d10 Mon Sep 17 00:00:00 2001 From: Shahrear Jahan Date: Fri, 25 Sep 2026 10:35:07 -0700 Subject: [PATCH 07/12] Make GPU-aware MPI a CMake option, on by default for Cray GPU builds --- src/CMakeLists.txt | 45 ++++++++++++++++++++++++++++++++------------- 1 file changed, 32 insertions(+), 13 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 32026a3e..fc9ba803 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -22,24 +22,43 @@ set(SOURCES add_library(polyMPO-core ${SOURCES}) -if(DEFINED ENV{PE_MPICH_GTL_DIR_nvidia80} +# Default ON only for Cray GPU builds (Kokkos CUDA + Cray MPICH GTL available). +# Can still be overridden with -DpolyMPO_ENABLE_GPU_AWARE_MPI=ON/OFF. +set(polyMPO_GPU_AWARE_MPI_DEFAULT OFF) +if(Kokkos_ENABLE_CUDA + AND DEFINED ENV{PE_MPICH_GTL_DIR_nvidia80} AND DEFINED ENV{PE_MPICH_GTL_LIBS_nvidia80}) + set(polyMPO_GPU_AWARE_MPI_DEFAULT ON) +endif() + +option(polyMPO_ENABLE_GPU_AWARE_MPI + "Use GPU-aware (CUDA-aware) MPI for halo exchange" + ${polyMPO_GPU_AWARE_MPI_DEFAULT}) +message(STATUS "polyMPO_ENABLE_GPU_AWARE_MPI: ${polyMPO_ENABLE_GPU_AWARE_MPI}") + +if(polyMPO_ENABLE_GPU_AWARE_MPI) + if(NOT Kokkos_ENABLE_CUDA) + message(FATAL_ERROR + "polyMPO_ENABLE_GPU_AWARE_MPI=ON requires Kokkos built with CUDA.") + endif() - target_link_options(polyMPO-core PUBLIC - "$ENV{PE_MPICH_GTL_DIR_nvidia80}" - "$ENV{PE_MPICH_GTL_LIBS_nvidia80}" - ) + if(DEFINED ENV{PE_MPICH_GTL_DIR_nvidia80} + AND DEFINED ENV{PE_MPICH_GTL_LIBS_nvidia80}) + target_link_options(polyMPO-core PUBLIC + "$ENV{PE_MPICH_GTL_DIR_nvidia80}" + "$ENV{PE_MPICH_GTL_LIBS_nvidia80}" + ) + message(STATUS "polyMPO: enabling Cray MPICH CUDA GTL linkage") + else() + message(FATAL_ERROR + "Cray MPICH CUDA GTL environment variables are unavailable. " + "Load cudatoolkit and craype-accel-nvidia80, or set " + "CRAY_ACCEL_TARGET=nvidia80 before configuring.") + endif() - message(STATUS - "polyMPO: enabling Cray MPICH CUDA GTL linkage") -else() - message(FATAL_ERROR - "Cray MPICH CUDA GTL environment variables are unavailable. " - "Load cudatoolkit and craype-accel-nvidia80, or set " - "CRAY_ACCEL_TARGET=nvidia80 before configuring.") + target_compile_definitions(polyMPO-core PUBLIC GPU_AWARE_MPI) endif() -target_compile_definitions(polyMPO-core PUBLIC GPU_AWARE_MPI) set_property(TARGET polyMPO-core PROPERTY CXX_STANDARD "17") set_property(TARGET polyMPO-core PROPERTY CXX_STANDARD_REQUIRED ON) set_property(TARGET polyMPO-core PROPERTY CXX_EXTENSIONS OFF) From 1d95e6ae3795c1046a6b74167332fdab119a31b7 Mon Sep 17 00:00:00 2001 From: Shahrear Jahan Date: Fri, 25 Sep 2026 11:13:59 -0700 Subject: [PATCH 08/12] CI: launch MPI tests with MPICH's mpiexec --- .github/workflows/cmake.yml | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/.github/workflows/cmake.yml b/.github/workflows/cmake.yml index 6c9ef06e..7cd32ca8 100644 --- a/.github/workflows/cmake.yml +++ b/.github/workflows/cmake.yml @@ -26,6 +26,11 @@ jobs: - name: Install mpi run: sudo apt-get install -yq mpich libmpich-dev + - name: Show MPI launchers + run: | + ls -l /usr/bin/mpirun /usr/bin/mpiexec* /etc/alternatives/mpirun || true + mpiexec.mpich --version || true + - name: Install Valgrind run: sudo apt-get install -yq valgrind @@ -238,6 +243,7 @@ jobs: -DIS_TESTING=on -DCMAKE_CXX_COMPILER=mpicxx -DCMAKE_INSTALL_PREFIX=${{ runner.temp }}/build-polyMPO/install + -DMPIRUN=/usr/bin/mpiexec.mpich - name: PolyMPO Build run: cmake --build ${{ runner.temp }}/build-polyMPO -j8 --target install @@ -247,4 +253,4 @@ jobs: - name: PolyMPO Print if: always() - run: cat ${{ runner.temp }}/build-polyMPO/Testing/Temporary/LastTest.log \ No newline at end of file + run: cat ${{ runner.temp }}/build-polyMPO/Testing/Temporary/LastTest.log From 95c18857ede551aea0751116300d6f4147219c35 Mon Sep 17 00:00:00 2001 From: Shahrear Jahan Date: Fri, 25 Sep 2026 11:39:43 -0700 Subject: [PATCH 09/12] CI: restrict UCX transports so MPICH initializes on GitHub runners --- .github/workflows/cmake.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.github/workflows/cmake.yml b/.github/workflows/cmake.yml index 7cd32ca8..d1d54adf 100644 --- a/.github/workflows/cmake.yml +++ b/.github/workflows/cmake.yml @@ -249,6 +249,8 @@ jobs: run: cmake --build ${{ runner.temp }}/build-polyMPO -j8 --target install - name: PolyMPO Test + env: + UCX_TLS: self,sm,tcp run: ctest --test-dir ${{ runner.temp }}/build-polyMPO - name: PolyMPO Print From 78452b8aa4df6e5175ee3e7765b8314ca8566234 Mon Sep 17 00:00:00 2001 From: Shahrear Jahan Date: Fri, 25 Sep 2026 11:56:15 -0700 Subject: [PATCH 10/12] Revert "CI: restrict UCX transports so MPICH initializes on GitHub runners" This reverts commit 95c18857ede551aea0751116300d6f4147219c35. --- .github/workflows/cmake.yml | 2 -- 1 file changed, 2 deletions(-) diff --git a/.github/workflows/cmake.yml b/.github/workflows/cmake.yml index d1d54adf..7cd32ca8 100644 --- a/.github/workflows/cmake.yml +++ b/.github/workflows/cmake.yml @@ -249,8 +249,6 @@ jobs: run: cmake --build ${{ runner.temp }}/build-polyMPO -j8 --target install - name: PolyMPO Test - env: - UCX_TLS: self,sm,tcp run: ctest --test-dir ${{ runner.temp }}/build-polyMPO - name: PolyMPO Print From 29031aac6ec23467474699254322b2ac9999a9cf Mon Sep 17 00:00:00 2001 From: Shahrear Jahan Date: Fri, 25 Sep 2026 11:56:15 -0700 Subject: [PATCH 11/12] Revert "CI: launch MPI tests with MPICH's mpiexec" This reverts commit 1d95e6ae3795c1046a6b74167332fdab119a31b7. --- .github/workflows/cmake.yml | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/.github/workflows/cmake.yml b/.github/workflows/cmake.yml index 7cd32ca8..6c9ef06e 100644 --- a/.github/workflows/cmake.yml +++ b/.github/workflows/cmake.yml @@ -26,11 +26,6 @@ jobs: - name: Install mpi run: sudo apt-get install -yq mpich libmpich-dev - - name: Show MPI launchers - run: | - ls -l /usr/bin/mpirun /usr/bin/mpiexec* /etc/alternatives/mpirun || true - mpiexec.mpich --version || true - - name: Install Valgrind run: sudo apt-get install -yq valgrind @@ -243,7 +238,6 @@ jobs: -DIS_TESTING=on -DCMAKE_CXX_COMPILER=mpicxx -DCMAKE_INSTALL_PREFIX=${{ runner.temp }}/build-polyMPO/install - -DMPIRUN=/usr/bin/mpiexec.mpich - name: PolyMPO Build run: cmake --build ${{ runner.temp }}/build-polyMPO -j8 --target install @@ -253,4 +247,4 @@ jobs: - name: PolyMPO Print if: always() - run: cat ${{ runner.temp }}/build-polyMPO/Testing/Temporary/LastTest.log + run: cat ${{ runner.temp }}/build-polyMPO/Testing/Temporary/LastTest.log \ No newline at end of file From fdb7efda17e4906e9fdf4d4ea07405361f750d5f Mon Sep 17 00:00:00 2001 From: Shahrear Jahan Date: Fri, 25 Sep 2026 12:14:43 -0700 Subject: [PATCH 12/12] Skip GPU-aware halo exchange test when not run with 4 MPI ranks --- test/CMakeLists.txt | 1 + test/testGPUAwareHaloExchange.cpp | 8 +++++--- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index b4af325a..346c01b4 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -93,6 +93,7 @@ pmpo_add_test(testFortran 1 ./testFortran) pmpo_add_test(testFortranMPMeshModule 1 ./testFortranMPMeshModule) pmpo_add_test(testFortranMPAppIDs 1 ./testFortranMPAppIDs) pmpo_add_test(test_gpu_aware_halo_exchange 4 ./testGPUAwareHaloExchange) +set_tests_properties(test_gpu_aware_halo_exchange PROPERTIES SKIP_RETURN_CODE 77) #set NC file for test set(TEST_NC_FILE_PLANAR "${CMAKE_SOURCE_DIR}/test/sample_mpas_meshes/planar_nonuniform_cvt_for_square_673elms.nc") diff --git a/test/testGPUAwareHaloExchange.cpp b/test/testGPUAwareHaloExchange.cpp index 9c4902fa..127fa303 100644 --- a/test/testGPUAwareHaloExchange.cpp +++ b/test/testGPUAwareHaloExchange.cpp @@ -45,11 +45,13 @@ int main(int argc, char** argv) if (size != 4) { if (rank == 0) { std::cerr - << "This test requires exactly 4 MPI ranks." + << "This test requires exactly 4 MPI ranks (got " + << size << "); skipping." << std::endl; } - - MPI_Abort(MPI_COMM_WORLD, 1); + Kokkos::finalize(); + MPI_Finalize(); + return 77; } // Create the existing polyMPO test mesh and MPMesh.