Improvements to integrator output

src/Distributed/allReduce.hh

@@ -27,22 +27,45 @@ namespace Spheral {
 #define SPHERAL_OP_PROD MPI_PROD
 #define SPHERAL_OP_LAND MPI_LAND
 #define SPHERAL_OP_LOR MPI_LOR
+#define SPHERAL_OP_MINLOC MPI_MINLOC
+#define SPHERAL_OP_MAXLOC MPI_MAXLOC
 
 template<typename Value>
 constexpr Value
 allReduce(const Value& value, const MPI_Op op,
           const MPI_Comm comm = Communicator::communicator()) {
+  CHECK(!(op == SPHERAL_OP_MINLOC || op == SPHERAL_OP_MAXLOC));
   Value tmp = value;
   Value result;
   MPI_Allreduce(&tmp, &result, 1,
                 DataTypeTraits<Value>::MpiDataType(), op, comm);
   return result;
 }
 
+template<typename Value>
+constexpr std::pair<Value, int>
+allReduceLoc(const Value value, const MPI_Op op,
+             const MPI_Comm comm = Communicator::communicator()) {
+  CHECK(op == SPHERAL_OP_MINLOC || op == SPHERAL_OP_MAXLOC);
+  struct {
+    Value val;
+    int rank;
+  } in, out;
+
+  MPI_Comm_rank(comm, &in.rank);
+  in.val = value;
+
+  MPI_Allreduce(&in, &out, 1, DataTypeTraits<Value>::MpiLocDataType(), op, comm);
+
+  return {out.val, out.rank};
+}
+
+
 template<typename Value>
 constexpr Value
 distScan(const Value& value, const MPI_Op op,
      const MPI_Comm comm = Communicator::communicator()) {
+  CHECK(!(op == SPHERAL_OP_MINLOC || op == SPHERAL_OP_MAXLOC));
   Value tmp = value;
   Value result;
   MPI_Scan(&tmp, &result, 1, DataTypeTraits<Value>::MpiDataType(), op, comm);
@@ -65,13 +88,22 @@ Barrier(const MPI_Comm comm = Communicator::communicator()) {
 #define SPHERAL_OP_PROD 4
 #define SPHERAL_OP_LAND 5
 #define SPHERAL_OP_LOR 6
+#define SPHERAL_OP_MINLOC 7
+#define SPHERAL_OP_MAXLOC 8
 
 template<typename Value>
 constexpr Value
 allReduce(const Value& value, const int /*op*/, const int = 0) {
   return value;
 }
 
+template<typename Value>
+inline std::pair<Value, int>
+allReduceLoc(const Value value, const int /*op*/,
+             const int = 0) {
+  return {value, 0};
+}
+
 template<typename Value>
 constexpr Value
 distScan(const Value& value, const int /*op*/, const int = 0) {

src/Integrator/Integrator.cc

@@ -48,6 +48,7 @@ Integrator(DataBase<Dimension>& dataBase,
   mCurrentTime(0.0),
   mCurrentCycle(0),
   mUpdateBoundaryFrequency(1),
+  mVerboseStep(1),
   mVerbose(false),
   mAllowDtCheck(false),
   mRequireConnectivity(true),
@@ -150,16 +151,16 @@ selectDt(const typename Dimension::Scalar dtMin,
 
   // In the parallel case we need to find the minimum timestep across all processors.
 #ifdef SPHERAL_ENABLE_GLOBALDT_REDUCTION
-  const auto globalDt = allReduce(dt.first, SPHERAL_OP_MIN);
+  const auto [globalDt, rank] = allReduceLoc(dt.first, SPHERAL_OP_MINLOC);
 #else
   const auto globalDt = dt.first;
+  const auto rank = Process::getRank();
 #endif
 
   // Are we verbose?
-  if (dt.first == globalDt and 
-      (verbose() or globalDt < mDtMin)) {
+  if (rank == Process::getRank() and (verbose() and currentCycle() % verboseStep() == 0)) {
     std::cout << "Selected timestep of "
-              << dt.first << std::endl
+              << dt.first << " on rank " << rank << std::endl
               << dt.second << std::endl;
   }
   std::cout.flush();

src/Integrator/Integrator.hh

@@ -192,6 +192,10 @@ public:
   bool verbose() const                                                              { return mVerbose; }
   void verbose(const bool x)                                                        { mVerbose = x; }
 
+  // Select whether the integrator is verbose or not during a cycle.
+  int verboseStep() const                                                              { return mVerboseStep; }
+  void verboseStep(const int x)                                                        { mVerboseStep = x; }
+
   // Should the integrator check interim timestep votes and abort steps?
   bool allowDtCheck() const                                                         { return mAllowDtCheck; }
   void allowDtCheck(const bool x)                                                   { mAllowDtCheck = x; }
@@ -235,7 +239,7 @@ protected:
 private:
   //--------------------------- Private Interface ---------------------------//
   Scalar mDtMin, mDtMax, mDtGrowth, mLastDt, mDtCheckFrac, mCurrentTime;
-  int mCurrentCycle, mUpdateBoundaryFrequency;
+  int mCurrentCycle, mUpdateBoundaryFrequency, mVerboseStep;
   bool mVerbose, mAllowDtCheck, mRequireConnectivity, mRequireGhostConnectivity, mRequireOverlapConnectivity, mRequireIntersectionConnectivity;
   std::reference_wrapper<DataBase<Dimension>> mDataBase;
   std::vector<Physics<Dimension>*> mPhysicsPackages;

src/PYB11/Integrator/Integrator.py

@@ -181,6 +181,7 @@ def copyGhostState(self,
     rigorousBoundaries = PYB11property("bool", "rigorousBoundaries", "rigorousBoundaries", doc="Toggle if ghost nodes should be recomputed every derivative estimate")
     updateBoundaryFrequency = PYB11property("int", "updateBoundaryFrequency", "updateBoundaryFrequency", doc="Optionally update the boundary ghost nodes only on this frequency of cycles")
     verbose = PYB11property("bool", "verbose", "verbose", doc="Verbose time step information every step")
+    verboseStep = PYB11property("int", "verboseStep", "verboseStep", doc="Select frequency of verbose time step output")
     allowDtCheck = PYB11property("bool", "allowDtCheck", "allowDtCheck", doc="Should the integrator check interim timestep votes and abort steps?")
     domainDecompositionIndependent = PYB11property("bool", "domainDecompositionIndependent", "domainDecompositionIndependent", doc="Order operations to be bit perfect reproducible regardless of domain decomposition")
     cullGhostNodes = PYB11property("bool", "cullGhostNodes", "cullGhostNodes", doc="Cull ghost nodes to just active set")

src/SimulationControl/SpheralController.py

@@ -438,12 +438,18 @@ def doPeriodicWork(self, force=False):
 
         # Do any periodic cycle work, as determined by the number of steps.
         for method, frequency in self._periodicWork:
-            if frequency is not None and (force or self.totalSteps % frequency == 0):
+            if frequency is None:
+                continue
+            f0 = frequency(t0) if callable(frequency) else frequency
+            if force or self.totalSteps % f0 == 0:
                 method(self.totalSteps, t1, dt)
 
         # Do any periodic time work, as determined by the current time.
         for method, frequency in self._periodicTimeWork:
-            if frequency is not None and (force or ((t0 // frequency) != (t1 // frequency))):
+            if frequency is None:
+                continue
+            f0 = frequency(t0) if callable(frequency) else frequency
+            if force or ((t0 // f0) != (t1 // f0)):
                 method(self.totalSteps, t1, dt)
 
         return

src/Utilities/DataTypeTraits.hh

@@ -110,6 +110,7 @@ struct DataTypeTraits<int> {
   SPHERAL_HOST_DEVICE static int zero() { return 0; }
 #ifdef SPHERAL_ENABLE_MPI
   static MPI_Datatype MpiDataType() { return MPI_INT; }
+  static MPI_Datatype MpiLocDataType() { return MPI_2INT; }
 #endif
   static axom::sidre::DataTypeId axomTypeID() { return axom::sidre::INT_ID; }
   using AxomType = int;
@@ -174,6 +175,7 @@ struct DataTypeTraits<float> {
   SPHERAL_HOST_DEVICE static float zero() { return 0.0; }
 #ifdef SPHERAL_ENABLE_MPI
   static MPI_Datatype MpiDataType() { return MPI_FLOAT; }
+  static MPI_Datatype MpiLocDataType() { return MPI_FLOAT_INT; }
 #endif
   static axom::sidre::DataTypeId axomTypeID() { return axom::sidre::FLOAT_ID; }
   using AxomType = float;
@@ -188,6 +190,7 @@ struct DataTypeTraits<double> {
   SPHERAL_HOST_DEVICE static double zero() { return 0.0; }
 #ifdef SPHERAL_ENABLE_MPI
   static MPI_Datatype MpiDataType() { return MPI_DOUBLE; }
+  static MPI_Datatype MpiLocDataType() { return MPI_DOUBLE_INT; }
 #endif
   static axom::sidre::DataTypeId axomTypeID() { return axom::sidre::DOUBLE_ID; }
   using AxomType = double;