Rewrite integrators script interface as Python code

src/core/EspressoSystemStandAlone.cpp

@@ -62,5 +62,5 @@ void EspressoSystemStandAlone::set_time_step(double time_step) const {
 }
 
 void EspressoSystemStandAlone::set_skin(double new_skin) const {
-  mpi_set_skin_local(new_skin);
+  ::set_skin(new_skin);
 }

src/core/integrate.cpp

@@ -93,14 +93,8 @@ static double verlet_reuse = 0.0;
 
 static int fluid_step = 0;
 
-bool set_py_interrupt = false;
 namespace {
 volatile std::sig_atomic_t ctrl_C = 0;
-
-void notify_sig_int() {
-  ctrl_C = 0;              // reset
-  set_py_interrupt = true; // global to notify Python
-}
 } // namespace
 
 namespace LeesEdwards {
@@ -259,10 +253,11 @@ int integrate(int n_steps, int reuse_forces) {
 
   // If any method vetoes (e.g. P3M not initialized), immediately bail out
   if (check_runtime_errors(comm_cart))
-    return 0;
+    return INTEG_ERROR_RUNTIME;
 
   // Additional preparations for the first integration step
-  if (reuse_forces == -1 || (recalc_forces && reuse_forces != 1)) {
+  if (reuse_forces == INTEG_REUSE_FORCES_NEVER or
+      (recalc_forces and reuse_forces != INTEG_REUSE_FORCES_ALWAYS)) {
     ESPRESSO_PROFILER_MARK_BEGIN("Initial Force Calculation");
     lb_lbcoupling_deactivate();
 
@@ -287,11 +282,17 @@ int integrate(int n_steps, int reuse_forces) {
   lb_lbcoupling_activate();
 
   if (check_runtime_errors(comm_cart))
-    return 0;
+    return INTEG_ERROR_RUNTIME;
 
   // Keep track of the number of Verlet updates (i.e. particle resorts)
   int n_verlet_updates = 0;
 
+  // Keep track of whether an interrupt signal was caught (only in singleton
+  // mode, since signal handlers are unreliable with more than 1 MPI rank)
+  auto const singleton_mode = comm_cart.size() == 1;
+  auto caught_sigint = false;
+  auto caught_error = false;
+
 #ifdef VALGRIND_MARKERS
   CALLGRIND_START_INSTRUMENTATION;
 #endif
@@ -384,12 +385,14 @@ int integrate(int n_steps, int reuse_forces) {
 
     integrated_steps++;
 
-    if (check_runtime_errors(comm_cart))
+    if (check_runtime_errors(comm_cart)) {
+      caught_error = true;
       break;
+    }
 
     // Check if SIGINT has been caught.
-    if (ctrl_C == 1) {
-      notify_sig_int();
+    if (singleton_mode and ctrl_C == 1) {
+      caught_sigint = true;
       break;
     }
 
@@ -416,40 +419,50 @@ int integrate(int n_steps, int reuse_forces) {
     synchronize_npt_state();
   }
 #endif
+  if (caught_sigint) {
+    ctrl_C = 0;
+    return INTEG_ERROR_SIGINT;
+  }
+  if (caught_error) {
+    return INTEG_ERROR_RUNTIME;
+  }
   return integrated_steps;
 }
 
-int python_integrate(int n_steps, bool recalc_forces_par,
-                     bool reuse_forces_par) {
+int integrate_with_signal_handler(int n_steps, int reuse_forces,
+                                  bool update_accumulators) {
 
   assert(n_steps >= 0);
 
   // Override the signal handler so that the integrator obeys Ctrl+C
   SignalHandler sa(SIGINT, [](int) { ctrl_C = 1; });
 
-  int reuse_forces = reuse_forces_par;
-
-  if (recalc_forces_par) {
-    if (reuse_forces) {
-      runtimeErrorMsg() << "cannot reuse old forces and recalculate forces";
-    }
-    reuse_forces = -1;
+  if (not update_accumulators or n_steps == 0) {
+    return integrate(n_steps, reuse_forces);
   }
 
+  auto const is_head_node = comm_cart.rank() == 0;
+
   /* if skin wasn't set, do an educated guess now */
   if (!skin_set) {
     auto const max_cut = maximal_cutoff(n_nodes);
     if (max_cut <= 0.0) {
-      runtimeErrorMsg()
-          << "cannot automatically determine skin, please set it manually";
-      return ES_ERROR;
+      if (is_head_node) {
+        throw std::runtime_error(
+            "cannot automatically determine skin, please set it manually");
+      }
+      return INTEG_ERROR_RUNTIME;
     }
     /* maximal skin that can be used without resorting is the maximal
      * range of the cell system minus what is needed for interactions. */
-    auto const new_skin =
-        std::min(0.4 * max_cut,
-                 *boost::min_element(cell_structure.max_cutoff()) - max_cut);
-    mpi_call_all(mpi_set_skin_local, new_skin);
+    auto const max_range = *boost::min_element(::cell_structure.max_cutoff());
+    auto const new_skin = std::min(0.4 * max_cut, max_range - max_cut);
+    ::set_skin(new_skin);
+  }
+
+  // re-acquire MpiCallbacks listener on worker nodes
+  if (not is_head_node) {
+    return 0;
   }
 
   using Accumulators::auto_update;
@@ -459,47 +472,23 @@ int python_integrate(int n_steps, bool recalc_forces_par,
     /* Integrate to either the next accumulator update, or the
      * end, depending on what comes first. */
     auto const steps = std::min((n_steps - i), auto_update_next_update());
-    if (mpi_integrate(steps, reuse_forces))
-      return ES_ERROR;
+    auto const retval = mpi_call(Communication::Result::main_rank, integrate,
+                                 steps, reuse_forces);
+    if (retval < 0) {
+      return retval; // propagate error code
+    }
 
-    reuse_forces = 1;
+    reuse_forces = INTEG_REUSE_FORCES_ALWAYS;
 
     auto_update(steps);
 
     i += steps;
   }
 
-  if (n_steps == 0) {
-    if (mpi_integrate(0, reuse_forces))
-      return ES_ERROR;
-  }
-
-  return ES_OK;
+  return 0;
 }
 
-static int mpi_steepest_descent_local(int steps) {
-  return integrate(steps, -1);
-}
-
-REGISTER_CALLBACK_MAIN_RANK(mpi_steepest_descent_local)
-
-int mpi_steepest_descent(int steps) {
-  return mpi_call(Communication::Result::main_rank, mpi_steepest_descent_local,
-                  steps);
-}
-
-static int mpi_integrate_local(int n_steps, int reuse_forces) {
-  integrate(n_steps, reuse_forces);
-
-  return check_runtime_errors_local();
-}
-
-REGISTER_CALLBACK_REDUCTION(mpi_integrate_local, std::plus<int>())
-
-int mpi_integrate(int n_steps, int reuse_forces) {
-  return mpi_call(Communication::Result::reduction, std::plus<int>(),
-                  mpi_integrate_local, n_steps, reuse_forces);
-}
+REGISTER_CALLBACK_MAIN_RANK(integrate)
 
 double interaction_range() {
   /* Consider skin only if there are actually interactions */
@@ -524,14 +513,12 @@ void set_time_step(double value) {
   on_timestep_change();
 }
 
-void mpi_set_skin_local(double value) {
+void set_skin(double value) {
   ::skin = value;
-  skin_set = true;
+  ::skin_set = true;
   on_skin_change();
 }
 
-REGISTER_CALLBACK(mpi_set_skin_local)
-
 void set_time(double value) {
   ::sim_time = value;
   ::recalc_forces = true;

src/core/integrate.hpp

@@ -36,6 +36,22 @@
 #define INTEG_METHOD_SD 7
 /**@}*/
 
+/** \name Integrator error codes */
+/**@{*/
+#define INTEG_ERROR_RUNTIME -1
+#define INTEG_ERROR_SIGINT -2
+/**@}*/
+
+/** \name Integrator flags */
+/**@{*/
+/// recalculate forces unconditionally (mostly used for timing)
+#define INTEG_REUSE_FORCES_NEVER -1
+/// recalculate forces if @ref recalc_forces is set
+#define INTEG_REUSE_FORCES_CONDITIONALLY 0
+/// do not recalculate forces (mostly when reading checkpoints with forces)
+#define INTEG_REUSE_FORCES_ALWAYS 1
+/**@}*/
+
 /** Switch determining which integrator to use. */
 extern int integ_switch;
 
@@ -45,9 +61,6 @@ extern double skin;
 /** If true, the forces will be recalculated before the next integration. */
 extern bool recalc_forces;
 
-/** Communicate signal handling to the Python interpreter */
-extern bool set_py_interrupt;
-
 double interaction_range();
 
 /** Check integrator parameters and incompatibilities between the integrator
@@ -57,13 +70,7 @@ void integrator_sanity_checks();
 
 /** Integrate equations of motion
  *  @param n_steps       Number of integration steps, can be zero
- *  @param reuse_forces  Decide when to re-calculate forces:
- *                       - -1: recalculate forces unconditionally
- *                         (mostly used for timing)
- *                       - 0: recalculate forces if @ref recalc_forces is set,
- *                         meaning it is probably necessary
- *                       - 1: do not recalculate forces (mostly when reading
- *                         checkpoints with forces)
+ *  @param reuse_forces  Decide when to re-calculate forces
  *
  *  @details This function calls two hooks for propagation kernels such as
  *  velocity verlet, velocity verlet + npt box changes, and steepest_descent.
@@ -90,37 +97,12 @@ void integrator_sanity_checks();
  *  High-level documentation of the integration and thermostatting schemes
  *  can be found in doc/sphinx/system_setup.rst and /doc/sphinx/running.rst
  *
- *  @return number of steps that have been integrated
+ *  @return number of steps that have been integrated, or a negative error code
  */
 int integrate(int n_steps, int reuse_forces);
 
-/** @brief Run the integration loop. Can be interrupted with Ctrl+C.
- *
- *  @param n_steps        Number of integration steps, can be zero
- *  @param recalc_forces  Whether to recalculate forces
- *  @param reuse_forces   Whether to re-use forces
- *  @retval ES_OK on success
- *  @retval ES_ERROR on error
- */
-int python_integrate(int n_steps, bool recalc_forces, bool reuse_forces);
-
-/** Start integrator.
- *  @param n_steps       how many steps to do.
- *  @param reuse_forces  whether to trust the old forces for the first half step
- *  @return nonzero on error
- */
-int mpi_integrate(int n_steps, int reuse_forces);
-
-/** Steepest descent main integration loop
- *
- *  Integration stops when the maximal force is lower than the user limit
- *  @ref SteepestDescentParameters::f_max "f_max" or when the maximal number
- *  of steps @p steps is reached.
- *
- *  @param steps Maximal number of integration steps
- *  @return number of integrated steps
- */
-int mpi_steepest_descent(int steps);
+int integrate_with_signal_handler(int n_steps, int reuse_forces,
+                                  bool update_accumulators);
 
 /** Get @c verlet_reuse */
 double get_verlet_reuse();
@@ -138,7 +120,7 @@ void increment_sim_time(double amount);
 void set_time_step(double value);
 
 /** @brief Set new skin. */
-void mpi_set_skin_local(double value);
+void set_skin(double value);
 
 /** @brief Set the simulation time. */
 void set_time(double value);

src/core/tuning.cpp

@@ -66,26 +66,25 @@ static void check_statistics(Utils::Statistics::RunningAverage<double> &acc) {
   }
 }
 
-static void throw_on_error() {
-  auto const n_errors = check_runtime_errors_local();
-  if (boost::mpi::all_reduce(::comm_cart, n_errors, std::plus<>()) != 0) {
+static void run_full_force_calc(int reuse_forces) {
+  auto const error_code = integrate(0, reuse_forces);
+  if (error_code == INTEG_ERROR_RUNTIME) {
     throw TuningFailed{};
   }
 }
 
 double benchmark_integration_step(int int_steps) {
   Utils::Statistics::RunningAverage<double> running_average;
 
-  integrate(0, 0);
-  throw_on_error();
+  // check if the system can be integrated with the current parameters
+  run_full_force_calc(INTEG_REUSE_FORCES_CONDITIONALLY);
 
-  /* perform force calculation test */
+  // measure force calculation time
   for (int i = 0; i < int_steps; i++) {
     auto const tick = MPI_Wtime();
-    integrate(0, -1);
+    run_full_force_calc(INTEG_REUSE_FORCES_NEVER);
     auto const tock = MPI_Wtime();
     running_average.add_sample((tock - tick));
-    throw_on_error();
   }
 
   if (this_node == 0) {
@@ -98,11 +97,6 @@ double benchmark_integration_step(int int_steps) {
   return retval;
 }
 
-static bool integration_failed() {
-  auto const count_local = check_runtime_errors_local();
-  return boost::mpi::all_reduce(::comm_cart, count_local, std::plus<>()) > 0;
-}
-
 /**
  * \brief Time the integration.
  * This times the integration and
@@ -112,16 +106,16 @@ static bool integration_failed() {
  * @return Time per integration in ms.
  */
 static double time_calc(int int_steps) {
-  integrate(0, 0);
-  if (integration_failed()) {
+  auto const error_code_init = integrate(0, INTEG_REUSE_FORCES_CONDITIONALLY);
+  if (error_code_init == INTEG_ERROR_RUNTIME) {
     return -1;
   }
 
   /* perform force calculation test */
   auto const tick = MPI_Wtime();
-  integrate(int_steps, -1);
+  auto const error_code = integrate(int_steps, INTEG_REUSE_FORCES_NEVER);
   auto const tock = MPI_Wtime();
-  if (integration_failed()) {
+  if (error_code == INTEG_ERROR_RUNTIME) {
     return -1;
   }
 
@@ -147,10 +141,10 @@ void tune_skin(double min_skin, double max_skin, double tol, int int_steps,
     b = max_permissible_skin;
 
   while (fabs(a - b) > tol) {
-    mpi_set_skin_local(a);
+    ::set_skin(a);
     auto const time_a = time_calc(int_steps);
 
-    mpi_set_skin_local(b);
+    ::set_skin(b);
     auto const time_b = time_calc(int_steps);
 
     if (time_a > time_b) {
@@ -160,5 +154,5 @@ void tune_skin(double min_skin, double max_skin, double tol, int int_steps,
     }
   }
   auto const new_skin = 0.5 * (a + b);
-  mpi_set_skin_local(new_skin);
+  ::set_skin(new_skin);
 }