Skip to content

Commit

Permalink
Encapsulate integrator
Browse files Browse the repository at this point in the history
  • Loading branch information
@jngrad
jngrad committed Oct 26, 2023
1 parent 882f20e commit d05263f
Show file tree
Hide file tree
Showing 24 changed files with 147 additions and 170 deletions.
70 changes: 33 additions & 37 deletions src/core/integrate.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,6 @@
#include "virtual_sites.hpp"

#include <boost/mpi/collectives/all_reduce.hpp>
#include <boost/mpi/collectives/reduce.hpp>

#ifdef CALIPER
#include <caliper/cali.h>
Expand Down Expand Up @@ -313,18 +312,17 @@ static void integrator_step_2(ParticleRange const &particles, double kT,
}
}

int integrate(System::System &system, int n_steps, int reuse_forces) {
namespace System {

int System::integrate(int n_steps, int reuse_forces) {
#ifdef CALIPER
CALI_CXX_MARK_FUNCTION;
#endif

auto &cell_structure = *system.cell_structure;
auto &box_geo = *system.box_geo;
auto &bond_breakage = *system.bond_breakage;
auto const time_step = system.get_time_step();
auto const time_step = get_time_step();

// Prepare particle structure and run sanity checks of all active algorithms
system.on_integration_start();
on_integration_start();

// If any method vetoes (e.g. P3M not initialized), immediately bail out
if (check_runtime_errors(comm_cart))
Expand All @@ -343,13 +341,13 @@ int integrate(System::System &system, int n_steps, int reuse_forces) {
#endif

// Communication step: distribute ghost positions
cells_update_ghosts(cell_structure, box_geo, global_ghost_flags());
cells_update_ghosts(*cell_structure, *box_geo, global_ghost_flags());

force_calc(system, time_step, temperature);
force_calc(*this, time_step, temperature);

if (integ_switch != INTEG_METHOD_STEEPEST_DESCENT) {
#ifdef ROTATION
convert_initial_torques(cell_structure.local_particles());
convert_initial_torques(cell_structure->local_particles());
#endif
}

Expand All @@ -375,8 +373,8 @@ int integrate(System::System &system, int n_steps, int reuse_forces) {
auto lb_active = false;
auto ek_active = false;
if (integ_switch != INTEG_METHOD_STEEPEST_DESCENT) {
lb_active = system.lb.is_solver_set();
ek_active = system.ek.is_ready_for_propagation();
lb_active = lb.is_solver_set();
ek_active = ek.is_ready_for_propagation();
}

#ifdef VALGRIND
Expand All @@ -392,19 +390,19 @@ int integrate(System::System &system, int n_steps, int reuse_forces) {
CALI_CXX_MARK_LOOP_ITERATION(integration_loop, step);
#endif

auto particles = cell_structure.local_particles();
auto particles = cell_structure->local_particles();

#ifdef BOND_CONSTRAINT
if (n_rigidbonds)
save_old_position(particles, cell_structure.ghost_particles());
save_old_position(particles, cell_structure->ghost_particles());
#endif

LeesEdwards::update_box_params(box_geo);
LeesEdwards::update_box_params(*box_geo);
bool early_exit = integrator_step_1(particles, time_step);
if (early_exit)
break;

LeesEdwards::run_kernel<LeesEdwards::Push>(box_geo);
LeesEdwards::run_kernel<LeesEdwards::Push>(*box_geo);

#ifdef NPT
if (integ_switch != INTEG_METHOD_NPT_ISO)
Expand All @@ -419,41 +417,39 @@ int integrate(System::System &system, int n_steps, int reuse_forces) {
#ifdef BOND_CONSTRAINT
// Correct particle positions that participate in a rigid/constrained bond
if (n_rigidbonds) {
correct_position_shake(cell_structure, box_geo);
correct_position_shake(*cell_structure, *box_geo);
}
#endif

#ifdef VIRTUAL_SITES
virtual_sites()->update();
#endif

if (cell_structure.get_resort_particles() >= Cells::RESORT_LOCAL)
if (cell_structure->get_resort_particles() >= Cells::RESORT_LOCAL)
n_verlet_updates++;

// Communication step: distribute ghost positions
cells_update_ghosts(cell_structure, box_geo, global_ghost_flags());
cells_update_ghosts(*cell_structure, *box_geo, global_ghost_flags());

particles = cell_structure.local_particles();
particles = cell_structure->local_particles();

force_calc(system, time_step, temperature);
force_calc(*this, time_step, temperature);

#ifdef VIRTUAL_SITES
virtual_sites()->after_force_calc(time_step);
#endif
integrator_step_2(particles, temperature, time_step);
LeesEdwards::run_kernel<LeesEdwards::UpdateOffset>(box_geo);
LeesEdwards::run_kernel<LeesEdwards::UpdateOffset>(*box_geo);
#ifdef BOND_CONSTRAINT
// SHAKE velocity updates
if (n_rigidbonds) {
correct_velocity_shake(cell_structure, box_geo);
correct_velocity_shake(*cell_structure, *box_geo);
}
#endif

// propagate one-step functionalities
if (integ_switch != INTEG_METHOD_STEEPEST_DESCENT) {
if (lb_active and ek_active) {
auto &lb = system.lb;
auto &ek = system.ek;
// assume that they are coupled, which is not necessarily true
auto const md_steps_per_lb_step = calc_md_steps_per_tau(lb.get_tau());
auto const md_steps_per_ek_step = calc_md_steps_per_tau(ek.get_tau());
Expand All @@ -475,7 +471,6 @@ int integrate(System::System &system, int n_steps, int reuse_forces) {
}
lb_lbcoupling_propagate();
} else if (lb_active) {
auto &lb = system.lb;
auto const md_steps_per_lb_step = calc_md_steps_per_tau(lb.get_tau());
lb_skipped_md_steps += 1;
if (lb_skipped_md_steps >= md_steps_per_lb_step) {
Expand All @@ -484,7 +479,6 @@ int integrate(System::System &system, int n_steps, int reuse_forces) {
}
lb_lbcoupling_propagate();
} else if (ek_active) {
auto &ek = system.ek;
auto const md_steps_per_ek_step = calc_md_steps_per_tau(ek.get_tau());
ek_skipped_md_steps += 1;
if (ek_skipped_md_steps >= md_steps_per_ek_step) {
Expand All @@ -498,9 +492,9 @@ int integrate(System::System &system, int n_steps, int reuse_forces) {
#endif

#ifdef COLLISION_DETECTION
handle_collisions(cell_structure);
handle_collisions(*cell_structure);
#endif
bond_breakage.process_queue(system);
bond_breakage->process_queue(*this);
}

integrated_steps++;
Expand All @@ -517,7 +511,7 @@ int integrate(System::System &system, int n_steps, int reuse_forces) {
}

} // for-loop over integration steps
LeesEdwards::update_box_params(box_geo);
LeesEdwards::update_box_params(*box_geo);
#ifdef CALIPER
CALI_CXX_MARK_LOOP_END(integration_loop);
#endif
Expand All @@ -531,7 +525,7 @@ int integrate(System::System &system, int n_steps, int reuse_forces) {
#endif

// Verlet list statistics
cell_structure.update_verlet_stats(n_steps, n_verlet_updates);
cell_structure->update_verlet_stats(n_steps, n_verlet_updates);

#ifdef NPT
if (integ_switch == INTEG_METHOD_NPT_ISO) {
Expand All @@ -548,17 +542,17 @@ int integrate(System::System &system, int n_steps, int reuse_forces) {
return integrated_steps;
}

int integrate_with_signal_handler(System::System &system, int n_steps,
int reuse_forces, bool update_accumulators) {
int System::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; });

/* if skin wasn't set, do an educated guess now */
if (not system.cell_structure->is_verlet_skin_set()) {
if (not cell_structure->is_verlet_skin_set()) {
try {
system.set_verlet_skin_heuristic();
set_verlet_skin_heuristic();
} catch (...) {
if (comm_cart.rank() == 0) {
throw;
Expand All @@ -568,7 +562,7 @@ int integrate_with_signal_handler(System::System &system, int n_steps,
}

if (not update_accumulators or n_steps == 0) {
return integrate(system, n_steps, reuse_forces);
return integrate(n_steps, reuse_forces);
}

using Accumulators::auto_update;
Expand All @@ -579,7 +573,7 @@ int integrate_with_signal_handler(System::System &system, int n_steps,
* end, depending on what comes first. */
auto const steps = std::min((n_steps - i), auto_update_next_update());

auto const local_retval = integrate(system, steps, reuse_forces);
auto const local_retval = integrate(steps, reuse_forces);

// make sure all ranks exit when one rank fails
std::remove_const_t<decltype(local_retval)> global_retval;
Expand All @@ -599,6 +593,8 @@ int integrate_with_signal_handler(System::System &system, int n_steps,
return 0;
}

} // namespace System

double get_time_step() { return System::get_system().get_time_step(); }

double get_sim_time() { return sim_time; }
Expand Down
39 changes: 0 additions & 39 deletions src/core/integrate.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,8 +29,6 @@

#include "config/config.hpp"

#include "system/System.hpp"

#ifdef WALBERLA
#include <string>

Expand Down Expand Up @@ -83,43 +81,6 @@ void walberla_agrid_sanity_checks(std::string method,
double agrid);
#endif // WALBERLA

/** Integrate equations of motion
* @param system system to propagate
* @param n_steps Number of integration steps, can be zero
* @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.
* One hook is called before and one after the force calculation.
* It is up to the propagation kernels to increment the simulation time.
*
* This function propagates the system according to the choice of integrator
* stored in @ref integ_switch. The general structure is:
* - if reuse_forces is zero, recalculate the forces based on the current
* state of the system
* - Loop over the number of simulation steps:
* -# initialization (e.g., RATTLE)
* -# First hook for propagation kernels
* -# Update dependent particles and properties (RATTLE, virtual sites)
* -# Calculate forces for the current state of the system. This includes
* forces added by the Langevin thermostat and the
* Lattice-Boltzmann-particle coupling
* -# Second hook for propagation kernels
* -# Update dependent properties (Virtual sites, RATTLE)
* -# Run single step algorithms (Lattice-Boltzmann propagation, collision
* detection, NpT update)
* - Final update of dependent properties and statistics/counters
*
* 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, or a negative error code
*/
int integrate(System::System &system, int n_steps, int reuse_forces);

int integrate_with_signal_handler(System::System &system, int n_steps,
int reuse_forces, bool update_accumulators);

/** Get time step */
double get_time_step();

Expand Down
37 changes: 37 additions & 0 deletions src/core/system/System.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -133,6 +133,43 @@ class System {
*/
double particle_short_range_energy_contribution(int pid);

/** Integrate equations of motion
* @param n_steps Number of integration steps, can be zero
* @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.
* One hook is called before and one after the force calculation.
* It is up to the propagation kernels to increment the simulation time.
*
* This function propagates the system according to the choice of integrator
* stored in @ref integ_switch. The general structure is:
* - if reuse_forces is zero, recalculate the forces based on the current
* state of the system
* - Loop over the number of simulation steps:
* -# initialization (e.g., RATTLE)
* -# First hook for propagation kernels
* -# Update dependent particles and properties (RATTLE, virtual sites)
* -# Calculate forces for the current state of the system. This includes
* forces added by the Langevin thermostat and the
* Lattice-Boltzmann-particle coupling
* -# Second hook for propagation kernels
* -# Update dependent properties (Virtual sites, RATTLE)
* -# Run single step algorithms (Lattice-Boltzmann propagation, collision
* detection, NpT update)
* - Final update of dependent properties and statistics/counters
*
* 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, or a negative error
* code
*/
int integrate(int n_steps, int reuse_forces);

int integrate_with_signal_handler(int n_steps, int reuse_forces,
bool update_accumulators);

/** \name Hook procedures
* These procedures are called if several significant changes to
* the system happen which may make a reinitialization of subsystems
Expand Down
7 changes: 3 additions & 4 deletions src/core/tuning.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,7 +66,7 @@ static void check_statistics(Utils::Statistics::RunningAverage<double> &acc) {
}

static void run_full_force_calc(System::System &system, int reuse_forces) {
auto const error_code = integrate(system, 0, reuse_forces);
auto const error_code = system.integrate(0, reuse_forces);
if (error_code == INTEG_ERROR_RUNTIME) {
throw TuningFailed{};
}
Expand Down Expand Up @@ -107,15 +107,14 @@ double benchmark_integration_step(System::System &system, int int_steps) {
*/
static double time_calc(System::System &system, int int_steps) {
auto const error_code_init =
integrate(system, 0, INTEG_REUSE_FORCES_CONDITIONALLY);
system.integrate(0, INTEG_REUSE_FORCES_CONDITIONALLY);
if (error_code_init == INTEG_ERROR_RUNTIME) {
return -1;
}

/* perform force calculation test */
auto const tick = MPI_Wtime();
auto const error_code =
integrate(system, int_steps, INTEG_REUSE_FORCES_NEVER);
auto const error_code = system.integrate(int_steps, INTEG_REUSE_FORCES_NEVER);
auto const tock = MPI_Wtime();
if (error_code == INTEG_ERROR_RUNTIME) {
return -1;
Expand Down
4 changes: 2 additions & 2 deletions src/core/tuning.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,8 +35,8 @@ class TuningFailed : public std::runtime_error {

/**
* @brief Benchmark the integration loop.
* Call @ref integrate() several times and measure the elapsed time
* without propagating the system. It therefore doesn't include e.g.
* Call @ref System::System::integrate() several times and measure the elapsed
* time without propagating the system. It therefore doesn't include e.g.
* Verlet list updates.
* @param system The system to tune.
* @param int_steps Number of integration steps.
Expand Down
4 changes: 2 additions & 2 deletions src/core/unit_tests/EspressoSystemStandAlone_test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -347,7 +347,7 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory) {
reset_particle_positions();

// recalculate forces without propagating the system
integrate(system, 0, INTEG_REUSE_FORCES_CONDITIONALLY);
system.integrate(0, INTEG_REUSE_FORCES_CONDITIONALLY);

// particles are arranged in a right triangle
auto const p1_opt = copy_particle_to_head_node(comm, pid1);
Expand Down Expand Up @@ -390,7 +390,7 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory) {
expected[pid] = p.pos();
}
}
integrate(system, 1, INTEG_REUSE_FORCES_CONDITIONALLY);
system.integrate(1, INTEG_REUSE_FORCES_CONDITIONALLY);
for (auto pid : pids) {
auto const p_opt = copy_particle_to_head_node(comm, pid);
if (rank == 0) {
Expand Down
Loading

0 comments on commit d05263f

Please sign in to comment.