From 62069c79f06912502b6848a6de4d61417492e66e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Mon, 24 Apr 2023 10:56:45 +0200 Subject: [PATCH 1/2] 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 | 2 +- src/core/integrate.cpp | 115 ++++++++---------- src/core/integrate.hpp | 60 ++++----- src/core/tuning.cpp | 34 +++--- .../EspressoSystemStandAlone_test.cpp | 4 +- src/core/unit_tests/Verlet_list_test.cpp | 6 +- src/python/espressomd/integrate.pxd | 33 ----- .../{integrate.pyx => integrate.py} | 41 +++---- src/python/espressomd/lb.pxd | 3 + src/python/espressomd/lb.pyx | 2 +- src/python/espressomd/script_interface.pxd | 1 + src/python/espressomd/script_interface.pyx | 19 ++- src/python/espressomd/utils.pyx | 3 +- src/script_interface/ObjectHandle.hpp | 2 +- .../cell_system/CellSystem.cpp | 2 +- .../integrators/CMakeLists.txt | 1 + .../integrators/Integrator.cpp | 64 ++++++++++ .../integrators/Integrator.hpp | 3 + .../integrators/SteepestDescent.cpp | 13 ++ .../integrators/SteepestDescent.hpp | 1 + .../integrators/VelocityVerlet.hpp | 2 +- testsuite/python/integrator_exceptions.py | 11 +- .../python/integrator_steepest_descent.py | 4 + testsuite/python/interactions_bonded.py | 2 +- testsuite/python/scafacos_interface.py | 1 + testsuite/python/sigint.py | 2 +- testsuite/python/test_checkpoint.py | 3 + 27 files changed, 231 insertions(+), 203 deletions(-) delete mode 100644 src/python/espressomd/integrate.pxd rename src/python/espressomd/{integrate.pyx => integrate.py} (87%) create mode 100644 src/script_interface/integrators/Integrator.cpp diff --git a/src/core/EspressoSystemStandAlone.cpp b/src/core/EspressoSystemStandAlone.cpp index 3dee8e90b87..024fd403b33 100644 --- a/src/core/EspressoSystemStandAlone.cpp +++ b/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); } diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index 575e30f53d1..a229ccecde8 100644 --- a/src/core/integrate.cpp +++ b/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 mpi_integrate(int n_steps, int reuse_forces) { - return mpi_call(Communication::Result::reduction, std::plus(), - 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; diff --git a/src/core/integrate.hpp b/src/core/integrate.hpp index adfc53f0347..50658357631 100644 --- a/src/core/integrate.hpp +++ b/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); diff --git a/src/core/tuning.cpp b/src/core/tuning.cpp index bf9e198b681..7f84be0497a 100644 --- a/src/core/tuning.cpp +++ b/src/core/tuning.cpp @@ -66,9 +66,9 @@ static void check_statistics(Utils::Statistics::RunningAverage &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{}; } } @@ -76,16 +76,15 @@ static void throw_on_error() { double benchmark_integration_step(int int_steps) { Utils::Statistics::RunningAverage 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); } diff --git a/src/core/unit_tests/EspressoSystemStandAlone_test.cpp b/src/core/unit_tests/EspressoSystemStandAlone_test.cpp index 7e1ec3083f5..f1e156175ab 100644 --- a/src/core/unit_tests/EspressoSystemStandAlone_test.cpp +++ b/src/core/unit_tests/EspressoSystemStandAlone_test.cpp @@ -312,7 +312,7 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory) { reset_particle_positions(); // recalculate forces without propagating the system - integrate(0, 0); + integrate(0, INTEG_REUSE_FORCES_CONDITIONALLY); // particles are arranged in a right triangle auto const p1_opt = copy_particle_to_head_node(comm, pid1); @@ -355,7 +355,7 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory) { expected[pid] = p.pos(); } } - integrate(1, 0); + integrate(1, INTEG_REUSE_FORCES_CONDITIONALLY); for (auto pid : pids) { auto const p_opt = copy_particle_to_head_node(comm, pid); if (rank == 0) { diff --git a/src/core/unit_tests/Verlet_list_test.cpp b/src/core/unit_tests/Verlet_list_test.cpp index 2d7866e3baf..2e8533f9ee0 100644 --- a/src/core/unit_tests/Verlet_list_test.cpp +++ b/src/core/unit_tests/Verlet_list_test.cpp @@ -224,7 +224,7 @@ BOOST_DATA_TEST_CASE_F(ParticleFactory, verlet_list_update, // integrate until both particles are closer than cutoff { - integrate(11, 0); + integrate(11, INTEG_REUSE_FORCES_CONDITIONALLY); auto const p1_opt = copy_particle_to_head_node(comm, pid1); auto const p2_opt = copy_particle_to_head_node(comm, pid2); if (rank == 0) { @@ -236,7 +236,7 @@ BOOST_DATA_TEST_CASE_F(ParticleFactory, verlet_list_update, // check forces and Verlet update { - integrate(1, 0); + integrate(1, INTEG_REUSE_FORCES_CONDITIONALLY); auto const p1_opt = copy_particle_to_head_node(comm, pid1); #ifdef EXTERNAL_FORCES auto const p2_opt = copy_particle_to_head_node(comm, pid2); @@ -273,7 +273,7 @@ BOOST_DATA_TEST_CASE_F(ParticleFactory, verlet_list_update, } } { - integrate(0, 0); + integrate(0, INTEG_REUSE_FORCES_CONDITIONALLY); auto const p3_opt = copy_particle_to_head_node(comm, pid3); if (rank == 0) { auto const &p3 = *p3_opt; diff --git a/src/python/espressomd/integrate.pxd b/src/python/espressomd/integrate.pxd deleted file mode 100644 index 91c372f321d..00000000000 --- a/src/python/espressomd/integrate.pxd +++ /dev/null @@ -1,33 +0,0 @@ -# -# Copyright (C) 2013-2022 The ESPResSo project -# -# This file is part of ESPResSo. -# -# ESPResSo is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# ESPResSo is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . -# -from libcpp cimport bool as cbool - -include "myconfig.pxi" - -cdef extern from "integrate.hpp": - double get_time_step() - -cdef extern from "integrate.hpp" nogil: - cdef int python_integrate(int n_steps, cbool recalc_forces, int reuse_forces) - cdef int mpi_steepest_descent(int max_steps) - cdef extern cbool set_py_interrupt - -cdef inline int _integrate(int nSteps, cbool recalc_forces, int reuse_forces): - with nogil: - return python_integrate(nSteps, recalc_forces, reuse_forces) diff --git a/src/python/espressomd/integrate.pyx b/src/python/espressomd/integrate.py similarity index 87% rename from src/python/espressomd/integrate.pyx rename to src/python/espressomd/integrate.py index 5c5b6e1ee7c..519dff9f702 100644 --- a/src/python/espressomd/integrate.pyx +++ b/src/python/espressomd/integrate.py @@ -16,12 +16,9 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . # -from cpython.exc cimport PyErr_CheckSignals, PyErr_SetInterrupt -include "myconfig.pxi" -from . import utils -from . cimport integrate from .script_interface import ScriptInterfaceHelper, script_interface_register from .code_features import assert_features +import signal @script_interface_register @@ -115,22 +112,17 @@ def run(self, steps=1, recalc_forces=False, reuse_forces=False): Reuse the forces from previous time step. """ - utils.check_type_or_throw_except(steps, 1, int, "steps must be an int") - utils.check_type_or_throw_except( - recalc_forces, 1, bool, "recalc_forces has to be a bool") - utils.check_type_or_throw_except( - reuse_forces, 1, bool, "reuse_forces has to be a bool") - if steps < 0: - raise ValueError("steps must be positive") + status = self.call_method( + "integrate", + with_nogil=True, + steps=steps, + recalc_forces=recalc_forces, + reuse_forces=reuse_forces) + self.handle_sigint(status) - _integrate(steps, recalc_forces, reuse_forces) - - if integrate.set_py_interrupt: - PyErr_SetInterrupt() - integrate.set_py_interrupt = False - PyErr_CheckSignals() - - utils.handle_errors("Encountered errors during integrate") + def handle_sigint(self, error_code): + if error_code == -2: + signal.raise_signal(signal.Signals.SIGINT) @script_interface_register @@ -161,7 +153,7 @@ class SteepestDescent(Integrator): _so_name = "Integrators::SteepestDescent" _so_creation_policy = "GLOBAL" - def run(self, steps=1, **kwargs): + def run(self, steps=1): """ Run the steepest descent. @@ -176,12 +168,9 @@ def run(self, steps=1, **kwargs): Number of integrated steps. """ - utils.check_type_or_throw_except(steps, 1, int, "steps must be an int") - assert steps >= 0, "steps has to be positive" - - integrated = mpi_steepest_descent(steps) - - utils.handle_errors("Encountered errors during integrate") + integrated = self.call_method( + "integrate", with_nogil=True, steps=steps) + self.handle_sigint(integrated) return integrated diff --git a/src/python/espressomd/lb.pxd b/src/python/espressomd/lb.pxd index 98ebec7a146..3ac79be7c11 100644 --- a/src/python/espressomd/lb.pxd +++ b/src/python/espressomd/lb.pxd @@ -119,6 +119,9 @@ cdef extern from "grid_based_algorithms/lb_interpolation.hpp" namespace "Interpo cdef InterpolationOrder linear cdef InterpolationOrder quadratic +cdef extern from "integrate.hpp": + double get_time_step() + ############################################## # # Wrapper-functions to handle unit conversions diff --git a/src/python/espressomd/lb.pyx b/src/python/espressomd/lb.pyx index dceef548081..e2875002e33 100644 --- a/src/python/espressomd/lb.pyx +++ b/src/python/espressomd/lb.pyx @@ -28,7 +28,7 @@ from . import highlander from . import utils from . cimport utils from .utils cimport Vector3i, Vector3d, Vector6d, Vector19d -from .integrate cimport get_time_step +from .lb cimport get_time_step cdef class FluidActor: diff --git a/src/python/espressomd/script_interface.pxd b/src/python/espressomd/script_interface.pxd index 112a05de941..0e10ea0ac5f 100644 --- a/src/python/espressomd/script_interface.pxd +++ b/src/python/espressomd/script_interface.pxd @@ -45,6 +45,7 @@ cdef extern from "script_interface/ScriptInterface.hpp" namespace "ScriptInterfa Variant get_parameter(const string & name) except + void set_parameter(const string & name, const Variant & value) except + Variant call_method(const string & name, const VariantMap & parameters) except + + Variant call_method_nogil "call_method"(const string & name, const VariantMap & parameters) nogil except + string_ref name() cdef extern from "script_interface/ContextManager.hpp" namespace "ScriptInterface::ContextManager": diff --git a/src/python/espressomd/script_interface.pyx b/src/python/espressomd/script_interface.pyx index b669b6e8abf..940c1f3e2cb 100644 --- a/src/python/espressomd/script_interface.pyx +++ b/src/python/espressomd/script_interface.pyx @@ -136,19 +136,23 @@ cdef class PScriptInterface: self.sip = sip - def call_method(self, method, handle_errors_message=None, **kwargs): + def call_method(self, method, handle_errors_message=None, + with_nogil=False, **kwargs): """ Call a method of the core class. Parameters ---------- - method : Creation policy. + method : :obj:`str` Name of the core method. handle_errors_message : :obj:`str`, optional Custom error message for runtime errors raised in a MPI context. + with_nogil : :obj:`bool`, optional + Run the core method without the GIL if ``True``. **kwargs Arguments for the method. """ + cdef ObjectHandle * handle = self.sip.get() cdef VariantMap parameters cdef Variant value @@ -156,7 +160,16 @@ cdef class PScriptInterface: parameters[utils.to_char_pointer(name)] = python_object_to_variant( kwargs[name]) - value = self.sip.get().call_method(utils.to_char_pointer(method), parameters) + # the internal buffer of a cython bytestring object can be accessed as + # a raw char pointer, but then the bytestring object must be kept alive + method_name_bytes_counted_reference = utils.to_char_pointer(method) + cdef char * method_name_char = method_name_bytes_counted_reference + + if with_nogil: + with nogil: + value = handle.call_method_nogil(method_name_char, parameters) + else: + value = handle.call_method(method_name_char, parameters) res = variant_to_python_object(value) if handle_errors_message is None: handle_errors_message = f"while calling method {method}()" diff --git a/src/python/espressomd/utils.pyx b/src/python/espressomd/utils.pyx index 987086c0baa..f1da1760064 100644 --- a/src/python/espressomd/utils.pyx +++ b/src/python/espressomd/utils.pyx @@ -70,7 +70,8 @@ cpdef check_type_or_throw_except(x, n, t, msg): def to_char_pointer(s): """ - Returns a char pointer which contains the information of the provided python string. + Returns a Cython bytes object which contains the information of the provided + Python string. Cython bytes objects implicitly cast to raw char pointers. Parameters ---------- diff --git a/src/script_interface/ObjectHandle.hpp b/src/script_interface/ObjectHandle.hpp index b28742d687c..a4325ee9129 100644 --- a/src/script_interface/ObjectHandle.hpp +++ b/src/script_interface/ObjectHandle.hpp @@ -133,7 +133,7 @@ class ObjectHandle { protected: /** - * @brief Local implementation of @c do_call_method. + * @brief Local implementation of @c call_method. * * If not overridden by the implementation, this does nothing. */ diff --git a/src/script_interface/cell_system/CellSystem.cpp b/src/script_interface/cell_system/CellSystem.cpp index 021b7479cb3..db39a9ad81b 100644 --- a/src/script_interface/cell_system/CellSystem.cpp +++ b/src/script_interface/cell_system/CellSystem.cpp @@ -96,7 +96,7 @@ CellSystem::CellSystem() { } throw Exception(""); } - mpi_set_skin_local(new_skin); + ::set_skin(new_skin); }, []() { return ::skin; }}, {"decomposition_type", AutoParameter::read_only, diff --git a/src/script_interface/integrators/CMakeLists.txt b/src/script_interface/integrators/CMakeLists.txt index a7d390ca37f..fe3c5359853 100644 --- a/src/script_interface/integrators/CMakeLists.txt +++ b/src/script_interface/integrators/CMakeLists.txt @@ -21,6 +21,7 @@ target_sources( espresso_script_interface PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/initialize.cpp ${CMAKE_CURRENT_SOURCE_DIR}/BrownianDynamics.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/Integrator.cpp ${CMAKE_CURRENT_SOURCE_DIR}/IntegratorHandle.cpp ${CMAKE_CURRENT_SOURCE_DIR}/SteepestDescent.cpp ${CMAKE_CURRENT_SOURCE_DIR}/StokesianDynamics.cpp diff --git a/src/script_interface/integrators/Integrator.cpp b/src/script_interface/integrators/Integrator.cpp new file mode 100644 index 00000000000..a4c469d1f7d --- /dev/null +++ b/src/script_interface/integrators/Integrator.cpp @@ -0,0 +1,64 @@ +/* + * Copyright (C) 2022-2023 The ESPResSo project + * + * This file is part of ESPResSo. + * + * ESPResSo is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * ESPResSo is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include "Integrator.hpp" + +#include "core/integrate.hpp" + +#include "script_interface/ScriptInterface.hpp" + +#include +#include + +namespace ScriptInterface { +namespace Integrators { + +Variant Integrator::integrate(VariantMap const ¶ms) { + auto const steps = get_value(params, "steps"); + auto const recalc_forces_flag = get_value(params, "recalc_forces"); + auto const reuse_forces_flag = get_value(params, "reuse_forces"); + + context()->parallel_try_catch([&]() { + if (steps < 0) { + throw std::domain_error("Parameter 'steps' must be positive"); + } + if (recalc_forces_flag and reuse_forces_flag) { + throw std::invalid_argument( + "cannot reuse old forces and recalculate forces"); + } + }); + + auto constexpr update_accumulators = true; + auto const reuse_forces = reuse_forces_flag + ? INTEG_REUSE_FORCES_ALWAYS + : INTEG_REUSE_FORCES_CONDITIONALLY; + return ::integrate_with_signal_handler(steps, reuse_forces, + update_accumulators); +} + +Variant Integrator::do_call_method(std::string const &name, + VariantMap const ¶ms) { + if (name == "integrate") { + return integrate(params); + } + return none; +} + +} // namespace Integrators +} // namespace ScriptInterface diff --git a/src/script_interface/integrators/Integrator.hpp b/src/script_interface/integrators/Integrator.hpp index 85f1f05af15..e2930901563 100644 --- a/src/script_interface/integrators/Integrator.hpp +++ b/src/script_interface/integrators/Integrator.hpp @@ -31,6 +31,9 @@ namespace Integrators { class Integrator : public ObjectHandle { public: + Variant do_call_method(std::string const &name, + VariantMap const ¶ms) override; + virtual Variant integrate(VariantMap const ¶ms); virtual void activate() const = 0; }; diff --git a/src/script_interface/integrators/SteepestDescent.cpp b/src/script_interface/integrators/SteepestDescent.cpp index 18f8bb180b3..5a243197ca4 100644 --- a/src/script_interface/integrators/SteepestDescent.cpp +++ b/src/script_interface/integrators/SteepestDescent.cpp @@ -32,6 +32,19 @@ namespace ScriptInterface { namespace Integrators { +Variant SteepestDescent::integrate(VariantMap const ¶ms) { + auto constexpr reuse_forces = INTEG_REUSE_FORCES_NEVER; + auto constexpr update_accumulators = false; + auto const steps = get_value(params, "steps"); + context()->parallel_try_catch([&]() { + if (steps < 0) { + throw std::domain_error("Parameter 'steps' must be positive"); + } + }); + return ::integrate_with_signal_handler(steps, reuse_forces, + update_accumulators); +} + SteepestDescent::SteepestDescent() { add_parameters({ {"f_max", AutoParameter::read_only, diff --git a/src/script_interface/integrators/SteepestDescent.hpp b/src/script_interface/integrators/SteepestDescent.hpp index 6e6aa6c80e3..73d80da6dfc 100644 --- a/src/script_interface/integrators/SteepestDescent.hpp +++ b/src/script_interface/integrators/SteepestDescent.hpp @@ -40,6 +40,7 @@ class SteepestDescent : public AutoParameters { SteepestDescent(); void do_construct(VariantMap const ¶ms) override; + Variant integrate(VariantMap const ¶ms) override; void activate() const override; ::SteepestDescentParameters const &get_instance() const { diff --git a/src/script_interface/integrators/VelocityVerlet.hpp b/src/script_interface/integrators/VelocityVerlet.hpp index 30ae3dc3d25..693b9039bdb 100644 --- a/src/script_interface/integrators/VelocityVerlet.hpp +++ b/src/script_interface/integrators/VelocityVerlet.hpp @@ -28,7 +28,7 @@ namespace ScriptInterface { namespace Integrators { -class VelocityVerlet : public Integrator { +class VelocityVerlet : public AutoParameters { void activate() const override; }; diff --git a/testsuite/python/integrator_exceptions.py b/testsuite/python/integrator_exceptions.py index 1e8a7c6b648..b5b5d8efa74 100644 --- a/testsuite/python/integrator_exceptions.py +++ b/testsuite/python/integrator_exceptions.py @@ -24,11 +24,11 @@ import unittest_decorators as utx -class Integrator_test(ut.TestCase): +class Test(ut.TestCase): system = espressomd.System(box_l=[1., 1., 1.]) system.time_step = 0.01 - msg = 'Encountered errors during integrate: ERROR: ' + msg = r'while calling method integrate\(\): ERROR: ' def setUp(self): self.system.part.add(pos=(0, 0, 0)) @@ -47,13 +47,14 @@ def test_00_common_interface(self): self.system.time_step = -2. with self.assertRaisesRegex(ValueError, 'time_step must be > 0.'): self.system.time_step = 0. - with self.assertRaisesRegex(ValueError, 'steps must be positive'): + with self.assertRaisesRegex(ValueError, "Parameter 'steps' must be positive"): self.system.integrator.run(steps=-1) - with self.assertRaisesRegex(Exception, self.msg + 'cannot automatically determine skin, please set it manually'): + with self.assertRaisesRegex(RuntimeError, 'cannot automatically determine skin, please set it manually'): self.system.integrator.run() self.system.cell_system.skin = 0.4 - with self.assertRaisesRegex(Exception, self.msg + 'cannot reuse old forces and recalculate forces'): + with self.assertRaisesRegex(ValueError, 'cannot reuse old forces and recalculate forces'): self.system.integrator.run(recalc_forces=True, reuse_forces=True) + self.assertIsNone(self.system.integrator.integrator.call_method("unk")) if espressomd.has_features("WCA"): wca = self.system.non_bonded_inter[0, 0].wca wca.set_params(epsilon=1., sigma=0.01) diff --git a/testsuite/python/integrator_steepest_descent.py b/testsuite/python/integrator_steepest_descent.py index 1574329d11f..dfb5cf91965 100644 --- a/testsuite/python/integrator_steepest_descent.py +++ b/testsuite/python/integrator_steepest_descent.py @@ -152,6 +152,10 @@ def test_integrator_exceptions(self): with self.assertRaises(RuntimeError): self.system.integrator.set_steepest_descent( f_max=0, gamma=1, max_displacement=-1) + with self.assertRaisesRegex(ValueError, "Parameter 'steps' must be positive"): + self.system.integrator.set_steepest_descent( + f_max=0, gamma=1, max_displacement=0.01) + self.system.integrator.run(steps=-1) def test_integrator_recovery(self): # the system is still in a valid state after a failure diff --git a/testsuite/python/interactions_bonded.py b/testsuite/python/interactions_bonded.py index c4ff9afc7da..3584e2e587b 100644 --- a/testsuite/python/interactions_bonded.py +++ b/testsuite/python/interactions_bonded.py @@ -253,7 +253,7 @@ def run_test(self, bond_instance, force_func, energy_func, min_dist, p_tensor_sim, p_tensor_expected, atol=1E-12) if test_breakage: p2.pos = p1.pos + self.axis * cutoff * 1.01 - with self.assertRaisesRegex(Exception, "Encountered errors during integrate"): + with self.assertRaisesRegex(Exception, r"while calling method integrate\(\)"): self.system.integrator.run(recalc_forces=True, steps=0) if test_same_pos_exception: p2.pos = p1.pos diff --git a/testsuite/python/scafacos_interface.py b/testsuite/python/scafacos_interface.py index 5ad5715b717..8e26a66a773 100644 --- a/testsuite/python/scafacos_interface.py +++ b/testsuite/python/scafacos_interface.py @@ -383,6 +383,7 @@ def test_electrostatics_plus_magnetostatics(self): gamma=0.001, max_displacement=0.01) system.integrator.run(100) + system.integrator.set_vv() # compute forces and energies p3m_E_coulomb, p3m_E_dipoles, p3m_forces, p3m_torques = self.p3m_data() diff --git a/testsuite/python/sigint.py b/testsuite/python/sigint.py index 71553f51479..9467bfe7d7d 100644 --- a/testsuite/python/sigint.py +++ b/testsuite/python/sigint.py @@ -31,7 +31,7 @@ def setUp(self): self.process = subprocess.Popen([sys.executable, script]) def test_signal_handling(self): - self.process.send_signal(signal.SIGINT) + self.process.send_signal(signal.Signals.SIGINT) # Wait for the signal to arrive and one integration step to be finished time.sleep(1) self.assertIsNotNone(self.process.poll()) diff --git a/testsuite/python/test_checkpoint.py b/testsuite/python/test_checkpoint.py index 2174d3b5a53..4c99af5bd01 100644 --- a/testsuite/python/test_checkpoint.py +++ b/testsuite/python/test_checkpoint.py @@ -275,7 +275,10 @@ def test_shape_based_constraints_serialization(self): p3, p4 = system.part.by_ids([3, 4]) old_force = np.copy(p3.f) system.constraints.remove(system.constraints[0]) + old_integrator = system.integrator.integrator + system.integrator.set_vv() system.integrator.run(0, recalc_forces=True) + system.integrator.integrator = old_integrator np.testing.assert_allclose( np.copy(p3.f), -np.copy(p4.f), rtol=1e-4) self.assertGreater(np.linalg.norm(np.copy(p3.f) - old_force), 1e6) From 752cb9738bf3f85844e36b2dcc9cda373f6ad693 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Mon, 24 Apr 2023 13:39:54 +0200 Subject: [PATCH 2/2] tests: Fix broken SIGINT test The old test would send the interrupt signal before the script had a chance to enter the integration loop, and thus the custom signal handler was never tested. The new version waits for the subprocess to write it's ready for the interrupt signal, and the traceback is checked to confirm the signal was re-raised by the script interface. --- testsuite/python/CMakeLists.txt | 20 +++++++++---- testsuite/python/sigint.py | 50 +++++++++++++++++++++++++++----- testsuite/python/sigint_child.py | 4 ++- 3 files changed, 61 insertions(+), 13 deletions(-) diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index 1072822819e..0c63b06706a 100644 --- a/testsuite/python/CMakeLists.txt +++ b/testsuite/python/CMakeLists.txt @@ -23,7 +23,11 @@ Configure a python test case. - Flag arguments: none. + Flag arguments: + + * ``NO_MPI``: run test case without MPI, i.e. in singleton mode. Only use + this option when relevant to the test, e.g. to check signal propagation + (MPI vendors introduce their own implementation-defined signal handlers). Single-value arguments: @@ -50,7 +54,7 @@ #]=======================================================================] function(python_test) - cmake_parse_arguments(TEST "" "FILE;MAX_NUM_PROC;GPU_SLOTS;SUFFIX" + cmake_parse_arguments(TEST "NO_MPI" "FILE;MAX_NUM_PROC;GPU_SLOTS;SUFFIX" "DEPENDS;DEPENDENCIES;LABELS;ARGUMENTS" ${ARGN}) get_filename_component(TEST_NAME ${TEST_FILE} NAME_WE) set(TEST_FILE_CONFIGURED "${CMAKE_CURRENT_BINARY_DIR}/${TEST_NAME}.py") @@ -65,7 +69,13 @@ function(python_test) set(TEST_FILE ${TEST_FILE_CONFIGURED}) if(NOT DEFINED TEST_MAX_NUM_PROC) - set(TEST_MAX_NUM_PROC 4) + if(${TEST_NO_MPI}) + set(TEST_MAX_NUM_PROC 1) + else() + set(TEST_MAX_NUM_PROC 4) + endif() + elseif(${TEST_NO_MPI} AND NOT ${TEST_MAX_NUM_PROC} EQUAL 1) + message(FATAL_ERROR "NO_MPI and MAX_NUM_PROC are mutually exclusive") endif() if(NOT DEFINED TEST_GPU_SLOTS) set(TEST_GPU_SLOTS 0) @@ -77,7 +87,7 @@ function(python_test) set(TEST_NUM_PROC ${TEST_MAX_NUM_PROC}) endif() - if(EXISTS ${MPIEXEC}) + if(EXISTS ${MPIEXEC} AND NOT ${TEST_NO_MPI}) espresso_set_mpiexec_tmpdir(${TEST_NAME}) add_test( NAME ${TEST_NAME} @@ -352,7 +362,7 @@ python_test(FILE lb_momentum_conservation.py MAX_NUM_PROC 1 GPU_SLOTS 1 SUFFIX 1_core) python_test(FILE p3m_electrostatic_pressure.py MAX_NUM_PROC 2 GPU_SLOTS 1) python_test(FILE p3m_madelung.py MAX_NUM_PROC 2 GPU_SLOTS 2 LABELS long) -python_test(FILE sigint.py DEPENDENCIES sigint_child.py MAX_NUM_PROC 1) +python_test(FILE sigint.py DEPENDENCIES sigint_child.py NO_MPI) python_test(FILE lb_density.py MAX_NUM_PROC 1) python_test(FILE observable_chain.py MAX_NUM_PROC 4) python_test(FILE mpiio.py MAX_NUM_PROC 4) diff --git a/testsuite/python/sigint.py b/testsuite/python/sigint.py index 9467bfe7d7d..738d8e9a2f5 100644 --- a/testsuite/python/sigint.py +++ b/testsuite/python/sigint.py @@ -22,19 +22,55 @@ import time import sys import pathlib +import os class SigintTest(ut.TestCase): - def setUp(self): - script = str(pathlib.Path(__file__).parent / 'sigint_child.py') - self.process = subprocess.Popen([sys.executable, script]) + script = str(pathlib.Path(__file__).parent / 'sigint_child.py') + + def check_signal_handling(self, process, sig): + # send signal + process.send_signal(sig) + # capture stderr and return code (negative of signum) + stdout, stderr = process.communicate(input=None, timeout=6.) + assert stdout is None + traceback = stderr.decode() + return_code = process.poll() + signum = -return_code + self.assertEqual(signum, sig.value) + if sig == signal.Signals.SIGTERM: + self.assertEqual(traceback, "") + elif sig == signal.Signals.SIGINT: + self.assertIn(" self.integrator.run(", traceback) + self.assertTrue(traceback.endswith( + " in handle_sigint\n signal.raise_signal(signal.Signals.SIGINT)\nKeyboardInterrupt\n")) def test_signal_handling(self): - self.process.send_signal(signal.Signals.SIGINT) - # Wait for the signal to arrive and one integration step to be finished - time.sleep(1) - self.assertIsNotNone(self.process.poll()) + signals = [signal.Signals.SIGINT, signal.Signals.SIGTERM] + processes = [] + # open asynchronous processes with non-blocking read access on stderr + for _ in range(len(signals)): + process = subprocess.Popen([sys.executable, self.script], + stderr=subprocess.PIPE) + os.set_blocking(process.stderr.fileno(), False) + processes.append(process) + + # wait for the script to reach the integration loop + time.sleep(0.5) + for process, sig in zip(processes, signals): + tick = time.time() + while True: + message = process.stderr.readline().decode() + if message == "start of integration loop\n": + # wait for the script to enter the integrator run method + time.sleep(0.1) + # send signal and check process behavior + self.check_signal_handling(process, sig) + break + tock = time.time() + assert tock - tick < 8., "subprocess timed out" + time.sleep(0.1) if __name__ == '__main__': diff --git a/testsuite/python/sigint_child.py b/testsuite/python/sigint_child.py index ed86ecec3cb..27a79e189a5 100644 --- a/testsuite/python/sigint_child.py +++ b/testsuite/python/sigint_child.py @@ -16,6 +16,7 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . # +import sys import numpy as np import espressomd @@ -26,5 +27,6 @@ for i in range(100): system.part.add(pos=np.random.random() * system.box_l) +print("start of integration loop", file=sys.stderr) while True: - system.integrator.run(1000) + system.integrator.run(100000)