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/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/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..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.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) 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)