From b98b69a9e36e740c91eb427b9100e5097122a61e Mon Sep 17 00:00:00 2001 From: Nath Date: Mon, 29 Jun 2026 19:40:25 -0400 Subject: [PATCH 01/11] First commit to port --- src/pmpo_c.cpp | 19 +++++++++++++++++++ src/pmpo_c.h | 2 ++ src/pmpo_fortran.f90 | 14 ++++++++++++++ src/pmpo_mesh.cpp | 12 ++++++++++++ src/pmpo_mesh.hpp | 2 ++ 5 files changed, 49 insertions(+) diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index 53c330e1..cc211190 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -1448,6 +1448,20 @@ void polympo_setSolveVelocityMesh_f(MPMesh_ptr p_mpmesh, const int nVertices, in Kokkos::deep_copy(solveVelocity, h_solveVelocity); } +void polympo_setIceAreaVertex_f(MPMesh_ptr p_mpmesh, const int nVertices, double* array){ + //chech validity + checkMPMeshValid(p_mpmesh); + auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; + + PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); + //copy the host array to the device + auto iceArea = p_mesh->getMeshField(); + auto h_iceArea = Kokkos::create_mirror_view(iceArea); + for(int i=0; ip_mesh; + p_mesh->calcOceanStressCoeff(); +} + void polympo_velocity_grid_solve_f(MPMesh_ptr p_mpmesh){ //Temporary grid Solve after calculateDivergence auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; diff --git a/src/pmpo_c.h b/src/pmpo_c.h index 23d78b70..9b759e64 100644 --- a/src/pmpo_c.h +++ b/src/pmpo_c.h @@ -106,6 +106,7 @@ void polympo_setElasticTimeStep_f(MPMesh_ptr p_mpmesh, const double elasticTimeS void polympo_setDynamicTimeStep_f(MPMesh_ptr p_mpmesh, const double dynamicTimeStep); void polympo_setSolveStressMesh_f(MPMesh_ptr p_mpmesh, const int nCells, int* array); void polympo_setSolveVelocityMesh_f(MPMesh_ptr p_mpmesh, const int nVertices, int* array); +void polympo_setIceAreaVertex_f(MPMesh_ptr p_mpmesh, const int nVertices, double* array); void polympo_calculateStressDivergence_f(MPMesh_ptr p_mpmesh); void polympo_getStressDivergence_f(MPMesh_ptr p_mpmesh, const int nVertices, double* uArray, double* vArray); void polympo_setTotalMassVtx_f(MPMesh_ptr p_mpmesh, const int nVertices, double* array); @@ -114,6 +115,7 @@ void polympo_set_surfaceTiltForce_f(MPMesh_ptr p_mpmesh, const int nVertices, do void polympo_set_totalMassVertexfVertex_f(MPMesh_ptr p_mpmesh, const int nVertices, double* array); void polympo_set_oceanStress_f(MPMesh_ptr p_mpmesh, const int nVertices, double* uArray, double* vArray); void polympo_set_oceanStressCoefficient_f(MPMesh_ptr p_mpmesh, const int nVertices, double* array); +void polympo_calculate_oceanStressCoefficient_f(MPMesh_ptr p_mpmesh); void polympo_velocity_grid_solve_f(MPMesh_ptr p_mpmesh); void polympo_set_boundary_normal_vertex_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, double* uArray, double* vArray); void polympo_set_free_slip_bc_f(MPMesh_ptr p_mpmesh); diff --git a/src/pmpo_fortran.f90 b/src/pmpo_fortran.f90 index d0bb903c..a69cb5b3 100644 --- a/src/pmpo_fortran.f90 +++ b/src/pmpo_fortran.f90 @@ -1011,6 +1011,14 @@ subroutine polympo_setSolveVelocityMesh(mpMesh, nVertices, array) & type(c_ptr), value :: array end subroutine + subroutine polympo_setIceAreaVertex(mpMesh, nVertices, array) & + bind(C, NAME='polympo_setIceAreaVertex_f') + use :: iso_c_binding + type(c_ptr), value :: mpMesh + integer(c_int), value :: nVertices + type(c_ptr), value :: array + end subroutine + subroutine polympo_calculateStressDivergence(mpMesh) & bind(C, NAME='polympo_calculateStressDivergence_f') use :: iso_c_binding @@ -1073,6 +1081,12 @@ subroutine polympo_set_oceanStressCoefficient(mpMesh, nVertices, array) & type(c_ptr), value :: array end subroutine + subroutine polympo_calculate_oceanStressCoefficient(mpMesh) & + bind(C, NAME='polympo_calculate_oceanStressCoefficient_f') + use :: iso_c_binding + type(c_ptr), value :: mpMesh + end subroutine + subroutine polympo_velocity_grid_solve(mpMesh) & bind(C, NAME='polympo_velocity_grid_solve_f') use :: iso_c_binding diff --git a/src/pmpo_mesh.cpp b/src/pmpo_mesh.cpp index 0e7849ca..cc89c2d4 100644 --- a/src/pmpo_mesh.cpp +++ b/src/pmpo_mesh.cpp @@ -157,6 +157,18 @@ namespace polyMPO{ }); } + void Mesh::calcOceanStressCoeff(){ + int numVerticesOwned = getNumVerticesOwned(); + auto iceAreaVtx = getMeshField(); + auto oceanStressCoeff = getMeshField(); + auto velocity = getMeshField(); + auto solve_velocity = getMeshField(); + + Kokkos::parallel_for("calcOceanStressCoeff", numVerticesOwned, KOKKOS_LAMBDA(const int vtx){ + if(solve_velocity(vtx) == 0) return; + oceanStressCoeff(vtx, 0) = 0.00536 * 1026.0 * iceAreaVtx(vtx, 0) * sqrt(velocity(vtx, 0)*velocity(vtx, 0) + velocity(vtx, 1)*velocity(vtx, 1)); + }); + } void Mesh::gridSolveGPU(){ //Mesh Fields diff --git a/src/pmpo_mesh.hpp b/src/pmpo_mesh.hpp index a851f114..d85d2a0f 100644 --- a/src/pmpo_mesh.hpp +++ b/src/pmpo_mesh.hpp @@ -271,6 +271,8 @@ class Mesh { double getDynamicTimeStep(){ return dynamicTimeStep_; } + + void calcOceanStressCoeff(); void gridSolveGPU(); void aggregateDeluDyn(); void applyFreeSlipBC(); From 12417d3fcc9a8e628639e1bf54114f9d727f2f81 Mon Sep 17 00:00:00 2001 From: Nath Date: Mon, 29 Jun 2026 21:59:24 -0400 Subject: [PATCH 02/11] Conditional for HALO exchange of nProcsTot=1 in polyMPO --- src/pmpo_c.cpp | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index cc211190..d199675a 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -1624,10 +1624,18 @@ void polympo_set_free_slip_bc_f(MPMesh_ptr p_mpmesh){ } void polympo_set_halo_vel_from_owner_f(MPMesh_ptr p_mpmesh){ + + int numProcsTot; + auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_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); } From 229c8cb0a8dbc95bd52666c51a295864df090dbd Mon Sep 17 00:00:00 2001 From: Nath Date: Mon, 13 Jul 2026 15:18:22 -0400 Subject: [PATCH 03/11] SetMPArea and modify in strain rate --- src/pmpo_MPMesh.cpp | 11 +++++- src/pmpo_c.cpp | 94 +++++++++++++++++++++++++++++++++++++++++++- src/pmpo_c.h | 4 ++ src/pmpo_fortran.f90 | 32 ++++++++++++++- src/pmpo_mesh.cpp | 7 +++- src/pmpo_mesh.hpp | 11 +++++- 6 files changed, 151 insertions(+), 8 deletions(-) diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index 3801dcd9..364fc982 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -14,11 +14,13 @@ void MPMesh::calculateStrain(){ auto MPsBasisGrads = p_MPs->getData(); auto MPsAppID = p_MPs->getData(); auto MPsStrainRate = p_MPs->getData(); + auto MPsArea = p_MPs->getData(); //Mesh Fields auto tanLatVertexRotatedOverRadius = p_mesh->getMeshField(); auto elm2VtxConn = p_mesh->getElm2VtxConn(); auto velField = p_mesh->getMeshField(); auto solveStress = p_mesh->getMeshField(); + auto elasticTimeStep = p_mesh->getElasticTimeStep(); auto setMPStrainRate = PS_LAMBDA(const int& elm, const int& mp, const int& mask){ if(mask){ @@ -29,7 +31,7 @@ void MPMesh::calculateStrain(){ MPsStrainRate(mp, 2) = 0.0; return; } - + int numVtx = elm2VtxConn(elm,0); double v11 = 0.0; @@ -52,6 +54,8 @@ void MPMesh::calculateStrain(){ MPsStrainRate(mp, 0) = v11 - vTanOverR; MPsStrainRate(mp, 1) = v22; MPsStrainRate(mp, 2) = 0.5*(v12 + v21 + uTanOverR); + + MPsArea(mp, 0) = MPsArea(mp, 0) * exp((v11+v22-vTanOverR)*elasticTimeStep); } }; p_MPs->parallel_for(setMPStrainRate, "setMPStrainRate"); @@ -75,13 +79,16 @@ void MPMesh::calculateStress(const int constitutive_relation){ if(mask){ Vec3d strain_rate (MPsStrainRate(mp, 0), MPsStrainRate(mp, 1), MPsStrainRate(mp, 2)); Vec3d stress(MPsStress(mp, 0), MPsStress(mp, 1), MPsStress(mp, 2)); + double rep_pressure=MPsRepPressure(mp,0); if (constitutive_relation == 1) - constitutive_evp(strain_rate, stress, MPsIcePressure(mp,0), MPsRepPressure(mp,0), MPsArea(mp,0), elasticTimeStep, dampingTimescale); + constitutive_evp(strain_rate, stress, MPsIcePressure(mp,0), rep_pressure, MPsArea(mp,0), elasticTimeStep, dampingTimescale); else if(constitutive_relation == 3) constitutive_linear(strain_rate, stress); + for (int m=0 ; m<3; m++) MPsStress(mp, m) = stress[m]*solveStress(elm); + MPsRepPressure(mp,0)=rep_pressure; } }; p_MPs->parallel_for(setMPStress, "setMPStress"); diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index d199675a..7c8d1585 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -728,7 +728,6 @@ void polympo_getMPStress_f(MPMesh_ptr p_mpmesh, const int nComps, const int numM void polympo_setAreaMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* areaMPHost){ Kokkos::Timer timer; checkMPMeshValid(p_mpmesh); - auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; //Rank information int self; @@ -754,10 +753,34 @@ void polympo_setAreaMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs pumipic::RecordTime("PolyMPO_setMPArea" + std::to_string(self), timer.seconds()); } +void polympo_getAreaMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* areaMPHost) { + Kokkos::Timer timer; + checkMPMeshValid(p_mpmesh); + auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; + + PMT_ALWAYS_ASSERT(nComps == 1); + PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount()); + + auto mpArea = p_MPs->getData(); + auto mpAppID = p_MPs->getData(); + + Kokkos::View mpAreaCopy("mpAreaCopy", nComps, numMPs); + auto getMPArea = PS_LAMBDA(const int& elm, const int& mp, const int& mask){ + if(mask){ + mpAreaCopy(0,mpAppID(mp)) = mpArea(mp,0); + } + }; + p_MPs->parallel_for(getMPArea, "getMPArea"); + kkDbl2dViewHostU arrayHost(areaMPHost, nComps, numMPs); + Kokkos::deep_copy(arrayHost, mpAreaCopy); + pumipic::RecordTime("PolyMPO_getMPArea", timer.seconds()); +} + + void polympo_setIcePressureMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* icePressureMPHost){ Kokkos::Timer timer; checkMPMeshValid(p_mpmesh); - + auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; //Rank information int self; @@ -783,6 +806,73 @@ void polympo_setIcePressureMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int pumipic::RecordTime("PolyMPO_setIcePressure" + std::to_string(self), timer.seconds()); } +void polympo_setOceanVelocity_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, const double* uArray, const double* vArray){ + checkMPMeshValid(p_mpmesh); + auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; + + PMT_ALWAYS_ASSERT(nComps == vec2d_nEntries); + PMT_ALWAYS_ASSERT(p_mesh->getNumVertices() == nVertices); + //copy the host array to the device + auto oceanVelocity = p_mesh->getMeshField(); + auto h_oceanVelocity = Kokkos::create_mirror_view(oceanVelocity); + for(int i=0; ip_mesh; + int numVerticesOwned = p_mesh->getNumVerticesOwned(); + + auto solveVelocity = p_mesh->getMeshField(); + auto velField = p_mesh->getMeshField(); + auto stress_divUV = p_mesh->getMeshField(); + auto oceanStress = p_mesh->getMeshField(); + + Kokkos::parallel_for("prep_arrays", numVerticesOwned, KOKKOS_LAMBDA(const int vtx){ + if(solveVelocity(vtx)==0){ + for (int k=0; k<2; k ++){ + velField(vtx, k) = 0.0; + stress_divUV(vtx, k) = 0.0; + oceanStress(vtx, k) = 0.0; + } + } + }); +} + +void polympo_setReplacementPressureMP_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; + MPI_Comm comm = p_MPs->getMPIComm(); + MPI_Comm_rank(comm, &self); + //Asserts + PMT_ALWAYS_ASSERT(nComps == 1); + PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount()); + //MP Data + auto mpReplacementPressure = p_MPs->getData(); + auto mpAppID = p_MPs->getData(); + + //Copy to device + kkViewHostU mpRepPressure_h(replacementPressureMPHost, nComps, numMPs); + Kokkos::View mpRepPressure_d("mpRepPressureDevice", nComps, numMPs); + Kokkos::deep_copy(mpRepPressure_d, mpRepPressure_h); + //Set in PS + auto setMPRepPressure = PS_LAMBDA(const int& elm, const int& mp, const int& mask){ + if(mask){ + mpReplacementPressure(mp,0) = mpRepPressure_d(0, mpAppID(mp)); + } + }; + p_MPs->parallel_for(setMPRepPressure, "setMPRepPressure"); + pumipic::RecordTime("PolyMPO_setReplacementPressure" + std::to_string(self), timer.seconds()); +} + void polympo_getReplacementPressureMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* replacementPressureMPHost){ Kokkos::Timer timer; checkMPMeshValid(p_mpmesh); diff --git a/src/pmpo_c.h b/src/pmpo_c.h index 9b759e64..99ba04f9 100644 --- a/src/pmpo_c.h +++ b/src/pmpo_c.h @@ -53,7 +53,11 @@ void polympo_calculateMPStress_f(MPMesh_ptr p_mpmesh, const int constitutive_mod void polympo_setMPStress_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpStressIn); void polympo_getMPStress_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpStressHost); void polympo_setAreaMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* areaMPHost); +void polympo_getAreaMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* areaMPHost); void polympo_setIcePressureMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* icePressureMPHost); +void polympo_setOceanVelocity_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, const double* uArray, const double* vArray); +void polympo_subcycle_prep_arrays_f(MPMesh_ptr p_mpmesh); +void polympo_setReplacementPressureMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* replacementPressureMPHost); void polympo_getReplacementPressureMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* replacementPressureMPHost); //Mesh info diff --git a/src/pmpo_fortran.f90 b/src/pmpo_fortran.f90 index a69cb5b3..522f3b22 100644 --- a/src/pmpo_fortran.f90 +++ b/src/pmpo_fortran.f90 @@ -419,6 +419,14 @@ subroutine polympo_setAreaMP(mpMesh, nComps, numMPs, array) & type(c_ptr), value :: array end subroutine + subroutine polympo_getAreaMP(mpMesh, nComps, numMPs, array) & + bind(C, NAME='polympo_getAreaMP_f') + use :: iso_c_binding + type(c_ptr), value :: mpMesh + integer(c_int), value :: nComps, numMPs + type(c_ptr), value :: array + end subroutine + subroutine polympo_setIcePressureMP(mpMesh, nComps, numMPs, array) & bind(C, NAME='polympo_setIcePressureMP_f') use :: iso_c_binding @@ -426,7 +434,29 @@ subroutine polympo_setIcePressureMP(mpMesh, nComps, numMPs, array) & integer(c_int), value :: nComps, numMPs type(c_ptr), value :: array end subroutine - + + subroutine polympo_setOceanVelocity(mpMesh, nComps, nVertices, uArray, vArray) & + bind(C, NAME='polympo_setOceanVelocity_f') + use :: iso_c_binding + type(c_ptr), value :: mpMesh + integer(c_int), value :: nComps, nVertices + type(c_ptr), value :: uArray, vArray + end subroutine + + subroutine polympo_subcycle_prep_arrays(mpMesh) & + bind(C, NAME='polympo_subcycle_prep_arrays_f') + use :: iso_c_binding + type(c_ptr), value :: mpMesh + end subroutine + + subroutine polympo_setReplacementPressureMP(mpMesh, nComps, numMPs, array) & + bind(C, NAME='polympo_setReplacementPressureMP_f') + use :: iso_c_binding + type(c_ptr), value :: mpMesh + integer(c_int), value :: nComps, numMPs + type(c_ptr), value :: array + end subroutine + subroutine polympo_getReplacementPressureMP(mpMesh, nComps, numMPs, array) & bind(C, NAME='polympo_getReplacementPressureMP_f') use :: iso_c_binding diff --git a/src/pmpo_mesh.cpp b/src/pmpo_mesh.cpp index cc89c2d4..5acb9454 100644 --- a/src/pmpo_mesh.cpp +++ b/src/pmpo_mesh.cpp @@ -78,6 +78,9 @@ namespace polyMPO{ oceanStress_ = MeshFView(meshFields2TypeAndString.at(MeshF_OceanStress).second, numVtxs_); oceanStressCoeff_ = MeshFView(meshFields2TypeAndString.at(MeshF_OceanStressCoeff).second, numVtxs_); + + oceanVelocity_ = MeshFView(meshFields2TypeAndString.at(MeshF_OceanVelocity).second, numVtxs_); + } void Mesh::setMeshElmBasedFieldSize(){ @@ -163,10 +166,12 @@ namespace polyMPO{ auto oceanStressCoeff = getMeshField(); auto velocity = getMeshField(); auto solve_velocity = getMeshField(); + auto oceanVelocity = getMeshField(); Kokkos::parallel_for("calcOceanStressCoeff", numVerticesOwned, KOKKOS_LAMBDA(const int vtx){ if(solve_velocity(vtx) == 0) return; - oceanStressCoeff(vtx, 0) = 0.00536 * 1026.0 * iceAreaVtx(vtx, 0) * sqrt(velocity(vtx, 0)*velocity(vtx, 0) + velocity(vtx, 1)*velocity(vtx, 1)); + auto relVelSq = pow(oceanVelocity(vtx, 0) - velocity(vtx, 0), 2) + pow(oceanVelocity(vtx, 1) - velocity(vtx, 1), 2); + oceanStressCoeff(vtx, 0) = 0.00536 * 1026.0 * iceAreaVtx(vtx, 0) * sqrt(relVelSq); }); } diff --git a/src/pmpo_mesh.hpp b/src/pmpo_mesh.hpp index d85d2a0f..40b7e03a 100644 --- a/src/pmpo_mesh.hpp +++ b/src/pmpo_mesh.hpp @@ -41,7 +41,8 @@ enum MeshFieldIndex{ MeshF_SurfaceTilt, MeshF_TotalMassFVtx, MeshF_OceanStress, - MeshF_OceanStressCoeff + MeshF_OceanStressCoeff, + MeshF_OceanVelocity }; enum MeshFieldType{ @@ -77,6 +78,7 @@ template <> struct meshFieldToType < MeshF_SurfaceTilt > { using type = Ko template <> struct meshFieldToType < MeshF_TotalMassFVtx > { using type = Kokkos::View; }; template <> struct meshFieldToType < MeshF_OceanStress > { using type = Kokkos::View; }; template <> struct meshFieldToType < MeshF_OceanStressCoeff > { using type = Kokkos::View; }; +template <> struct meshFieldToType < MeshF_OceanVelocity > { using type = Kokkos::View; }; template using MeshFView = typename meshFieldToType::type; @@ -108,7 +110,8 @@ const std::map> meshFields {MeshF_SurfaceTilt, {MeshFType_VtxBased,"MeshField_SurfaceTilt"}}, {MeshF_TotalMassFVtx, {MeshFType_VtxBased,"MeshField_TotalMassFVtx"}}, {MeshF_OceanStress, {MeshFType_VtxBased,"MeshField_OceanStress"}}, - {MeshF_OceanStressCoeff, {MeshFType_VtxBased,"MeshField_OceanStressCoeff"}} + {MeshF_OceanStressCoeff, {MeshFType_VtxBased,"MeshField_OceanStressCoeff"}}, + {MeshF_OceanVelocity, {MeshFType_VtxBased,"MeshField_OceanVelocity"}} }; enum mesh_type {mesh_unrecognized_lower = -1, @@ -167,6 +170,7 @@ class Mesh { MeshFView totalMassFVtx_; MeshFView oceanStress_; MeshFView oceanStressCoeff_; + MeshFView oceanVelocity_; bool isRotatedFlag = false; double elasticTimeStep_; @@ -363,6 +367,9 @@ auto Mesh::getMeshField(){ else if constexpr (index==MeshF_OceanStressCoeff){ return oceanStressCoeff_; } + else if constexpr (index==MeshF_OceanVelocity){ + return oceanVelocity_; + } fprintf(stderr,"Mesh Field Index error!\n"); exit(1); } From f0c25f71d87d2897cf208af78e26df5621dcc0c2 Mon Sep 17 00:00:00 2001 From: Nath Date: Tue, 14 Jul 2026 16:33:51 -0400 Subject: [PATCH 04/11] Passing oceanDragCoeff and using same vel corrections as rebase 29 June --- src/pmpo_MPMesh_assembly.hpp | 27 +++++++++++++++++++++++---- src/pmpo_c.cpp | 23 +++++++++++++++++++++-- src/pmpo_c.h | 3 ++- src/pmpo_fortran.f90 | 11 ++++++++++- src/pmpo_mesh.cpp | 9 +++++++-- src/pmpo_mesh.hpp | 11 +++++++++-- src/pmpo_utils.hpp | 21 +++++++++++++++++++++ src/pmpo_wachspressBasis.hpp | 33 ++++++++++++++++++++++++++++----- 8 files changed, 121 insertions(+), 17 deletions(-) diff --git a/src/pmpo_MPMesh_assembly.hpp b/src/pmpo_MPMesh_assembly.hpp index 9f184e05..e4099d54 100644 --- a/src/pmpo_MPMesh_assembly.hpp +++ b/src/pmpo_MPMesh_assembly.hpp @@ -342,16 +342,24 @@ void MPMesh::assemblyVtx1(){ p_mesh->fillMeshField(numVtx, numEntries, 0.0); auto meshField = p_mesh->getMeshField(); + auto vtxRotLon = p_mesh->getMeshField(); + //Material Points auto mpData = p_MPs->getData(); auto weight = p_MPs->getData(); auto mpPositions = p_MPs->getData(); - + auto curPosRotLatLon = p_MPs->getData(); + auto MPsAppID = p_MPs->getData(); //Earth Radius double radius = 1.0; if(p_mesh->getGeomType() == geom_spherical_surf) radius=p_mesh->getSphereRadius(); + bool use_correction_term = false; + if constexpr (meshFieldIndex == MeshF_Vel) { + use_correction_term = true; + } + //Reconstruct auto reconstruct = PS_LAMBDA(const int& elm, const int& mp, const int& mask) { if(mask) { //if material point is 'active'/'enabled' @@ -367,9 +375,20 @@ void MPMesh::assemblyVtx1(){ VtxCoeffs_new(vID,0, 2)*CoordDiffs[2] + VtxCoeffs_new(vID,0, 3)*CoordDiffs[3]); - for (int k=0; kp_mesh; + + //check the size + PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); + + //copy the host array to the device + auto coordsArray = p_mesh->getMeshField(); + auto h_coordsArray = Kokkos::create_mirror_view(coordsArray); + for(int i=0; ip_mesh; - p_mesh->calcOceanStressCoeff(); + p_mesh->calcOceanStressCoeff(configIceOceanDragCoeff); } void polympo_velocity_grid_solve_f(MPMesh_ptr p_mpmesh){ diff --git a/src/pmpo_c.h b/src/pmpo_c.h index 99ba04f9..d19b62b5 100644 --- a/src/pmpo_c.h +++ b/src/pmpo_c.h @@ -89,6 +89,7 @@ int polympo_getMeshFElmType_f(); void polympo_setMeshVtxCoords_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* xArray, const double* yArray, const double* zArray); void polympo_getMeshVtxCoords_f(MPMesh_ptr p_mpmesh, const int nVertices, double* xArray, double* yArray, double* zArray); void polympo_setMeshVtxRotLat_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* latitude); +void polympo_setMeshVtxRotLon_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* longitude); void polympo_getMeshVtxRotLat_f(MPMesh_ptr p_mpmesh, const int nVertices, double* latitude); void polympo_setMeshVtxVel_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* uVelocity, const double* vVelocity); void polympo_getMeshVtxVel_f(MPMesh_ptr p_mpmesh, const int nVertices, double* uVelocity, double* vVelocity); @@ -119,7 +120,7 @@ void polympo_set_surfaceTiltForce_f(MPMesh_ptr p_mpmesh, const int nVertices, do void polympo_set_totalMassVertexfVertex_f(MPMesh_ptr p_mpmesh, const int nVertices, double* array); void polympo_set_oceanStress_f(MPMesh_ptr p_mpmesh, const int nVertices, double* uArray, double* vArray); void polympo_set_oceanStressCoefficient_f(MPMesh_ptr p_mpmesh, const int nVertices, double* array); -void polympo_calculate_oceanStressCoefficient_f(MPMesh_ptr p_mpmesh); +void polympo_calculate_oceanStressCoefficient_f(MPMesh_ptr p_mpmesh, const double configIceOceanDragCoeff); void polympo_velocity_grid_solve_f(MPMesh_ptr p_mpmesh); void polympo_set_boundary_normal_vertex_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, double* uArray, double* vArray); void polympo_set_free_slip_bc_f(MPMesh_ptr p_mpmesh); diff --git a/src/pmpo_fortran.f90 b/src/pmpo_fortran.f90 index 522f3b22..709b895e 100644 --- a/src/pmpo_fortran.f90 +++ b/src/pmpo_fortran.f90 @@ -786,6 +786,14 @@ subroutine polympo_getMeshVtxRotLat(mpMesh, nVertices, latitude) & type(c_ptr), value :: latitude end subroutine + subroutine polympo_setMeshVtxRotLon(mpMesh, nVertices, longitude) & + bind(C, NAME='polympo_setMeshVtxRotLon_f') + use :: iso_c_binding + type(c_ptr), value :: mpMesh + integer(c_int), value :: nVertices + type(c_ptr), intent(in), value :: longitude + end subroutine + !--------------------------------------------------------------------------- !> @brief set the vertices velocity from a host array !> @param mpmesh(in/out) MPMesh object @@ -1111,10 +1119,11 @@ subroutine polympo_set_oceanStressCoefficient(mpMesh, nVertices, array) & type(c_ptr), value :: array end subroutine - subroutine polympo_calculate_oceanStressCoefficient(mpMesh) & + subroutine polympo_calculate_oceanStressCoefficient(mpMesh, configIceOceanDragCoeff) & bind(C, NAME='polympo_calculate_oceanStressCoefficient_f') use :: iso_c_binding type(c_ptr), value :: mpMesh + real(c_double), value::configIceOceanDragCoeff end subroutine subroutine polympo_velocity_grid_solve(mpMesh) & diff --git a/src/pmpo_mesh.cpp b/src/pmpo_mesh.cpp index 5acb9454..4d7b2eae 100644 --- a/src/pmpo_mesh.cpp +++ b/src/pmpo_mesh.cpp @@ -20,6 +20,10 @@ namespace polyMPO{ PMT_ALWAYS_ASSERT(vtxRotLatMapEntry.first == MeshFType_VtxBased); vtxRotLat_ = MeshFView(vtxRotLatMapEntry.second,numVtxs_); + auto vtxRotLonMapEntry = meshFields2TypeAndString.at(MeshF_VtxRotLon); + PMT_ALWAYS_ASSERT(vtxRotLonMapEntry.first == MeshFType_VtxBased); + vtxRotLon_ = MeshFView(vtxRotLonMapEntry.second,numVtxs_); + auto vtxVelMapEntry = meshFields2TypeAndString.at(MeshF_Vel); PMT_ALWAYS_ASSERT(vtxVelMapEntry.first == MeshFType_VtxBased); vtxVel_ = MeshFView(vtxVelMapEntry.second,numVtxs_); @@ -160,18 +164,19 @@ namespace polyMPO{ }); } - void Mesh::calcOceanStressCoeff(){ + void Mesh::calcOceanStressCoeff(const double configIceOceanDragCoeff){ int numVerticesOwned = getNumVerticesOwned(); auto iceAreaVtx = getMeshField(); auto oceanStressCoeff = getMeshField(); auto velocity = getMeshField(); auto solve_velocity = getMeshField(); auto oceanVelocity = getMeshField(); + auto seaiceDensitySeaWater_ = polyMPO::seaiceDensitySeaWater; Kokkos::parallel_for("calcOceanStressCoeff", numVerticesOwned, KOKKOS_LAMBDA(const int vtx){ if(solve_velocity(vtx) == 0) return; auto relVelSq = pow(oceanVelocity(vtx, 0) - velocity(vtx, 0), 2) + pow(oceanVelocity(vtx, 1) - velocity(vtx, 1), 2); - oceanStressCoeff(vtx, 0) = 0.00536 * 1026.0 * iceAreaVtx(vtx, 0) * sqrt(relVelSq); + oceanStressCoeff(vtx, 0) = configIceOceanDragCoeff * seaiceDensitySeaWater_ * iceAreaVtx(vtx, 0) * sqrt(relVelSq); }); } diff --git a/src/pmpo_mesh.hpp b/src/pmpo_mesh.hpp index 40b7e03a..c0e3ffa4 100644 --- a/src/pmpo_mesh.hpp +++ b/src/pmpo_mesh.hpp @@ -19,6 +19,7 @@ enum MeshFieldIndex{ MeshF_Unsupported, MeshF_VtxCoords, MeshF_VtxRotLat, + MeshF_VtxRotLon, MeshF_ElmCenterXYZ, MeshF_DualTriangleArea, MeshF_Vel, @@ -55,6 +56,7 @@ enum MeshFieldType{ template struct meshFieldToType; template <> struct meshFieldToType < MeshF_VtxCoords > { using type = Kokkos::View; }; template <> struct meshFieldToType < MeshF_VtxRotLat > { using type = DoubleView; }; +template <> struct meshFieldToType < MeshF_VtxRotLon > { using type = DoubleView; }; template <> struct meshFieldToType < MeshF_ElmCenterXYZ > { using type = Kokkos::View; }; template <> struct meshFieldToType < MeshF_DualTriangleArea > { using type = Kokkos::View; }; template <> struct meshFieldToType < MeshF_Vel > { using type = Kokkos::View; }; @@ -88,6 +90,7 @@ const std::map> meshFields {MeshF_Unsupported, {MeshFType_Unsupported,"MeshField_Unsupported"}}, {MeshF_VtxCoords, {MeshFType_VtxBased,"MeshField_VerticesCoords"}}, {MeshF_VtxRotLat, {MeshFType_VtxBased,"MeshField_VerticesLatitude"}}, + {MeshF_VtxRotLon, {MeshFType_VtxBased,"MeshField_VerticesLongitude"}}, {MeshF_ElmCenterXYZ, {MeshFType_ElmBased,"MeshField_ElementCenterXYZ"}}, {MeshF_DualTriangleArea, {MeshFType_VtxBased,"MeshField_DualTriangleArea"}}, {MeshF_Vel, {MeshFType_VtxBased,"MeshField_Velocity"}}, @@ -144,6 +147,7 @@ class Mesh { //start of meshFields MeshFView vtxCoords_; MeshFView vtxRotLat_; + MeshFView vtxRotLon_; MeshFView elmCenterXYZ_; MeshFView dualTriangleArea_; @@ -159,7 +163,7 @@ class Mesh { //GnomonicProjection MeshFView vtxGnomProj_; MeshFView elmCenterGnomProj_; - + MeshFView tanLatVertexRotatedOverRadius_; MeshFView solveStress_; MeshFView solveVelocity_; @@ -276,7 +280,7 @@ class Mesh { return dynamicTimeStep_; } - void calcOceanStressCoeff(); + void calcOceanStressCoeff(const double configIceOceanDragCoeff); void gridSolveGPU(); void aggregateDeluDyn(); void applyFreeSlipBC(); @@ -298,6 +302,9 @@ auto Mesh::getMeshField(){ else if constexpr (index==MeshF_VtxRotLat){ return vtxRotLat_; } + else if constexpr (index==MeshF_VtxRotLon){ + return vtxRotLon_; + } else if constexpr (index==MeshF_ElmCenterXYZ){ return elmCenterXYZ_; } diff --git a/src/pmpo_utils.hpp b/src/pmpo_utils.hpp index 84fa98e8..c7ff689d 100644 --- a/src/pmpo_utils.hpp +++ b/src/pmpo_utils.hpp @@ -45,6 +45,9 @@ using DoubleView = Kokkos::View; using IntView = Kokkos::View; using BoolView = Kokkos::View; +//CONSTANTS +inline constexpr double seaiceDensitySeaWater = 1026.0; + class Vec2d { private: vec2d_t coords_; @@ -554,6 +557,24 @@ void lat_lon_from_xyz(double& lat, double& lon, Vec3d& xyz, double r){ lat = Kokkos::asin(xyz[2]/r); } +KOKKOS_INLINE_FUNCTION +double seaice_mpm_wrap_longitude(const double longitude){ + const auto PI = 3.141592653589; + const auto TAU = 2*PI; + const double wrapped = longitude - Kokkos::floor((longitude + PI) / TAU) * TAU; + return wrapped; +} + +KOKKOS_INLINE_FUNCTION +void seaice_mpm_coord_parallel_transport(const double long_start, const double long_end, + const double latitude, double transport[2]){ + double dLon = long_end - long_start; + dLon = seaice_mpm_wrap_longitude(dLon); + const double psi = dLon * Kokkos::sin(latitude); + transport[0] = Kokkos::cos(psi); + transport[1] = Kokkos::sin(psi); +} + }//namespace polyMPO end #endif diff --git a/src/pmpo_wachspressBasis.hpp b/src/pmpo_wachspressBasis.hpp index 051537df..fc37f7a4 100644 --- a/src/pmpo_wachspressBasis.hpp +++ b/src/pmpo_wachspressBasis.hpp @@ -15,10 +15,13 @@ void sphericalInterpolation(MPMesh& mpMesh){ auto vtxCoords = p_mesh->getMeshField(); int numVtxs = p_mesh->getNumVertices(); auto elm2VtxConn = p_mesh->getElm2VtxConn(); + auto vtxRotLon = p_mesh->getMeshField(); auto p_MPs = mpMesh.p_MPs; auto MPsPosition = p_MPs->getPositions(); auto MPsBasis = p_MPs->getData(); + auto curPosRotLatLon = p_MPs->getData(); + auto MPsAppID = p_MPs->getData(); constexpr MaterialPointSlice mpfIndex = meshFieldIndexToMPSlice; auto mpField = p_MPs->getData(); @@ -26,14 +29,34 @@ void sphericalInterpolation(MPMesh& mpMesh){ const int numEntries = mpSliceToNumEntries(); auto meshField = p_mesh->getMeshField(); + bool use_correction_term = false; + if constexpr (meshFieldIndex == MeshF_OnSurfVeloIncr) { + use_correction_term = true; + } + auto interpolation = PS_LAMBDA(const int& elm, const int& mp, const int& mask) { if(mask) { //if material point is 'active'/'enabled' int numVtx = elm2VtxConn(elm,0); - for(int entry=0; entry Date: Tue, 14 Jul 2026 16:55:50 -0400 Subject: [PATCH 05/11] Reconstrucntion of velocity test commented out as now in src involves spherical corrections --- test/testFortranMPReconstruction.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/testFortranMPReconstruction.f90 b/test/testFortranMPReconstruction.f90 index c98f1b88..76b4b80f 100644 --- a/test/testFortranMPReconstruction.f90 +++ b/test/testFortranMPReconstruction.f90 @@ -127,14 +127,14 @@ program main !Test vtx order 1 reconstruction call polympo_reconstruct_coeff_with_MPI(mpmesh) call polympo_setReconstructionOfMass(mpMesh,1,polympo_getMeshFVtxType()) - call polympo_setReconstructionOfVel(mpMesh, 1, polympo_getMeshFVtxType()) + !call polympo_setReconstructionOfVel(mpMesh, 1, polympo_getMeshFVtxType()) call polympo_applyReconstruction(mpMesh) call polympo_getMeshVtxMass(mpMesh,nVertices,c_loc(meshVtxMass1)) call polympo_getMeshVtxVel(mpMesh, nVertices, c_loc(meshVtxVelu), c_loc(meshVtxVelv)) do i = 1, nVertices call assert(meshVtxMass1(i) < TEST_VAL+TOLERANCE1 .and. meshVtxMass1(i) > TEST_VAL-TOLERANCE1, "Error: wrong vtx mass order 1") - call assert(meshVtxVelu(i) < TEST_VAL+TOLERANCE1 .and. meshVtxVelu(i) > TEST_VAL-TOLERANCE1, "Error: wrong vtx velU order 1") - call assert(meshVtxVelv(i) < TEST_VAL+TOLERANCE1 .and. meshVtxVelv(i) > TEST_VAL-TOLERANCE1, "Error: wrong vtx velV order 1") + !call assert(meshVtxVelu(i) < TEST_VAL+TOLERANCE1 .and. meshVtxVelu(i) > TEST_VAL-TOLERANCE1, "Error: wrong vtx velU order 1") + !call assert(meshVtxVelv(i) < TEST_VAL+TOLERANCE1 .and. meshVtxVelv(i) > TEST_VAL-TOLERANCE1, "Error: wrong vtx velV order 1") end do ! Test vtx order 0 reconstruction From 215764a6ca65e6975a2ef96c80b28f2b980c981e Mon Sep 17 00:00:00 2001 From: Nath Date: Wed, 15 Jul 2026 18:02:29 -0400 Subject: [PATCH 06/11] More formatting --- src/pmpo_MPMesh.cpp | 22 +++++++++++----------- src/pmpo_c.cpp | 4 ++-- src/pmpo_fortran.f90 | 16 ++++++++-------- src/pmpo_materialPoints.hpp | 2 +- 4 files changed, 22 insertions(+), 22 deletions(-) diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index 364fc982..53a5e343 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -31,7 +31,7 @@ void MPMesh::calculateStrain(){ MPsStrainRate(mp, 2) = 0.0; return; } - + int numVtx = elm2VtxConn(elm,0); double v11 = 0.0; @@ -40,7 +40,7 @@ void MPMesh::calculateStrain(){ double v22 = 0.0; double uTanOverR = 0.0; double vTanOverR = 0.0; - + for (int i = 0; i < numVtx; i++){ int iVertex = elm2VtxConn(elm, i+1)-1; v11 = v11 + MPsBasisGrads(mp, i*2 + 0) * velField(iVertex, 0); @@ -85,7 +85,7 @@ void MPMesh::calculateStress(const int constitutive_relation){ constitutive_evp(strain_rate, stress, MPsIcePressure(mp,0), rep_pressure, MPsArea(mp,0), elasticTimeStep, dampingTimescale); else if(constitutive_relation == 3) constitutive_linear(strain_rate, stress); - + for (int m=0 ; m<3; m++) MPsStress(mp, m) = stress[m]*solveStress(elm); MPsRepPressure(mp,0)=rep_pressure; @@ -152,7 +152,7 @@ void MPMesh::calculateStressDivergence(){ (1.0 - ramp) * invM * w_vtx; factor = factor * tanLatVertexRotatedOverRadius(vID, 0); - + auto factor1 = ramp * (w_vtx/radius) * (VtxCoeffs_new(vID, 1, 0) + VtxCoeffs_new(vID, 1, 1)*CoordDiffs[1] + VtxCoeffs_new(vID, 1, 2)*CoordDiffs[2] + VtxCoeffs_new(vID, 1, 3)*CoordDiffs[3]) - @@ -287,7 +287,7 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){ Vec3d dx = MPnew-MP; while(true){ int numConnElms = elm2ElmConn(iElm,0); - + Vec3d center(elmCenter(iElm, 0), elmCenter(iElm, 1), elmCenter(iElm, 2)); Vec3d delta = MPnew - center; @@ -749,8 +749,8 @@ void MPMesh::T2LTracking(Vec2dView dx){ Vec2d MP(mpPositions(mp,0),mpPositions(mp,1));//XXX:the input is XYZ, but we only support 2d vector if(mask){ int iElm = elm; - Vec2d MPnew = MP + dx(mp); - + Vec2d MPnew = MP + dx(mp); + while(true){ int numVtx = elm2VtxConn(iElm,0); bool goToNeighbour = false; @@ -760,7 +760,7 @@ void MPMesh::T2LTracking(Vec2dView dx){ v[i] = elm2VtxConn(iElm,i+1)-1; //get edges and perpendiculardx Vec2d e[maxVtxsPerElm]; - double pdx[maxVtxsPerElm]; + double pdx[maxVtxsPerElm]; for(int i=0; i< numVtx; i++){ int idx_ip1 = (i+1)%numVtx; Vec2d v_i(vtxCoords(v[i],0),vtxCoords(v[i],1)); @@ -768,17 +768,17 @@ void MPMesh::T2LTracking(Vec2dView dx){ e[i] = v_ip1 - v_i; pdx[i] = (v_i - MP).cross(dx(mp)); } - + for(int i=0; ip_mesh; @@ -1873,7 +1873,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/src/pmpo_fortran.f90 b/src/pmpo_fortran.f90 index 709b895e..98b5f646 100644 --- a/src/pmpo_fortran.f90 +++ b/src/pmpo_fortran.f90 @@ -369,7 +369,7 @@ subroutine polympo_calculateMPStrainRate(mpMesh) & use :: iso_c_binding type(c_ptr), value :: mpMesh end subroutine - + subroutine polympo_setMPStrainRate(mpMesh, nComps, numMPs, array) & bind(C, NAME='polympo_setMPStrainRate_f') use :: iso_c_binding @@ -386,7 +386,7 @@ subroutine polympo_getMPStrainRate(mpMesh, nComps, numMPs, array) & type(c_ptr), value :: array end subroutine - + !MP Stress subroutine polympo_calculateMPStress(mpMesh, constitutive_model) & bind(C, NAME='polympo_calculateMPStress_f') @@ -1032,7 +1032,7 @@ subroutine polympo_setDynamicTimeStep(mpMesh, dynamicTimeStep) & type(c_ptr), value :: mpMesh real(c_double), value :: dynamicTimeStep end subroutine - + subroutine polympo_setSolveStressMesh(mpMesh, nCells, array) & bind(C, NAME='polympo_setSolveStressMesh_f') use :: iso_c_binding @@ -1086,7 +1086,7 @@ subroutine polympo_set_airStress(mpMesh, nVertices, uArray, vArray) & integer(c_int), value :: nVertices type(c_ptr), value :: uArray, vArray end subroutine - + subroutine polympo_set_surfaceTiltForce(mpMesh, nVertices, uArray, vArray) & bind(C, NAME='polympo_set_surfaceTiltForce_f') use :: iso_c_binding @@ -1144,13 +1144,13 @@ subroutine polympo_set_free_slip_bc(mpMesh) & bind(C, NAME='polympo_set_free_slip_bc_f') use :: iso_c_binding type(c_ptr), value :: mpMesh - end subroutine + end subroutine subroutine polympo_set_halo_vel_from_owner(mpMesh) & bind(C, NAME='polympo_set_halo_vel_from_owner_f') use :: iso_c_binding type(c_ptr), value :: mpMesh - end subroutine + end subroutine !--------------------------------------------------------------------------- !> @brief calculate the MPs from given mesh vertices rotational latitude @@ -1163,13 +1163,13 @@ subroutine polympo_push(mpMesh) & use :: iso_c_binding type(c_ptr), value :: mpMesh end subroutine - + subroutine polympo_push_ahead(mpMesh) & bind(C, NAME='polympo_push_ahead_f') use :: iso_c_binding type(c_ptr), value :: mpMesh end subroutine - + !--------------------------------------------------------------------------- !> @brief calculate the MPs from given mesh vertices rotational latitude !--------------------------------------------------------------------------- diff --git a/src/pmpo_materialPoints.hpp b/src/pmpo_materialPoints.hpp index 6a48615f..c47f911b 100644 --- a/src/pmpo_materialPoints.hpp +++ b/src/pmpo_materialPoints.hpp @@ -148,7 +148,7 @@ class MaterialPoints { void rebuild(IntView addedMP2elm, IntView addedMPAppID); void startRebuild(IntView tgtElm, int addedNumMPs, IntView addedMP2elm, IntView addedMPAppID, Kokkos::View addedMPMask); void startRebuild(IntView tgtElm, int addedNumMPs, IntView addedMP2elm, IntView addedMPAppID); - + void finishRebuild(); bool rebuildOngoing(); From ab4fa2d6ca8e35f2ecc4b40ebc96c754b82be457 Mon Sep 17 00:00:00 2001 From: Shahrear Jahan Santho Date: Mon, 29 Jun 2026 06:56:17 -0700 Subject: [PATCH 07/11] 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 | 4 +- 6 files changed, 737 insertions(+), 98 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 53a5e343..41fcb8c8 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -178,7 +178,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 e4099d54..62f98912 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()); @@ -399,7 +399,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 ae914492..4fbfd4d1 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -1744,8 +1744,8 @@ 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); + + mpMesh->communicate_and_take_halo_contributions1_improved(vtxFieldVel, numVertices, 2, 1, 1); } //Advection Calcualtions From 476b5809c8af8b366bbc6d76cc13ad76825cbaa3 Mon Sep 17 00:00:00 2001 From: Shahrear Jahan Santho Date: Sat, 11 Jul 2026 21:59:58 -0700 Subject: [PATCH 08/11] Add CUDA-aware communication path to MPMesh --- src/CMakeLists.txt | 18 ++ src/pmpo_MPMesh.hpp | 366 ++++++++++++++++++++++++----------- src/pmpo_MPMesh_assembly.hpp | 15 +- 3 files changed, 284 insertions(+), 115 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index d0fb8e74..5eee8a31 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -21,6 +21,24 @@ set(SOURCES ) add_library(polyMPO-core ${SOURCES}) + +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 CUDA_AWARE_MPI) set_property(TARGET polyMPO-core PROPERTY CXX_STANDARD "17") set_property(TARGET polyMPO-core PROPERTY CXX_STANDARD_REQUIRED ON) diff --git a/src/pmpo_MPMesh.hpp b/src/pmpo_MPMesh.hpp index dfd922e5..094289de 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{ @@ -484,10 +489,106 @@ class MPMesh{ #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. + // Use explicit CUDA memory space for MPI device buffers when CUDA is + // enabled. This avoids ambiguity in the default Kokkos::View memory space + // and gives Cray MPICH/GTL plain CUDA allocations to register/export. +#ifdef KOKKOS_ENABLE_CUDA + // Plain cudaMalloc'd device buffer, exposed as an unmanaged Kokkos::View + // via .view(). Used only for the 4 buffers below that get handed + // directly to MPI_Isend/Irecv under CUDA-aware MPI. + // + // Why not just a Kokkos::View: when Kokkos is + // built with Kokkos_ENABLE_IMPL_CUDA_MALLOC_ASYNC=ON (the default since + // Kokkos 4.2), View allocations use cudaMallocAsync/memory pools. + // cuIpcGetMemHandle (which Cray MPICH/GTL uses for intra-node GPU-to-GPU + // sends) rejects pool allocations with CUDA_ERROR_INVALID_VALUE. Since + // polyMPO isn't allowed to touch the Kokkos build config, these 4 + // buffers bypass Kokkos's allocator entirely via a direct cudaMalloc, + // which cuIpcGetMemHandle always accepts, regardless of how the rest of + // Kokkos (or the rest of the app's Views) is configured. + 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; + } + + 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 + // No CUDA backend: no CUDA IPC/pool-allocation concern, so just wrap a + // normal Kokkos::View with the same .allocate()/.data()/.view() + // interface as RawCudaMPIBuffer above, so the cache struct and its call + // sites below don't need to branch on KOKKOS_ENABLE_CUDA. + 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 batched GPU buffers. + // Every neighbor's data lives in one shared allocation (see + // CudaAwareMPIFieldCache below) and MPI is given "base pointer + byte + // offset" per neighbor rather than a separate allocation per neighbor. + // If you ever need to fall back to one allocation per neighbor (e.g. an + // MPI/GPU stack that mishandles offset device pointers for CUDA IPC), + // restore the per-proc-buffer version from version control. bool cudaAwareMPICacheValid = true; bool cudaAwareMPIDisabled = false; bool cudaAwareMPIEnvChecked = false; @@ -500,12 +601,23 @@ 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; + + // Single batched GPU buffers (one allocation each, instead of one + // Kokkos::View per neighbor proc). Per-proc slices are + // [offset, offset + count) for the ID buffers, and + // [offset * numEntries, (offset + count) * numEntries) for the data + // buffers. MPI is given "buffer base pointer + offset", not a + // separate allocation per proc. + CudaAwareMPIIntBuffer sendEntityGPU; + CudaAwareMPIIntBuffer recvIDGPU; + + CudaAwareMPIDoubleBuffer sendDataGPU; + CudaAwareMPIDoubleBuffer recvDataGPU; }; std::map, CudaAwareMPIFieldCache> cudaAwareMPICaches; @@ -521,10 +633,30 @@ class MPMesh{ return cudaAwareMPIForceCPU; } - // Fully CUDA-aware MPI version: + // Fully CUDA-aware MPI version, batched buffer variant: // 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. // + // Every neighbor's send/recv entity-ID list and data live in ONE big + // GPU buffer each (laid out back-to-back in proc order), instead of one + // Kokkos::View allocation per neighbor. Packing/unpacking is a single + // kernel launch over all neighbors' entities at once instead of one + // launch per neighbor, and MPI_Isend/Irecv use "buffer base pointer + + // offset" into that single buffer per proc. This is what actually + // shrinks MPI_Wait time: fewer, larger, more uniform in-flight + // transfers instead of many small independent ones. + // + // Note: an earlier version of this cache used one Kokkos::View + // allocation per neighbor specifically to avoid handing MPI a + // "base pointer + offset" GPU address, out of concern for CUDA IPC + // issues on Cray MPICH/GTL. That failure mode was root-caused to + // Kokkos allocating device Views via cudaMallocAsync (invalid for + // cuIpcGetMemHandle), not to offset pointers themselves, and is fixed + // by building Kokkos with -DKokkos_ENABLE_IMPL_CUDA_MALLOC_ASYNC=OFF. + // If you ever do hit IPC trouble that tracks back to offset pointers + // specifically, the per-proc-buffer version can be restored from + // version control. + // // Important: // This function caches communication metadata and GPU buffers per // (mode, numEntries). If the communication pattern changes, clear @@ -561,7 +693,7 @@ class MPMesh{ #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()" + << "[CUDA_AWARE_MPI] Using batched single-buffer GPU-aware MPI path in communicate_and_take_halo_contributions1_improved()" << "\n"; cudaAwareMPILogged = true; } @@ -587,16 +719,8 @@ class MPMesh{ 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; @@ -611,109 +735,124 @@ class MPMesh{ } } + int totalSend = 0; + int totalRecv = 0; + for(int proc = 0; proc < numProcsTot; proc++){ - if(proc == self) continue; + cudaAwareCache.sendOffsets[proc] = totalSend; + totalSend += cudaAwareCache.sendCounts[proc]; + + cudaAwareCache.recvOffsets[proc] = totalRecv; + totalRecv += cudaAwareCache.recvCounts[proc]; + } - const int sendCount = cudaAwareCache.sendCounts[proc]; - const int recvCount = cudaAwareCache.recvCounts[proc]; + cudaAwareCache.totalSendCount = totalSend; + cudaAwareCache.totalRecvCount = totalRecv; - if(sendCount > 0){ - cudaAwareCache.sendEntityGPUPerProc[proc] = - Kokkos::View( - "cudaAwareMPISendEntityGPUPerProc", - sendCount); + cudaAwareCache.sendEntityGPU.allocate(totalSend); + cudaAwareCache.sendDataGPU.allocate(totalSend * numEntries); + cudaAwareCache.recvIDGPU.allocate(totalRecv); + cudaAwareCache.recvDataGPU.allocate(totalRecv * numEntries); - cudaAwareCache.sendDataGPUPerProc[proc] = - Kokkos::View( - "cudaAwareMPISendDataGPUPerProc", - sendCount * numEntries); + // ---- Build the flattened send-entity list (host, then one deep_copy) ---- + if(totalSend > 0){ + auto sendEntityCPU = + Kokkos::View( + "sendEntityCPU_batched", totalSend); - auto sendEntityCPU = - Kokkos::View( - "sendEntityCPU", - sendCount); + if(mode == 0){ + for(int proc = 0; proc < numProcsTot; proc++){ + if(proc == self) continue; + if(cudaAwareCache.sendCounts[proc] <= 0) continue; - if(mode == 0){ assert(haloOwnerLocalIDs[proc].size() == - static_cast(sendCount)); + static_cast(cudaAwareCache.sendCounts[proc])); + } - int localIndex = 0; + std::vector cursor(cudaAwareCache.sendOffsets); - for(int iEnt = 0; iEnt < numHalosTot; iEnt++){ - int ownerProc = haloOwnerProcs[iEnt]; + for(int iEnt = 0; iEnt < numHalosTot; iEnt++){ + int ownerProc = haloOwnerProcs[iEnt]; + if(ownerProc == self) continue; - if(ownerProc != proc) continue; + sendEntityCPU(cursor[ownerProc]) = numOwnersTot + iEnt; + cursor[ownerProc]++; + } - assert(localIndex < sendCount); + for(int proc = 0; proc < numProcsTot; proc++){ + if(proc == self) continue; - sendEntityCPU(localIndex) = numOwnersTot + iEnt; + assert(cursor[proc] == + cudaAwareCache.sendOffsets[proc] + + cudaAwareCache.sendCounts[proc]); + } + } + else{ + for(int proc = 0; proc < numProcsTot; proc++){ + if(proc == self) continue; - localIndex++; - } + int sendCount = cudaAwareCache.sendCounts[proc]; + if(sendCount <= 0) continue; - assert(localIndex == sendCount); - } - else{ assert(ownerOwnerLocalIDs[proc].size() == static_cast(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); } - if(recvCount > 0){ - cudaAwareCache.recvIDGPUPerProc[proc] = - Kokkos::View( - "cudaAwareMPIRecvIDGPUPerProc", - recvCount); + Kokkos::deep_copy(cudaAwareCache.sendEntityGPU.view(), sendEntityCPU); + } - cudaAwareCache.recvDataGPUPerProc[proc] = - Kokkos::View( - "cudaAwareMPIRecvDataGPUPerProc", - recvCount * numEntries); + // ---- Build the flattened recv-ID list (host, then one deep_copy) ---- + if(totalRecv > 0){ + auto recvIDCPU = + Kokkos::View( + "recvIDCPU_batched", totalRecv); - auto recvIDCPU = - Kokkos::View( - "recvIDCPU", - recvCount); + if(mode == 0){ + for(int proc = 0; proc < numProcsTot; proc++){ + if(proc == self) continue; + + int recvCount = cudaAwareCache.recvCounts[proc]; + if(recvCount <= 0) continue; - if(mode == 0){ assert(ownerOwnerLocalIDs[proc].size() == static_cast(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; + } + else{ + std::vector cursor(cudaAwareCache.recvOffsets); - assert(localIndex < recvCount); + for(int iEnt = 0; iEnt < numHalosTot; iEnt++){ + int ownerProc = haloOwnerProcs[iEnt]; + if(ownerProc == self) continue; - recvIDCPU(localIndex) = numOwnersTot + iEnt; + recvIDCPU(cursor[ownerProc]) = numOwnersTot + iEnt; + cursor[ownerProc]++; + } - localIndex++; - } + for(int proc = 0; proc < numProcsTot; proc++){ + if(proc == self) continue; - assert(localIndex == recvCount); + assert(cursor[proc] == + cudaAwareCache.recvOffsets[proc] + + cudaAwareCache.recvCounts[proc]); } - - Kokkos::deep_copy( - cudaAwareCache.recvIDGPUPerProc[proc], - recvIDCPU); } - } - Kokkos::fence(); + Kokkos::deep_copy(cudaAwareCache.recvIDGPU.view(), recvIDCPU); + } cudaAwareCache.valid = true; @@ -725,17 +864,14 @@ class MPMesh{ 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]; + // ---- Pack: ONE kernel over all neighbors' send entities at once ---- + if(cudaAwareCache.totalSendCount > 0){ + auto sendEntityGPU = cudaAwareCache.sendEntityGPU.view(); + auto sendDataGPU = cudaAwareCache.sendDataGPU.view(); Kokkos::parallel_for( - "pack cached cuda-aware mpi send buffer per proc", - sendCount, + "pack cached cuda-aware mpi send buffer batched", + cudaAwareCache.totalSendCount, KOKKOS_LAMBDA(const int i){ int entity = sendEntityGPU(i); @@ -766,8 +902,13 @@ class MPMesh{ if(cudaAwareCache.recvCounts[proc] > 0){ MPI_Request reqData; + double* recvPtr = + cudaAwareCache.recvDataGPU.data() + + static_cast(cudaAwareCache.recvOffsets[proc]) * + numEntries; + mpiError = MPI_Irecv( - cudaAwareCache.recvDataGPUPerProc[proc].data(), + recvPtr, cudaAwareCache.recvCounts[proc] * numEntries, MPI_DOUBLE, proc, @@ -781,8 +922,13 @@ class MPMesh{ if(cudaAwareCache.sendCounts[proc] > 0){ MPI_Request reqData; + double* sendPtr = + cudaAwareCache.sendDataGPU.data() + + static_cast(cudaAwareCache.sendOffsets[proc]) * + numEntries; + mpiError = MPI_Isend( - cudaAwareCache.sendDataGPUPerProc[proc].data(), + sendPtr, cudaAwareCache.sendCounts[proc] * numEntries, MPI_DOUBLE, proc, @@ -818,7 +964,7 @@ class MPMesh{ if(self == 0){ std::cout - << "[CUDA_AWARE_MPI] Device-pointer MPI failed." + << "[CUDA_AWARE_MPI] Batched device-pointer MPI failed." << std::endl; } @@ -838,6 +984,8 @@ class MPMesh{ return; } + timer.reset(); + if(self == 0){ std::cout << "[CUDA_AWARE_MPI] Failure happened after MPI requests were posted. " @@ -856,18 +1004,15 @@ class MPMesh{ 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: ONE kernel over all neighbors' recv entities at once ---- + 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, + "halo add cached cuda-aware mpi batched", + cudaAwareCache.totalRecvCount, KOKKOS_LAMBDA(const int i){ const int vertex = recvIDGPU(i); @@ -885,8 +1030,8 @@ class MPMesh{ } else{ Kokkos::parallel_for( - "halo assign cached cuda-aware mpi per proc", - recvCount, + "halo assign cached cuda-aware mpi batched", + cudaAwareCache.totalRecvCount, KOKKOS_LAMBDA(const int i){ const int vertex = recvIDGPU(i); @@ -933,4 +1078,3 @@ class MPMesh{ }//namespace polyMPO end #endif - diff --git a/src/pmpo_MPMesh_assembly.hpp b/src/pmpo_MPMesh_assembly.hpp index 62f98912..7c4120ce 100644 --- a/src/pmpo_MPMesh_assembly.hpp +++ b/src/pmpo_MPMesh_assembly.hpp @@ -130,6 +130,7 @@ void MPMesh::reconstruct_coeff_full(){ radius=p_mesh->getSphereRadius(); //Assemble matrix for each vertex + timer.reset(); auto assemble = PS_LAMBDA(const int& elm, const int& mp, const int& mask) { if(mask) { //if material point is 'active'/'enabled' int nVtxE = elm2VtxConn(elm,0); //number of vertices bounding the element @@ -152,7 +153,7 @@ void MPMesh::reconstruct_coeff_full(){ }; p_MPs->parallel_for(assemble, "assembly"); Kokkos::fence(); - pumipic::RecordTime("Assemble Matrix Per Process" + std::to_string(self), timer.seconds()); + pumipic::RecordTime("VR 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 @@ -162,20 +163,26 @@ void MPMesh::reconstruct_coeff_full(){ int op = 0; if (numProcsTot >1){ communicate_and_take_halo_contributions1_improved(vtxMatrices, numVertices, numEntriesMatrix, mode, op); + pumipic::RecordTime("VR Matrix Gather Halo/MPI " + std::to_string(self), timer.seconds()); + + timer.reset(); mode=1; op=1; communicate_and_take_halo_contributions1_improved(vtxMatrices, numVertices, numEntriesMatrix, mode, op); + pumipic::RecordTime("VR Matrix Scatter Halo/MPI " + std::to_string(self), timer.seconds()); } - 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; + timer.reset(); invertMatrix(vtxMatrices, radius); + pumipic::RecordTime("VR Invert Matrix " + std::to_string(self), timer.seconds()); } void MPMesh::invertMatrix(const Kokkos::View& vtxMatrices, const double& radius){ From 4d6dcda95defeff52aa03e92d487080f95ec9419 Mon Sep 17 00:00:00 2001 From: Shahrear Jahan Santho Date: Thu, 30 Jul 2026 15:35:41 -0500 Subject: [PATCH 09/11] Cuda Aware MPI Halo Exchange-Diagnostics lables --- src/pmpo_MPMesh.cpp | 42 ++++++++++-- src/pmpo_MPMesh.hpp | 102 +++++++++++++++++++--------- src/pmpo_MPMesh_assembly.hpp | 124 ++++++++++++++++++++++++++++++----- src/pmpo_c.cpp | 35 +++++++++- 4 files changed, 245 insertions(+), 58 deletions(-) diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index 41fcb8c8..3f9451e9 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -96,12 +96,17 @@ 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); MPI_Comm_size(comm, &numProcsTot); + Kokkos::fence(); //B0: drain any device work left from the previous phase + MPI_Barrier(comm); //B0: align all ranks so this timed region starts at the same instant + + Kokkos::Timer totalTimer; + Kokkos::Timer computeTimer; + //Mesh Information auto elm2VtxConn = p_mesh->getElm2VtxConn(); int numVtxOwned = p_mesh->getNumVerticesOwned(); @@ -172,18 +177,41 @@ 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()); + Kokkos::fence(); //drain device work -> compute really is done on this rank + + const double computeTime = computeTimer.seconds(); //T_before: pure local compute time, no waiting + + MPI_Barrier(comm); //B1: fast ranks wait here for the slowest rank + + const double computeTimeSync = computeTimer.seconds(); //time until every rank reached the barrier + const double computeImbalance = computeTimeSync - computeTime; //this rank's wait time = compute load imbalance + + Kokkos::Timer communicationTimer; - timer.reset(); 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_contributions1_improved(stress_divUV, numVertices, 2, 0, 0, "Stress_Divergence"); //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()); + Kokkos::fence(); //drain device work from the communication step + + const double communicationTime = communicationTimer.seconds(); //pure local communication time, no waiting + + MPI_Barrier(comm); //B2: fast ranks wait here for the slowest rank + + const double communicationTimeSync = communicationTimer.seconds(); + const double communicationImbalance = communicationTimeSync - communicationTime; //communication load imbalance + + const double totalTime = totalTimer.seconds(); + + pumipic::RecordTime("Stress_Divergence_Compute_" + std::to_string(self),computeTime); + pumipic::RecordTime("Stress_Divergence_Compute_Imbalance_" + std::to_string(self),computeImbalance); + + pumipic::RecordTime("Stress_Divergence_Communication_" + std::to_string(self),communicationTime); + pumipic::RecordTime("Stress_Divergence_Communication_Imbalance_" + std::to_string(self),communicationImbalance); + + pumipic::RecordTime("Stress_Divergence_Total_" + std::to_string(self),totalTime); } void MPMesh::calcBasis() { diff --git a/src/pmpo_MPMesh.hpp b/src/pmpo_MPMesh.hpp index 094289de..537ea1f4 100644 --- a/src/pmpo_MPMesh.hpp +++ b/src/pmpo_MPMesh.hpp @@ -667,7 +667,8 @@ class MPMesh{ int nEntities, int numEntries, int mode, - int op){ + int op, + const std::string& label){ int self, numProcsTot; @@ -676,6 +677,14 @@ class MPMesh{ MPI_Comm_rank(comm, &self); MPI_Comm_size(comm, &numProcsTot); + const char* diagnosticsEnv = + + std::getenv("POLYMPO_MPI_DIAGNOSTICS"); + + const bool mpiDiagnostics = + diagnosticsEnv != nullptr && + std::atoi(diagnosticsEnv) != 0; + assert(mode == 0 || mode == 1); assert(op == 0 || op == 1); assert(nEntities == numOwnersTot + numHalosTot); @@ -856,10 +865,10 @@ class MPMesh{ 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()); + if(mpiDiagnostics){ + pumipic::RecordTime(label + "_MPI_Diagnostics_CacheBuild_m" + std::to_string(mode) + "_e" + std::to_string(numEntries) + + "_rank" + std::to_string(self), timer.seconds()); + } timer.reset(); } @@ -884,20 +893,41 @@ class MPMesh{ 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()); + if(mpiDiagnostics){ + pumipic::RecordTime( + label + "_MPI_Diagnostics_Pack_m" + std::to_string(mode) + "_e" + std::to_string(numEntries) + "_rank" + std::to_string(self), timer.seconds()); + } timer.reset(); - Kokkos::Timer mpiTotalTimer; + double postTime = 0.0; + double waitallTime = 0.0; + std::vector requests; requests.reserve(2 * numProcsTot); int mpiError = MPI_SUCCESS; + // Data volume exchanged by this rank in this call (recorded once per + // call, independent of the post/wait timing below), tagged by the + // caller (label) so SD, VR, and Reconstruction can be told apart. + const double bytesSent = + static_cast(cudaAwareCache.totalSendCount) * numEntries * sizeof(double); + const double bytesRecv = + static_cast(cudaAwareCache.totalRecvCount) * numEntries * sizeof(double); + + pumipic::RecordTime(label + "_MPI_BytesSent_" + std::to_string(self), bytesSent); + pumipic::RecordTime(label + "_MPI_BytesRecv_" + std::to_string(self), bytesRecv); + + //Post both the Irecv and the matching Isend for a proc together, in + //the same loop iteration and under their own counts (recvCounts for + //Irecv, sendCounts for Isend). This replaces the previous two-pass + //version (a first loop that posted Irecv only, plus a second loop + //that posted Isend only) which existed only because of a leftover, + //commented-out duplicate of this same block. + int numNeighbors = 0; for(int proc = 0; proc < numProcsTot; proc++){ if(proc == self) continue; + bool hasComm = false; if(cudaAwareCache.recvCounts[proc] > 0){ MPI_Request reqData; @@ -917,6 +947,7 @@ class MPMesh{ &reqData); if(mpiError != MPI_SUCCESS) break; requests.push_back(reqData); + hasComm = true; } if(cudaAwareCache.sendCounts[proc] > 0){ @@ -937,13 +968,28 @@ class MPMesh{ &reqData); if(mpiError != MPI_SUCCESS) break; requests.push_back(reqData); + hasComm = true; } + if(hasComm) numNeighbors++; } + pumipic::RecordTime(label + "_MPI_NumNeighbors_" + std::to_string(self), static_cast(numNeighbors)); + + postTime = timer.seconds(); - pumipic::RecordTime( - "SD: CUDA-aware MPI Post m" + std::to_string(mode) + - " e" + std::to_string(numEntries) + "-" + std::to_string(self), - timer.seconds()); + pumipic::RecordTime(label + "_MPI_Post_" + std::to_string(self), postTime); + + //Barrier here, not before posting: Isend/Irecv are non-blocking and + //their cost is local (looping + building MPI_Request objects), so a + //barrier before posting would only be measuring how skewed ranks + //were on entry to this function, which the calling function's own + //barriers already capture. Placed here, right before Waitall, it + //makes every rank enter Waitall at the same instant, so the + //waitallTime below reflects real message-arrival/network imbalance + //instead of being contaminated by skew left over from posting. + Kokkos::Timer barrierTimer; + MPI_Barrier(comm); + const double barrierWaitTime = barrierTimer.seconds(); + pumipic::RecordTime(label + "_MPI_BarrierWait_" + std::to_string(self), barrierWaitTime); timer.reset(); @@ -954,10 +1000,9 @@ class MPMesh{ 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()); + waitallTime = timer.seconds(); + + pumipic::RecordTime(label + "_MPI_Waitall_" + std::to_string(self), waitallTime); if(mpiError != MPI_SUCCESS){ cudaAwareMPIDisabled = true; @@ -984,8 +1029,6 @@ class MPMesh{ return; } - timer.reset(); - if(self == 0){ std::cout << "[CUDA_AWARE_MPI] Failure happened after MPI requests were posted. " @@ -997,11 +1040,6 @@ class MPMesh{ 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(); // ---- Unpack: ONE kernel over all neighbors' recv entities at once ---- @@ -1045,10 +1083,12 @@ class MPMesh{ 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()); + if(mpiDiagnostics){ + pumipic::RecordTime(label + "_MPI_Diagnostics_Contribution_m" + std::to_string(mode) + "_e" + + std::to_string(numEntries) + "_rank" + std::to_string(self), timer.seconds()); + } + + } #else @@ -1061,7 +1101,9 @@ class MPMesh{ int nEntities, int numEntries, int mode, - int op){ + int op, + const std::string& label){ + (void)label; // no per-call diagnostics on the CPU-staged fallback path communicate_and_take_halo_contributions1( meshField, diff --git a/src/pmpo_MPMesh_assembly.hpp b/src/pmpo_MPMesh_assembly.hpp index 7c4120ce..8e6189a6 100644 --- a/src/pmpo_MPMesh_assembly.hpp +++ b/src/pmpo_MPMesh_assembly.hpp @@ -96,12 +96,27 @@ 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); - + + Kokkos::fence(); //B0: drain any device work left from the previous phase + MPI_Barrier(comm); //B0: align all ranks so this timed region starts at the same instant + + Kokkos::Timer totalTimer; + Kokkos::Timer phaseTimer; + + double preCommunicationComputeTime = 0.0; + double preCommunicationComputeImbalance = 0.0; + double gatherCommunicationTime = 0.0; + double gatherCommunicationImbalance = 0.0; + double scatterCommunicationTime = 0.0; + double scatterCommunicationImbalance = 0.0; + double postCommunicationComputeTime = 0.0; + double postCommunicationComputeImbalance = 0.0; + static int coeff_count=0; if(!self) std::cout<<"===="<<__FUNCTION__<<" "<getMeshField(); //Material Points + phaseTimer.reset(); + calcBasis(); auto weight = p_MPs->getData(); @@ -130,7 +147,6 @@ void MPMesh::reconstruct_coeff_full(){ radius=p_mesh->getSphereRadius(); //Assemble matrix for each vertex - timer.reset(); auto assemble = PS_LAMBDA(const int& elm, const int& mp, const int& mask) { if(mask) { //if material point is 'active'/'enabled' int nVtxE = elm2VtxConn(elm,0); //number of vertices bounding the element @@ -153,25 +169,40 @@ void MPMesh::reconstruct_coeff_full(){ }; p_MPs->parallel_for(assemble, "assembly"); Kokkos::fence(); - pumipic::RecordTime("VR Assemble Matrix Per Process" + std::to_string(self), timer.seconds()); + + preCommunicationComputeTime = phaseTimer.seconds(); //T_before: pure local compute time, no waiting + MPI_Barrier(comm); //B1: fast ranks wait here for the slowest rank + preCommunicationComputeImbalance = phaseTimer.seconds() - preCommunicationComputeTime; + //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); - pumipic::RecordTime("VR Matrix Gather Halo/MPI " + std::to_string(self), timer.seconds()); - timer.reset(); + phaseTimer.reset(); + communicate_and_take_halo_contributions1_improved(vtxMatrices, numVertices, numEntriesMatrix, mode, op, "Reconstruction_Gather"); + + Kokkos::fence(); + gatherCommunicationTime = phaseTimer.seconds(); //T_before + MPI_Barrier(comm); //B2: fast ranks wait here for the slowest rank + gatherCommunicationImbalance = phaseTimer.seconds() - gatherCommunicationTime; + mode=1; op=1; - communicate_and_take_halo_contributions1_improved(vtxMatrices, numVertices, numEntriesMatrix, mode, op); - pumipic::RecordTime("VR Matrix Scatter Halo/MPI " + std::to_string(self), timer.seconds()); + + phaseTimer.reset(); + communicate_and_take_halo_contributions1_improved(vtxMatrices, numVertices, numEntriesMatrix, mode, op, "Reconstruction_Scatter"); + Kokkos::fence(); + scatterCommunicationTime = phaseTimer.seconds(); //T_before + MPI_Barrier(comm); //B3: fast ranks wait here for the slowest rank + scatterCommunicationImbalance = phaseTimer.seconds() - scatterCommunicationTime; } + phaseTimer.reset(); //Store the 1st matrix element Kokkos::ViewvtxMatrixMass_l("vtxMass", numVertices); Kokkos::parallel_for("storeMatrixMass", numVertices, KOKKOS_LAMBDA(const int vtx){ @@ -180,9 +211,27 @@ void MPMesh::reconstruct_coeff_full(){ Kokkos::fence(); this->vtxMatrixMass = vtxMatrixMass_l; - timer.reset(); invertMatrix(vtxMatrices, radius); - pumipic::RecordTime("VR Invert Matrix " + std::to_string(self), timer.seconds()); + Kokkos::fence(); + postCommunicationComputeTime = phaseTimer.seconds(); //T_before + MPI_Barrier(comm); //B4: fast ranks wait here for the slowest rank + postCommunicationComputeImbalance = phaseTimer.seconds() - postCommunicationComputeTime; + + const double computeTime = preCommunicationComputeTime + postCommunicationComputeTime; + const double communicationTime = gatherCommunicationTime + scatterCommunicationTime; + const double totalTime = totalTimer.seconds(); + + pumipic::RecordTime("Reconstruction_PreComm_Compute_" + std::to_string(self), preCommunicationComputeTime); + pumipic::RecordTime("Reconstruction_PreComm_Compute_Imbalance_" + std::to_string(self), preCommunicationComputeImbalance); + pumipic::RecordTime("Reconstruction_Gather_Communication_" + std::to_string(self),gatherCommunicationTime); + pumipic::RecordTime("Reconstruction_Gather_Communication_Imbalance_" + std::to_string(self),gatherCommunicationImbalance); + pumipic::RecordTime("Reconstruction_Scatter_Communication_" + std::to_string(self), scatterCommunicationTime); + pumipic::RecordTime("Reconstruction_Scatter_Communication_Imbalance_" + std::to_string(self), scatterCommunicationImbalance); + pumipic::RecordTime("Reconstruction_PostComm_Compute_" + std::to_string(self), postCommunicationComputeTime); + pumipic::RecordTime("Reconstruction_PostComm_Compute_Imbalance_" + std::to_string(self), postCommunicationComputeImbalance); + pumipic::RecordTime("Reconstruction_Compute_Total_" + std::to_string(self), computeTime); + pumipic::RecordTime("Reconstruction_Communication_Total_" + std::to_string(self), communicationTime); + pumipic::RecordTime("Reconstruction_Total_" + std::to_string(self), totalTime); } void MPMesh::invertMatrix(const Kokkos::View& vtxMatrices, const double& radius){ @@ -328,13 +377,17 @@ 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(); MPI_Comm_rank(comm, &self); MPI_Comm_size(comm, &numProcsTot); + Kokkos::fence(); //B0: drain any device work left from the previous phase + MPI_Barrier(comm); //B0: align all ranks so this timed region starts at the same instant + + Kokkos::Timer totalTimer; + Kokkos::Timer computeTimer; + auto VtxCoeffs_new=this->precomputedVtxCoeffs_new; //Mesh Information @@ -402,13 +455,48 @@ void MPMesh::assemblyVtx1(){ }; p_MPs->parallel_for(reconstruct, "reconstruct"); Kokkos::fence(); - pumipic::RecordTime("Assemble Field per process" + std::to_string(self), timer.seconds()); + const double computeTime = computeTimer.seconds(); //T_before: pure local compute time, no waiting + + MPI_Barrier(comm); //B1: fast ranks wait here for the slowest rank + const double computeTimeSync = computeTimer.seconds(); + const double computeImbalance = computeTimeSync - computeTime; + + Kokkos::Timer communicationTimer; - timer.reset(); if(numProcsTot>1){ - communicate_and_take_halo_contributions1_improved(meshField, numVertices, numEntries, 0, 0); + communicate_and_take_halo_contributions1_improved(meshField, numVertices, numEntries, 0, 0, "Velocity_Reconstruction"); } - pumipic::RecordTime("Communicate Field Values" + std::to_string(self), timer.seconds()); + Kokkos::fence(); + + const double communicationTime = communicationTimer.seconds(); //T_before + + MPI_Barrier(comm); //B2: fast ranks wait here for the slowest rank + const double communicationTimeSync = communicationTimer.seconds(); + const double communicationImbalance = communicationTimeSync - communicationTime; + + const double totalTime = totalTimer.seconds(); + + if constexpr (meshFieldIndex == MeshF_Vel) { + pumipic::RecordTime( + "Velocity_Reconstruction_Compute_" + std::to_string(self), + computeTime); + pumipic::RecordTime( + "Velocity_Reconstruction_Compute_Imbalance_" + std::to_string(self), + computeImbalance); + + pumipic::RecordTime( + "Velocity_Reconstruction_Communication_" + std::to_string(self), + communicationTime); + pumipic::RecordTime( + "Velocity_Reconstruction_Communication_Imbalance_" + std::to_string(self), + communicationImbalance); + + pumipic::RecordTime( + "Velocity_Reconstruction_Total_" + std::to_string(self), + totalTime); + +} + } template diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index 4fbfd4d1..d0f4f928 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -571,6 +571,9 @@ void polympo_setMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, Kokkos::Timer timer; checkMPMeshValid(p_mpmesh); auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; + int self; + MPI_Comm comm = p_MPs->getMPIComm(); + MPI_Comm_rank(comm, &self); PMT_ALWAYS_ASSERT(nComps == vec2d_nEntries); PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount()); //PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID()); @@ -587,7 +590,7 @@ void polympo_setMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, } }; p_MPs->parallel_for(setMPVel, "setMPVel"); - pumipic::RecordTime("PolyMPO_setMPVel", timer.seconds()); + pumipic::RecordTime("PolyMPO_setMPVel" + std::to_string(self),timer.seconds()); } void polympo_getMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpVelHost) { @@ -1249,6 +1252,10 @@ void polympo_getMeshVtxVel_f(MPMesh_ptr p_mpmesh, const int nVertices, double* u //check mpMesh is valid checkMPMeshValid(p_mpmesh); auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; + auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; + int self; + MPI_Comm comm = p_MPs->getMPIComm(); + MPI_Comm_rank(comm, &self); //check the size PMT_ALWAYS_ASSERT(p_mesh->getNumVertices() == nVertices); @@ -1260,7 +1267,7 @@ void polympo_getMeshVtxVel_f(MPMesh_ptr p_mpmesh, const int nVertices, double* u uVelOut[i] = h_coordsArray(i,0); vVelOut[i] = h_coordsArray(i,1); } - pumipic::RecordTime("PolyMPO_getMeshVtxVel", timer.seconds()); + pumipic::RecordTime("PolyMPO_getMeshVtxVel" + std::to_string(self), timer.seconds()); } void polympo_setMeshVtxMass_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* vtxMass){ @@ -1574,8 +1581,16 @@ void polympo_setIceAreaVertex_f(MPMesh_ptr p_mpmesh, const int nVertices, double void polympo_calculateStressDivergence_f(MPMesh_ptr p_mpmesh){ //chech validity checkMPMeshValid(p_mpmesh); + Kokkos::Timer wrapperTimer; + auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; + int self; + MPI_Comm comm = p_MPs->getMPIComm(); + MPI_Comm_rank(comm, &self); + ((polyMPO::MPMesh*)p_mpmesh) -> calculateStressDivergence(); + pumipic::RecordTime("Wrapper_StressDivergence_Total_" + std::to_string(self), wrapperTimer.seconds()); + } void polympo_getStressDivergence_f(MPMesh_ptr p_mpmesh, const int nVertices, double* uArray, double* vArray){ @@ -1745,7 +1760,7 @@ void polympo_set_halo_vel_from_owner_f(MPMesh_ptr p_mpmesh){ int numVertices = p_mesh->getNumVertices(); auto vtxFieldVel = p_mesh->getMeshField(); - mpMesh->communicate_and_take_halo_contributions1_improved(vtxFieldVel, numVertices, 2, 1, 1); + mpMesh->communicate_and_take_halo_contributions1_improved(vtxFieldVel, numVertices, 2, 1, 1, "halo_vel_from_owner"); } //Advection Calcualtions @@ -1824,7 +1839,14 @@ void polympo_applyReconstruction_f(MPMesh_ptr p_mpmesh){ void polympo_reconstruct_coeff_with_MPI_f(MPMesh_ptr p_mpmesh){ checkMPMeshValid(p_mpmesh); auto mpmesh = ((polyMPO::MPMesh*)p_mpmesh); + Kokkos::Timer wrapperTimer; + int self; + MPI_Comm comm = mpmesh->p_MPs->getMPIComm(); + MPI_Comm_rank(comm, &self); + mpmesh->reconstruct_coeff_full(); + + pumipic::RecordTime("Wrapper_ReconstructCoeff_Total_" + std::to_string(self), wrapperTimer.seconds()); } void polympo_reconstruct_iceArea_with_MPI_f(MPMesh_ptr p_mpmesh){ @@ -1836,7 +1858,14 @@ void polympo_reconstruct_iceArea_with_MPI_f(MPMesh_ptr p_mpmesh){ void polympo_reconstruct_velocity_with_MPI_f(MPMesh_ptr p_mpmesh){ checkMPMeshValid(p_mpmesh); auto mpmesh = ((polyMPO::MPMesh*)p_mpmesh); + Kokkos::Timer wrapperTimer; + int self; + MPI_Comm comm = mpmesh->p_MPs->getMPIComm(); + MPI_Comm_rank(comm, &self); + mpmesh->assemblyVtx1(); + + pumipic::RecordTime("Wrapper_ReconstructVelocity_Total_" + std::to_string(self), wrapperTimer.seconds()); } void polympo_init_deludelvDyn_f(MPMesh_ptr p_mpmesh){ From cac3f1e08ac37c08ba2c891ef0fd01d7ddede0b6 Mon Sep 17 00:00:00 2001 From: Shahrear Jahan Santho Date: Wed, 12 Aug 2026 15:42:52 -0700 Subject: [PATCH 10/11] Cuda Aware MPI Halo Exchange - Optimize: derived-datatype send (pack path removed) + MPI_Waitany incremental unpack --- src/pmpo_MPMesh.cpp | 28 ++- src/pmpo_MPMesh.hpp | 495 ++++++++++++++++++++++++++++---------------- 2 files changed, 332 insertions(+), 191 deletions(-) diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index 3f9451e9..6eb635c3 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -100,9 +100,14 @@ void MPMesh::calculateStressDivergence(){ MPI_Comm comm = p_MPs->getMPIComm(); MPI_Comm_rank(comm, &self); MPI_Comm_size(comm, &numProcsTot); + + //Kokkos::Timer b0Timer; + //Kokkos::fence(); - Kokkos::fence(); //B0: drain any device work left from the previous phase - MPI_Barrier(comm); //B0: align all ranks so this timed region starts at the same instant + //const double b0DrainTime = b0Timer.seconds(); //B0: drain any device work left from the previous phase + //MPI_Barrier(comm); + //const double b0Sync = b0Timer.seconds(); //B0: align all ranks so this timed region starts at the same instant + //const double b0BarrierWait = b0Sync - b0DrainTime; Kokkos::Timer totalTimer; Kokkos::Timer computeTimer; @@ -181,10 +186,10 @@ void MPMesh::calculateStressDivergence(){ const double computeTime = computeTimer.seconds(); //T_before: pure local compute time, no waiting - MPI_Barrier(comm); //B1: fast ranks wait here for the slowest rank + //MPI_Barrier(comm); //B1: fast ranks wait here for the slowest rank - const double computeTimeSync = computeTimer.seconds(); //time until every rank reached the barrier - const double computeImbalance = computeTimeSync - computeTime; //this rank's wait time = compute load imbalance + //const double computeTimeSync = computeTimer.seconds(); //time until every rank reached the barrier + //const double b1BarrierWait = computeTimeSync - computeTime; //this rank's wait time = compute load imbalance Kokkos::Timer communicationTimer; @@ -198,18 +203,21 @@ void MPMesh::calculateStressDivergence(){ const double communicationTime = communicationTimer.seconds(); //pure local communication time, no waiting - MPI_Barrier(comm); //B2: fast ranks wait here for the slowest rank + //MPI_Barrier(comm); //B2: fast ranks wait here for the slowest rank - const double communicationTimeSync = communicationTimer.seconds(); - const double communicationImbalance = communicationTimeSync - communicationTime; //communication load imbalance + //const double communicationTimeSync = communicationTimer.seconds(); + //const double b2BarrierWait = communicationTimeSync - communicationTime; //communication load imbalance const double totalTime = totalTimer.seconds(); + //pumipic::RecordTime("Stress_Divergence_B0_Drain_" + std::to_string(self), b0DrainTime); + //pumipic::RecordTime("Stress_Divergence_B0_BarrierWait_" + std::to_string(self), b0BarrierWait); + pumipic::RecordTime("Stress_Divergence_Compute_" + std::to_string(self),computeTime); - pumipic::RecordTime("Stress_Divergence_Compute_Imbalance_" + std::to_string(self),computeImbalance); + //pumipic::RecordTime("Stress_Divergence_Compute_B1_Barrier_Wait_Time_" + std::to_string(self),b1BarrierWait); pumipic::RecordTime("Stress_Divergence_Communication_" + std::to_string(self),communicationTime); - pumipic::RecordTime("Stress_Divergence_Communication_Imbalance_" + std::to_string(self),communicationImbalance); + //pumipic::RecordTime("Stress_Divergence_Communication_B2_Wait_Time_" + std::to_string(self), b2BarrierWait); pumipic::RecordTime("Stress_Divergence_Total_" + std::to_string(self),totalTime); } diff --git a/src/pmpo_MPMesh.hpp b/src/pmpo_MPMesh.hpp index 537ea1f4..11849ece 100644 --- a/src/pmpo_MPMesh.hpp +++ b/src/pmpo_MPMesh.hpp @@ -13,6 +13,54 @@ #include #endif +// ---------------------------------------------------------------------- +// Halo-exchange optimization: pure derived-datatype send (pack path fully +// removed) + MPI_Waitany incremental unpack. +// +// communicate_and_take_halo_contributions1_improved() (CUDA_AWARE_MPI path +// only - the CPU-staged fallback below is untouched): +// +// 1. Eliminate the pack step / sendDataGPU buffer, unconditionally. Each +// neighbor proc gets its own committed MPI derived datatype +// (MPI_Type_create_hindexed_block over a "rowType" built from +// meshField's own strides) that gathers that proc's rows directly out +// of meshField's GPU memory. MPI_Isend reads straight from the field - +// no pack kernel, no extra buffer, no fallback to a pack path. +// +// 2. Hide the unpack step behind MPI_Waitall latency instead of paying +// for it serially afterward. Recv requests are drained with +// MPI_Waitany instead of MPI_Waitall, and each neighbor's contribution +// is unpacked (async kernel launch, no intermediate fence) the instant +// that neighbor's message lands, so unpacking early arrivals overlaps +// with waiting on stragglers. A single Kokkos::fence() at the end +// guarantees every launched unpack kernel has completed before +// meshField is used downstream. +// +// The receive side still uses a batched contiguous GPU buffer + a real +// unpack kernel (not a derived receive-datatype): for op==0 the unpack is a +// scatter-ADD (Kokkos::atomic_add), because a single owned vertex can +// receive contributions from multiple different remote ranks across +// separate messages, and a derived datatype can only place bytes at an +// offset, not reduce concurrently-arriving values. Only the gather (send) +// side is a pure "pick these rows" operation, which is what derived +// datatypes are actually good for here. +// +// Known risk, confirmed relevant on this codebase's field layout +// (LayoutLeft): sending directly from meshField's memory means MPI's +// datatype engine has to gather scattered elements (MPI_Type_vector nested +// in MPI_Type_create_hindexed_block) instead of a few large contiguous +// blocks, which is expensive for many MPI implementations' datatype +// engines to execute well relative to a hand-written parallel GPU pack +// kernel. It also means MPI is handed a pointer straight into meshField's +// own (Kokkos-managed, possibly CudaMallocAsync pool-allocated) memory, +// rather than a dedicated raw-cudaMalloc'd staging buffer - see +// RawCudaMPIBuffer's comment below for why that distinction matters on +// Cray MPICH/GTL. A pre-change, fully protected copy of this file's +// send-side logic is preserved in pmpo_MPMesh_backup.hpp and in git history +// (the "Add CUDA-aware MPI halo exchange" / "Add CUDA-aware communication +// path to MPMesh" commits) for comparison/revert. +// ---------------------------------------------------------------------- + namespace polyMPO{ template @@ -42,14 +90,14 @@ 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 template void communicate_and_take_halo_contributions1( @@ -74,7 +122,7 @@ class MPMesh{ std::vector> recvDataVec; pumipic::RecordTime("SD: Recv Vec Allocation-" + std::to_string(self), timer.seconds()); - + timer.reset(); communicateFields1( @@ -154,7 +202,7 @@ class MPMesh{ } pumipic::RecordTime("SD: Copy CPU-GPU2-" + std::to_string(self), timer.seconds()); - + timer.reset(); if(op == 0){ @@ -198,9 +246,9 @@ class MPMesh{ int mode, std::vector>& recvIDVec, std::vector>& recvDataVec); - - template + + template void communicateFields1( const ViewType& fieldData, const int numEntities, @@ -208,7 +256,7 @@ class MPMesh{ int mode, std::vector>& recvIDVec, std::vector>& recvDataVec){ - + int self, numProcsTot; MPI_Comm comm = p_MPs->getMPIComm(); @@ -394,7 +442,7 @@ class MPMesh{ MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE); } - + MPMesh(Mesh* inMesh, MaterialPoints* inMPs): p_mesh(inMesh), @@ -494,18 +542,20 @@ class MPMesh{ // and gives Cray MPICH/GTL plain CUDA allocations to register/export. #ifdef KOKKOS_ENABLE_CUDA // Plain cudaMalloc'd device buffer, exposed as an unmanaged Kokkos::View - // via .view(). Used only for the 4 buffers below that get handed - // directly to MPI_Isend/Irecv under CUDA-aware MPI. + // via .view(). Used only for the recv-side buffers below (recvIDGPU / + // recvDataGPU) that get handed directly to MPI_Irecv under CUDA-aware + // MPI. The send side no longer needs a plain buffer here - it sends + // directly out of meshField via a derived datatype. // // Why not just a Kokkos::View: when Kokkos is // built with Kokkos_ENABLE_IMPL_CUDA_MALLOC_ASYNC=ON (the default since // Kokkos 4.2), View allocations use cudaMallocAsync/memory pools. // cuIpcGetMemHandle (which Cray MPICH/GTL uses for intra-node GPU-to-GPU // sends) rejects pool allocations with CUDA_ERROR_INVALID_VALUE. Since - // polyMPO isn't allowed to touch the Kokkos build config, these 4 - // buffers bypass Kokkos's allocator entirely via a direct cudaMalloc, - // which cuIpcGetMemHandle always accepts, regardless of how the rest of - // Kokkos (or the rest of the app's Views) is configured. + // polyMPO isn't allowed to touch the Kokkos build config, these buffers + // bypass Kokkos's allocator entirely via a direct cudaMalloc, which + // cuIpcGetMemHandle always accepts, regardless of how the rest of Kokkos + // (or the rest of the app's Views) is configured. template struct RawCudaMPIBuffer{ T* ptr = nullptr; @@ -583,12 +633,6 @@ class MPMesh{ #endif // Cached CUDA-aware MPI communication metadata and batched GPU buffers. - // Every neighbor's data lives in one shared allocation (see - // CudaAwareMPIFieldCache below) and MPI is given "base pointer + byte - // offset" per neighbor rather than a separate allocation per neighbor. - // If you ever need to fall back to one allocation per neighbor (e.g. an - // MPI/GPU stack that mishandles offset device pointers for CUDA IPC), - // restore the per-proc-buffer version from version control. bool cudaAwareMPICacheValid = true; bool cudaAwareMPIDisabled = false; bool cudaAwareMPIEnvChecked = false; @@ -607,17 +651,51 @@ class MPMesh{ int totalSendCount = 0; int totalRecvCount = 0; - // Single batched GPU buffers (one allocation each, instead of one - // Kokkos::View per neighbor proc). Per-proc slices are - // [offset, offset + count) for the ID buffers, and - // [offset * numEntries, (offset + count) * numEntries) for the data - // buffers. MPI is given "buffer base pointer + offset", not a - // separate allocation per proc. - CudaAwareMPIIntBuffer sendEntityGPU; + // Receive side keeps a single batched GPU buffer (one allocation each, + // instead of one Kokkos::View per neighbor proc). Per-proc slices are + // [offset, offset + count) for recvIDGPU, and + // [offset * numEntries, (offset + count) * numEntries) for + // recvDataGPU. Still needed because op==0 unpacking is a scatter-ADD + // (Kokkos::atomic_add, since a single owned vertex can receive + // contributions from several different remote ranks across separate + // messages), which a derived receive-datatype can't do on its own. CudaAwareMPIIntBuffer recvIDGPU; - - CudaAwareMPIDoubleBuffer sendDataGPU; CudaAwareMPIDoubleBuffer recvDataGPU; + + // ---- Send side ---- + // No pack buffer/kernel anymore, unconditionally. rowType describes + // one entity's numEntries doubles exactly as they sit in meshField's + // own memory (built from meshField's actual strides, so it's correct + // whether the view is LayoutRight, LayoutLeft, or padded). + // sendProcTypes[proc] wraps rowType with that proc's list of entity + // byte-displacements via MPI_Type_create_hindexed_block, so MPI_Isend + // can gather directly out of meshField's GPU memory with no separate + // pack kernel/buffer, always. + MPI_Datatype rowType = MPI_DATATYPE_NULL; + std::vector sendProcTypes; // size numProcsTot; MPI_DATATYPE_NULL where unused + + void freeTypes(){ + if(rowType != MPI_DATATYPE_NULL){ + MPI_Type_free(&rowType); + rowType = MPI_DATATYPE_NULL; + } + + for(auto& t : sendProcTypes){ + if(t != MPI_DATATYPE_NULL){ + MPI_Type_free(&t); + t = MPI_DATATYPE_NULL; + } + } + + sendProcTypes.clear(); + } + + // Assumes this cache (and therefore the MPMesh that owns it) is torn + // down before MPI_Finalize - same assumption the rest of this class + // already makes about the MPI resources it holds. + ~CudaAwareMPIFieldCache(){ + freeTypes(); + } }; std::map, CudaAwareMPIFieldCache> cudaAwareMPICaches; @@ -633,18 +711,22 @@ class MPMesh{ return cudaAwareMPIForceCPU; } - // Fully CUDA-aware MPI version, batched buffer variant: - // 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. + // Fully CUDA-aware MPI version, batched buffer + derived-datatype send + // variant: + // + // Field data is sent using GPU pointers gathered directly by a per- + // neighbor MPI derived datatype (no pack kernel/buffer, unconditionally), + // and received into a single batched GPU buffer per cache entry, same + // as before. Receive IDs are cached once from the fixed halo/owner + // mapping and are not sent every call. // - // Every neighbor's send/recv entity-ID list and data live in ONE big - // GPU buffer each (laid out back-to-back in proc order), instead of one - // Kokkos::View allocation per neighbor. Packing/unpacking is a single - // kernel launch over all neighbors' entities at once instead of one - // launch per neighbor, and MPI_Isend/Irecv use "buffer base pointer + - // offset" into that single buffer per proc. This is what actually - // shrinks MPI_Wait time: fewer, larger, more uniform in-flight - // transfers instead of many small independent ones. + // Receive completion is drained with MPI_Waitany instead of + // MPI_Waitall: each neighbor's contribution is unpacked (async kernel + // launch) the moment that neighbor's message arrives, so unpacking + // already-arrived neighbors overlaps with waiting on the remaining + // (slower) ones, instead of the previous "wait for everyone, then + // unpack everyone" ordering. A single Kokkos::fence() at the end + // guarantees all launched unpack kernels have completed. // // Note: an earlier version of this cache used one Kokkos::View // allocation per neighbor specifically to avoid handing MPI a @@ -653,14 +735,16 @@ class MPMesh{ // Kokkos allocating device Views via cudaMallocAsync (invalid for // cuIpcGetMemHandle), not to offset pointers themselves, and is fixed // by building Kokkos with -DKokkos_ENABLE_IMPL_CUDA_MALLOC_ASYNC=OFF. - // If you ever do hit IPC trouble that tracks back to offset pointers - // specifically, the per-proc-buffer version can be restored from - // version control. + // The same offset-pointer reasoning applies to the derived-datatype + // send below (MPI is given meshField's base pointer plus per-entity + // byte displacements); if you hit IPC trouble that tracks back to that, + // the previous pack-buffer send path can be restored from + // pmpo_MPMesh_backup.hpp / version control. // // Important: - // This function caches communication metadata and GPU buffers per - // (mode, numEntries). If the communication pattern changes, clear - // cudaAwareMPICaches before the next call. + // This function caches communication metadata, GPU buffers, and MPI + // derived datatypes 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, @@ -678,7 +762,7 @@ class MPMesh{ MPI_Comm_size(comm, &numProcsTot); const char* diagnosticsEnv = - + std::getenv("POLYMPO_MPI_DIAGNOSTICS"); const bool mpiDiagnostics = @@ -702,7 +786,7 @@ class MPMesh{ #ifdef POLYMPO_VERBOSE_MPI if(self == 0 && !cudaAwareMPILogged){ std::cout - << "[CUDA_AWARE_MPI] Using batched single-buffer GPU-aware MPI path in communicate_and_take_halo_contributions1_improved()" + << "[CUDA_AWARE_MPI] Using batched single-buffer GPU-aware MPI path (derived-datatype send + MPI_Waitany incremental unpack) in communicate_and_take_halo_contributions1_improved()" << "\n"; cudaAwareMPILogged = true; } @@ -724,6 +808,11 @@ class MPMesh{ if(needRebuild){ + // Drop any datatypes committed for a previous topology before + // rebuilding (no-op the first time, when nothing has been built + // yet). + cudaAwareCache.freeTypes(); + cudaAwareCache.cachedNumProcs = numProcsTot; cudaAwareCache.sendCounts.assign(numProcsTot, 0); @@ -758,63 +847,89 @@ class MPMesh{ cudaAwareCache.totalSendCount = totalSend; cudaAwareCache.totalRecvCount = totalRecv; - cudaAwareCache.sendEntityGPU.allocate(totalSend); - cudaAwareCache.sendDataGPU.allocate(totalSend * numEntries); cudaAwareCache.recvIDGPU.allocate(totalRecv); cudaAwareCache.recvDataGPU.allocate(totalRecv * numEntries); - // ---- Build the flattened send-entity list (host, then one deep_copy) ---- + // ---- Build the send-side derived datatypes, unconditionally ---- + // Instead of a flattened host entity list + GPU pack buffer, build + // one MPI_Datatype per neighbor proc that gathers that proc's rows + // directly out of meshField's own memory. 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() == - static_cast(cudaAwareCache.sendCounts[proc])); - } + // rowType: one entity's numEntries doubles, as they actually sit + // in meshField's memory. Query the real strides rather than + // assuming a layout, so this is correct for LayoutRight (entries + // within an entity contiguous - typical host default) and + // LayoutLeft (entries within an entity strided by the entity + // count - typical Kokkos CUDA default) alike. Note: for + // LayoutLeft this produces an MPI_Type_vector describing + // widely-scattered elements, which is the expensive case flagged + // in the file-level comment above. + const size_t strideEntity = meshField.stride_0(); + const size_t strideEntry = meshField.stride_1(); + + if(strideEntry == 1){ + MPI_Type_contiguous(numEntries, MPI_DOUBLE, &cudaAwareCache.rowType); + } + else{ + MPI_Type_vector( + numEntries, + 1, + static_cast(strideEntry), + MPI_DOUBLE, + &cudaAwareCache.rowType); + } + MPI_Type_commit(&cudaAwareCache.rowType); - std::vector cursor(cudaAwareCache.sendOffsets); + std::vector> sendEntityIDsByProc(numProcsTot); + if(mode == 0){ for(int iEnt = 0; iEnt < numHalosTot; iEnt++){ int ownerProc = haloOwnerProcs[iEnt]; if(ownerProc == self) continue; - sendEntityCPU(cursor[ownerProc]) = numOwnersTot + iEnt; - cursor[ownerProc]++; + sendEntityIDsByProc[ownerProc].push_back(numOwnersTot + iEnt); } - + } + else{ for(int proc = 0; proc < numProcsTot; proc++){ if(proc == self) continue; - assert(cursor[proc] == - cudaAwareCache.sendOffsets[proc] + - cudaAwareCache.sendCounts[proc]); + for(auto& ownerID : ownerOwnerLocalIDs[proc]){ + sendEntityIDsByProc[proc].push_back(ownerID); + } } } - else{ - for(int proc = 0; proc < numProcsTot; proc++){ - if(proc == self) continue; - int sendCount = cudaAwareCache.sendCounts[proc]; - if(sendCount <= 0) continue; + cudaAwareCache.sendProcTypes.assign(numProcsTot, MPI_DATATYPE_NULL); - assert(ownerOwnerLocalIDs[proc].size() == - static_cast(sendCount)); + for(int proc = 0; proc < numProcsTot; proc++){ + if(proc == self) continue; - int base = cudaAwareCache.sendOffsets[proc]; + const int sendCount = cudaAwareCache.sendCounts[proc]; + if(sendCount <= 0) continue; - for(int i = 0; i < sendCount; i++){ - sendEntityCPU(base + i) = ownerOwnerLocalIDs[proc][i]; - } + assert(sendEntityIDsByProc[proc].size() == + static_cast(sendCount)); + + std::vector displacements(sendCount); + + for(int i = 0; i < sendCount; i++){ + displacements[i] = + static_cast(sendEntityIDsByProc[proc][i]) * + static_cast(strideEntity) * + static_cast(sizeof(double)); } - } - Kokkos::deep_copy(cudaAwareCache.sendEntityGPU.view(), sendEntityCPU); + MPI_Type_create_hindexed_block( + sendCount, + 1, + displacements.data(), + cudaAwareCache.rowType, + &cudaAwareCache.sendProcTypes[proc]); + + MPI_Type_commit(&cudaAwareCache.sendProcTypes[proc]); + } } // ---- Build the flattened recv-ID list (host, then one deep_copy) ---- @@ -873,38 +988,34 @@ class MPMesh{ timer.reset(); } - // ---- Pack: ONE kernel over all neighbors' send entities at once ---- - if(cudaAwareCache.totalSendCount > 0){ - auto sendEntityGPU = cudaAwareCache.sendEntityGPU.view(); - auto sendDataGPU = cudaAwareCache.sendDataGPU.view(); - - Kokkos::parallel_for( - "pack cached cuda-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); - } - }); - } - + // ---- No explicit pack step, unconditionally ---- + // MPI_Isend below reads meshField directly through the derived + // datatypes built above, so there is no separate pack kernel or + // sendDataGPU buffer to launch/fence on here. We still fence once so + // that any of the caller's kernels which wrote meshField complete + // before MPI starts reading it directly (mirrors the fence the old + // pack step used to provide, just guarding meshField itself now + // instead of a pack buffer). Kokkos::fence(); if(mpiDiagnostics){ pumipic::RecordTime( - label + "_MPI_Diagnostics_Pack_m" + std::to_string(mode) + "_e" + std::to_string(numEntries) + "_rank" + std::to_string(self), timer.seconds()); + label + "_MPI_Diagnostics_Pack_m" + std::to_string(mode) + "_e" + std::to_string(numEntries) + "_rank" + std::to_string(self), 0.0); } timer.reset(); - double postTime = 0.0; - double waitallTime = 0.0; + double postTime = 0.0; + double waitTime = 0.0; + + std::vector recvRequests; + std::vector recvReqProc; + std::vector sendRequests; + + recvRequests.reserve(numProcsTot); + recvReqProc.reserve(numProcsTot); + sendRequests.reserve(numProcsTot); - std::vector requests; - requests.reserve(2 * numProcsTot); int mpiError = MPI_SUCCESS; // Data volume exchanged by this rank in this call (recorded once per @@ -918,12 +1029,16 @@ class MPMesh{ pumipic::RecordTime(label + "_MPI_BytesSent_" + std::to_string(self), bytesSent); pumipic::RecordTime(label + "_MPI_BytesRecv_" + std::to_string(self), bytesRecv); - //Post both the Irecv and the matching Isend for a proc together, in - //the same loop iteration and under their own counts (recvCounts for - //Irecv, sendCounts for Isend). This replaces the previous two-pass - //version (a first loop that posted Irecv only, plus a second loop - //that posted Isend only) which existed only because of a leftover, - //commented-out duplicate of this same block. + // Post both the Irecv and the matching Isend for a proc together, in + // the same loop iteration and under their own counts (recvCounts for + // Irecv, sendCounts for Isend). + // + // The send is posted directly against meshField using that proc's + // derived datatype (no sendDataGPU pointer/offset), unconditionally, + // and recv requests are tracked separately from send requests (with a + // parallel recvReqProc[] telling us which proc each recv request + // belongs to) so the recvs can be drained incrementally with + // MPI_Waitany below instead of all being blocked on together. int numNeighbors = 0; for(int proc = 0; proc < numProcsTot; proc++){ if(proc == self) continue; @@ -946,28 +1061,24 @@ class MPMesh{ comm, &reqData); if(mpiError != MPI_SUCCESS) break; - requests.push_back(reqData); - hasComm = true; + recvRequests.push_back(reqData); + recvReqProc.push_back(proc); + hasComm = true; } if(cudaAwareCache.sendCounts[proc] > 0){ MPI_Request reqData; - double* sendPtr = - cudaAwareCache.sendDataGPU.data() + - static_cast(cudaAwareCache.sendOffsets[proc]) * - numEntries; - mpiError = MPI_Isend( - sendPtr, - cudaAwareCache.sendCounts[proc] * numEntries, - MPI_DOUBLE, + meshField.data(), + 1, + cudaAwareCache.sendProcTypes[proc], proc, 2, comm, &reqData); if(mpiError != MPI_SUCCESS) break; - requests.push_back(reqData); + sendRequests.push_back(reqData); hasComm = true; } if(hasComm) numNeighbors++; @@ -978,31 +1089,84 @@ class MPMesh{ pumipic::RecordTime(label + "_MPI_Post_" + std::to_string(self), postTime); - //Barrier here, not before posting: Isend/Irecv are non-blocking and - //their cost is local (looping + building MPI_Request objects), so a - //barrier before posting would only be measuring how skewed ranks - //were on entry to this function, which the calling function's own - //barriers already capture. Placed here, right before Waitall, it - //makes every rank enter Waitall at the same instant, so the - //waitallTime below reflects real message-arrival/network imbalance - //instead of being contaminated by skew left over from posting. - Kokkos::Timer barrierTimer; - MPI_Barrier(comm); - const double barrierWaitTime = barrierTimer.seconds(); - pumipic::RecordTime(label + "_MPI_BarrierWait_" + std::to_string(self), barrierWaitTime); - timer.reset(); - if(mpiError == MPI_SUCCESS && !requests.empty()){ + // ---- Drain recvs with MPI_Waitany, unpacking each neighbor's + // contribution the moment it lands instead of waiting for every + // neighbor before unpacking any of them. Each unpack kernel launch is + // asynchronous (no fence in the loop), so while neighbor i's data is + // being scattered on the GPU, the CPU is already back inside + // MPI_Waitany waiting on the rest - hiding unpack behind the wait for + // stragglers rather than paying for it serially afterward. ---- + if(mpiError == MPI_SUCCESS && !recvRequests.empty()){ + auto recvIDGPUView = cudaAwareCache.recvIDGPU.view(); + auto recvDataGPUView = cudaAwareCache.recvDataGPU.view(); + + for(size_t reqIdx = 0; reqIdx < recvRequests.size(); reqIdx++){ + int idx = MPI_UNDEFINED; + + mpiError = MPI_Waitany( + static_cast(recvRequests.size()), + recvRequests.data(), + &idx, + MPI_STATUS_IGNORE); + + if(mpiError != MPI_SUCCESS || idx == MPI_UNDEFINED) break; + + const int proc = recvReqProc[idx]; + const int base = cudaAwareCache.recvOffsets[proc]; + const int count = cudaAwareCache.recvCounts[proc]; + + if(op == 0){ + Kokkos::parallel_for( + "halo add cached cuda-aware mpi incremental", + Kokkos::RangePolicy<>(base, base + count), + KOKKOS_LAMBDA(const int i){ + const int vertex = recvIDGPUView(i); + + for(int k = 0; k < numEntries; k++){ +#ifdef POLYMPO_ASSUME_UNIQUE_HALO_CONTRIBS + meshField(vertex, k) += + recvDataGPUView(i * numEntries + k); +#else + Kokkos::atomic_add( + &meshField(vertex, k), + recvDataGPUView(i * numEntries + k)); +#endif + } + }); + } + else{ + Kokkos::parallel_for( + "halo assign cached cuda-aware mpi incremental", + Kokkos::RangePolicy<>(base, base + count), + KOKKOS_LAMBDA(const int i){ + const int vertex = recvIDGPUView(i); + + for(int k = 0; k < numEntries; k++){ + meshField(vertex, k) = + recvDataGPUView(i * numEntries + k); + } + }); + } + // Deliberately no fence here - see comment above the loop. + } + } + + // Sends don't feed any unpack step, but we still need to know they + // completed before treating this call as done, so wait on them once, + // after the recv/unpack loop rather than before it (so posting the + // sends can't itself delay draining the recvs). + if(mpiError == MPI_SUCCESS && !sendRequests.empty()){ mpiError = MPI_Waitall( - static_cast(requests.size()), - requests.data(), + static_cast(sendRequests.size()), + sendRequests.data(), MPI_STATUSES_IGNORE); } - waitallTime = timer.seconds(); + waitTime = timer.seconds(); - pumipic::RecordTime(label + "_MPI_Waitall_" + std::to_string(self), waitallTime); + pumipic::RecordTime(label + "_MPI_Waitall_" + std::to_string(self), waitTime); if(mpiError != MPI_SUCCESS){ cudaAwareMPIDisabled = true; @@ -1013,7 +1177,7 @@ class MPMesh{ << std::endl; } - if(requests.empty()){ + if(recvRequests.empty() && sendRequests.empty()){ if(self == 0){ std::cout << "[CUDA_AWARE_MPI] Falling back to CPU-staged communication." @@ -1042,52 +1206,21 @@ class MPMesh{ timer.reset(); - // ---- Unpack: ONE kernel over all neighbors' recv entities at once ---- - 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 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); -#else - Kokkos::atomic_add( - &meshField(vertex, k), - recvDataGPU(i * numEntries + k)); -#endif - } - }); - } - else{ - Kokkos::parallel_for( - "halo assign cached cuda-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); - } - }); - } - } - + // Final fence: guarantees every unpack kernel launched inside the + // MPI_Waitany loop above has actually finished before meshField is + // used downstream. When neighbor arrivals are staggered, most of that + // unpack work already finished while we were still waiting on + // stragglers, so in that case this fence's cost is close to just the + // last-arriving neighbor's unpack kernel rather than the sum of all + // of them. If arrivals are tightly bunched instead, expect this to + // look a lot like the old post-Waitall unpack cost. Kokkos::fence(); if(mpiDiagnostics){ pumipic::RecordTime(label + "_MPI_Diagnostics_Contribution_m" + std::to_string(mode) + "_e" + std::to_string(numEntries) + "_rank" + std::to_string(self), timer.seconds()); } - + } From cb5fceae651805ccb568b6d3c71017824a3fa9c4 Mon Sep 17 00:00:00 2001 From: Shahrear Jahan Santho Date: Wed, 12 Aug 2026 21:21:06 -0700 Subject: [PATCH 11/11] Cuda Aware Halo Exchange - Diagnostic labels with Barrier-off --- src/pmpo_MPMesh.cpp | 8 +- src/pmpo_MPMesh.hpp | 495 ++++++++++----------------- src/pmpo_MPMesh_assembly.hpp | 72 ++-- src/pmpo_c.cpp | 16 +- src/pmpo_c.h | 4 +- src/pmpo_fortran.f90 | 6 +- test/testFortran.f90 | 4 +- test/testFortranMPAdvection.f90 | 6 +- test/testFortranMPReconstruction.f90 | 4 +- 9 files changed, 243 insertions(+), 372 deletions(-) diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index 6eb635c3..10aa52aa 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -484,7 +484,13 @@ void MPMesh::startCommunication(){ int self, numProcsTot; MPI_Comm comm = p_MPs->getMPIComm(); MPI_Comm_rank(comm, &self); - MPI_Comm_size(comm, &numProcsTot); + MPI_Comm_size(comm, &numProcsTot); + + std::cout << "[RankSummary] Rank=" << self + << " Vertices(total)=" << p_mesh->getNumVertices() + << " Vertices(owned)=" << p_mesh->getNumVerticesOwned() + << " Elements=" << p_mesh->getNumElements() + << std::endl; //The routine should work for elements too, although currently the communication //is done for vertices. For elements, the follwoing three variables should correspond diff --git a/src/pmpo_MPMesh.hpp b/src/pmpo_MPMesh.hpp index 11849ece..3ddb1032 100644 --- a/src/pmpo_MPMesh.hpp +++ b/src/pmpo_MPMesh.hpp @@ -13,54 +13,6 @@ #include #endif -// ---------------------------------------------------------------------- -// Halo-exchange optimization: pure derived-datatype send (pack path fully -// removed) + MPI_Waitany incremental unpack. -// -// communicate_and_take_halo_contributions1_improved() (CUDA_AWARE_MPI path -// only - the CPU-staged fallback below is untouched): -// -// 1. Eliminate the pack step / sendDataGPU buffer, unconditionally. Each -// neighbor proc gets its own committed MPI derived datatype -// (MPI_Type_create_hindexed_block over a "rowType" built from -// meshField's own strides) that gathers that proc's rows directly out -// of meshField's GPU memory. MPI_Isend reads straight from the field - -// no pack kernel, no extra buffer, no fallback to a pack path. -// -// 2. Hide the unpack step behind MPI_Waitall latency instead of paying -// for it serially afterward. Recv requests are drained with -// MPI_Waitany instead of MPI_Waitall, and each neighbor's contribution -// is unpacked (async kernel launch, no intermediate fence) the instant -// that neighbor's message lands, so unpacking early arrivals overlaps -// with waiting on stragglers. A single Kokkos::fence() at the end -// guarantees every launched unpack kernel has completed before -// meshField is used downstream. -// -// The receive side still uses a batched contiguous GPU buffer + a real -// unpack kernel (not a derived receive-datatype): for op==0 the unpack is a -// scatter-ADD (Kokkos::atomic_add), because a single owned vertex can -// receive contributions from multiple different remote ranks across -// separate messages, and a derived datatype can only place bytes at an -// offset, not reduce concurrently-arriving values. Only the gather (send) -// side is a pure "pick these rows" operation, which is what derived -// datatypes are actually good for here. -// -// Known risk, confirmed relevant on this codebase's field layout -// (LayoutLeft): sending directly from meshField's memory means MPI's -// datatype engine has to gather scattered elements (MPI_Type_vector nested -// in MPI_Type_create_hindexed_block) instead of a few large contiguous -// blocks, which is expensive for many MPI implementations' datatype -// engines to execute well relative to a hand-written parallel GPU pack -// kernel. It also means MPI is handed a pointer straight into meshField's -// own (Kokkos-managed, possibly CudaMallocAsync pool-allocated) memory, -// rather than a dedicated raw-cudaMalloc'd staging buffer - see -// RawCudaMPIBuffer's comment below for why that distinction matters on -// Cray MPICH/GTL. A pre-change, fully protected copy of this file's -// send-side logic is preserved in pmpo_MPMesh_backup.hpp and in git history -// (the "Add CUDA-aware MPI halo exchange" / "Add CUDA-aware communication -// path to MPMesh" commits) for comparison/revert. -// ---------------------------------------------------------------------- - namespace polyMPO{ template @@ -90,14 +42,14 @@ 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 template void communicate_and_take_halo_contributions1( @@ -122,7 +74,7 @@ class MPMesh{ std::vector> recvDataVec; pumipic::RecordTime("SD: Recv Vec Allocation-" + std::to_string(self), timer.seconds()); - + timer.reset(); communicateFields1( @@ -202,7 +154,7 @@ class MPMesh{ } pumipic::RecordTime("SD: Copy CPU-GPU2-" + std::to_string(self), timer.seconds()); - + timer.reset(); if(op == 0){ @@ -246,9 +198,9 @@ class MPMesh{ int mode, std::vector>& recvIDVec, std::vector>& recvDataVec); + - - template + template void communicateFields1( const ViewType& fieldData, const int numEntities, @@ -256,7 +208,7 @@ class MPMesh{ int mode, std::vector>& recvIDVec, std::vector>& recvDataVec){ - + int self, numProcsTot; MPI_Comm comm = p_MPs->getMPIComm(); @@ -442,7 +394,7 @@ class MPMesh{ MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE); } - + MPMesh(Mesh* inMesh, MaterialPoints* inMPs): p_mesh(inMesh), @@ -542,20 +494,18 @@ class MPMesh{ // and gives Cray MPICH/GTL plain CUDA allocations to register/export. #ifdef KOKKOS_ENABLE_CUDA // Plain cudaMalloc'd device buffer, exposed as an unmanaged Kokkos::View - // via .view(). Used only for the recv-side buffers below (recvIDGPU / - // recvDataGPU) that get handed directly to MPI_Irecv under CUDA-aware - // MPI. The send side no longer needs a plain buffer here - it sends - // directly out of meshField via a derived datatype. + // via .view(). Used only for the 4 buffers below that get handed + // directly to MPI_Isend/Irecv under CUDA-aware MPI. // // Why not just a Kokkos::View: when Kokkos is // built with Kokkos_ENABLE_IMPL_CUDA_MALLOC_ASYNC=ON (the default since // Kokkos 4.2), View allocations use cudaMallocAsync/memory pools. // cuIpcGetMemHandle (which Cray MPICH/GTL uses for intra-node GPU-to-GPU // sends) rejects pool allocations with CUDA_ERROR_INVALID_VALUE. Since - // polyMPO isn't allowed to touch the Kokkos build config, these buffers - // bypass Kokkos's allocator entirely via a direct cudaMalloc, which - // cuIpcGetMemHandle always accepts, regardless of how the rest of Kokkos - // (or the rest of the app's Views) is configured. + // polyMPO isn't allowed to touch the Kokkos build config, these 4 + // buffers bypass Kokkos's allocator entirely via a direct cudaMalloc, + // which cuIpcGetMemHandle always accepts, regardless of how the rest of + // Kokkos (or the rest of the app's Views) is configured. template struct RawCudaMPIBuffer{ T* ptr = nullptr; @@ -633,6 +583,12 @@ class MPMesh{ #endif // Cached CUDA-aware MPI communication metadata and batched GPU buffers. + // Every neighbor's data lives in one shared allocation (see + // CudaAwareMPIFieldCache below) and MPI is given "base pointer + byte + // offset" per neighbor rather than a separate allocation per neighbor. + // If you ever need to fall back to one allocation per neighbor (e.g. an + // MPI/GPU stack that mishandles offset device pointers for CUDA IPC), + // restore the per-proc-buffer version from version control. bool cudaAwareMPICacheValid = true; bool cudaAwareMPIDisabled = false; bool cudaAwareMPIEnvChecked = false; @@ -651,51 +607,17 @@ class MPMesh{ int totalSendCount = 0; int totalRecvCount = 0; - // Receive side keeps a single batched GPU buffer (one allocation each, - // instead of one Kokkos::View per neighbor proc). Per-proc slices are - // [offset, offset + count) for recvIDGPU, and - // [offset * numEntries, (offset + count) * numEntries) for - // recvDataGPU. Still needed because op==0 unpacking is a scatter-ADD - // (Kokkos::atomic_add, since a single owned vertex can receive - // contributions from several different remote ranks across separate - // messages), which a derived receive-datatype can't do on its own. + // Single batched GPU buffers (one allocation each, instead of one + // Kokkos::View per neighbor proc). Per-proc slices are + // [offset, offset + count) for the ID buffers, and + // [offset * numEntries, (offset + count) * numEntries) for the data + // buffers. MPI is given "buffer base pointer + offset", not a + // separate allocation per proc. + CudaAwareMPIIntBuffer sendEntityGPU; CudaAwareMPIIntBuffer recvIDGPU; - CudaAwareMPIDoubleBuffer recvDataGPU; - - // ---- Send side ---- - // No pack buffer/kernel anymore, unconditionally. rowType describes - // one entity's numEntries doubles exactly as they sit in meshField's - // own memory (built from meshField's actual strides, so it's correct - // whether the view is LayoutRight, LayoutLeft, or padded). - // sendProcTypes[proc] wraps rowType with that proc's list of entity - // byte-displacements via MPI_Type_create_hindexed_block, so MPI_Isend - // can gather directly out of meshField's GPU memory with no separate - // pack kernel/buffer, always. - MPI_Datatype rowType = MPI_DATATYPE_NULL; - std::vector sendProcTypes; // size numProcsTot; MPI_DATATYPE_NULL where unused - - void freeTypes(){ - if(rowType != MPI_DATATYPE_NULL){ - MPI_Type_free(&rowType); - rowType = MPI_DATATYPE_NULL; - } - for(auto& t : sendProcTypes){ - if(t != MPI_DATATYPE_NULL){ - MPI_Type_free(&t); - t = MPI_DATATYPE_NULL; - } - } - - sendProcTypes.clear(); - } - - // Assumes this cache (and therefore the MPMesh that owns it) is torn - // down before MPI_Finalize - same assumption the rest of this class - // already makes about the MPI resources it holds. - ~CudaAwareMPIFieldCache(){ - freeTypes(); - } + CudaAwareMPIDoubleBuffer sendDataGPU; + CudaAwareMPIDoubleBuffer recvDataGPU; }; std::map, CudaAwareMPIFieldCache> cudaAwareMPICaches; @@ -711,22 +633,18 @@ class MPMesh{ return cudaAwareMPIForceCPU; } - // Fully CUDA-aware MPI version, batched buffer + derived-datatype send - // variant: - // - // Field data is sent using GPU pointers gathered directly by a per- - // neighbor MPI derived datatype (no pack kernel/buffer, unconditionally), - // and received into a single batched GPU buffer per cache entry, same - // as before. Receive IDs are cached once from the fixed halo/owner - // mapping and are not sent every call. + // Fully CUDA-aware MPI version, batched buffer variant: + // 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. // - // Receive completion is drained with MPI_Waitany instead of - // MPI_Waitall: each neighbor's contribution is unpacked (async kernel - // launch) the moment that neighbor's message arrives, so unpacking - // already-arrived neighbors overlaps with waiting on the remaining - // (slower) ones, instead of the previous "wait for everyone, then - // unpack everyone" ordering. A single Kokkos::fence() at the end - // guarantees all launched unpack kernels have completed. + // Every neighbor's send/recv entity-ID list and data live in ONE big + // GPU buffer each (laid out back-to-back in proc order), instead of one + // Kokkos::View allocation per neighbor. Packing/unpacking is a single + // kernel launch over all neighbors' entities at once instead of one + // launch per neighbor, and MPI_Isend/Irecv use "buffer base pointer + + // offset" into that single buffer per proc. This is what actually + // shrinks MPI_Wait time: fewer, larger, more uniform in-flight + // transfers instead of many small independent ones. // // Note: an earlier version of this cache used one Kokkos::View // allocation per neighbor specifically to avoid handing MPI a @@ -735,16 +653,14 @@ class MPMesh{ // Kokkos allocating device Views via cudaMallocAsync (invalid for // cuIpcGetMemHandle), not to offset pointers themselves, and is fixed // by building Kokkos with -DKokkos_ENABLE_IMPL_CUDA_MALLOC_ASYNC=OFF. - // The same offset-pointer reasoning applies to the derived-datatype - // send below (MPI is given meshField's base pointer plus per-entity - // byte displacements); if you hit IPC trouble that tracks back to that, - // the previous pack-buffer send path can be restored from - // pmpo_MPMesh_backup.hpp / version control. + // If you ever do hit IPC trouble that tracks back to offset pointers + // specifically, the per-proc-buffer version can be restored from + // version control. // // Important: - // This function caches communication metadata, GPU buffers, and MPI - // derived datatypes per (mode, numEntries). If the communication - // pattern changes, clear cudaAwareMPICaches before the next call. + // 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, @@ -762,7 +678,7 @@ class MPMesh{ MPI_Comm_size(comm, &numProcsTot); const char* diagnosticsEnv = - + std::getenv("POLYMPO_MPI_DIAGNOSTICS"); const bool mpiDiagnostics = @@ -786,7 +702,7 @@ class MPMesh{ #ifdef POLYMPO_VERBOSE_MPI if(self == 0 && !cudaAwareMPILogged){ std::cout - << "[CUDA_AWARE_MPI] Using batched single-buffer GPU-aware MPI path (derived-datatype send + MPI_Waitany incremental unpack) in communicate_and_take_halo_contributions1_improved()" + << "[CUDA_AWARE_MPI] Using batched single-buffer GPU-aware MPI path in communicate_and_take_halo_contributions1_improved()" << "\n"; cudaAwareMPILogged = true; } @@ -808,11 +724,6 @@ class MPMesh{ if(needRebuild){ - // Drop any datatypes committed for a previous topology before - // rebuilding (no-op the first time, when nothing has been built - // yet). - cudaAwareCache.freeTypes(); - cudaAwareCache.cachedNumProcs = numProcsTot; cudaAwareCache.sendCounts.assign(numProcsTot, 0); @@ -847,89 +758,63 @@ class MPMesh{ cudaAwareCache.totalSendCount = totalSend; cudaAwareCache.totalRecvCount = totalRecv; + cudaAwareCache.sendEntityGPU.allocate(totalSend); + cudaAwareCache.sendDataGPU.allocate(totalSend * numEntries); cudaAwareCache.recvIDGPU.allocate(totalRecv); cudaAwareCache.recvDataGPU.allocate(totalRecv * numEntries); - // ---- Build the send-side derived datatypes, unconditionally ---- - // Instead of a flattened host entity list + GPU pack buffer, build - // one MPI_Datatype per neighbor proc that gathers that proc's rows - // directly out of meshField's own memory. + // ---- Build the flattened send-entity list (host, then one deep_copy) ---- if(totalSend > 0){ + auto sendEntityCPU = + Kokkos::View( + "sendEntityCPU_batched", totalSend); - // rowType: one entity's numEntries doubles, as they actually sit - // in meshField's memory. Query the real strides rather than - // assuming a layout, so this is correct for LayoutRight (entries - // within an entity contiguous - typical host default) and - // LayoutLeft (entries within an entity strided by the entity - // count - typical Kokkos CUDA default) alike. Note: for - // LayoutLeft this produces an MPI_Type_vector describing - // widely-scattered elements, which is the expensive case flagged - // in the file-level comment above. - const size_t strideEntity = meshField.stride_0(); - const size_t strideEntry = meshField.stride_1(); - - if(strideEntry == 1){ - MPI_Type_contiguous(numEntries, MPI_DOUBLE, &cudaAwareCache.rowType); - } - else{ - MPI_Type_vector( - numEntries, - 1, - static_cast(strideEntry), - MPI_DOUBLE, - &cudaAwareCache.rowType); - } - MPI_Type_commit(&cudaAwareCache.rowType); + if(mode == 0){ + for(int proc = 0; proc < numProcsTot; proc++){ + if(proc == self) continue; + if(cudaAwareCache.sendCounts[proc] <= 0) continue; - std::vector> sendEntityIDsByProc(numProcsTot); + assert(haloOwnerLocalIDs[proc].size() == + static_cast(cudaAwareCache.sendCounts[proc])); + } + + std::vector cursor(cudaAwareCache.sendOffsets); - if(mode == 0){ for(int iEnt = 0; iEnt < numHalosTot; iEnt++){ int ownerProc = haloOwnerProcs[iEnt]; if(ownerProc == self) continue; - sendEntityIDsByProc[ownerProc].push_back(numOwnersTot + iEnt); + sendEntityCPU(cursor[ownerProc]) = numOwnersTot + iEnt; + cursor[ownerProc]++; } - } - else{ + for(int proc = 0; proc < numProcsTot; proc++){ if(proc == self) continue; - for(auto& ownerID : ownerOwnerLocalIDs[proc]){ - sendEntityIDsByProc[proc].push_back(ownerID); - } + assert(cursor[proc] == + cudaAwareCache.sendOffsets[proc] + + cudaAwareCache.sendCounts[proc]); } } + else{ + for(int proc = 0; proc < numProcsTot; proc++){ + if(proc == self) continue; - cudaAwareCache.sendProcTypes.assign(numProcsTot, MPI_DATATYPE_NULL); - - for(int proc = 0; proc < numProcsTot; proc++){ - if(proc == self) continue; - - const int sendCount = cudaAwareCache.sendCounts[proc]; - if(sendCount <= 0) continue; + int sendCount = cudaAwareCache.sendCounts[proc]; + if(sendCount <= 0) continue; - assert(sendEntityIDsByProc[proc].size() == - static_cast(sendCount)); + assert(ownerOwnerLocalIDs[proc].size() == + static_cast(sendCount)); - std::vector displacements(sendCount); + int base = cudaAwareCache.sendOffsets[proc]; - for(int i = 0; i < sendCount; i++){ - displacements[i] = - static_cast(sendEntityIDsByProc[proc][i]) * - static_cast(strideEntity) * - static_cast(sizeof(double)); + for(int i = 0; i < sendCount; i++){ + sendEntityCPU(base + i) = ownerOwnerLocalIDs[proc][i]; + } } - - MPI_Type_create_hindexed_block( - sendCount, - 1, - displacements.data(), - cudaAwareCache.rowType, - &cudaAwareCache.sendProcTypes[proc]); - - MPI_Type_commit(&cudaAwareCache.sendProcTypes[proc]); } + + Kokkos::deep_copy(cudaAwareCache.sendEntityGPU.view(), sendEntityCPU); } // ---- Build the flattened recv-ID list (host, then one deep_copy) ---- @@ -988,34 +873,38 @@ class MPMesh{ timer.reset(); } - // ---- No explicit pack step, unconditionally ---- - // MPI_Isend below reads meshField directly through the derived - // datatypes built above, so there is no separate pack kernel or - // sendDataGPU buffer to launch/fence on here. We still fence once so - // that any of the caller's kernels which wrote meshField complete - // before MPI starts reading it directly (mirrors the fence the old - // pack step used to provide, just guarding meshField itself now - // instead of a pack buffer). + // ---- Pack: ONE kernel over all neighbors' send entities at once ---- + if(cudaAwareCache.totalSendCount > 0){ + auto sendEntityGPU = cudaAwareCache.sendEntityGPU.view(); + auto sendDataGPU = cudaAwareCache.sendDataGPU.view(); + + Kokkos::parallel_for( + "pack cached cuda-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); + } + }); + } + Kokkos::fence(); if(mpiDiagnostics){ pumipic::RecordTime( - label + "_MPI_Diagnostics_Pack_m" + std::to_string(mode) + "_e" + std::to_string(numEntries) + "_rank" + std::to_string(self), 0.0); + label + "_MPI_Diagnostics_Pack_m" + std::to_string(mode) + "_e" + std::to_string(numEntries) + "_rank" + std::to_string(self), timer.seconds()); } timer.reset(); - double postTime = 0.0; - double waitTime = 0.0; - - std::vector recvRequests; - std::vector recvReqProc; - std::vector sendRequests; - - recvRequests.reserve(numProcsTot); - recvReqProc.reserve(numProcsTot); - sendRequests.reserve(numProcsTot); + double postTime = 0.0; + double waitallTime = 0.0; + std::vector requests; + requests.reserve(2 * numProcsTot); int mpiError = MPI_SUCCESS; // Data volume exchanged by this rank in this call (recorded once per @@ -1029,16 +918,12 @@ class MPMesh{ pumipic::RecordTime(label + "_MPI_BytesSent_" + std::to_string(self), bytesSent); pumipic::RecordTime(label + "_MPI_BytesRecv_" + std::to_string(self), bytesRecv); - // Post both the Irecv and the matching Isend for a proc together, in - // the same loop iteration and under their own counts (recvCounts for - // Irecv, sendCounts for Isend). - // - // The send is posted directly against meshField using that proc's - // derived datatype (no sendDataGPU pointer/offset), unconditionally, - // and recv requests are tracked separately from send requests (with a - // parallel recvReqProc[] telling us which proc each recv request - // belongs to) so the recvs can be drained incrementally with - // MPI_Waitany below instead of all being blocked on together. + //Post both the Irecv and the matching Isend for a proc together, in + //the same loop iteration and under their own counts (recvCounts for + //Irecv, sendCounts for Isend). This replaces the previous two-pass + //version (a first loop that posted Irecv only, plus a second loop + //that posted Isend only) which existed only because of a leftover, + //commented-out duplicate of this same block. int numNeighbors = 0; for(int proc = 0; proc < numProcsTot; proc++){ if(proc == self) continue; @@ -1061,24 +946,28 @@ class MPMesh{ comm, &reqData); if(mpiError != MPI_SUCCESS) break; - recvRequests.push_back(reqData); - recvReqProc.push_back(proc); - hasComm = true; + requests.push_back(reqData); + hasComm = true; } if(cudaAwareCache.sendCounts[proc] > 0){ MPI_Request reqData; + double* sendPtr = + cudaAwareCache.sendDataGPU.data() + + static_cast(cudaAwareCache.sendOffsets[proc]) * + numEntries; + mpiError = MPI_Isend( - meshField.data(), - 1, - cudaAwareCache.sendProcTypes[proc], + sendPtr, + cudaAwareCache.sendCounts[proc] * numEntries, + MPI_DOUBLE, proc, 2, comm, &reqData); if(mpiError != MPI_SUCCESS) break; - sendRequests.push_back(reqData); + requests.push_back(reqData); hasComm = true; } if(hasComm) numNeighbors++; @@ -1089,84 +978,31 @@ class MPMesh{ pumipic::RecordTime(label + "_MPI_Post_" + std::to_string(self), postTime); - timer.reset(); + //Barrier here, not before posting: Isend/Irecv are non-blocking and + //their cost is local (looping + building MPI_Request objects), so a + //barrier before posting would only be measuring how skewed ranks + //were on entry to this function, which the calling function's own + //barriers already capture. Placed here, right before Waitall, it + //makes every rank enter Waitall at the same instant, so the + //waitallTime below reflects real message-arrival/network imbalance + //instead of being contaminated by skew left over from posting. + //Kokkos::Timer barrierTimer; + //MPI_Barrier(comm); + //const double InternalBarrierWait = barrierTimer.seconds(); + //pumipic::RecordTime(label + "_MPI_BarrierWait_" + std::to_string(self), InternalBarrierWait); - // ---- Drain recvs with MPI_Waitany, unpacking each neighbor's - // contribution the moment it lands instead of waiting for every - // neighbor before unpacking any of them. Each unpack kernel launch is - // asynchronous (no fence in the loop), so while neighbor i's data is - // being scattered on the GPU, the CPU is already back inside - // MPI_Waitany waiting on the rest - hiding unpack behind the wait for - // stragglers rather than paying for it serially afterward. ---- - if(mpiError == MPI_SUCCESS && !recvRequests.empty()){ - auto recvIDGPUView = cudaAwareCache.recvIDGPU.view(); - auto recvDataGPUView = cudaAwareCache.recvDataGPU.view(); - - for(size_t reqIdx = 0; reqIdx < recvRequests.size(); reqIdx++){ - int idx = MPI_UNDEFINED; - - mpiError = MPI_Waitany( - static_cast(recvRequests.size()), - recvRequests.data(), - &idx, - MPI_STATUS_IGNORE); - - if(mpiError != MPI_SUCCESS || idx == MPI_UNDEFINED) break; - - const int proc = recvReqProc[idx]; - const int base = cudaAwareCache.recvOffsets[proc]; - const int count = cudaAwareCache.recvCounts[proc]; - - if(op == 0){ - Kokkos::parallel_for( - "halo add cached cuda-aware mpi incremental", - Kokkos::RangePolicy<>(base, base + count), - KOKKOS_LAMBDA(const int i){ - const int vertex = recvIDGPUView(i); - - for(int k = 0; k < numEntries; k++){ -#ifdef POLYMPO_ASSUME_UNIQUE_HALO_CONTRIBS - meshField(vertex, k) += - recvDataGPUView(i * numEntries + k); -#else - Kokkos::atomic_add( - &meshField(vertex, k), - recvDataGPUView(i * numEntries + k)); -#endif - } - }); - } - else{ - Kokkos::parallel_for( - "halo assign cached cuda-aware mpi incremental", - Kokkos::RangePolicy<>(base, base + count), - KOKKOS_LAMBDA(const int i){ - const int vertex = recvIDGPUView(i); - - for(int k = 0; k < numEntries; k++){ - meshField(vertex, k) = - recvDataGPUView(i * numEntries + k); - } - }); - } - // Deliberately no fence here - see comment above the loop. - } - } + timer.reset(); - // Sends don't feed any unpack step, but we still need to know they - // completed before treating this call as done, so wait on them once, - // after the recv/unpack loop rather than before it (so posting the - // sends can't itself delay draining the recvs). - if(mpiError == MPI_SUCCESS && !sendRequests.empty()){ + if(mpiError == MPI_SUCCESS && !requests.empty()){ mpiError = MPI_Waitall( - static_cast(sendRequests.size()), - sendRequests.data(), + static_cast(requests.size()), + requests.data(), MPI_STATUSES_IGNORE); } - waitTime = timer.seconds(); + waitallTime = timer.seconds(); - pumipic::RecordTime(label + "_MPI_Waitall_" + std::to_string(self), waitTime); + pumipic::RecordTime(label + "_MPI_Waitall_" + std::to_string(self), waitallTime); if(mpiError != MPI_SUCCESS){ cudaAwareMPIDisabled = true; @@ -1177,7 +1013,7 @@ class MPMesh{ << std::endl; } - if(recvRequests.empty() && sendRequests.empty()){ + if(requests.empty()){ if(self == 0){ std::cout << "[CUDA_AWARE_MPI] Falling back to CPU-staged communication." @@ -1206,21 +1042,52 @@ class MPMesh{ timer.reset(); - // Final fence: guarantees every unpack kernel launched inside the - // MPI_Waitany loop above has actually finished before meshField is - // used downstream. When neighbor arrivals are staggered, most of that - // unpack work already finished while we were still waiting on - // stragglers, so in that case this fence's cost is close to just the - // last-arriving neighbor's unpack kernel rather than the sum of all - // of them. If arrivals are tightly bunched instead, expect this to - // look a lot like the old post-Waitall unpack cost. + // ---- Unpack: ONE kernel over all neighbors' recv entities at once ---- + 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 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); +#else + Kokkos::atomic_add( + &meshField(vertex, k), + recvDataGPU(i * numEntries + k)); +#endif + } + }); + } + else{ + Kokkos::parallel_for( + "halo assign cached cuda-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(); if(mpiDiagnostics){ pumipic::RecordTime(label + "_MPI_Diagnostics_Contribution_m" + std::to_string(mode) + "_e" + std::to_string(numEntries) + "_rank" + std::to_string(self), timer.seconds()); } - + } diff --git a/src/pmpo_MPMesh_assembly.hpp b/src/pmpo_MPMesh_assembly.hpp index 8e6189a6..d57bf4c7 100644 --- a/src/pmpo_MPMesh_assembly.hpp +++ b/src/pmpo_MPMesh_assembly.hpp @@ -103,19 +103,19 @@ void MPMesh::reconstruct_coeff_full(){ MPI_Comm_size(comm, &numProcsTot); Kokkos::fence(); //B0: drain any device work left from the previous phase - MPI_Barrier(comm); //B0: align all ranks so this timed region starts at the same instant + //MPI_Barrier(comm); //B0: align all ranks so this timed region starts at the same instant Kokkos::Timer totalTimer; Kokkos::Timer phaseTimer; double preCommunicationComputeTime = 0.0; - double preCommunicationComputeImbalance = 0.0; + //double preCommunicationComputeImbalance = 0.0; double gatherCommunicationTime = 0.0; - double gatherCommunicationImbalance = 0.0; + //double gatherCommunicationImbalance = 0.0; double scatterCommunicationTime = 0.0; - double scatterCommunicationImbalance = 0.0; + //double scatterCommunicationImbalance = 0.0; double postCommunicationComputeTime = 0.0; - double postCommunicationComputeImbalance = 0.0; + //double postCommunicationComputeImbalance = 0.0; static int coeff_count=0; if(!self) std::cout<<"===="<<__FUNCTION__<<" "<p_MPs; - p_MPs->setElmIDoffset(offset); + p_MPs->setElmIDoffset(offset); + + int self; + MPI_Comm_rank(p_MPs->getMPIComm(), &self); + std::cout << "[RankSummary] Rank=" << self + << " MaterialPoints=" << p_MPs->getCount() + << std::endl; } void polympo_startRebuildMPs_f(MPMesh_ptr p_mpmesh, @@ -567,7 +573,7 @@ void polympo_getMPMass_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs pumipic::RecordTime("PolyMPO_getMPMass", timer.seconds()); } -void polympo_setMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpVelIn) { +void polympo_setMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpVelIn, const int callSiteId) { Kokkos::Timer timer; checkMPMeshValid(p_mpmesh); auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; @@ -590,7 +596,7 @@ void polympo_setMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, } }; p_MPs->parallel_for(setMPVel, "setMPVel"); - pumipic::RecordTime("PolyMPO_setMPVel" + std::to_string(self),timer.seconds()); + pumipic::RecordTime("PolyMPO_setMPVel_site" + std::to_string(callSiteId) + "_" + std::to_string(self), timer.seconds()); } void polympo_getMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpVelHost) { @@ -1247,7 +1253,7 @@ void polympo_setMeshVtxVel_f(MPMesh_ptr p_mpmesh, const int nVertices, const dou Kokkos::deep_copy(coordsArray, h_coordsArray); } -void polympo_getMeshVtxVel_f(MPMesh_ptr p_mpmesh, const int nVertices, double* uVelOut, double* vVelOut){ +void polympo_getMeshVtxVel_f(MPMesh_ptr p_mpmesh, const int nVertices, double* uVelOut, double* vVelOut, const int callSiteId){ Kokkos::Timer timer; //check mpMesh is valid checkMPMeshValid(p_mpmesh); @@ -1267,7 +1273,7 @@ void polympo_getMeshVtxVel_f(MPMesh_ptr p_mpmesh, const int nVertices, double* u uVelOut[i] = h_coordsArray(i,0); vVelOut[i] = h_coordsArray(i,1); } - pumipic::RecordTime("PolyMPO_getMeshVtxVel" + std::to_string(self), timer.seconds()); + pumipic::RecordTime("PolyMPO_getMeshVtxVel_site" + std::to_string(callSiteId) + "_" + std::to_string(self), timer.seconds()); } void polympo_setMeshVtxMass_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* vtxMass){ diff --git a/src/pmpo_c.h b/src/pmpo_c.h index d19b62b5..b938570f 100644 --- a/src/pmpo_c.h +++ b/src/pmpo_c.h @@ -44,7 +44,7 @@ void polympo_getMPTgtRotLatLon_f(MPMesh_ptr p_mpmesh, const int nComps, const in //MP fields void polympo_setMPMass_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpMassIn); void polympo_getMPMass_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpMassHost); -void polympo_setMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpVelIn); +void polympo_setMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpVelIn, const int callSiteId); void polympo_getMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpVelHost); void polympo_calculateMPStrainRate_f(MPMesh_ptr p_mpmesh); void polympo_setMPStrainRate_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpStrainRateIn); @@ -92,7 +92,7 @@ void polympo_setMeshVtxRotLat_f(MPMesh_ptr p_mpmesh, const int nVertices, const void polympo_setMeshVtxRotLon_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* longitude); void polympo_getMeshVtxRotLat_f(MPMesh_ptr p_mpmesh, const int nVertices, double* latitude); void polympo_setMeshVtxVel_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* uVelocity, const double* vVelocity); -void polympo_getMeshVtxVel_f(MPMesh_ptr p_mpmesh, const int nVertices, double* uVelocity, double* vVelocity); +void polympo_getMeshVtxVel_f(MPMesh_ptr p_mpmesh, const int nVertices, double* uVelocity, double* vVelocity, const int callSiteId); void polympo_setMeshVtxMass_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* vtxMass); void polympo_getMeshVtxMass_f(MPMesh_ptr p_mpmesh, const int nVertices, double* vtxMass); void polympo_setMeshElmMass_f(MPMesh_ptr p_mpmesh, const int nCells, const double* elmMass); diff --git a/src/pmpo_fortran.f90 b/src/pmpo_fortran.f90 index 98b5f646..e07db516 100644 --- a/src/pmpo_fortran.f90 +++ b/src/pmpo_fortran.f90 @@ -340,12 +340,13 @@ subroutine polympo_getMPMass(mpMesh, nComps, numMPs, array) & !> @param numMPs(in) number of the MPs !> @param array(in) input MP velocity 1D array (numMPs*2) !--------------------------------------------------------------------------- - subroutine polympo_setMPVel(mpMesh, nComps, numMPs, array) & + subroutine polympo_setMPVel(mpMesh, nComps, numMPs, array, callSiteId) & bind(C, NAME='polympo_setMPVel_f') use :: iso_c_binding type(c_ptr), value :: mpMesh integer(c_int), value :: nComps, numMPs type(c_ptr), intent(in), value :: array + integer(c_int), value :: callSiteId end subroutine !--------------------------------------------------------------------------- @@ -818,12 +819,13 @@ subroutine polympo_setMeshVtxVel(mpMesh, nVertices, uVel, vVel) & !> @param vVel(in/out) output vertices v-component velocity !> 1D array (numVtx), allocated by user !--------------------------------------------------------------------------- - subroutine polympo_getMeshVtxVel(mpMesh, nVertices, uVel, vVel) & + subroutine polympo_getMeshVtxVel(mpMesh, nVertices, uVel, vVel, callSiteId) & bind(C, NAME='polympo_getMeshVtxVel_f') use :: iso_c_binding type(c_ptr), value :: mpMesh integer(c_int), value :: nVertices type(c_ptr), value :: uVel, vVel + integer(c_int), value :: callSiteId end subroutine !--------------------------------------------------------------------------- diff --git a/test/testFortran.f90 b/test/testFortran.f90 index 2eff4d9f..5d8538dd 100644 --- a/test/testFortran.f90 +++ b/test/testFortran.f90 @@ -76,7 +76,7 @@ program main MParray(i,j) = (i-1)*numMPs + j end do end do - call polympo_setMPVel(mpMesh, numCompsVel, numMPs, c_loc(MParray)) + call polympo_setMPVel(mpMesh, numCompsVel, numMPs, c_loc(MParray), 5) ! check MP Fields MParray = -1 @@ -134,7 +134,7 @@ program main call polympo_setMeshVtxVel(mpMesh, nverts, c_loc(xArray),c_loc(yArray)) xArray = -1 yArray = -1 - call polympo_getMeshVtxVel(mpMesh, nverts, c_loc(xArray),c_loc(yArray)) + call polympo_getMeshVtxVel(mpMesh, nverts, c_loc(xArray),c_loc(yArray), 4) do i = 1, nverts call assert((xArray(i) .eq. i+value1), "Assert MeshVel u-component Velocity Fail") call assert((yArray(i) .eq. value2-i), "Assert MeshVel v-component Velocity Fail") diff --git a/test/testFortranMPAdvection.f90 b/test/testFortranMPAdvection.f90 index 4ac1c6c9..1f48729c 100644 --- a/test/testFortranMPAdvection.f90 +++ b/test/testFortranMPAdvection.f90 @@ -105,7 +105,7 @@ subroutine runReconstructionTest(mpMesh, numMPs, numPush, nCells, nVertices, mp2 mpVel = TEST_VAL call polympo_setMPMass(mpMesh,1,numMPs,c_loc(mpMass)) - call polympo_setMPVel(mpMesh,2,numMPs,c_loc(mpVel)) + call polympo_setMPVel(mpMesh,2,numMPs,c_loc(mpVel), 7) ! Although this test just does 0th order reconstruction testing, and just needs the BasisSlice, ! calculating the coefficeints too as that will involve calculating the Basis Slice @@ -169,7 +169,7 @@ subroutine runApiTest(mpMesh, numMPs, nVertices, nCells, numPush, mpLatLon, mpPo call polympo_setMPPositions(mpMesh,3,numMPs,c_loc(mpPosition)) call polympo_setMeshVtxOnSurfDispIncr(mpMesh,nCompsDisp,nVertices,c_loc(dispIncr)) call polympo_setMPMass(mpMesh,1,numMPs,c_loc(mpMass)) - call polympo_setMPVel(mpMesh,2,numMPs,c_loc(mpVel)) + call polympo_setMPVel(mpMesh,2,numMPs,c_loc(mpVel), 8) call polympo_setMeshVtxCoords(mpMesh, nVertices, c_loc(xVertex), c_loc(yVertex), c_loc(zVertex)) call polympo_setMeshVtxRotLat(mpMesh,nVertices,c_loc(latVertex)) @@ -177,7 +177,7 @@ subroutine runApiTest(mpMesh, numMPs, nVertices, nCells, numPush, mpLatLon, mpPo call polympo_getMPVel(mpMesh, 2, numMPs, c_loc(mpVel)) call polympo_getMeshElmMass(mpMesh,nCells,c_loc(meshElmMass)) call polympo_getMeshVtxMass(mpMesh,nVertices,c_loc(meshVtxMass)) - call polympo_getMeshVtxVel(mpMesh,nVertices, c_loc(xVertex),c_loc(yVertex)) + call polympo_getMeshVtxVel(mpMesh,nVertices, c_loc(xVertex),c_loc(yVertex), 6) end do deallocate(dispIncr) diff --git a/test/testFortranMPReconstruction.f90 b/test/testFortranMPReconstruction.f90 index 76b4b80f..86f88238 100644 --- a/test/testFortranMPReconstruction.f90 +++ b/test/testFortranMPReconstruction.f90 @@ -120,7 +120,7 @@ program main call polympo_setMPPositions(mpMesh,3,numMPs,c_loc(mpPosition)) call polympo_setMPMass(mpMesh,1,numMPs,c_loc(mpMass)) - call polympo_setMPVel(mpMesh,2,numMPs,c_loc(mpVel)) + call polympo_setMPVel(mpMesh,2,numMPs,c_loc(mpVel), 6) !First Order reconstruction done before 0th order reconstruction as calculation of coefficients will !fill the MpBasis slice @@ -130,7 +130,7 @@ program main !call polympo_setReconstructionOfVel(mpMesh, 1, polympo_getMeshFVtxType()) call polympo_applyReconstruction(mpMesh) call polympo_getMeshVtxMass(mpMesh,nVertices,c_loc(meshVtxMass1)) - call polympo_getMeshVtxVel(mpMesh, nVertices, c_loc(meshVtxVelu), c_loc(meshVtxVelv)) + call polympo_getMeshVtxVel(mpMesh, nVertices, c_loc(meshVtxVelu), c_loc(meshVtxVelv), 5) do i = 1, nVertices call assert(meshVtxMass1(i) < TEST_VAL+TOLERANCE1 .and. meshVtxMass1(i) > TEST_VAL-TOLERANCE1, "Error: wrong vtx mass order 1") !call assert(meshVtxVelu(i) < TEST_VAL+TOLERANCE1 .and. meshVtxVelu(i) > TEST_VAL-TOLERANCE1, "Error: wrong vtx velU order 1")