From bde47e88dbcb6cabedf2171b311315ee06327ceb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 13 Sep 2022 18:07:20 +0200 Subject: [PATCH 1/3] core: Remove MPI dependency from OIF code --- src/core/forces.cpp | 8 +++-- .../object-in-fluid/oif_global_forces.cpp | 33 ++++--------------- .../object-in-fluid/oif_global_forces.hpp | 1 - 3 files changed, 13 insertions(+), 29 deletions(-) diff --git a/src/core/forces.cpp b/src/core/forces.cpp index dfbb17019ed..2c4b2ef85ee 100644 --- a/src/core/forces.cpp +++ b/src/core/forces.cpp @@ -211,9 +211,13 @@ void force_calc(CellStructure &cell_structure, double time_step, double kT) { // There are two global quantities that need to be evaluated: // object's surface and object's volume. for (int i = 0; i < max_oif_objects; i++) { - auto const area_volume = calc_oif_global(i, cell_structure); - if (fabs(area_volume[0]) < 1e-100 && fabs(area_volume[1]) < 1e-100) + auto const area_volume = boost::mpi::all_reduce( + comm_cart, calc_oif_global(i, cell_structure), std::plus<>()); + auto const oif_part_area = std::abs(area_volume[0]); + auto const oif_part_vol = std::abs(area_volume[1]); + if (oif_part_area < 1e-100 and oif_part_vol < 1e-100) { break; + } add_oif_global_forces(area_volume, i, cell_structure); } } diff --git a/src/core/object-in-fluid/oif_global_forces.cpp b/src/core/object-in-fluid/oif_global_forces.cpp index ef09ccedfc9..6cdc35d3a13 100644 --- a/src/core/object-in-fluid/oif_global_forces.cpp +++ b/src/core/object-in-fluid/oif_global_forces.cpp @@ -22,7 +22,6 @@ #include "BoxGeometry.hpp" #include "Particle.hpp" #include "cell_system/CellStructure.hpp" -#include "communication.hpp" #include "grid.hpp" #include "bonded_interactions/bonded_interaction_data.hpp" @@ -32,12 +31,7 @@ #include #include -#include - -#include - -using Utils::area_triangle; -using Utils::get_n_triangle; +int max_oif_objects = 0; Utils::Vector2d calc_oif_global(int molType, CellStructure &cs) { // first-fold-then-the-same approach @@ -60,20 +54,19 @@ Utils::Vector2d calc_oif_global(int molType, CellStructure &cs) { auto const p33 = p11 + box_geo.get_mi_vector(partners[1]->pos(), p11); // unfolded positions correct - auto const VOL_A = area_triangle(p11, p22, p33); + auto const VOL_A = Utils::area_triangle(p11, p22, p33); partArea += VOL_A; - auto const VOL_norm = get_n_triangle(p11, p22, p33); + auto const VOL_norm = Utils::get_n_triangle(p11, p22, p33); auto const VOL_dn = VOL_norm.norm(); auto const VOL_hz = 1.0 / 3.0 * (p11[2] + p22[2] + p33[2]); - VOL_partVol += VOL_A * -1 * VOL_norm[2] / VOL_dn * VOL_hz; + VOL_partVol -= VOL_A * VOL_norm[2] / VOL_dn * VOL_hz; } return false; }); - auto const area_volume_local = Utils::Vector2d{{partArea, VOL_partVol}}; - return boost::mpi::all_reduce(comm_cart, area_volume_local, std::plus<>()); + return {{partArea, VOL_partVol}}; } void add_oif_global_forces(Utils::Vector2d const &area_volume, int molType, @@ -96,8 +89,8 @@ void add_oif_global_forces(Utils::Vector2d const &area_volume, int molType, // unfolded positions correct // starting code from volume force - auto const VOL_norm = get_n_triangle(p11, p22, p33).normalize(); - auto const VOL_A = area_triangle(p11, p22, p33); + auto const VOL_norm = Utils::get_n_triangle(p11, p22, p33).normalize(); + auto const VOL_A = Utils::area_triangle(p11, p22, p33); auto const VOL_vv = (VOL_volume - iaparams->V0) / iaparams->V0; auto const VOL_force = @@ -127,15 +120,3 @@ void add_oif_global_forces(Utils::Vector2d const &area_volume, int molType, return false; }); } - -int max_oif_objects = 0; - -void mpi_set_max_oif_objects_local(int max_oif_objects) { - ::max_oif_objects = max_oif_objects; -} - -REGISTER_CALLBACK(mpi_set_max_oif_objects_local) - -void mpi_set_max_oif_objects(int max_oif_objects) { - mpi_call_all(mpi_set_max_oif_objects_local, max_oif_objects); -} diff --git a/src/core/object-in-fluid/oif_global_forces.hpp b/src/core/object-in-fluid/oif_global_forces.hpp index 685c457b4be..a0bd8924350 100644 --- a/src/core/object-in-fluid/oif_global_forces.hpp +++ b/src/core/object-in-fluid/oif_global_forces.hpp @@ -43,5 +43,4 @@ void add_oif_global_forces(Utils::Vector2d const &area_volume, int molType, extern int max_oif_objects; -void mpi_set_max_oif_objects(int max_oif_objects); #endif From 67eeeadf0341d7c8fb895809f2832ca84b2144a7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 14 Sep 2022 22:58:01 +0200 Subject: [PATCH 2/3] script_interface: Make MPI communicator public The MPI communicator is now accessible through the object context. --- src/script_interface/Context.hpp | 7 ++++ src/script_interface/GlobalContext.hpp | 5 ++- src/script_interface/LocalContext.hpp | 3 ++ src/script_interface/ObjectList.hpp | 3 +- src/script_interface/ObjectMap.hpp | 3 +- .../ParallelExceptionHandler.cpp | 2 - .../cell_system/CellSystem.cpp | 10 ++--- .../object_container_mpi_guard.cpp | 7 ++-- .../object_container_mpi_guard.hpp | 4 +- .../tests/ObjectHandle_test.cpp | 14 ++++++- .../tests/serialization_mpi_guard_test.cpp | 42 ++++++++++++------- 11 files changed, 69 insertions(+), 31 deletions(-) diff --git a/src/script_interface/Context.hpp b/src/script_interface/Context.hpp index 5aace50f5a3..6a979517c1c 100644 --- a/src/script_interface/Context.hpp +++ b/src/script_interface/Context.hpp @@ -35,6 +35,12 @@ #include #include +namespace boost { +namespace mpi { +class communicator; +} // namespace mpi +} // namespace boost + namespace ScriptInterface { /** * @brief Context of an object handle. @@ -99,6 +105,7 @@ class Context : public std::enable_shared_from_this { virtual bool is_head_node() const = 0; virtual void parallel_try_catch(std::function const &cb) const = 0; + virtual boost::mpi::communicator const &get_comm() const = 0; virtual ~Context() = default; }; diff --git a/src/script_interface/GlobalContext.hpp b/src/script_interface/GlobalContext.hpp index d4de6ebbc02..03f2d0dc329 100644 --- a/src/script_interface/GlobalContext.hpp +++ b/src/script_interface/GlobalContext.hpp @@ -37,6 +37,7 @@ #include +#include #include #include @@ -71,6 +72,7 @@ class GlobalContext : public Context { std::shared_ptr m_node_local_context; bool m_is_head_node; + boost::mpi::communicator const &m_comm; ParallelExceptionHandler m_parallel_exception_handler; @@ -89,7 +91,7 @@ class GlobalContext : public Context { GlobalContext(Communication::MpiCallbacks &callbacks, std::shared_ptr node_local_context) : m_local_objects(), m_node_local_context(std::move(node_local_context)), - m_is_head_node(callbacks.comm().rank() == 0), + m_is_head_node(callbacks.comm().rank() == 0), m_comm(callbacks.comm()), // NOLINTNEXTLINE(bugprone-throw-keyword-missing) m_parallel_exception_handler(callbacks.comm()), cb_make_handle(&callbacks, @@ -170,6 +172,7 @@ class GlobalContext : public Context { void parallel_try_catch(std::function const &cb) const override { m_parallel_exception_handler.parallel_try_catch(cb); } + boost::mpi::communicator const &get_comm() const override { return m_comm; } }; } // namespace ScriptInterface diff --git a/src/script_interface/LocalContext.hpp b/src/script_interface/LocalContext.hpp index dcb81be5edc..11a79a7a2fb 100644 --- a/src/script_interface/LocalContext.hpp +++ b/src/script_interface/LocalContext.hpp @@ -43,12 +43,14 @@ namespace ScriptInterface { class LocalContext : public Context { Utils::Factory m_factory; bool m_is_head_node; + boost::mpi::communicator const &m_comm; ParallelExceptionHandler m_parallel_exception_handler; public: LocalContext(Utils::Factory factory, boost::mpi::communicator const &comm) : m_factory(std::move(factory)), m_is_head_node(comm.rank() == 0), + m_comm(comm), // NOLINTNEXTLINE(bugprone-throw-keyword-missing) m_parallel_exception_handler(comm) {} @@ -79,6 +81,7 @@ class LocalContext : public Context { void parallel_try_catch(std::function const &cb) const override { m_parallel_exception_handler.parallel_try_catch(cb); } + boost::mpi::communicator const &get_comm() const override { return m_comm; } }; } // namespace ScriptInterface diff --git a/src/script_interface/ObjectList.hpp b/src/script_interface/ObjectList.hpp index e43506afd15..3bcad1caa57 100644 --- a/src/script_interface/ObjectList.hpp +++ b/src/script_interface/ObjectList.hpp @@ -133,7 +133,8 @@ class ObjectList : public BaseType { private: std::string get_internal_state() const override { - object_container_mpi_guard(BaseType::name(), m_elements.size()); + object_container_mpi_guard(BaseType::name(), m_elements.size(), + BaseType::context()->get_comm().size()); std::vector object_states(m_elements.size()); diff --git a/src/script_interface/ObjectMap.hpp b/src/script_interface/ObjectMap.hpp index 2556a4e198e..cfd222f6241 100644 --- a/src/script_interface/ObjectMap.hpp +++ b/src/script_interface/ObjectMap.hpp @@ -164,7 +164,8 @@ class ObjectMap : public BaseType { private: std::string get_internal_state() const override { - object_container_mpi_guard(BaseType::name(), m_elements.size()); + object_container_mpi_guard(BaseType::name(), m_elements.size(), + BaseType::context()->get_comm().size()); using packed_type = std::pair; std::vector object_states(m_elements.size()); diff --git a/src/script_interface/ParallelExceptionHandler.cpp b/src/script_interface/ParallelExceptionHandler.cpp index e70feff9da4..b694cccea6f 100644 --- a/src/script_interface/ParallelExceptionHandler.cpp +++ b/src/script_interface/ParallelExceptionHandler.cpp @@ -21,8 +21,6 @@ #include "Exception.hpp" -#include "core/MpiCallbacks.hpp" -#include "core/communication.hpp" #include "core/error_handling/RuntimeError.hpp" #include "core/errorhandling.hpp" diff --git a/src/script_interface/cell_system/CellSystem.cpp b/src/script_interface/cell_system/CellSystem.cpp index 553af6321cb..688da7bc598 100644 --- a/src/script_interface/cell_system/CellSystem.cpp +++ b/src/script_interface/cell_system/CellSystem.cpp @@ -25,7 +25,6 @@ #include "core/cell_system/HybridDecomposition.hpp" #include "core/cell_system/RegularDecomposition.hpp" #include "core/cells.hpp" -#include "core/communication.hpp" #include "core/event.hpp" #include "core/grid.hpp" #include "core/integrate.hpp" @@ -159,7 +158,7 @@ Variant CellSystem::do_call_method(std::string const &name, {"n_square", hd.count_particles_in_n_square()}}}; } state["verlet_reuse"] = get_verlet_reuse(); - state["n_nodes"] = ::n_nodes; + state["n_nodes"] = context()->get_comm().size(); return state; } if (name == "get_pairs") { @@ -187,7 +186,7 @@ Variant CellSystem::do_call_method(std::string const &name, } if (name == "get_neighbors") { std::vector> neighbors_global; - context()->parallel_try_catch([&neighbors_global, ¶ms]() { + context()->parallel_try_catch([this, &neighbors_global, ¶ms]() { auto const dist = get_value(params, "distance"); auto const pid = get_value(params, "pid"); auto const ret = mpi_get_short_range_neighbors_local(pid, dist, true); @@ -195,7 +194,8 @@ Variant CellSystem::do_call_method(std::string const &name, if (ret) { neighbors_local = *ret; } - boost::mpi::gather(::comm_cart, neighbors_local, neighbors_global, 0); + boost::mpi::gather(context()->get_comm(), neighbors_local, + neighbors_global, 0); }); std::vector neighbors; for (auto const &neighbors_local : neighbors_global) { @@ -235,7 +235,7 @@ std::vector CellSystem::mpi_resort_particles(bool global_flag) const { } auto const size = static_cast(::cell_structure.local_particles().size()); std::vector n_part_per_node; - boost::mpi::gather(::comm_cart, size, n_part_per_node, 0); + boost::mpi::gather(context()->get_comm(), size, n_part_per_node, 0); return n_part_per_node; } diff --git a/src/script_interface/object_container_mpi_guard.cpp b/src/script_interface/object_container_mpi_guard.cpp index 4f1e680e280..da62287b3f2 100644 --- a/src/script_interface/object_container_mpi_guard.cpp +++ b/src/script_interface/object_container_mpi_guard.cpp @@ -19,8 +19,7 @@ #include "script_interface/object_container_mpi_guard.hpp" -#include "core/communication.hpp" - +#include #include #include @@ -28,8 +27,8 @@ #include void object_container_mpi_guard(boost::string_ref const &name, - std::size_t n_elements) { - if (comm_cart.size() > 1 and n_elements) { + std::size_t n_elements, int world_size) { + if (world_size > 1 and n_elements) { std::stringstream error_msg; error_msg << "Non-empty object containers do not support checkpointing in " << "MPI environments. Container " << name << " contains " diff --git a/src/script_interface/object_container_mpi_guard.hpp b/src/script_interface/object_container_mpi_guard.hpp index 90eededb143..bb9723c5e6f 100644 --- a/src/script_interface/object_container_mpi_guard.hpp +++ b/src/script_interface/object_container_mpi_guard.hpp @@ -19,6 +19,7 @@ #ifndef ESPRESSO_OBJECT_CONTAINER_MPI_GUARD_HPP #define ESPRESSO_OBJECT_CONTAINER_MPI_GUARD_HPP +#include #include #include @@ -39,8 +40,9 @@ * * @param name Name of the object container * @param n_elements Number of elements in the container + * @param world_size MPI world size */ void object_container_mpi_guard(boost::string_ref const &name, - std::size_t n_elements); + std::size_t n_elements, int world_size); #endif diff --git a/src/script_interface/tests/ObjectHandle_test.cpp b/src/script_interface/tests/ObjectHandle_test.cpp index ed24a0ec4a2..5ea001e084a 100644 --- a/src/script_interface/tests/ObjectHandle_test.cpp +++ b/src/script_interface/tests/ObjectHandle_test.cpp @@ -34,6 +34,7 @@ #include #include +#include #include #include #include @@ -165,6 +166,12 @@ BOOST_AUTO_TEST_CASE(do_call_method_) { MockCall::CallMethod{&name, ¶ms})); } +namespace boost { +namespace mpi { +class communicator {}; +} // namespace mpi +} // namespace boost + namespace Testing { /** * Logging mock for Context. @@ -195,6 +202,10 @@ struct LogContext : public Context { bool is_head_node() const override { return true; } void parallel_try_catch(std::function const &) const override {} + boost::mpi::communicator const &get_comm() const override { return m_comm; } + +private: + boost::mpi::communicator m_comm; }; } // namespace Testing @@ -249,5 +260,6 @@ BOOST_AUTO_TEST_CASE(interface_) { auto o = log_ctx->make_shared({}, {}); BOOST_CHECK(log_ctx->is_head_node()); BOOST_CHECK_EQUAL(log_ctx->name(o.get()), "Dummy"); - static_cast(log_ctx->parallel_try_catch([]() {})); + log_ctx->parallel_try_catch([]() {}); + std::ignore = log_ctx->get_comm(); } diff --git a/src/script_interface/tests/serialization_mpi_guard_test.cpp b/src/script_interface/tests/serialization_mpi_guard_test.cpp index 88fea91e851..72749046f97 100644 --- a/src/script_interface/tests/serialization_mpi_guard_test.cpp +++ b/src/script_interface/tests/serialization_mpi_guard_test.cpp @@ -22,6 +22,7 @@ #define BOOST_TEST_DYN_LINK #include +#include "script_interface/GlobalContext.hpp" #include "script_interface/ObjectList.hpp" #include @@ -51,24 +52,35 @@ struct ObjectContainer : ScriptInterface::ObjectList { BOOST_AUTO_TEST_CASE(parallel_exception) { boost::mpi::communicator world; - auto const obj_ptr = std::make_shared(); - auto const predicate = [](std::exception const &ex) { - std::string message = - "Non-empty object containers do not support checkpointing " - "in MPI environments. Container contains 1 elements."; - return ex.what() == message; - }; + Utils::Factory factory; + factory.register_new("ObjectContainer"); + Communication::MpiCallbacks cb{world}; + auto ctx = std::make_shared( + cb, std::make_shared(factory, world)); - Testing::ObjectContainer list; - BOOST_CHECK_NO_THROW(list.serialize()); + if (world.rank() == 0) { + auto const obj_ptr = std::make_shared(); + auto const predicate = [](std::exception const &ex) { + std::string message = + "Non-empty object containers do not support checkpointing in MPI " + "environments. Container ObjectContainer contains 1 elements."; + return ex.what() == message; + }; - list.add(obj_ptr); - if (world.size() > 1) { - BOOST_CHECK_EXCEPTION(list.serialize(), std::runtime_error, predicate); - } + auto list_so = ctx->make_shared("ObjectContainer", {}); + auto &list = dynamic_cast(*list_so); + BOOST_CHECK_NO_THROW(list.serialize()); + + list.add(obj_ptr); + if (world.size() > 1) { + BOOST_CHECK_EXCEPTION(list.serialize(), std::runtime_error, predicate); + } - list.remove(obj_ptr); - BOOST_CHECK_NO_THROW(list.serialize()); + list.remove(obj_ptr); + BOOST_CHECK_NO_THROW(list.serialize()); + } else { + cb.loop(); + } } int main(int argc, char **argv) { From 59b5bea019a956b34e2f14b863596c5bd18cfbeb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 16 Sep 2022 14:51:17 +0200 Subject: [PATCH 3/3] script_interface: Make main analysis methods parallel --- src/core/analysis/statistics.cpp | 29 ++++--------- src/core/energy.cpp | 30 ++++--------- src/core/energy.hpp | 6 +-- src/core/grid_based_algorithms/lb.cpp | 13 +++--- src/core/grid_based_algorithms/lb.hpp | 7 ++-- .../grid_based_algorithms/lb_interface.cpp | 17 +++----- src/core/observables/EnergyObservable.hpp | 2 +- src/core/observables/PressureObservable.hpp | 2 +- src/core/observables/PressureTensor.hpp | 2 +- src/core/pressure.cpp | 15 +++---- src/core/pressure.hpp | 2 +- .../reaction_methods/ReactionAlgorithm.cpp | 8 ++-- src/core/reaction_methods/WidomInsertion.cpp | 4 +- .../EspressoSystemStandAlone_test.cpp | 14 ++++--- src/python/espressomd/analyze.py | 4 +- src/script_interface/analysis/Analysis.cpp | 26 +++++++----- src/script_interface/communication.hpp | 42 +++++++++++++++++++ 17 files changed, 123 insertions(+), 100 deletions(-) create mode 100644 src/script_interface/communication.hpp diff --git a/src/core/analysis/statistics.cpp b/src/core/analysis/statistics.cpp index 3db0b719511..9ad044ca7f7 100644 --- a/src/core/analysis/statistics.cpp +++ b/src/core/analysis/statistics.cpp @@ -73,32 +73,21 @@ double mindist(PartCfg &partCfg, std::vector const &set1, return std::sqrt(mindist_sq); } -static Utils::Vector3d mpi_particle_momentum_local() { - auto const particles = cell_structure.local_particles(); - auto const momentum = - std::accumulate(particles.begin(), particles.end(), Utils::Vector3d{}, - [](Utils::Vector3d &m, Particle const &p) { - return m + p.mass() * p.v(); - }); - - return momentum; -} - -REGISTER_CALLBACK_REDUCTION(mpi_particle_momentum_local, - std::plus()) - Utils::Vector3d calc_linear_momentum(bool include_particles, bool include_lbfluid) { - Utils::Vector3d linear_momentum{}; + Utils::Vector3d momentum{}; if (include_particles) { - linear_momentum += - mpi_call(::Communication::Result::reduction, - std::plus(), mpi_particle_momentum_local); + auto const particles = cell_structure.local_particles(); + momentum = + std::accumulate(particles.begin(), particles.end(), Utils::Vector3d{}, + [](Utils::Vector3d &m, Particle const &p) { + return m + p.mass() * p.v(); + }); } if (include_lbfluid) { - linear_momentum += lb_lbfluid_calc_fluid_momentum(); + momentum += lb_lbfluid_calc_fluid_momentum(); } - return linear_momentum; + return momentum; } Utils::Vector3d center_of_mass(PartCfg &partCfg, int p_type) { diff --git a/src/core/energy.cpp b/src/core/energy.cpp index aceabcfc33a..1bec058ccb1 100644 --- a/src/core/energy.cpp +++ b/src/core/energy.cpp @@ -43,7 +43,7 @@ #include -static std::shared_ptr calculate_energy_local() { +std::shared_ptr calculate_energy() { auto obs_energy_ptr = std::make_shared(1); @@ -117,23 +117,19 @@ static std::shared_ptr calculate_energy_local() { // NOLINTNEXTLINE(clang-analyzer-cplusplus.NewDeleteLeaks) } -REGISTER_CALLBACK_MAIN_RANK(calculate_energy_local) +REGISTER_CALLBACK_MAIN_RANK(calculate_energy) -std::shared_ptr calculate_energy() { - return mpi_call(Communication::Result::main_rank, calculate_energy_local); +double mpi_calculate_potential_energy() { + auto const obs = mpi_call(Communication::Result::main_rank, calculate_energy); + return obs->accumulate(-obs->kinetic[0]); } -double calculate_current_potential_energy_of_system() { - auto const obs_energy = calculate_energy(); - return obs_energy->accumulate(-obs_energy->kinetic[0]); +double mpi_observable_compute_energy() { + auto const obs = mpi_call(Communication::Result::main_rank, calculate_energy); + return obs->accumulate(0); } -double observable_compute_energy() { - auto const obs_energy = calculate_energy(); - return obs_energy->accumulate(0); -} - -double particle_short_range_energy_contribution_local(int pid) { +double particle_short_range_energy_contribution(int pid) { double ret = 0.0; if (cell_structure.get_resort_particles()) { @@ -158,11 +154,3 @@ double particle_short_range_energy_contribution_local(int pid) { } return ret; } - -REGISTER_CALLBACK_REDUCTION(particle_short_range_energy_contribution_local, - std::plus()) - -double particle_short_range_energy_contribution(int pid) { - return mpi_call(Communication::Result::reduction, std::plus(), - particle_short_range_energy_contribution_local, pid); -} diff --git a/src/core/energy.hpp b/src/core/energy.hpp index 63380c6c200..91f5e4176ef 100644 --- a/src/core/energy.hpp +++ b/src/core/energy.hpp @@ -34,16 +34,16 @@ std::shared_ptr calculate_energy(); /** Calculate the total energy of the system. */ -double calculate_current_potential_energy_of_system(); +double mpi_calculate_potential_energy(); /** Helper function for @ref Observables::Energy. */ -double observable_compute_energy(); +double mpi_observable_compute_energy(); /** * @brief Compute short-range energy of a particle. * * Iterates through particles inside cell and neighboring cells and computes - * energy contribution for a specificparticle. + * energy contribution for a specific particle. * * @param pid Particle id * @return Non-bonded energy of the particle. diff --git a/src/core/grid_based_algorithms/lb.cpp b/src/core/grid_based_algorithms/lb.cpp index 6de37da6819..b9c60828073 100644 --- a/src/core/grid_based_algorithms/lb.cpp +++ b/src/core/grid_based_algorithms/lb.cpp @@ -1316,15 +1316,15 @@ Utils::Vector3d lb_calc_local_momentum_density(Lattice::index_t index, } /** Calculate momentum of the LB fluid. - * @param[out] result Fluid momentum in MD units * @param[in] lb_parameters LB parameters * @param[in] lb_fields Hydrodynamic fields of the fluid * @param[in] lb_lattice The underlying lattice */ -void lb_calc_fluid_momentum(double *result, const LB_Parameters &lb_parameters, - const std::vector &lb_fields, - const Lattice &lb_lattice) { - Utils::Vector3d momentum_density{}, momentum{}; +Utils::Vector3d +mpi_lb_calc_fluid_momentum_local(LB_Parameters const &lb_parameters, + std::vector const &lb_fields, + Lattice const &lb_lattice) { + Utils::Vector3d momentum_density{}, momentum{}, result{}; for (int x = 1; x <= lb_lattice.grid[0]; x++) { for (int y = 1; y <= lb_lattice.grid[1]; y++) { @@ -1338,7 +1338,8 @@ void lb_calc_fluid_momentum(double *result, const LB_Parameters &lb_parameters, } momentum *= lb_parameters.agrid / lb_parameters.tau; - boost::mpi::reduce(comm_cart, momentum.data(), 3, result, std::plus<>(), 0); + boost::mpi::reduce(::comm_cart, momentum, result, std::plus<>(), 0); + return result; } void lb_collect_boundary_forces(double *result) { diff --git a/src/core/grid_based_algorithms/lb.hpp b/src/core/grid_based_algorithms/lb.hpp index dac4508134f..f7dc8eae44d 100644 --- a/src/core/grid_based_algorithms/lb.hpp +++ b/src/core/grid_based_algorithms/lb.hpp @@ -250,9 +250,10 @@ void lb_bounce_back(LB_Fluid &lbfluid, const LB_Parameters &lb_parameters, #endif /* LB_BOUNDARIES */ -void lb_calc_fluid_momentum(double *result, const LB_Parameters &lb_parameters, - const std::vector &lb_fields, - const Lattice &lb_lattice); +Utils::Vector3d +mpi_lb_calc_fluid_momentum_local(LB_Parameters const &lb_parameters, + std::vector const &lb_fields, + Lattice const &lb_lattice); void lb_collect_boundary_forces(double *result); void lb_initialize_fields(std::vector &fields, LB_Parameters const &lb_parameters, diff --git a/src/core/grid_based_algorithms/lb_interface.cpp b/src/core/grid_based_algorithms/lb_interface.cpp index 3630d31d47c..662a1c59b93 100644 --- a/src/core/grid_based_algorithms/lb_interface.cpp +++ b/src/core/grid_based_algorithms/lb_interface.cpp @@ -1099,23 +1099,18 @@ void lb_lbnode_set_pop(const Utils::Vector3i &ind, } } -static void mpi_lb_lbfluid_calc_fluid_momentum_local() { - lb_calc_fluid_momentum(nullptr, lbpar, lbfields, lblattice); -} - -REGISTER_CALLBACK(mpi_lb_lbfluid_calc_fluid_momentum_local) - Utils::Vector3d lb_lbfluid_calc_fluid_momentum() { - Utils::Vector3d fluid_momentum{}; + Utils::Vector3d momentum{}; if (lattice_switch == ActiveLB::GPU) { #ifdef CUDA - lb_calc_fluid_momentum_GPU(fluid_momentum.data()); + if (::comm_cart.rank() == 0) { + lb_calc_fluid_momentum_GPU(momentum.data()); + } #endif } else if (lattice_switch == ActiveLB::CPU) { - mpi_call(mpi_lb_lbfluid_calc_fluid_momentum_local); - lb_calc_fluid_momentum(fluid_momentum.data(), lbpar, lbfields, lblattice); + momentum = mpi_lb_calc_fluid_momentum_local(lbpar, lbfields, lblattice); } - return fluid_momentum; + return momentum; } const Utils::Vector3d diff --git a/src/core/observables/EnergyObservable.hpp b/src/core/observables/EnergyObservable.hpp index 8f0a906971c..40132153c01 100644 --- a/src/core/observables/EnergyObservable.hpp +++ b/src/core/observables/EnergyObservable.hpp @@ -32,7 +32,7 @@ class Energy : public Observable { std::vector shape() const override { return {1}; } std::vector operator()() const override { std::vector res{1}; - res[0] = observable_compute_energy(); + res[0] = mpi_observable_compute_energy(); return res; } }; diff --git a/src/core/observables/PressureObservable.hpp b/src/core/observables/PressureObservable.hpp index f9d50f09b45..1f9ec9cc8eb 100644 --- a/src/core/observables/PressureObservable.hpp +++ b/src/core/observables/PressureObservable.hpp @@ -30,7 +30,7 @@ class Pressure : public Observable { public: std::vector shape() const override { return {1}; } std::vector operator()() const override { - auto const ptensor = observable_compute_pressure_tensor(); + auto const ptensor = mpi_observable_compute_pressure_tensor(); std::vector res{1}; res[0] = (ptensor[0] + ptensor[4] + ptensor[8]) / 3.; return res; diff --git a/src/core/observables/PressureTensor.hpp b/src/core/observables/PressureTensor.hpp index 485dd4b7a95..d9fc5eabcde 100644 --- a/src/core/observables/PressureTensor.hpp +++ b/src/core/observables/PressureTensor.hpp @@ -30,7 +30,7 @@ class PressureTensor : public Observable { public: std::vector shape() const override { return {3, 3}; } std::vector operator()() const override { - return observable_compute_pressure_tensor().as_vector(); + return mpi_observable_compute_pressure_tensor().as_vector(); } }; diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp index 3fa348ea420..5b770bf8ae1 100644 --- a/src/core/pressure.cpp +++ b/src/core/pressure.cpp @@ -52,7 +52,7 @@ #include #include -static std::shared_ptr calculate_pressure_local() { +std::shared_ptr calculate_pressure() { auto obs_pressure_ptr = std::make_shared(9); @@ -125,17 +125,14 @@ static std::shared_ptr calculate_pressure_local() { return obs_pressure_ptr; } -REGISTER_CALLBACK_MAIN_RANK(calculate_pressure_local) - -std::shared_ptr calculate_pressure() { - return mpi_call(Communication::Result::main_rank, calculate_pressure_local); -} +REGISTER_CALLBACK_MAIN_RANK(calculate_pressure) -Utils::Vector9d observable_compute_pressure_tensor() { - auto const obs_pressure = calculate_pressure(); +Utils::Vector9d mpi_observable_compute_pressure_tensor() { + auto const obs = + mpi_call(Communication::Result::main_rank, calculate_pressure); Utils::Vector9d pressure_tensor{}; for (std::size_t j = 0; j < 9; j++) { - pressure_tensor[j] = obs_pressure->accumulate(0, j); + pressure_tensor[j] = obs->accumulate(0, j); } return pressure_tensor; } diff --git a/src/core/pressure.hpp b/src/core/pressure.hpp index 9ab357126b0..7fd55d37841 100644 --- a/src/core/pressure.hpp +++ b/src/core/pressure.hpp @@ -35,6 +35,6 @@ std::shared_ptr calculate_pressure(); /** Helper function for @ref Observables::PressureTensor. */ -Utils::Vector9d observable_compute_pressure_tensor(); +Utils::Vector9d mpi_observable_compute_pressure_tensor(); #endif diff --git a/src/core/reaction_methods/ReactionAlgorithm.cpp b/src/core/reaction_methods/ReactionAlgorithm.cpp index f5feb49a82c..1f8330318be 100644 --- a/src/core/reaction_methods/ReactionAlgorithm.cpp +++ b/src/core/reaction_methods/ReactionAlgorithm.cpp @@ -54,7 +54,7 @@ int ReactionAlgorithm::do_reaction(int reaction_steps) { // calculate potential energy; only consider potential energy since we // assume that the kinetic part drops out in the process of calculating // ensemble averages (kinetic part may be separated and crossed out) - auto current_E_pot = calculate_current_potential_energy_of_system(); + auto current_E_pot = mpi_calculate_potential_energy(); for (int i = 0; i < reaction_steps; i++) { int reaction_id = i_random(static_cast(reactions.size())); generic_oneway_reaction(*reactions[reaction_id], current_E_pot); @@ -295,7 +295,7 @@ void ReactionAlgorithm::generic_oneway_reaction( auto const E_pot_new = (particle_inside_exclusion_range_touched) ? std::numeric_limits::max() - : calculate_current_potential_energy_of_system(); + : mpi_calculate_potential_energy(); auto const bf = calculate_acceptance_probability( current_reaction, E_pot_old, E_pot_new, old_particle_numbers); @@ -592,14 +592,14 @@ bool ReactionAlgorithm::do_global_mc_move_for_particles_of_type( return false; } - auto const E_pot_old = calculate_current_potential_energy_of_system(); + auto const E_pot_old = mpi_calculate_potential_energy(); auto const original_positions = generate_new_particle_positions( type, particle_number_of_type_to_be_changed); auto const E_pot_new = (particle_inside_exclusion_range_touched) ? std::numeric_limits::max() - : calculate_current_potential_energy_of_system(); + : mpi_calculate_potential_energy(); auto const beta = 1.0 / kT; diff --git a/src/core/reaction_methods/WidomInsertion.cpp b/src/core/reaction_methods/WidomInsertion.cpp index 783557ca75a..0945b735716 100644 --- a/src/core/reaction_methods/WidomInsertion.cpp +++ b/src/core/reaction_methods/WidomInsertion.cpp @@ -36,7 +36,7 @@ double WidomInsertion::calculate_particle_insertion_potential_energy( throw std::runtime_error("Trying to remove some non-existing particles " "from the system via the inverse Widom scheme."); - auto const E_pot_old = calculate_current_potential_energy_of_system(); + auto const E_pot_old = mpi_calculate_potential_energy(); // make reaction attempt std::vector p_ids_created_particles; @@ -46,7 +46,7 @@ double WidomInsertion::calculate_particle_insertion_potential_energy( hidden_particles_properties) = make_reaction_attempt(current_reaction); - auto const E_pot_new = calculate_current_potential_energy_of_system(); + auto const E_pot_new = mpi_calculate_potential_energy(); // reverse reaction attempt // reverse reaction // 1) delete created product particles diff --git a/src/core/unit_tests/EspressoSystemStandAlone_test.cpp b/src/core/unit_tests/EspressoSystemStandAlone_test.cpp index ddf19a283ed..0366ac86be9 100644 --- a/src/core/unit_tests/EspressoSystemStandAlone_test.cpp +++ b/src/core/unit_tests/EspressoSystemStandAlone_test.cpp @@ -104,6 +104,10 @@ static void mpi_remove_translational_motion() { mpi_call_all(mpi_remove_translational_motion_local); } +static auto mpi_calculate_energy() { + return mpi_call(Communication::Result::main_rank, calculate_energy); +} + #ifdef P3M static void mpi_set_tuned_p3m_local(double prefactor) { auto p3m = P3MParameters{false, @@ -201,9 +205,9 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory, set_particle_v(pid2, {static_cast(i), 0., 0.}); auto const &p = get_particle_data(pid2); auto const kinetic_energy = 0.5 * p.mass() * p.v().norm2(); - auto const obs_energy = calculate_energy(); + auto const obs_energy = mpi_calculate_energy(); BOOST_CHECK_CLOSE(obs_energy->kinetic[0], kinetic_energy, tol); - BOOST_CHECK_CLOSE(observable_compute_energy(), kinetic_energy, tol); + BOOST_CHECK_CLOSE(mpi_observable_compute_energy(), kinetic_energy, tol); } } @@ -232,7 +236,7 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory, auto const lj_energy = 4.0 * eps * (Utils::sqr(frac6) - frac6 + shift); // measure energies - auto const obs_energy = calculate_energy(); + auto const obs_energy = mpi_calculate_energy(); for (int i = 0; i < n_pairs; ++i) { auto const ref_inter = (i == lj_pair_ab) ? lj_energy : 0.; auto const ref_intra = (i == lj_pair_bb) ? lj_energy : 0.; @@ -259,7 +263,7 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory, add_particle_bond(pid2, std::vector{fene_bond_id, pid3}); // measure energies - auto const obs_energy = calculate_energy(); + auto const obs_energy = mpi_calculate_energy(); auto const none_energy = 0.0; auto const harm_energy = 0.5 * harm_bond.k * Utils::sqr(harm_bond.r - dist); auto const fene_energy = @@ -291,7 +295,7 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory, place_particle(pid2, pos2); auto const r = (pos2 - pos1).norm(); // check P3M energy - auto const obs_energy = calculate_energy(); + auto const obs_energy = mpi_calculate_energy(); // at very short distances, the real-space contribution to // the energy is much larger than the k-space contribution auto const energy_ref = -prefactor / r; diff --git a/src/python/espressomd/analyze.py b/src/python/espressomd/analyze.py index 69168368145..83f9c140e09 100644 --- a/src/python/espressomd/analyze.py +++ b/src/python/espressomd/analyze.py @@ -83,7 +83,7 @@ def acf_1d(signal, n_with_padding, n): @script_interface_register class _ObservableStat(ScriptInterfaceHelper): _so_name = "ScriptInterface::Analysis::ObservableStat" - _so_creation_policy = "LOCAL" + _so_creation_policy = "GLOBAL" def _generate_summary(self, obj, dim, calc_sp): """ @@ -346,7 +346,7 @@ class Analysis(ScriptInterfaceHelper): """ _so_name = "ScriptInterface::Analysis::Analysis" - _so_creation_policy = "LOCAL" + _so_creation_policy = "GLOBAL" _so_bind_methods = ( "linear_momentum", "center_of_mass", diff --git a/src/script_interface/analysis/Analysis.cpp b/src/script_interface/analysis/Analysis.cpp index 5ebebd25802..14c1b0bdf7d 100644 --- a/src/script_interface/analysis/Analysis.cpp +++ b/src/script_interface/analysis/Analysis.cpp @@ -28,6 +28,8 @@ #include "core/partCfg_global.hpp" #include "core/particle_node.hpp" +#include "script_interface/communication.hpp" + #include #include @@ -68,6 +70,20 @@ static void check_particle_type(int p_type) { Variant Analysis::do_call_method(std::string const &name, VariantMap const ¶meters) { + if (name == "linear_momentum") { + auto const local = calc_linear_momentum( + get_value_or(parameters, "include_particles", true), + get_value_or(parameters, "include_lbfluid", true)); + return mpi_reduce_sum(context()->get_comm(), local).as_vector(); + } + if (name == "particle_energy") { + auto const pid = get_value(parameters, "pid"); + auto const local = particle_short_range_energy_contribution(pid); + return mpi_reduce_sum(context()->get_comm(), local); + } + if (not context()->is_head_node()) { + return {}; + } if (name == "min_dist") { auto const p_types1 = get_value>(parameters, "p_types1"); auto const p_types2 = get_value>(parameters, "p_types2"); @@ -89,12 +105,6 @@ Variant Analysis::do_call_method(std::string const &name, auto const result = angular_momentum(partCfg(), p_type); return result.as_vector(); } - if (name == "linear_momentum") { - auto const result = calc_linear_momentum( - get_value_or(parameters, "include_particles", true), - get_value_or(parameters, "include_lbfluid", true)); - return result.as_vector(); - } if (name == "nbhood") { auto const pos = get_value(parameters, "pos"); auto const radius = get_value(parameters, "r_catch"); @@ -107,10 +117,6 @@ Variant Analysis::do_call_method(std::string const &name, return result.as_vector(); } #endif // DPD - if (name == "particle_energy") { - auto const pid = get_value(parameters, "pid"); - return particle_short_range_energy_contribution(pid); - } if (name == "calc_re") { auto const chain_start = get_value(parameters, "chain_start"); auto const chain_length = get_value(parameters, "chain_length"); diff --git a/src/script_interface/communication.hpp b/src/script_interface/communication.hpp new file mode 100644 index 00000000000..d66dbeef1b3 --- /dev/null +++ b/src/script_interface/communication.hpp @@ -0,0 +1,42 @@ +/* + * Copyright (C) 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 . + */ +#ifndef ESPRESSO_SRC_SCRIPT_INTERFACE_COMMUNICATION_HPP +#define ESPRESSO_SRC_SCRIPT_INTERFACE_COMMUNICATION_HPP + +#include +#include + +#include + +namespace ScriptInterface { + +/** + * @brief Reduce object by sum on the head node. + * Worker nodes get a default-constructed object. + */ +template +T mpi_reduce_sum(boost::mpi::communicator const &comm, T const &result) { + T out{}; + boost::mpi::reduce(comm, result, out, std::plus{}, 0); + return out; +} + +} // namespace ScriptInterface + +#endif