Merge pull request #278 from LLNL/bugfix/central_mpi_calls

Cleanup MPI calls

RELEASE_NOTES.md

@@ -1,3 +1,22 @@
+Version vYYYY.MM.p -- Release date YYYY-MM-DD
+==============================================
+  * Important Notes:
+
+Notable changes include:
+
+  * New features/ API changes:
+    * MPI variables are now wrapped as
+      ```
+      SPHERAL_OP_SUM, SPHERAL_OP_MAX, SPHERAL_OP_MIN
+      ```
+
+  * Build changes / improvements:
+    * Distributed source directory must always be built now
+
+  * Bug Fixes / improvements:
+    * Wrappers for MPI calls are simplified and improved
+
+
 Version v2024.06.1 -- Release date 2024-07-09
 ==============================================
 

src/Boundary/InflowOutflowBoundary.cc

@@ -11,7 +11,7 @@
 #include "Field/FieldBase.hh"
 #include "Hydro/HydroFieldNames.hh"
 #include "Geometry/GeometryRegistrar.hh"
-#include "Utilities/allReduce.hh"
+#include "Distributed/allReduce.hh"
 #include "Utilities/planarReflectingOperator.hh"
 #include "Utilities/DBC.hh"
 
@@ -133,7 +133,7 @@ updateGhostNodes(NodeList<Dimension>& nodeList) {
       const auto xd = mPlane.signedDistance(pos[i]);
       xmin = std::min(xmin, xd);
     }
-    xmin = allReduce(xmin, MPI_MIN, Communicator::communicator());
+    xmin = allReduce(xmin, SPHERAL_OP_MIN);
     // CHECK(xmin >= 0.0);
 
     // Offset the current ghost points appropriately.
@@ -309,8 +309,8 @@ InflowOutflowBoundary<Dimension>::initializeProblemStartup(const bool /*final*/)
     for (const auto i: nodeIDs) {
       vinflow += vel[i].dot(nhat);
     }
-    vinflow = (allReduce(vinflow, MPI_SUM, Communicator::communicator())/
-               std::max(1.0e-30, allReduce(double(nodeIDs.size()), MPI_SUM, Communicator::communicator())));  // Negative implies outflow
+    vinflow = (allReduce(vinflow, SPHERAL_OP_SUM)/
+               std::max(1.0e-30, allReduce(double(nodeIDs.size()), SPHERAL_OP_SUM)));  // Negative implies outflow
 
     // Figure out a timestep limit such that we don't move more than the ghost
     // node thickness.
@@ -320,8 +320,8 @@ InflowOutflowBoundary<Dimension>::initializeProblemStartup(const bool /*final*/)
       xmin = std::min(xmin, xd);
       xmax = std::max(xmax, xd);
     }
-    xmin = allReduce(xmin, MPI_MIN, Communicator::communicator());
-    xmax = allReduce(xmax, MPI_MAX, Communicator::communicator());
+    xmin = allReduce(xmin, SPHERAL_OP_MIN);
+    xmax = allReduce(xmax, SPHERAL_OP_MAX);
     mXmin[nodeList.name()] = xmin;
     mDT = std::min(mDT, std::abs(xmax - xmin)/std::max(1e-30, std::abs(vinflow)));   // Protect from negative outflow velocity
     // cerr << "Timestep constraint: " << mDT << endl;
@@ -458,7 +458,7 @@ InflowOutflowBoundary<Dimension>::finalize(const Scalar /*time*/,
       nodeList.neighbor().updateNodes();
     }
   }
-  altered = (allReduce((altered ? 1 : 0), MPI_MAX, Communicator::communicator()) == 1);
+  altered = (allReduce((altered ? 1 : 0), SPHERAL_OP_MAX) == 1);
 
   // If any NodeLists were altered, recompute the boundary conditions.
   if (altered) {

src/Boundary/PlanarBoundary.cc

@@ -12,7 +12,6 @@
 #include "NodeList/FluidNodeList.hh"
 #include "Mesh/Mesh.hh"
 #include "Utilities/DBC.hh"
-#include "Utilities/allReduce.hh"
 
 #include "PlanarBoundary.hh"
 

src/Boundary/findNodesTouchingThroughPlanes.cc

@@ -2,7 +2,7 @@
 // Find the set of nodes that see through a pair of planes.
 //------------------------------------------------------------------------------
 #include "findNodesTouchingThroughPlanes.hh"
-#include "Utilities/allReduce.hh"
+#include "Distributed/allReduce.hh"
 
 namespace Spheral {
 
@@ -50,7 +50,7 @@ findNodesTouchingThroughPlanes(const NodeList<Dimension>& nodeList,
       const auto  hmaxi = 1.0/Hi.eigenValues().minElement();
       if (hmaxi > hmax and std::min(exitPlane.minimumDistance(ri), enterPlane.minimumDistance(ri)) < kernelExtent*hmaxi) hmax = hmaxi;
     }
-    hmax = allReduce(hmax, MPI_MAX, Communicator::communicator());
+    hmax = allReduce(hmax, SPHERAL_OP_MAX);
 
     // Now find all points within this range of the exit plane.
     if (hmax > 0.0) {

src/Damage/DamageModel.cc

@@ -196,7 +196,7 @@ finalize(const Scalar /*time*/,
         nD += nD_thread;
       }
     }
-    nD = allReduce(nD, MPI_SUM, Communicator::communicator());
+    nD = allReduce(nD, SPHERAL_OP_SUM);
     const auto ntot = std::max(1, dataBase.globalNumInternalNodes());
     const auto dfrac = double(nD)/double(ntot);
     mComputeIntersectConnectivity = (dfrac > 0.2);  // Should tune this number...

src/Damage/IvanoviSALEDamageModel.cc

@@ -38,7 +38,6 @@
 #include "Boundary/Boundary.hh"
 #include "Neighbor/Neighbor.hh"
 #include "Utilities/mortonOrderIndices.hh"
-#include "Utilities/allReduce.hh"
 #include "Utilities/uniform_random.hh"
 
 #include <boost/functional/hash.hpp>  // hash_combine

src/Damage/ProbabilisticDamageModel.cc

@@ -29,7 +29,7 @@
 #include "Boundary/Boundary.hh"
 #include "Neighbor/Neighbor.hh"
 #include "Utilities/mortonOrderIndices.hh"
-#include "Utilities/allReduce.hh"
+#include "Distributed/allReduce.hh"
 #include "Utilities/uniform_random.hh"
 
 #include <boost/functional/hash.hpp>  // hash_combine
@@ -121,7 +121,7 @@ initializeProblemStartupDependencies(DataBase<Dimension>& dataBase,
   for (auto i = 0u; i < mMask.numInternalElements(); ++i) {
     if (mMask[i] == 1) ++nused_local;
   }
-  const size_t nused_global = allReduce(nused_local, MPI_SUM, Communicator::communicator());
+  const size_t nused_global = allReduce(nused_local, SPHERAL_OP_SUM);
 
   // Compute the Morton-ordering for hashing with the global seed to seed each
   // point-wise random number generator.
@@ -158,8 +158,8 @@ initializeProblemStartupDependencies(DataBase<Dimension>& dataBase,
       randomGenerators[i]();                // Recommended to discard first value in sequence
     }
   }
-  mVmin = allReduce(mVmin, MPI_MIN, Communicator::communicator());
-  mVmax = allReduce(mVmax, MPI_MAX, Communicator::communicator());
+  mVmin = allReduce(mVmin, SPHERAL_OP_MIN);
+  mVmax = allReduce(mVmax, SPHERAL_OP_MAX);
 
   // Generate min/max ranges of flaws for each point.
   const auto mInv = 1.0/mmWeibull;
@@ -200,12 +200,12 @@ initializeProblemStartupDependencies(DataBase<Dimension>& dataBase,
 
   // Some diagnostic output.
   if (nused_global > 0) {
-    minNumFlaws = allReduce(minNumFlaws, MPI_MIN, Communicator::communicator());
-    maxNumFlaws = allReduce(maxNumFlaws, MPI_MAX, Communicator::communicator());
-    totalNumFlaws = allReduce(totalNumFlaws, MPI_SUM, Communicator::communicator());
-    epsMin = allReduce(epsMin, MPI_MIN, Communicator::communicator());
-    epsMax = allReduce(epsMax, MPI_MAX, Communicator::communicator());
-    numFlawsRatio = allReduce(numFlawsRatio, MPI_SUM, Communicator::communicator())/nused_global;
+    minNumFlaws = allReduce(minNumFlaws, SPHERAL_OP_MIN);
+    maxNumFlaws = allReduce(maxNumFlaws, SPHERAL_OP_MAX);
+    totalNumFlaws = allReduce(totalNumFlaws, SPHERAL_OP_SUM);
+    epsMin = allReduce(epsMin, SPHERAL_OP_MIN);
+    epsMax = allReduce(epsMax, SPHERAL_OP_MAX);
+    numFlawsRatio = allReduce(numFlawsRatio, SPHERAL_OP_SUM)/nused_global;
     if (Process::getRank() == 0) {
       cerr << "ProbabilisticDamageModel for " << nodes.name() << ":" << endl
            << " Min, max, max/min volumes: " << mVmin << " " << mVmax << " " << mVmax*safeInv(mVmin) << endl

src/Damage/computeFragmentField.cc

@@ -68,14 +68,7 @@ globalReduceToUniqueElements(vector<int>& x) {
     copy(otherX.begin(), otherX.end(), back_inserter(x));
   }
   reduceToUniqueElements(x);
-  BEGIN_CONTRACT_SCOPE
-  {
-    int tmp = x.size();
-    int sum;
-    MPI_Allreduce(&tmp, &sum, 1, MPI_INT, MPI_SUM, Communicator::communicator());
-    ENSURE(sum == (int)x.size()*numProcs);
-  }
-  END_CONTRACT_SCOPE
+  ENSURE(allReduce((int)x.size(), SPHERAL_OP_SUM) == (int)x.size()*numProcs);
 #endif
 }
 
@@ -132,12 +125,7 @@ computeFragmentField(const NodeList<Dimension>& nodes,
                                                          gIDs.end());
   int maxGlobalID = 0;
   if (maxGlobalItr != gIDs.end()) maxGlobalID = *maxGlobalItr;
-#ifdef USE_MPI
-  {
-    int tmp = maxGlobalID;
-    MPI_Allreduce(&tmp, &maxGlobalID, 1, MPI_INT, MPI_MAX, Communicator::communicator());
-  }
-#endif
+  maxGlobalID = allReduce(maxGlobalID, SPHERAL_OP_MAX);
   maxGlobalID += 1;
   CHECK(maxGlobalID >= numGlobalNodesRemaining);
 
@@ -161,12 +149,7 @@ computeFragmentField(const NodeList<Dimension>& nodes,
   }
 
   // Reduce the count of remaining nodes by the number of dust nodes.
-#ifdef USE_MPI
-  {
-    int tmp = numDustNodes;
-    MPI_Allreduce(&tmp, &numDustNodes, 1, MPI_INT, MPI_SUM, Communicator::communicator());
-  }
-#endif
+  numDustNodes = allReduce(numDustNodes, SPHERAL_OP_SUM);
   CHECK(numDustNodes >= 0 && numDustNodes <= numGlobalNodesRemaining);
   numGlobalNodesRemaining -= numDustNodes;
   CHECK(numGlobalNodesRemaining >= 0);
@@ -184,12 +167,7 @@ computeFragmentField(const NodeList<Dimension>& nodes,
                                                            globalNodesRemaining.end());
     int globalMinID = maxGlobalID;
     if (globalMinItr != globalNodesRemaining.end()) globalMinID = *globalMinItr;
-#ifdef USE_MPI
-    {
-      int tmp = globalMinID;
-      MPI_Allreduce(&tmp, &globalMinID, 1, MPI_INT, MPI_MIN, Communicator::communicator());
-    }
-#endif
+    globalMinID = allReduce(globalMinID, SPHERAL_OP_MIN);
     CHECK(globalMinID < maxGlobalID);
 
     // Is this node on this domain?
@@ -201,30 +179,17 @@ computeFragmentField(const NodeList<Dimension>& nodes,
                                                  gIDs.end(),
                                                  globalMinID);
     localNode = (ilocalItr != gIDs.end());
-    BEGIN_CONTRACT_SCOPE
-    {
-      int tmp = localNode ? 1 : 0;
-      int sum;
-      MPI_Allreduce(&tmp, &sum, 1, MPI_INT, MPI_SUM, Communicator::communicator());
-      CHECK(sum == 1);
-    }
-    END_CONTRACT_SCOPE
+    CHECK(allReduce(localNode ? 1 : 0, SPHERAL_OP_SUM) == 1);
     int tmp = numProcs;
     if (localNode) {
       CHECK(ilocalItr != gIDs.end());
       ilocal = distance(gIDs.begin(), ilocalItr);
       tmp = procID;
       CHECK(result(ilocal) == maxGlobalID);
     }
-    MPI_Allreduce(&tmp, &nodeDomain, 1, MPI_INT, MPI_MIN, Communicator::communicator());
+    nodeDomain = allReduce(tmp, SPHERAL_OP_MIN);
     CHECK(nodeDomain >= 0 && nodeDomain < numProcs);
-    BEGIN_CONTRACT_SCOPE
-    {
-      int tmp;
-      MPI_Allreduce(&nodeDomain, &tmp, 1, MPI_INT, MPI_SUM, Communicator::communicator());
-      CHECK(tmp == numProcs*nodeDomain);
-    }
-    END_CONTRACT_SCOPE
+    CHECK(allReduce(nodeDomain, SPHERAL_OP_SUM) == numProcs*nodeDomain);
 #endif
 
     // Get the position and H for this node.
@@ -273,15 +238,8 @@ computeFragmentField(const NodeList<Dimension>& nodes,
     }
     CHECK(fragID >= 0 && fragID < numFragments);
 #ifdef USE_MPI
-    BEGIN_CONTRACT_SCOPE
-    {
-      int tmp;
-      MPI_Allreduce(&fragID, &tmp, 1, MPI_INT, MPI_SUM, Communicator::communicator());
-      CHECK(tmp == numProcs*fragID);
-      MPI_Allreduce(&numFragments, &tmp, 1, MPI_INT, MPI_SUM, Communicator::communicator());
-      CHECK(tmp == numProcs*numFragments);
-    }
-    END_CONTRACT_SCOPE
+    CHECK(allReduce(fragID, SPHERAL_OP_SUM) == numProcs*fragID);
+    CHECK(allReduce(numFragments, SPHERAL_OP_SUM) == numProcs*numFragments);
 #endif
 
     // Remove the known maxGlobalID from the stack of fragment IDs.
@@ -317,13 +275,7 @@ computeFragmentField(const NodeList<Dimension>& nodes,
       if (removeItr != globalNodesRemaining.end())
         globalNodesRemaining.erase(removeItr);
     }
-    numGlobalNodesRemaining = globalNodesRemaining.size();
-#ifdef USE_MPI
-    {
-      int tmp = numGlobalNodesRemaining;
-      MPI_Allreduce(&tmp, &numGlobalNodesRemaining, 1, MPI_INT, MPI_SUM, Communicator::communicator());
-    }
-#endif
+    numGlobalNodesRemaining = allReduce(globalNodesRemaining.size(), SPHERAL_OP_SUM);
 
     BEGIN_CONTRACT_SCOPE
     {

src/Damage/weibullFlawDistributionBenzAsphaug.cc

@@ -19,7 +19,7 @@
 #include "Strength/SolidFieldNames.hh"
 #include "DataBase/State.hh"
 #include "Distributed/Communicator.hh"
-#include "Utilities/allReduce.hh"
+#include "Distributed/allReduce.hh"
 
 using std::unordered_map;
 using std::vector;
@@ -99,7 +99,7 @@ weibullFlawDistributionBenzAsphaug(double volume,
         CHECK(rho(i) > 0.0);
         volume += mass(i)/rho(i);
       }
-      volume = allReduce(volume, MPI_SUM, Communicator::communicator());
+      volume = allReduce(volume, SPHERAL_OP_SUM);
     }
     volume = std::max(volume, 1e-100);
     CHECK(volume > 0.0);
@@ -168,11 +168,11 @@ weibullFlawDistributionBenzAsphaug(double volume,
 
     // Prepare some diagnostic output.
     const auto nused = std::max(1, mask.sumElements());
-    minNumFlaws = allReduce(minNumFlaws, MPI_MIN, Communicator::communicator());
-    maxNumFlaws = allReduce(maxNumFlaws, MPI_MAX, Communicator::communicator());
-    totalNumFlaws = allReduce(totalNumFlaws, MPI_SUM, Communicator::communicator());
-    epsMax = allReduce(epsMax, MPI_MAX, Communicator::communicator());
-    sumFlaws = allReduce(sumFlaws, MPI_SUM, Communicator::communicator());
+    minNumFlaws = allReduce(minNumFlaws, SPHERAL_OP_MIN);
+    maxNumFlaws = allReduce(maxNumFlaws, SPHERAL_OP_MAX);
+    totalNumFlaws = allReduce(totalNumFlaws, SPHERAL_OP_SUM);
+    epsMax = allReduce(epsMax, SPHERAL_OP_MAX);
+    sumFlaws = allReduce(sumFlaws, SPHERAL_OP_SUM);
     if (procID == 0) {
       cerr << "weibullFlawDistributionBenzAsphaug: Min num flaws per node: " << minNumFlaws << endl
            << "                                    Max num flaws per node: " << maxNumFlaws << endl

src/Damage/weibullFlawDistributionOwen.cc

@@ -12,7 +12,7 @@
 #include "Strength/SolidFieldNames.hh"
 #include "DataBase/State.hh"
 #include "Distributed/Communicator.hh"
-#include "Utilities/allReduce.hh"
+#include "Distributed/allReduce.hh"
 
 #include <boost/functional/hash.hpp>  // hash_combine
 
@@ -108,8 +108,8 @@ weibullFlawDistributionOwen(const unsigned seed,
         Vmax = max(Vmax, Vi);
       }
     }
-    Vmin = allReduce(Vmin*volumeMultiplier, MPI_MIN, Communicator::communicator());
-    Vmax = allReduce(Vmax*volumeMultiplier, MPI_MAX, Communicator::communicator());
+    Vmin = allReduce(Vmin*volumeMultiplier, SPHERAL_OP_MIN);
+    Vmax = allReduce(Vmax*volumeMultiplier, SPHERAL_OP_MAX);
     CHECK(Vmin > 0.0);
     CHECK(Vmax >= Vmin);
 
@@ -157,12 +157,12 @@ weibullFlawDistributionOwen(const unsigned seed,
     // Some diagnostic output.
     const auto nused = std::max(1, mask.sumElements());
     if (nglobal > 0) {
-      minNumFlaws = allReduce(minNumFlaws, MPI_MIN, Communicator::communicator());
-      maxNumFlaws = allReduce(maxNumFlaws, MPI_MAX, Communicator::communicator());
-      totalNumFlaws = allReduce(totalNumFlaws, MPI_SUM, Communicator::communicator());
-      epsMin = allReduce(epsMin, MPI_MIN, Communicator::communicator());
-      epsMax = allReduce(epsMax, MPI_MAX, Communicator::communicator());
-      sumFlaws = allReduce(sumFlaws, MPI_SUM, Communicator::communicator());
+      minNumFlaws = allReduce(minNumFlaws, SPHERAL_OP_MIN);
+      maxNumFlaws = allReduce(maxNumFlaws, SPHERAL_OP_MAX);
+      totalNumFlaws = allReduce(totalNumFlaws, SPHERAL_OP_SUM);
+      epsMin = allReduce(epsMin, SPHERAL_OP_MIN);
+      epsMax = allReduce(epsMax, SPHERAL_OP_MAX);
+      sumFlaws = allReduce(sumFlaws, SPHERAL_OP_SUM);
     }
     if (procID == 0) {
       cerr << "weibullFlawDistributionOwen: Min num flaws per node: " << minNumFlaws << endl