Rewrite integrators script interface as Python code

The nogil attribute is required to make the integrator run smoothly
with threads (e.g. from the visualizer). Functions marked with the
nogil attribute can throw C++ exceptions since Cython 0.19. The
Cython functions that handle signals can be replaced with the
`signal.raise_signal()` function introduced in Python 3.8.

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,7 +253,7 @@ 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)) {
@@ -287,11 +281,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 +384,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 +418,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,8 +471,11 @@ 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;
 
@@ -469,37 +484,10 @@ int python_integrate(int n_steps, bool recalc_forces_par,
     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 +512,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,12 @@
 #define INTEG_METHOD_SD 7
 /**@}*/
 
+/** \name Integrator error codes */
+/**@{*/
+#define INTEG_ERROR_RUNTIME -1
+#define INTEG_ERROR_SIGINT -2
+/**@}*/
+
 /** Switch determining which integrator to use. */
 extern int integ_switch;
 
@@ -45,9 +51,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
@@ -90,37 +93,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 +116,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,22 @@ 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) {
-    throw TuningFailed{};
-  }
-}
-
 double benchmark_integration_step(int int_steps) {
   Utils::Statistics::RunningAverage<double> running_average;
 
-  integrate(0, 0);
-  throw_on_error();
+  if (integrate(0, 0) == INTEG_ERROR_RUNTIME) {
+    throw TuningFailed{};
+  }
 
   /* perform force calculation test */
   for (int i = 0; i < int_steps; i++) {
     auto const tick = MPI_Wtime();
-    integrate(0, -1);
+    auto const error_code = integrate(0, -1);
     auto const tock = MPI_Wtime();
+    if (error_code == INTEG_ERROR_RUNTIME) {
+      throw TuningFailed{};
+    }
     running_average.add_sample((tock - tick));
-    throw_on_error();
   }
 
   if (this_node == 0) {
@@ -147,10 +143,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 +156,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);
 }