Skip to content

Commit

Permalink
Rewrite integrators script interface as Python code (#4713)
Browse files Browse the repository at this point in the history 
Fixes #4712

Description of changes:
- rewrite integrator as Python (partial fix for #4542)
- remove two MPI callbacks (partial fix for #4614)
- remove a global variable (partial fix for #2628)
  • Loading branch information
@kodiakhq
kodiakhq[bot] authored Apr 24, 2023
2 parents e709c1a + 752cb97 commit 8b0c0bf
Show file tree
Hide file tree
Showing 29 changed files with 291 additions and 215 deletions.
2 changes: 1 addition & 1 deletion src/core/EspressoSystemStandAlone.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}
115 changes: 51 additions & 64 deletions src/core/integrate.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 {
Expand Down Expand Up @@ -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();

Expand All @@ -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
Expand Down Expand Up @@ -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;
}

Expand All @@ -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;
Expand All @@ -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 */
Expand All @@ -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;
Expand Down
60 changes: 21 additions & 39 deletions src/core/integrate.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand All @@ -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
Expand All @@ -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.
Expand All @@ -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();
Expand All @@ -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);
Expand Down
34 changes: 14 additions & 20 deletions src/core/tuning.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand All @@ -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
Expand All @@ -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;
}

Expand All @@ -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) {
Expand All @@ -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);
}
Loading

0 comments on commit 8b0c0bf

Please sign in to comment.