From 4942938ced78260d671ee04b609146f623e9d9da Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Mon, 24 Jul 2023 20:01:19 +0200 Subject: [PATCH] Implement parallel observables Co-authored-by: Alexander Reinauer --- src/core/accumulators.cpp | 4 +- src/core/accumulators.hpp | 2 +- src/core/accumulators/AccumulatorBase.hpp | 7 +- src/core/accumulators/Correlator.cpp | 24 +- src/core/accumulators/Correlator.hpp | 4 +- .../accumulators/MeanVarianceCalculator.cpp | 8 +- .../accumulators/MeanVarianceCalculator.hpp | 2 +- src/core/accumulators/TimeSeries.cpp | 8 +- src/core/accumulators/TimeSeries.hpp | 2 +- src/core/dpd.cpp | 16 +- src/core/dpd.hpp | 7 +- src/core/energy.cpp | 17 - src/core/energy.hpp | 10 +- .../grid_based_algorithms/lb_interface.cpp | 58 +--- .../grid_based_algorithms/lb_interface.hpp | 25 +- src/core/integrate.cpp | 23 +- src/core/observables/BondAngles.hpp | 17 +- src/core/observables/BondDihedrals.hpp | 21 +- src/core/observables/CosPersistenceAngles.hpp | 16 +- .../observables/CylindricalDensityProfile.hpp | 35 +- .../CylindricalFluxDensityProfile.hpp | 38 ++- ...BFluxDensityProfileAtParticlePositions.cpp | 62 ++-- ...BFluxDensityProfileAtParticlePositions.hpp | 6 +- .../CylindricalLBProfileObservable.hpp | 1 - .../CylindricalLBVelocityProfile.cpp | 59 +++- .../CylindricalLBVelocityProfile.hpp | 3 +- ...alLBVelocityProfileAtParticlePositions.cpp | 61 ++-- ...alLBVelocityProfileAtParticlePositions.hpp | 6 +- .../CylindricalProfileObservable.hpp | 7 - .../CylindricalVelocityProfile.hpp | 44 ++- src/core/observables/DPDStress.hpp | 5 +- src/core/observables/DensityProfile.hpp | 32 +- src/core/observables/EnergyObservable.hpp | 7 +- src/core/observables/FluxDensityProfile.hpp | 39 ++- src/core/observables/ForceDensityProfile.hpp | 43 ++- .../observables/LBFluidPressureTensor.hpp | 5 +- src/core/observables/LBVelocityProfile.cpp | 58 ++-- src/core/observables/LBVelocityProfile.hpp | 3 +- src/core/observables/Observable.hpp | 5 +- src/core/observables/ParticleDipoleFields.hpp | 7 +- src/core/observables/ParticleDistances.hpp | 14 +- src/core/observables/ParticleTraits.hpp | 8 +- src/core/observables/PidObservable.cpp | 13 +- src/core/observables/PidObservable.hpp | 172 +++++++++- src/core/observables/PressureObservable.hpp | 13 +- src/core/observables/PressureTensor.hpp | 13 +- src/core/observables/RDF.cpp | 67 ++-- src/core/observables/RDF.hpp | 14 +- src/core/observables/TotalForce.hpp | 17 +- src/core/observables/fetch_particles.hpp | 42 +-- src/core/observables/utils_histogram.hpp | 61 ++++ src/core/pressure.cpp | 16 - src/core/pressure.hpp | 3 - .../EspressoSystemStandAlone_test.cpp | 47 ++- .../unit_tests/lb_particle_coupling_test.cpp | 3 +- .../particle_observables/algorithms.hpp | 30 +- .../particle_observables/observable.hpp | 1 + .../particle_observables/properties.hpp | 8 + src/particle_observables/tests/algorithms.cpp | 19 +- .../tests/observables.cpp | 6 +- src/python/espressomd/accumulators.py | 8 +- src/python/espressomd/observables.py | 2 +- .../accumulators/Correlator.hpp | 49 ++- .../accumulators/MeanVarianceCalculator.hpp | 6 +- .../accumulators/TimeSeries.hpp | 3 +- src/script_interface/analysis/Analysis.cpp | 12 +- src/script_interface/communication.hpp | 26 -- .../CylindricalLBProfileObservable.hpp | 29 +- .../CylindricalPidProfileObservable.hpp | 29 +- .../observables/LBProfileObservable.hpp | 18 +- .../observables/Observable.hpp | 6 +- .../observables/PidObservable.hpp | 6 +- .../observables/PidProfileObservable.hpp | 12 +- src/script_interface/observables/RDF.hpp | 9 +- .../particle_data/ParticleHandle.cpp | 4 +- .../tests/Accumulators_test.cpp | 314 +++++++++++------- src/script_interface/tests/CMakeLists.txt | 3 +- .../tests/GlobalContext_test.cpp | 8 +- src/script_interface/tests/reduction_test.cpp | 19 -- src/script_interface/walberla/EKReaction.hpp | 5 +- .../walberla/EKSpeciesNode.cpp | 14 +- src/script_interface/walberla/LBFluid.cpp | 4 +- src/script_interface/walberla/LBFluidNode.cpp | 23 +- .../include/utils/mpi/reduce_optional.hpp | 57 ++++ src/utils/tests/CMakeLists.txt | 2 + src/utils/tests/reduce_optional_test.cpp | 52 +++ testsuite/python/CMakeLists.txt | 2 +- testsuite/python/observable_profile.py | 7 +- 88 files changed, 1419 insertions(+), 674 deletions(-) create mode 100644 src/core/observables/utils_histogram.hpp create mode 100644 src/utils/include/utils/mpi/reduce_optional.hpp create mode 100644 src/utils/tests/reduce_optional_test.cpp diff --git a/src/core/accumulators.cpp b/src/core/accumulators.cpp index e54d91a3bf5..6f707afd0a3 100644 --- a/src/core/accumulators.cpp +++ b/src/core/accumulators.cpp @@ -39,12 +39,12 @@ struct AutoUpdateAccumulator { std::vector auto_update_accumulators; } // namespace -void auto_update(int steps) { +void auto_update(boost::mpi::communicator const &comm, int steps) { for (auto &acc : auto_update_accumulators) { assert(steps <= acc.frequency); acc.counter -= steps; if (acc.counter <= 0) { - acc.acc->update(); + acc.acc->update(comm); acc.counter = acc.frequency; } diff --git a/src/core/accumulators.hpp b/src/core/accumulators.hpp index 3082b7b9a2c..68a9da09a05 100644 --- a/src/core/accumulators.hpp +++ b/src/core/accumulators.hpp @@ -29,7 +29,7 @@ namespace Accumulators { * they need to be updated and if so does. * */ -void auto_update(int steps); +void auto_update(boost::mpi::communicator const &comm, int steps); int auto_update_next_update(); void auto_update_add(AccumulatorBase *); void auto_update_remove(AccumulatorBase *); diff --git a/src/core/accumulators/AccumulatorBase.hpp b/src/core/accumulators/AccumulatorBase.hpp index 807fc41841f..7bf17d4a799 100644 --- a/src/core/accumulators/AccumulatorBase.hpp +++ b/src/core/accumulators/AccumulatorBase.hpp @@ -22,6 +22,11 @@ #include #include +// Forward declaration +namespace boost::mpi { +class communicator; +} + namespace Accumulators { class AccumulatorBase { @@ -31,7 +36,7 @@ class AccumulatorBase { int &delta_N() { return m_delta_N; } - virtual void update() = 0; + virtual void update(boost::mpi::communicator const &comm) = 0; /** Dimensions needed to reshape the flat array returned by the accumulator */ virtual std::vector shape() const = 0; diff --git a/src/core/accumulators/Correlator.cpp b/src/core/accumulators/Correlator.cpp index c0a1349ae43..cbc1c453a07 100644 --- a/src/core/accumulators/Correlator.cpp +++ b/src/core/accumulators/Correlator.cpp @@ -310,11 +310,22 @@ void Correlator::initialize() { } } -void Correlator::update() { +void Correlator::update(boost::mpi::communicator const &comm) { if (finalized) { throw std::runtime_error( "No data can be added after finalize() was called."); } + + if (comm.rank() != 0) { + // worker nodes just need to update the observables and exit + A_obs->operator()(comm); + if (A_obs != B_obs) { + B_obs->operator()(comm); + } + + return; + } + // We must now go through the hierarchy and make sure there is space for the // new datapoint. For every hierarchy level we have to decide if it is // necessary to move something @@ -360,9 +371,9 @@ void Correlator::update() { newest[0] = (newest[0] + 1) % (m_tau_lin + 1); n_vals[0]++; - A[0][newest[0]] = A_obs->operator()(); + A[0][newest[0]] = A_obs->operator()(comm); if (A_obs != B_obs) { - B[0][newest[0]] = B_obs->operator()(); + B[0][newest[0]] = B_obs->operator()(comm); } else { B[0][newest[0]] = A[0][newest[0]]; } @@ -411,7 +422,7 @@ void Correlator::update() { } } -int Correlator::finalize() { +int Correlator::finalize(boost::mpi::communicator const &comm) { using index_type = decltype(result)::index; if (finalized) { throw std::runtime_error("Correlator::finalize() can only be called once."); @@ -423,6 +434,11 @@ int Correlator::finalize() { // mark the correlation as finalized finalized = true; + // worker nodes don't need to do anything + if (comm.rank() != 0) { + return 0; + } + for (int ll = 0; ll < m_hierarchy_depth - 1; ll++) { long vals_ll; // number of values remaining in the lowest level if (n_vals[ll] > m_tau_lin + 1) diff --git a/src/core/accumulators/Correlator.hpp b/src/core/accumulators/Correlator.hpp index 690c1a91c6b..42588719223 100644 --- a/src/core/accumulators/Correlator.hpp +++ b/src/core/accumulators/Correlator.hpp @@ -177,12 +177,12 @@ class Correlator : public AccumulatorBase { * the correlation estimate is updated. * TODO: Not all correlation estimates have to be updated. */ - void update() override; + void update(boost::mpi::communicator const &comm) override; /** At the end of data collection, go through the whole hierarchy and * correlate data left there. */ - int finalize(); + int finalize(boost::mpi::communicator const &comm); /** Return correlation result */ std::vector get_correlation(); diff --git a/src/core/accumulators/MeanVarianceCalculator.cpp b/src/core/accumulators/MeanVarianceCalculator.cpp index 89e04180574..06373928ed7 100644 --- a/src/core/accumulators/MeanVarianceCalculator.cpp +++ b/src/core/accumulators/MeanVarianceCalculator.cpp @@ -32,7 +32,13 @@ #include namespace Accumulators { -void MeanVarianceCalculator::update() { m_acc(m_obs->operator()()); } +void MeanVarianceCalculator::update(boost::mpi::communicator const &comm) { + if (comm.rank() == 0) { + m_acc(m_obs->operator()(comm)); + } else { + m_obs->operator()(comm); + } +} std::vector MeanVarianceCalculator::mean() { return m_acc.mean(); } diff --git a/src/core/accumulators/MeanVarianceCalculator.hpp b/src/core/accumulators/MeanVarianceCalculator.hpp index f4a23acdcd9..f2aa5ec6852 100644 --- a/src/core/accumulators/MeanVarianceCalculator.hpp +++ b/src/core/accumulators/MeanVarianceCalculator.hpp @@ -38,7 +38,7 @@ class MeanVarianceCalculator : public AccumulatorBase { int delta_N) : AccumulatorBase(delta_N), m_obs(obs), m_acc(obs->n_values()) {} - void update() override; + void update(boost::mpi::communicator const &comm) override; std::vector mean(); std::vector variance(); std::vector std_error(); diff --git a/src/core/accumulators/TimeSeries.cpp b/src/core/accumulators/TimeSeries.cpp index e30fd154a1a..206b6a0160d 100644 --- a/src/core/accumulators/TimeSeries.cpp +++ b/src/core/accumulators/TimeSeries.cpp @@ -28,7 +28,13 @@ #include namespace Accumulators { -void TimeSeries::update() { m_data.emplace_back(m_obs->operator()()); } +void TimeSeries::update(boost::mpi::communicator const &comm) { + if (comm.rank() == 0) { + m_data.emplace_back(m_obs->operator()(comm)); + } else { + m_obs->operator()(comm); + } +} std::string TimeSeries::get_internal_state() const { std::stringstream ss; diff --git a/src/core/accumulators/TimeSeries.hpp b/src/core/accumulators/TimeSeries.hpp index b5b5fd847dd..e7d78c49de6 100644 --- a/src/core/accumulators/TimeSeries.hpp +++ b/src/core/accumulators/TimeSeries.hpp @@ -43,7 +43,7 @@ class TimeSeries : public AccumulatorBase { TimeSeries(std::shared_ptr obs, int delta_N) : AccumulatorBase(delta_N), m_obs(std::move(obs)) {} - void update() override; + void update(boost::mpi::communicator const &comm) override; std::string get_internal_state() const; void set_internal_state(std::string const &); diff --git a/src/core/dpd.cpp b/src/core/dpd.cpp index 52a427b743e..369b513c19e 100644 --- a/src/core/dpd.cpp +++ b/src/core/dpd.cpp @@ -27,9 +27,7 @@ #include "dpd.hpp" -#include "MpiCallbacks.hpp" #include "cells.hpp" -#include "communication.hpp" #include "event.hpp" #include "grid.hpp" #include "integrate.hpp" @@ -44,6 +42,8 @@ #include #include +#include + #include #include #include @@ -151,7 +151,6 @@ static auto dpd_viscous_stress_local() { return stress; } -REGISTER_CALLBACK_REDUCTION(dpd_viscous_stress_local, std::plus<>()) /** * @brief Viscous stress tensor of the DPD interaction. @@ -168,12 +167,13 @@ REGISTER_CALLBACK_REDUCTION(dpd_viscous_stress_local, std::plus<>()) * * @return Stress tensor contribution. */ -Utils::Vector9d dpd_stress() { - auto const stress = mpi_call(Communication::Result::reduction, std::plus<>(), - dpd_viscous_stress_local); - auto const volume = box_geo.volume(); +Utils::Vector9d dpd_stress(boost::mpi::communicator const &comm) { + auto const local_stress = dpd_viscous_stress_local(); + std::remove_const_t global_stress; + + boost::mpi::reduce(comm, local_stress, global_stress, std::plus<>(), 0); - return Utils::flatten(stress) / volume; + return Utils::flatten(global_stress) / box_geo.volume(); } #endif // DPD diff --git a/src/core/dpd.hpp b/src/core/dpd.hpp index e921e67da9f..d464e3aa98b 100644 --- a/src/core/dpd.hpp +++ b/src/core/dpd.hpp @@ -34,6 +34,11 @@ #include +// Forward declaration +namespace boost::mpi { +class communicator; +} + struct IA_parameters; void dpd_init(double kT, double time_step); @@ -42,7 +47,7 @@ Utils::Vector3d dpd_pair_force(Particle const &p1, Particle const &p2, IA_parameters const &ia_params, Utils::Vector3d const &d, double dist, double dist2); -Utils::Vector9d dpd_stress(); +Utils::Vector9d dpd_stress(boost::mpi::communicator const &comm); #endif // DPD #endif diff --git a/src/core/energy.cpp b/src/core/energy.cpp index f45174017bd..8b1645e2108 100644 --- a/src/core/energy.cpp +++ b/src/core/energy.cpp @@ -38,7 +38,6 @@ #include "magnetostatics/dipoles.hpp" #include -#include #include @@ -114,18 +113,6 @@ std::shared_ptr calculate_energy() { // NOLINTNEXTLINE(clang-analyzer-cplusplus.NewDeleteLeaks) } -REGISTER_CALLBACK_MAIN_RANK(calculate_energy) - -double mpi_calculate_potential_energy() { - auto const obs = mpi_call(Communication::Result::main_rank, calculate_energy); - return obs->accumulate(-obs->kinetic[0]); -} - -double mpi_observable_compute_energy() { - auto const obs = mpi_call(Communication::Result::main_rank, calculate_energy); - return obs->accumulate(0); -} - double particle_short_range_energy_contribution(int pid) { double ret = 0.0; @@ -159,8 +146,4 @@ void calc_long_range_fields() { auto const &dipoles = System::get_system().dipoles; dipoles.calc_long_range_field(particles); } - -REGISTER_CALLBACK(calc_long_range_fields) - -void mpi_calc_long_range_fields() { mpi_call_all(calc_long_range_fields); } #endif diff --git a/src/core/energy.hpp b/src/core/energy.hpp index b20825a7270..76ba28a2fff 100644 --- a/src/core/energy.hpp +++ b/src/core/energy.hpp @@ -33,12 +33,6 @@ /** Parallel energy calculation. */ std::shared_ptr calculate_energy(); -/** Calculate the total energy of the system. */ -double mpi_calculate_potential_energy(); - -/** Helper function for @ref Observables::Energy. */ -double mpi_observable_compute_energy(); - /** * @brief Compute short-range energy of a particle. * @@ -50,9 +44,9 @@ double mpi_observable_compute_energy(); */ double particle_short_range_energy_contribution(int pid); -#endif #ifdef DIPOLE_FIELD_TRACKING /** Calculate dipole fields. */ void calc_long_range_fields(); -void mpi_calc_long_range_fields(); +#endif // DIPOLE_FIELD_TRACKING + #endif diff --git a/src/core/grid_based_algorithms/lb_interface.cpp b/src/core/grid_based_algorithms/lb_interface.cpp index c573b88ee1e..3c183b170c3 100644 --- a/src/core/grid_based_algorithms/lb_interface.cpp +++ b/src/core/grid_based_algorithms/lb_interface.cpp @@ -28,6 +28,7 @@ #include "grid.hpp" #include +#include #include #include @@ -126,40 +127,15 @@ double get_kT() { double get_lattice_speed() { return get_agrid() / get_tau(); } -#ifdef WALBERLA -namespace Walberla { - -static Utils::Vector3d get_momentum() { return lb_walberla()->get_momentum(); } - -static boost::optional -get_velocity_at_pos(Utils::Vector3d pos) { - return lb_walberla()->get_velocity_at_pos(pos); -} - -REGISTER_CALLBACK_ONE_RANK(get_velocity_at_pos) - -static boost::optional -get_interpolated_density_at_pos(Utils::Vector3d pos) { - return lb_walberla()->get_interpolated_density_at_pos(pos); -} - -REGISTER_CALLBACK_ONE_RANK(get_interpolated_density_at_pos) - -static Utils::VectorXd<9> get_pressure_tensor() { - return lb_walberla()->get_pressure_tensor(); -} - -REGISTER_CALLBACK_REDUCTION(get_pressure_tensor, std::plus<>()) - -} // namespace Walberla -#endif // WALBERLA - -Utils::VectorXd<9> const get_pressure_tensor() { +Utils::VectorXd<9> const +get_pressure_tensor(boost::mpi::communicator const &comm) { if (lattice_switch == ActiveLB::WALBERLA_LB) { #ifdef WALBERLA - return ::Communication::mpiCallbacks().call( - ::Communication::Result::reduction, std::plus<>(), - Walberla::get_pressure_tensor); + auto const local_pressure_tensor = lb_walberla()->get_pressure_tensor(); + std::remove_const_t pressure_tensor; + boost::mpi::reduce(comm, local_pressure_tensor, pressure_tensor, + std::plus<>(), 0); + return pressure_tensor; #endif } throw NoLBActive(); @@ -168,30 +144,28 @@ Utils::VectorXd<9> const get_pressure_tensor() { Utils::Vector3d calc_fluid_momentum() { if (lattice_switch == ActiveLB::WALBERLA_LB) { #ifdef WALBERLA - return Walberla::get_momentum(); + return lb_walberla()->get_momentum(); #endif } throw NoLBActive(); } -Utils::Vector3d const get_interpolated_velocity(Utils::Vector3d const &pos) { +boost::optional +get_interpolated_velocity(Utils::Vector3d const &pos) { if (lattice_switch == ActiveLB::WALBERLA_LB) { #ifdef WALBERLA - auto const folded_pos = folded_position(pos, box_geo); - return mpi_call(::Communication::Result::one_rank, - Walberla::get_velocity_at_pos, folded_pos / get_agrid()); + auto const folded_pos = folded_position(pos, box_geo) / get_agrid(); + return lb_walberla()->get_velocity_at_pos(folded_pos); #endif } throw NoLBActive(); } -double get_interpolated_density(Utils::Vector3d const &pos) { +boost::optional get_interpolated_density(Utils::Vector3d const &pos) { if (lattice_switch == ActiveLB::WALBERLA_LB) { #ifdef WALBERLA - auto const folded_pos = folded_position(pos, box_geo); - return mpi_call(::Communication::Result::one_rank, - Walberla::get_interpolated_density_at_pos, - folded_pos / get_agrid()); + auto const folded_pos = folded_position(pos, box_geo) / get_agrid(); + return lb_walberla()->get_interpolated_density_at_pos(folded_pos); #endif } throw NoLBActive(); diff --git a/src/core/grid_based_algorithms/lb_interface.hpp b/src/core/grid_based_algorithms/lb_interface.hpp index ad8a811f322..6e8236da2af 100644 --- a/src/core/grid_based_algorithms/lb_interface.hpp +++ b/src/core/grid_based_algorithms/lb_interface.hpp @@ -23,10 +23,17 @@ #include +#include + #include #include #include +// Forward Declarations +namespace boost::mpi { +class communicator; +} // namespace boost::mpi + /** @brief LB implementation currently active. */ enum class ActiveLB : int { NONE, WALBERLA_LB }; @@ -102,23 +109,27 @@ double get_lattice_speed(); * over all nodes and dividing by the number of nodes. * Returns the lower triangle of the LB pressure tensor. */ -Utils::VectorXd<9> const get_pressure_tensor(); +Utils::VectorXd<9> const +get_pressure_tensor(boost::mpi::communicator const &comm); Utils::Vector3d calc_fluid_momentum(); /** - * @brief Calculates the interpolated fluid velocity on the head node process. + * @brief Calculates the interpolated fluid velocity on the responsible node + * process. * @param pos Position at which the velocity is to be calculated. - * @retval interpolated fluid velocity. + * @retval optional interpolated fluid velocity. */ -Utils::Vector3d const get_interpolated_velocity(Utils::Vector3d const &pos); +boost::optional +get_interpolated_velocity(Utils::Vector3d const &pos); /** - * @brief Calculates the interpolated fluid density on the head node process. + * @brief Calculates the interpolated fluid density on the responsible node + * process. * @param pos Position at which the density is to be calculated. - * @retval interpolated fluid density. + * @retval optional interpolated fluid density. */ -double get_interpolated_density(Utils::Vector3d const &pos); +boost::optional get_interpolated_density(Utils::Vector3d const &pos); } // namespace LB diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index 21e8e60fe28..23b7f323263 100644 --- a/src/core/integrate.cpp +++ b/src/core/integrate.cpp @@ -57,6 +57,7 @@ #include "thermostat.hpp" #include "virtual_sites.hpp" +#include #include #include @@ -516,11 +517,6 @@ int integrate_with_signal_handler(int n_steps, int reuse_forces, ::set_skin(new_skin); } - // re-acquire MpiCallbacks listener on worker nodes - if (not is_head_node) { - return 0; - } - using Accumulators::auto_update; using Accumulators::auto_update_next_update; @@ -528,15 +524,20 @@ int integrate_with_signal_handler(int n_steps, int reuse_forces, /* 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()); - auto const retval = mpi_call(Communication::Result::main_rank, integrate, - steps, reuse_forces); - if (retval < 0) { - return retval; // propagate error code + + auto const local_retval = integrate(steps, reuse_forces); + + // make sure all ranks exit when one rank fails + std::remove_const_t global_retval; + boost::mpi::all_reduce(comm_cart, local_retval, global_retval, + std::plus()); + if (global_retval < 0) { + return global_retval; // propagate error code } reuse_forces = INTEG_REUSE_FORCES_ALWAYS; - auto_update(steps); + auto_update(comm_cart, steps); i += steps; } @@ -544,8 +545,6 @@ int integrate_with_signal_handler(int n_steps, int reuse_forces, return 0; } -REGISTER_CALLBACK_MAIN_RANK(integrate) - double interaction_range() { /* Consider skin only if there are actually interactions */ auto const max_cut = maximal_cutoff(n_nodes == 1); diff --git a/src/core/observables/BondAngles.hpp b/src/core/observables/BondAngles.hpp index 2effe5c9f82..56500e1a378 100644 --- a/src/core/observables/BondAngles.hpp +++ b/src/core/observables/BondAngles.hpp @@ -48,15 +48,22 @@ class BondAngles : public PidObservable { } std::vector - evaluate(ParticleReferenceRange particles, + evaluate(boost::mpi::communicator const &comm, + ParticleReferenceRange const &local_particles, const ParticleObservables::traits &traits) const override { + auto const positions_sorted = detail::get_all_particle_positions( + comm, local_particles, ids(), traits, false); + + if (comm.rank() != 0) { + return {}; + } + std::vector res(n_values()); - auto v1 = box_geo.get_mi_vector(traits.position(particles[1]), - traits.position(particles[0])); + auto v1 = box_geo.get_mi_vector(positions_sorted[1], positions_sorted[0]); auto n1 = v1.norm(); for (std::size_t i = 0, end = n_values(); i < end; i++) { - auto v2 = box_geo.get_mi_vector(traits.position(particles[i + 2]), - traits.position(particles[i + 1])); + auto v2 = box_geo.get_mi_vector(positions_sorted[i + 2], + positions_sorted[i + 1]); auto const n2 = v2.norm(); auto const cosine = std::clamp((v1 * v2) / (n1 * n2), -TINY_COS_VALUE, TINY_COS_VALUE); diff --git a/src/core/observables/BondDihedrals.hpp b/src/core/observables/BondDihedrals.hpp index 2df648b7bdf..864628cd915 100644 --- a/src/core/observables/BondDihedrals.hpp +++ b/src/core/observables/BondDihedrals.hpp @@ -23,7 +23,6 @@ #include "PidObservable.hpp" #include "grid.hpp" -#include #include #include @@ -53,17 +52,23 @@ class BondDihedrals : public PidObservable { } std::vector - evaluate(ParticleReferenceRange particles, + evaluate(boost::mpi::communicator const &comm, + ParticleReferenceRange const &local_particles, const ParticleObservables::traits &traits) const override { + auto const positions_sorted = detail::get_all_particle_positions( + comm, local_particles, ids(), traits, false); + + if (comm.rank() != 0) { + return {}; + } + std::vector res(n_values()); - auto v1 = box_geo.get_mi_vector(traits.position(particles[1]), - traits.position(particles[0])); - auto v2 = box_geo.get_mi_vector(traits.position(particles[2]), - traits.position(particles[1])); + auto v1 = box_geo.get_mi_vector(positions_sorted[1], positions_sorted[0]); + auto v2 = box_geo.get_mi_vector(positions_sorted[2], positions_sorted[1]); auto c1 = Utils::vector_product(v1, v2); for (std::size_t i = 0, end = n_values(); i < end; i++) { - auto v3 = box_geo.get_mi_vector(traits.position(particles[i + 3]), - traits.position(particles[i + 2])); + auto v3 = box_geo.get_mi_vector(positions_sorted[i + 3], + positions_sorted[i + 2]); auto c2 = vector_product(v2, v3); /* the 2-argument arctangent returns an angle in the range [-pi, pi] that * allows for an unambiguous determination of the 4th particle position */ diff --git a/src/core/observables/CosPersistenceAngles.hpp b/src/core/observables/CosPersistenceAngles.hpp index 8faf5662f53..10c7f892c1d 100644 --- a/src/core/observables/CosPersistenceAngles.hpp +++ b/src/core/observables/CosPersistenceAngles.hpp @@ -23,7 +23,6 @@ #include "PidObservable.hpp" #include "grid.hpp" -#include #include #include @@ -49,15 +48,24 @@ class CosPersistenceAngles : public PidObservable { } std::vector - evaluate(ParticleReferenceRange particles, + evaluate(boost::mpi::communicator const &comm, + ParticleReferenceRange const &local_particles, const ParticleObservables::traits &traits) const override { + + auto const positions_sorted = detail::get_all_particle_positions( + comm, local_particles, ids(), traits, false); + + if (comm.rank() != 0) { + return {}; + } + auto const no_of_angles = n_values(); auto const no_of_bonds = no_of_angles + 1; std::vector angles(no_of_angles); std::vector bond_vectors(no_of_bonds); auto get_bond_vector = [&](auto index) { - return box_geo.get_mi_vector(traits.position(particles[index + 1]), - traits.position(particles[index])); + return box_geo.get_mi_vector(positions_sorted[index + 1], + positions_sorted[index]); }; for (std::size_t i = 0; i < no_of_bonds; ++i) { auto const tmp = get_bond_vector(i); diff --git a/src/core/observables/CylindricalDensityProfile.hpp b/src/core/observables/CylindricalDensityProfile.hpp index c439d304071..01e387e6ba6 100644 --- a/src/core/observables/CylindricalDensityProfile.hpp +++ b/src/core/observables/CylindricalDensityProfile.hpp @@ -35,15 +35,38 @@ class CylindricalDensityProfile : public CylindricalPidProfileObservable { public: using CylindricalPidProfileObservable::CylindricalPidProfileObservable; std::vector - evaluate(ParticleReferenceRange particles, + evaluate(boost::mpi::communicator const &comm, + ParticleReferenceRange const &local_particles, const ParticleObservables::traits &traits) const override { + using pos_type = Utils::Vector3d; + + std::vector local_folded_positions{}; + local_folded_positions.reserve(local_particles.size()); + + for (auto const &p : local_particles) { + auto const pos = folded_position(traits.position(p), box_geo); + auto const pos_shifted = pos - transform_params->center(); + + local_folded_positions.emplace_back( + Utils::transform_coordinate_cartesian_to_cylinder( + pos_shifted, transform_params->axis(), + transform_params->orientation())); + } + + std::vector global_folded_positions{}; + boost::mpi::gather(comm, local_folded_positions, global_folded_positions, + 0); + + if (comm.rank() != 0) { + return {}; + } + Utils::CylindricalHistogram histogram(n_bins(), limits()); - for (auto const &p : particles) { - histogram.update(Utils::transform_coordinate_cartesian_to_cylinder( - folded_position(traits.position(p), box_geo) - - transform_params->center(), - transform_params->axis(), transform_params->orientation())); + for (auto const &vec : global_folded_positions) { + for (auto const &p : vec) { + histogram.update(p); + } } histogram.normalize(); diff --git a/src/core/observables/CylindricalFluxDensityProfile.hpp b/src/core/observables/CylindricalFluxDensityProfile.hpp index 4a05729b168..5796705a601 100644 --- a/src/core/observables/CylindricalFluxDensityProfile.hpp +++ b/src/core/observables/CylindricalFluxDensityProfile.hpp @@ -22,11 +22,11 @@ #include "BoxGeometry.hpp" #include "CylindricalPidProfileObservable.hpp" #include "grid.hpp" +#include "utils_histogram.hpp" #include #include -#include #include #include #include @@ -37,20 +37,44 @@ class CylindricalFluxDensityProfile : public CylindricalPidProfileObservable { using CylindricalPidProfileObservable::CylindricalPidProfileObservable; std::vector - evaluate(Utils::Span> particles, + evaluate(boost::mpi::communicator const &comm, + ParticleReferenceRange const &local_particles, const ParticleObservables::traits &traits) const override { - Utils::CylindricalHistogram histogram(n_bins(), limits()); + using pos_type = decltype(traits.position(std::declval())); + using vel_type = decltype(traits.velocity(std::declval())); + + std::vector local_folded_positions{}; + local_folded_positions.reserve(local_particles.size()); + std::vector local_velocities{}; + local_velocities.reserve(local_particles.size()); - // Write data to the histogram - for (auto p : particles) { + for (auto const &p : local_particles) { auto const pos = folded_position(traits.position(p), box_geo) - transform_params->center(); - histogram.update( + + local_folded_positions.emplace_back( Utils::transform_coordinate_cartesian_to_cylinder( - pos, transform_params->axis(), transform_params->orientation()), + pos, transform_params->axis(), transform_params->orientation())); + local_velocities.emplace_back( Utils::transform_vector_cartesian_to_cylinder( traits.velocity(p), transform_params->axis(), pos)); } + + auto const world_size = comm.size(); + std::vector global_folded_positions{}; + std::vector global_velocities{}; + global_folded_positions.reserve(world_size); + global_velocities.reserve(world_size); + boost::mpi::gather(comm, local_folded_positions, global_folded_positions, + 0); + boost::mpi::gather(comm, local_velocities, global_velocities, 0); + + if (comm.rank() != 0) { + return {}; + } + + Utils::CylindricalHistogram histogram(n_bins(), limits()); + accumulate(histogram, global_folded_positions, global_velocities); histogram.normalize(); return histogram.get_histogram(); } diff --git a/src/core/observables/CylindricalLBFluxDensityProfileAtParticlePositions.cpp b/src/core/observables/CylindricalLBFluxDensityProfileAtParticlePositions.cpp index 0d8cea95aaa..d881e82dcb9 100644 --- a/src/core/observables/CylindricalLBFluxDensityProfileAtParticlePositions.cpp +++ b/src/core/observables/CylindricalLBFluxDensityProfileAtParticlePositions.cpp @@ -21,43 +21,57 @@ #include "BoxGeometry.hpp" #include "grid.hpp" #include "grid_based_algorithms/lb_interface.hpp" +#include "utils_histogram.hpp" #include -#include #include +#include +#include + +#include #include namespace Observables { std::vector CylindricalLBFluxDensityProfileAtParticlePositions::evaluate( - Utils::Span> particles, + boost::mpi::communicator const &comm, + ParticleReferenceRange const &local_particles, const ParticleObservables::traits &traits) const { - Utils::CylindricalHistogram histogram(n_bins(), limits()); - // First collect all positions (since we want to call the LB function to - // get the fluid velocities only once). + using pos_type = decltype(traits.position(std::declval())); + using flux_type = Utils::Vector3d; + + std::vector local_folded_positions{}; + std::vector local_flux_densities{}; + local_folded_positions.reserve(local_particles.size()); + local_flux_densities.reserve(local_particles.size()); + + auto const vel_conv = LB::get_lattice_speed(); - for (auto p : particles) { + for (auto const &p : local_particles) { auto const pos = folded_position(traits.position(p), box_geo); - auto const v = LB::get_interpolated_velocity(pos) * LB::get_lattice_speed(); - auto const flux_dens = LB::get_interpolated_density(pos) * v; - - histogram.update(Utils::transform_coordinate_cartesian_to_cylinder( - pos - transform_params->center(), - transform_params->axis(), - transform_params->orientation()), - Utils::transform_vector_cartesian_to_cylinder( - flux_dens, transform_params->axis(), - pos - transform_params->center())); + auto const pos_shifted = pos - transform_params->center(); + auto const vel = *(LB::get_interpolated_velocity(pos)); + auto const dens = *(LB::get_interpolated_density(pos)); + auto const pos_cyl = Utils::transform_coordinate_cartesian_to_cylinder( + pos_shifted, transform_params->axis(), transform_params->orientation()); + auto const flux_cyl = Utils::transform_vector_cartesian_to_cylinder( + vel * vel_conv * dens, transform_params->axis(), pos_shifted); + local_folded_positions.emplace_back(pos_cyl); + local_flux_densities.emplace_back(flux_cyl); } - // normalize by number of hits per bin - auto hist_tmp = histogram.get_histogram(); - auto tot_count = histogram.get_tot_count(); - std::transform(hist_tmp.begin(), hist_tmp.end(), tot_count.begin(), - hist_tmp.begin(), [](auto hi, auto ci) { - return ci > 0 ? hi / static_cast(ci) : 0.; - }); - return hist_tmp; + std::vector global_folded_positions{}; + std::vector global_flux_densities{}; + boost::mpi::gather(comm, local_folded_positions, global_folded_positions, 0); + boost::mpi::gather(comm, local_flux_densities, global_flux_densities, 0); + + if (comm.rank() != 0) { + return {}; + } + + Utils::CylindricalHistogram histogram(n_bins(), limits()); + accumulate(histogram, global_folded_positions, global_flux_densities); + return normalize_by_bin_size(histogram); } } // namespace Observables diff --git a/src/core/observables/CylindricalLBFluxDensityProfileAtParticlePositions.hpp b/src/core/observables/CylindricalLBFluxDensityProfileAtParticlePositions.hpp index 7ac255b75c5..1bad6c6cf77 100644 --- a/src/core/observables/CylindricalLBFluxDensityProfileAtParticlePositions.hpp +++ b/src/core/observables/CylindricalLBFluxDensityProfileAtParticlePositions.hpp @@ -24,10 +24,7 @@ #include "Particle.hpp" #include "observables/ParticleTraits.hpp" -#include - #include -#include #include namespace Observables { @@ -37,7 +34,8 @@ class CylindricalLBFluxDensityProfileAtParticlePositions using CylindricalPidProfileObservable::CylindricalPidProfileObservable; std::vector - evaluate(Utils::Span> particles, + evaluate(boost::mpi::communicator const &comm, + ParticleReferenceRange const &local_particles, const ParticleObservables::traits &traits) const override; std::vector shape() const override { diff --git a/src/core/observables/CylindricalLBProfileObservable.hpp b/src/core/observables/CylindricalLBProfileObservable.hpp index 8eda18957aa..1a34b0966cf 100644 --- a/src/core/observables/CylindricalLBProfileObservable.hpp +++ b/src/core/observables/CylindricalLBProfileObservable.hpp @@ -21,7 +21,6 @@ #include "CylindricalProfileObservable.hpp" -#include #include #include #include diff --git a/src/core/observables/CylindricalLBVelocityProfile.cpp b/src/core/observables/CylindricalLBVelocityProfile.cpp index b0b98bc7be3..2f0ae68ce59 100644 --- a/src/core/observables/CylindricalLBVelocityProfile.cpp +++ b/src/core/observables/CylindricalLBVelocityProfile.cpp @@ -20,33 +20,56 @@ #include "CylindricalLBVelocityProfile.hpp" #include "grid_based_algorithms/lb_interface.hpp" +#include "utils_histogram.hpp" #include #include -#include -#include +#include +#include + #include namespace Observables { -std::vector CylindricalLBVelocityProfile::operator()() const { - Utils::CylindricalHistogram histogram(n_bins(), limits()); - for (auto const &p : sampling_positions) { - auto const velocity = - LB::get_interpolated_velocity(p) * LB::get_lattice_speed(); - auto const pos_shifted = p - transform_params->center(); - auto const pos_cyl = Utils::transform_coordinate_cartesian_to_cylinder( - pos_shifted, transform_params->axis(), transform_params->orientation()); - histogram.update(pos_cyl, - Utils::transform_vector_cartesian_to_cylinder( - velocity, transform_params->axis(), pos_shifted)); +std::vector CylindricalLBVelocityProfile::operator()( + boost::mpi::communicator const &comm) const { + using vel_type = Utils::Vector3d; + + decltype(sampling_positions) local_positions{}; + std::vector local_velocities{}; + + auto const vel_conv = LB::get_lattice_speed(); + + for (auto const &pos : sampling_positions) { + if (auto const vel = LB::get_interpolated_velocity(pos)) { + auto const pos_shifted = pos - transform_params->center(); + auto const pos_cyl = Utils::transform_coordinate_cartesian_to_cylinder( + pos_shifted, transform_params->axis(), + transform_params->orientation()); + auto const vel_cyl = Utils::transform_vector_cartesian_to_cylinder( + (*vel) * vel_conv, transform_params->axis(), pos_shifted); + + local_positions.emplace_back(pos_cyl); + local_velocities.emplace_back(vel_cyl); + } } - auto hist_data = histogram.get_histogram(); - auto const tot_count = histogram.get_tot_count(); - std::transform(hist_data.begin(), hist_data.end(), tot_count.begin(), - hist_data.begin(), std::divides()); - return hist_data; + + auto const world_size = comm.size(); + std::vector global_positions{}; + std::vector global_velocities{}; + global_positions.reserve(world_size); + global_velocities.reserve(world_size); + boost::mpi::gather(comm, local_positions, global_positions, 0); + boost::mpi::gather(comm, local_velocities, global_velocities, 0); + + if (comm.rank() != 0) { + return {}; + } + + Utils::CylindricalHistogram histogram(n_bins(), limits()); + accumulate(histogram, global_positions, global_velocities); + return normalize_by_bin_size(histogram); } } // namespace Observables diff --git a/src/core/observables/CylindricalLBVelocityProfile.hpp b/src/core/observables/CylindricalLBVelocityProfile.hpp index 5a7cdf925d3..4b1ffaf9e69 100644 --- a/src/core/observables/CylindricalLBVelocityProfile.hpp +++ b/src/core/observables/CylindricalLBVelocityProfile.hpp @@ -28,7 +28,8 @@ namespace Observables { class CylindricalLBVelocityProfile : public CylindricalLBProfileObservable { public: using CylindricalLBProfileObservable::CylindricalLBProfileObservable; - std::vector operator()() const override; + std::vector + operator()(boost::mpi::communicator const &comm) const override; std::vector shape() const override { auto const b = n_bins(); return {b[0], b[1], b[2], 3}; diff --git a/src/core/observables/CylindricalLBVelocityProfileAtParticlePositions.cpp b/src/core/observables/CylindricalLBVelocityProfileAtParticlePositions.cpp index d7ec93b4c0c..a139561a164 100644 --- a/src/core/observables/CylindricalLBVelocityProfileAtParticlePositions.cpp +++ b/src/core/observables/CylindricalLBVelocityProfileAtParticlePositions.cpp @@ -21,40 +21,59 @@ #include "BoxGeometry.hpp" #include "grid.hpp" #include "grid_based_algorithms/lb_interface.hpp" +#include "utils_histogram.hpp" #include -#include #include -#include +#include +#include + +#include #include namespace Observables { std::vector CylindricalLBVelocityProfileAtParticlePositions::evaluate( - Utils::Span> particles, + boost::mpi::communicator const &comm, + ParticleReferenceRange const &local_particles, const ParticleObservables::traits &traits) const { - Utils::CylindricalHistogram histogram(n_bins(), limits()); + using pos_type = decltype(traits.position(std::declval())); + using vel_type = Utils::Vector3d; + + std::vector local_folded_positions{}; + std::vector local_velocities{}; + local_folded_positions.reserve(local_particles.size()); + local_velocities.reserve(local_particles.size()); + + auto const vel_conv = LB::get_lattice_speed(); - for (auto const &p : particles) { + for (auto const &p : local_particles) { auto const pos = folded_position(traits.position(p), box_geo); - auto const v = LB::get_interpolated_velocity(pos) * LB::get_lattice_speed(); - - histogram.update( - Utils::transform_coordinate_cartesian_to_cylinder( - pos - transform_params->center(), transform_params->axis(), - transform_params->orientation()), - Utils::transform_vector_cartesian_to_cylinder( - v, transform_params->axis(), pos - transform_params->center())); + auto const pos_shifted = pos - transform_params->center(); + auto const vel = *(LB::get_interpolated_velocity(pos)); + auto const pos_cyl = Utils::transform_coordinate_cartesian_to_cylinder( + pos_shifted, transform_params->axis(), transform_params->orientation()); + auto const vel_cyl = Utils::transform_vector_cartesian_to_cylinder( + vel * vel_conv, transform_params->axis(), pos_shifted); + local_folded_positions.emplace_back(pos_cyl); + local_velocities.emplace_back(vel_cyl); } - // normalize by number of hits per bin - auto hist_tmp = histogram.get_histogram(); - auto tot_count = histogram.get_tot_count(); - std::transform(hist_tmp.begin(), hist_tmp.end(), tot_count.begin(), - hist_tmp.begin(), [](auto hi, auto ci) { - return ci > 0 ? hi / static_cast(ci) : 0.; - }); - return hist_tmp; + auto const world_size = comm.size(); + std::vector global_folded_positions{}; + std::vector global_velocities{}; + global_folded_positions.reserve(world_size); + global_velocities.reserve(world_size); + boost::mpi::gather(comm, local_folded_positions, global_folded_positions, 0); + boost::mpi::gather(comm, local_velocities, global_velocities, 0); + + if (comm.rank() != 0) { + return {}; + } + + Utils::CylindricalHistogram histogram(n_bins(), limits()); + accumulate(histogram, global_folded_positions, global_velocities); + return normalize_by_bin_size(histogram); } } // namespace Observables diff --git a/src/core/observables/CylindricalLBVelocityProfileAtParticlePositions.hpp b/src/core/observables/CylindricalLBVelocityProfileAtParticlePositions.hpp index 4698c93c428..7f16b046e94 100644 --- a/src/core/observables/CylindricalLBVelocityProfileAtParticlePositions.hpp +++ b/src/core/observables/CylindricalLBVelocityProfileAtParticlePositions.hpp @@ -24,10 +24,7 @@ #include "Particle.hpp" #include "observables/ParticleTraits.hpp" -#include - #include -#include #include namespace Observables { @@ -37,7 +34,8 @@ class CylindricalLBVelocityProfileAtParticlePositions using CylindricalPidProfileObservable::CylindricalPidProfileObservable; std::vector - evaluate(ParticleReferenceRange particles, + evaluate(boost::mpi::communicator const &comm, + ParticleReferenceRange const &particles, const ParticleObservables::traits &) const override; std::vector shape() const override { diff --git a/src/core/observables/CylindricalProfileObservable.hpp b/src/core/observables/CylindricalProfileObservable.hpp index b36664bd507..988343f7ac3 100644 --- a/src/core/observables/CylindricalProfileObservable.hpp +++ b/src/core/observables/CylindricalProfileObservable.hpp @@ -21,19 +21,12 @@ #include "ProfileObservable.hpp" -#include -#include #include -#include #include -#include -#include -#include #include #include -#include namespace Observables { diff --git a/src/core/observables/CylindricalVelocityProfile.hpp b/src/core/observables/CylindricalVelocityProfile.hpp index 58cb0f8ab34..ab54d4610d1 100644 --- a/src/core/observables/CylindricalVelocityProfile.hpp +++ b/src/core/observables/CylindricalVelocityProfile.hpp @@ -22,12 +22,11 @@ #include "BoxGeometry.hpp" #include "CylindricalPidProfileObservable.hpp" #include "grid.hpp" +#include "utils_histogram.hpp" #include -#include #include -#include #include #include #include @@ -38,29 +37,46 @@ class CylindricalVelocityProfile : public CylindricalPidProfileObservable { using CylindricalPidProfileObservable::CylindricalPidProfileObservable; std::vector - evaluate(ParticleReferenceRange particles, + evaluate(boost::mpi::communicator const &comm, + ParticleReferenceRange const &local_particles, const ParticleObservables::traits &traits) const override { - Utils::CylindricalHistogram histogram(n_bins(), limits()); + using pos_type = Utils::Vector3d; + using vel_type = Utils::Vector3d; + + std::vector local_folded_positions{}; + std::vector local_velocities{}; + local_folded_positions.reserve(local_particles.size()); + local_velocities.reserve(local_particles.size()); - for (auto p : particles) { + for (auto const &p : local_particles) { auto const pos = folded_position(traits.position(p), box_geo) - transform_params->center(); - histogram.update( + local_folded_positions.emplace_back( Utils::transform_coordinate_cartesian_to_cylinder( - pos, transform_params->axis(), transform_params->orientation()), + pos, transform_params->axis(), transform_params->orientation())); + local_velocities.emplace_back( Utils::transform_vector_cartesian_to_cylinder( traits.velocity(p), transform_params->axis(), pos)); } - auto hist_tmp = histogram.get_histogram(); - auto tot_count = histogram.get_tot_count(); - for (std::size_t ind = 0; ind < hist_tmp.size(); ++ind) { - if (tot_count[ind] > 0) { - hist_tmp[ind] /= static_cast(tot_count[ind]); - } + auto const world_size = comm.size(); + std::vector global_folded_positions{}; + std::vector global_velocities{}; + global_folded_positions.reserve(world_size); + global_velocities.reserve(world_size); + boost::mpi::gather(comm, local_folded_positions, global_folded_positions, + 0); + boost::mpi::gather(comm, local_velocities, global_velocities, 0); + + if (comm.rank() != 0) { + return {}; } - return hist_tmp; + + Utils::CylindricalHistogram histogram(n_bins(), limits()); + accumulate(histogram, global_folded_positions, global_velocities); + return normalize_by_bin_size(histogram); } + std::vector shape() const override { auto const b = n_bins(); return {b[0], b[1], b[2], 3}; diff --git a/src/core/observables/DPDStress.hpp b/src/core/observables/DPDStress.hpp index 2170c3368ab..fa537f08ed3 100644 --- a/src/core/observables/DPDStress.hpp +++ b/src/core/observables/DPDStress.hpp @@ -30,7 +30,10 @@ namespace Observables { class DPDStress : public Observable { public: std::vector shape() const override { return {3, 3}; } - std::vector operator()() const override { return dpd_stress(); } + std::vector + operator()(boost::mpi::communicator const &comm) const override { + return dpd_stress(comm); + } }; } // Namespace Observables diff --git a/src/core/observables/DensityProfile.hpp b/src/core/observables/DensityProfile.hpp index 1dfd46ed995..325d38908c9 100644 --- a/src/core/observables/DensityProfile.hpp +++ b/src/core/observables/DensityProfile.hpp @@ -25,8 +25,8 @@ #include "grid.hpp" #include -#include +#include #include namespace Observables { @@ -36,13 +36,37 @@ class DensityProfile : public PidProfileObservable { using PidProfileObservable::PidProfileObservable; std::vector - evaluate(Utils::Span> particles, + evaluate(boost::mpi::communicator const &comm, + ParticleReferenceRange const &local_particles, const ParticleObservables::traits &traits) const override { + using pos_type = decltype(traits.position(std::declval())); + + std::vector local_folded_positions{}; + local_folded_positions.reserve(local_particles.size()); + + for (auto const &p : local_particles) { + local_folded_positions.emplace_back( + folded_position(traits.position(p), box_geo)); + } + + auto const world_size = comm.size(); + std::vector global_folded_positions{}; + global_folded_positions.reserve(world_size); + boost::mpi::gather(comm, local_folded_positions, global_folded_positions, + 0); + + if (comm.rank() != 0) { + return {}; + } + Utils::Histogram histogram(n_bins(), limits()); - for (auto p : particles) { - histogram.update(folded_position(traits.position(p), box_geo)); + for (auto const &vec : global_folded_positions) { + for (auto const &p : vec) { + histogram.update(p); + } } + histogram.normalize(); return histogram.get_histogram(); } diff --git a/src/core/observables/EnergyObservable.hpp b/src/core/observables/EnergyObservable.hpp index 40132153c01..9b09a33582e 100644 --- a/src/core/observables/EnergyObservable.hpp +++ b/src/core/observables/EnergyObservable.hpp @@ -30,10 +30,9 @@ namespace Observables { class Energy : public Observable { public: std::vector shape() const override { return {1}; } - std::vector operator()() const override { - std::vector res{1}; - res[0] = mpi_observable_compute_energy(); - return res; + std::vector + operator()(boost::mpi::communicator const &comm) const override { + return {calculate_energy()->accumulate(0)}; } }; diff --git a/src/core/observables/FluxDensityProfile.hpp b/src/core/observables/FluxDensityProfile.hpp index 690d17211bd..1304ce263eb 100644 --- a/src/core/observables/FluxDensityProfile.hpp +++ b/src/core/observables/FluxDensityProfile.hpp @@ -23,11 +23,14 @@ #include "Particle.hpp" #include "PidProfileObservable.hpp" #include "grid.hpp" +#include "utils_histogram.hpp" #include -#include + +#include #include +#include #include namespace Observables { @@ -40,14 +43,38 @@ class FluxDensityProfile : public PidProfileObservable { } std::vector - evaluate(Utils::Span> particles, + evaluate(boost::mpi::communicator const &comm, + ParticleReferenceRange const &local_particles, const ParticleObservables::traits &traits) const override { - Utils::Histogram histogram(n_bins(), limits()); + using pos_type = decltype(traits.position(std::declval())); + using vel_type = decltype(traits.velocity(std::declval())); + + std::vector local_folded_positions{}; + std::vector local_velocities{}; + local_folded_positions.reserve(local_particles.size()); + local_velocities.reserve(local_particles.size()); + + for (auto const &p : local_particles) { + local_folded_positions.emplace_back( + folded_position(traits.position(p), box_geo)); + local_velocities.emplace_back(traits.velocity(p)); + } + + auto const world_size = comm.size(); + std::vector global_folded_positions{}; + std::vector global_velocities{}; + global_folded_positions.reserve(world_size); + global_velocities.reserve(world_size); + boost::mpi::gather(comm, local_folded_positions, global_folded_positions, + 0); + boost::mpi::gather(comm, local_velocities, global_velocities, 0); - for (auto p : particles) { - auto const ppos = folded_position(traits.position(p), box_geo); - histogram.update(ppos, traits.velocity(p)); + if (comm.rank() != 0) { + return {}; } + + Utils::Histogram histogram(n_bins(), limits()); + accumulate(histogram, global_folded_positions, global_velocities); histogram.normalize(); return histogram.get_histogram(); } diff --git a/src/core/observables/ForceDensityProfile.hpp b/src/core/observables/ForceDensityProfile.hpp index 348e2f3ab44..3bdda403ca0 100644 --- a/src/core/observables/ForceDensityProfile.hpp +++ b/src/core/observables/ForceDensityProfile.hpp @@ -23,11 +23,15 @@ #include "Particle.hpp" #include "PidProfileObservable.hpp" #include "grid.hpp" +#include "utils_histogram.hpp" #include -#include +#include +#include + #include +#include #include namespace Observables { @@ -41,13 +45,38 @@ class ForceDensityProfile : public PidProfileObservable { } std::vector - evaluate(ParticleReferenceRange particles, - const ParticleObservables::traits &) const override { - Utils::Histogram histogram(n_bins(), limits()); - for (auto const &p : particles) { - histogram.update(folded_position(p.get().pos(), box_geo), - p.get().force()); + evaluate(boost::mpi::communicator const &comm, + ParticleReferenceRange const &local_particles, + const ParticleObservables::traits &traits) const override { + using pos_type = decltype(traits.position(std::declval())); + using force_type = decltype(traits.force(std::declval())); + + std::vector local_folded_positions{}; + local_folded_positions.reserve(local_particles.size()); + std::vector local_forces{}; + local_forces.reserve(local_particles.size()); + + for (auto const &p : local_particles) { + local_folded_positions.emplace_back( + folded_position(traits.position(p), box_geo)); + local_forces.emplace_back(traits.force(p)); } + + auto const world_size = comm.size(); + std::vector global_folded_positions{}; + std::vector global_forces{}; + global_folded_positions.reserve(world_size); + global_forces.reserve(world_size); + boost::mpi::gather(comm, local_folded_positions, global_folded_positions, + 0); + boost::mpi::gather(comm, local_forces, global_forces, 0); + + if (comm.rank() != 0) { + return {}; + } + + Utils::Histogram histogram(n_bins(), limits()); + accumulate(histogram, global_folded_positions, global_forces); histogram.normalize(); return histogram.get_histogram(); } diff --git a/src/core/observables/LBFluidPressureTensor.hpp b/src/core/observables/LBFluidPressureTensor.hpp index 933d4c17c11..cafeaa643d2 100644 --- a/src/core/observables/LBFluidPressureTensor.hpp +++ b/src/core/observables/LBFluidPressureTensor.hpp @@ -31,10 +31,11 @@ namespace Observables { class LBFluidPressureTensor : public Observable { public: std::vector shape() const override { return {3, 3}; } - std::vector operator()() const override { + std::vector + operator()(boost::mpi::communicator const &comm) const override { auto const unit_conversion = 1. / (LB::get_agrid() * Utils::sqr(LB::get_tau())); - auto const tensor = LB::get_pressure_tensor() * unit_conversion; + auto const tensor = LB::get_pressure_tensor(comm) * unit_conversion; return tensor.as_vector(); } }; diff --git a/src/core/observables/LBVelocityProfile.cpp b/src/core/observables/LBVelocityProfile.cpp index 4c8b2740d37..c15c6c63100 100644 --- a/src/core/observables/LBVelocityProfile.cpp +++ b/src/core/observables/LBVelocityProfile.cpp @@ -19,35 +19,55 @@ #include "LBVelocityProfile.hpp" #include "grid_based_algorithms/lb_interface.hpp" +#include "utils_histogram.hpp" #include +#include + +#include +#include -#include #include -#include #include namespace Observables { -std::vector LBVelocityProfile::operator()() const { - Utils::Histogram histogram(n_bins(), limits()); - for (auto const &p : sampling_positions) { - const auto v = LB::get_interpolated_velocity(p) * LB::get_lattice_speed(); - histogram.update(p, v); - } - auto hist_tmp = histogram.get_histogram(); - auto const tot_count = histogram.get_tot_count(); - for (std::size_t ind = 0; ind < hist_tmp.size(); ++ind) { - if (tot_count[ind] == 0 and not allow_empty_bins) { - auto const error = "Decrease sampling delta(s), bin " + - std::to_string(ind) + " has no hit"; - throw std::runtime_error(error); - } - if (tot_count[ind] > 0) { - hist_tmp[ind] /= static_cast(tot_count[ind]); +std::vector +LBVelocityProfile::operator()(boost::mpi::communicator const &comm) const { + using vel_type = Utils::Vector3d; + + decltype(sampling_positions) local_positions{}; + std::vector local_velocities{}; + + auto const vel_conv = LB::get_lattice_speed(); + + for (auto const &pos : sampling_positions) { + if (auto const vel = LB::get_interpolated_velocity(pos)) { + local_positions.emplace_back(pos); + local_velocities.emplace_back((*vel) * vel_conv); } } - return hist_tmp; + + auto const world_size = comm.size(); + std::vector global_positions{}; + std::vector global_velocities{}; + global_positions.reserve(world_size); + global_velocities.reserve(world_size); + boost::mpi::gather(comm, local_positions, global_positions, 0); + boost::mpi::gather(comm, local_velocities, global_velocities, 0); + + if (comm.rank() != 0) { + return {}; + } + + Utils::Histogram histogram(n_bins(), limits()); + accumulate(histogram, global_positions, global_velocities); + try { + return normalize_by_bin_size(histogram, allow_empty_bins); + } catch (empty_bin_exception const &) { + throw std::runtime_error( + "Decrease sampling delta(s), some bins have no hit"); + } } } // namespace Observables diff --git a/src/core/observables/LBVelocityProfile.hpp b/src/core/observables/LBVelocityProfile.hpp index 53a25d012f0..2d0e8e7d472 100644 --- a/src/core/observables/LBVelocityProfile.hpp +++ b/src/core/observables/LBVelocityProfile.hpp @@ -33,7 +33,8 @@ class LBVelocityProfile : public LBProfileObservable { auto const b = n_bins(); return {b[0], b[1], b[2], 3}; } - std::vector operator()() const override; + std::vector + operator()(boost::mpi::communicator const &comm) const override; }; } // Namespace Observables diff --git a/src/core/observables/Observable.hpp b/src/core/observables/Observable.hpp index a3ab98dca9a..ac2c7b6409c 100644 --- a/src/core/observables/Observable.hpp +++ b/src/core/observables/Observable.hpp @@ -25,6 +25,8 @@ #include #include +#include + namespace Observables { /** Base class for observables. @@ -42,7 +44,8 @@ class Observable { Observable() = default; virtual ~Observable() = default; /** Calculate the set of values measured by the observable */ - virtual std::vector operator()() const = 0; + virtual std::vector + operator()(boost::mpi::communicator const &comm) const = 0; /** Size of the flat array returned by the observable */ std::size_t n_values() const { diff --git a/src/core/observables/ParticleDipoleFields.hpp b/src/core/observables/ParticleDipoleFields.hpp index 19c42904408..f58899f34be 100644 --- a/src/core/observables/ParticleDipoleFields.hpp +++ b/src/core/observables/ParticleDipoleFields.hpp @@ -38,13 +38,14 @@ class ParticleDipoleFields using ParticleObservable< ParticleObservables::DipoleFields>::ParticleObservable; std::vector - evaluate(ParticleReferenceRange particles, + evaluate(boost::mpi::communicator const &comm, + ParticleReferenceRange const &local_particles, const ParticleObservables::traits &traits) const override { #ifdef DIPOLE_FIELD_TRACKING - mpi_calc_long_range_fields(); + calc_long_range_fields(); #endif return ParticleObservable::evaluate( - particles, traits); + comm, local_particles, traits); } }; diff --git a/src/core/observables/ParticleDistances.hpp b/src/core/observables/ParticleDistances.hpp index 24d43973a29..1db61a1a690 100644 --- a/src/core/observables/ParticleDistances.hpp +++ b/src/core/observables/ParticleDistances.hpp @@ -45,13 +45,21 @@ class ParticleDistances : public PidObservable { } std::vector - evaluate(ParticleReferenceRange particles, + evaluate(boost::mpi::communicator const &comm, + ParticleReferenceRange const &local_particles, const ParticleObservables::traits &traits) const override { + auto const positions_sorted = detail::get_all_particle_positions( + comm, local_particles, ids(), traits, false); + + if (comm.rank() != 0) { + return {}; + } + std::vector res(n_values()); for (std::size_t i = 0, end = n_values(); i < end; i++) { - auto const v = box_geo.get_mi_vector(traits.position(particles[i]), - traits.position(particles[i + 1])); + auto const v = + box_geo.get_mi_vector(positions_sorted[i], positions_sorted[i + 1]); res[i] = v.norm(); } return res; diff --git a/src/core/observables/ParticleTraits.hpp b/src/core/observables/ParticleTraits.hpp index bc84f6f3e67..fe57ed0d419 100644 --- a/src/core/observables/ParticleTraits.hpp +++ b/src/core/observables/ParticleTraits.hpp @@ -19,8 +19,10 @@ #ifndef OBSERVABLES_PARTICLE_TRAITS #define OBSERVABLES_PARTICLE_TRAITS +#include "BoxGeometry.hpp" #include "Particle.hpp" #include "config/config.hpp" +#include "grid.hpp" #include "rotation.hpp" namespace ParticleObservables { @@ -30,7 +32,11 @@ namespace ParticleObservables { * of observables independent of the particle type. */ template <> struct traits { - auto position(Particle const &p) const { return p.pos(); } + auto id(Particle const &p) const { return p.id(); } + auto position(Particle const &p) const { + return unfolded_position(p.pos(), p.image_box(), box_geo.length()); + } + auto position_folded(Particle const &p) const { return p.pos(); } auto velocity(Particle const &p) const { return p.v(); } auto force(Particle const &p) const { return p.force(); } auto mass(Particle const &p) const { diff --git a/src/core/observables/PidObservable.cpp b/src/core/observables/PidObservable.cpp index 6773eae7f40..2932f91c64e 100644 --- a/src/core/observables/PidObservable.cpp +++ b/src/core/observables/PidObservable.cpp @@ -21,16 +21,15 @@ #include "ParticleTraits.hpp" #include "fetch_particles.hpp" -#include +#include + #include namespace Observables { -std::vector PidObservable::operator()() const { - std::vector particles = fetch_particles(ids()); - - std::vector> particle_refs( - particles.begin(), particles.end()); - return this->evaluate(ParticleReferenceRange(particle_refs), +std::vector +PidObservable::operator()(boost::mpi::communicator const &comm) const { + auto const &local_particles = fetch_particles(ids()); + return this->evaluate(comm, local_particles, ParticleObservables::traits{}); } } // namespace Observables diff --git a/src/core/observables/PidObservable.hpp b/src/core/observables/PidObservable.hpp index c19ccbad433..ee2e0ba3f0f 100644 --- a/src/core/observables/PidObservable.hpp +++ b/src/core/observables/PidObservable.hpp @@ -26,14 +26,17 @@ #include "Particle.hpp" #include "ParticleTraits.hpp" -#include #include #include +#include +#include + #include +#include +#include #include -#include #include #include #include @@ -42,7 +45,7 @@ namespace Observables { using ParticleReferenceRange = - Utils::Span>; + std::vector>; /** %Particle-based observable. * @@ -54,12 +57,14 @@ class PidObservable : virtual public Observable { std::vector m_ids; virtual std::vector - evaluate(ParticleReferenceRange particles, + evaluate(boost::mpi::communicator const &comm, + ParticleReferenceRange const &local_particles, const ParticleObservables::traits &traits) const = 0; public: explicit PidObservable(std::vector ids) : m_ids(std::move(ids)) {} - std::vector operator()() const final; + std::vector + operator()(boost::mpi::communicator const &comm) const final; std::vector const &ids() const { return m_ids; } }; @@ -85,6 +90,97 @@ template struct shape_impl> { return ret; } }; +template struct shape_impl> { + static std::vector eval(std::size_t n_part) { + return shape_impl::eval(n_part); + } +}; + +inline auto get_argsort(boost::mpi::communicator const &comm, + std::vector const &local_pids, + std::vector const &sorted_pids) { + std::vector argsort{}; + + std::vector> global_pids; + boost::mpi::gather(comm, local_pids, global_pids, 0); + if (comm.rank() == 0) { + auto const n_part = sorted_pids.size(); + std::vector unsorted_pids; + unsorted_pids.reserve(n_part); + for (auto const &vec : global_pids) { + for (auto const pid : vec) { + unsorted_pids.emplace_back(pid); + } + } + // get vector of indices that sorts the data vectors + std::vector iota(n_part); + std::iota(iota.begin(), iota.end(), 0u); + argsort.reserve(n_part); + auto const pid_begin = std::begin(unsorted_pids); + auto const pid_end = std::end(unsorted_pids); + for (auto const pid : sorted_pids) { + auto const pid_pos = std::find(pid_begin, pid_end, pid); + auto const i = + static_cast(std::distance(pid_begin, pid_pos)); + argsort.emplace_back(iota[i]); + } + } + return argsort; +} + +/** Get the positions of all the particles in the system in the order + * specified by \a sorted_pids. Only the head node returns a vector + * of positions, all other nodes return an empty vector. + * This requires several MPI communications to construct. + */ +inline auto +get_all_particle_positions(boost::mpi::communicator const &comm, + ParticleReferenceRange const &local_particles, + std::vector const &sorted_pids, + ParticleObservables::traits const &traits, + bool use_folded_positions = false) { + using pos_type = decltype(traits.position(std::declval())); + std::vector local_positions{}; + std::vector local_pids{}; + local_positions.reserve(local_particles.size()); + local_pids.reserve(local_particles.size()); + + for (auto const &particle : local_particles) { + if (use_folded_positions) { + local_positions.emplace_back(traits.position_folded(particle)); + } else { + local_positions.emplace_back(traits.position(particle)); + } + local_pids.emplace_back(traits.id(particle)); + } + + auto const argsort = detail::get_argsort(comm, local_pids, sorted_pids); + + std::vector> global_positions{}; + global_positions.reserve(comm.size()); + boost::mpi::gather(comm, local_positions, global_positions, 0); + + if (comm.rank() != 0) { + return std::vector(); + } + + std::vector global_positions_flattened{}; + global_positions_flattened.reserve(sorted_pids.size()); + for (auto const &vec : global_positions) { + for (auto const &pos : vec) { + global_positions_flattened.emplace_back(pos); + } + } + assert(global_positions_flattened.size() == sorted_pids.size()); + + std::vector positions_sorted{}; + positions_sorted.reserve(sorted_pids.size()); + for (auto const i : argsort) { + positions_sorted.emplace_back(global_positions_flattened[i]); + } + + return positions_sorted; +} } // namespace detail /** @@ -108,15 +204,69 @@ template class ParticleObservable : public PidObservable { using std::declval; return detail::shape_impl()( - declval()))>::eval(ids().size()); + declval()))>::eval(ids().size()); } + template struct is_map : std::false_type {}; + template + struct is_map> : std::true_type {}; + std::vector - evaluate(ParticleReferenceRange particles, - const ParticleObservables::traits &) const override { - std::vector res; - Utils::flatten(ObsType{}(particles), std::back_inserter(res)); - return res; + evaluate(boost::mpi::communicator const &comm, + ParticleReferenceRange const &local_particles, + ParticleObservables::traits const &traits) const override { + if constexpr (is_map::value) { + std::vector local_traits{}; + local_traits.reserve(local_particles.size()); + Utils::flatten(ObsType{}(local_particles), + std::back_inserter(local_traits)); + std::vector local_pids{}; + local_pids.reserve(local_particles.size()); + Utils::flatten(ParticleObservables::Identities{}(local_particles), + std::back_inserter(local_pids)); + + std::vector> global_traits{}; + boost::mpi::gather(comm, local_traits, global_traits, 0); + + auto const argsort = detail::get_argsort(comm, local_pids, ids()); + + if (comm.rank() != 0) { + return {}; + } + + // get total size of the global traits vector + auto const size = std::accumulate( + global_traits.begin(), global_traits.end(), 0u, + [](auto const acc, auto const &vec) { return acc + vec.size(); }); + + auto const n_dims = size / ids().size(); + + std::vector global_traits_flattened{}; + global_traits_flattened.reserve(size); + + for (auto const &vec : global_traits) { + for (auto const val : vec) { + global_traits_flattened.emplace_back(val); + } + } + + std::vector output{}; + output.reserve(size); + + for (auto const i : argsort) { + for (std::size_t j = 0; j < n_dims; ++j) { + output.emplace_back(global_traits_flattened[i * n_dims + j]); + } + } + return output; + } else { + auto observable = ObsType{}; + auto const local_result = observable(local_particles); + std::remove_const_t result{}; + boost::mpi::reduce(comm, local_result, result, observable, 0); + + return result.first; + } } }; diff --git a/src/core/observables/PressureObservable.hpp b/src/core/observables/PressureObservable.hpp index 1f9ec9cc8eb..a75580d7e99 100644 --- a/src/core/observables/PressureObservable.hpp +++ b/src/core/observables/PressureObservable.hpp @@ -21,6 +21,7 @@ #include "Observable.hpp" #include "pressure.hpp" + #include #include @@ -29,11 +30,13 @@ namespace Observables { class Pressure : public Observable { public: std::vector shape() const override { return {1}; } - std::vector operator()() const override { - auto const ptensor = mpi_observable_compute_pressure_tensor(); - std::vector res{1}; - res[0] = (ptensor[0] + ptensor[4] + ptensor[8]) / 3.; - return res; + std::vector + operator()(boost::mpi::communicator const &comm) const override { + auto const obs = calculate_pressure(); + + return {(obs->accumulate(0., 0u) + obs->accumulate(0., 4u) + + obs->accumulate(0., 8u)) / + 3.}; } }; diff --git a/src/core/observables/PressureTensor.hpp b/src/core/observables/PressureTensor.hpp index d9fc5eabcde..e15ec57b0a8 100644 --- a/src/core/observables/PressureTensor.hpp +++ b/src/core/observables/PressureTensor.hpp @@ -21,6 +21,7 @@ #include "Observable.hpp" #include "pressure.hpp" + #include #include @@ -29,8 +30,16 @@ namespace Observables { class PressureTensor : public Observable { public: std::vector shape() const override { return {3, 3}; } - std::vector operator()() const override { - return mpi_observable_compute_pressure_tensor().as_vector(); + std::vector + operator()(boost::mpi::communicator const &comm) const override { + auto const obs = calculate_pressure(); + + std::vector result; + result.reserve(9); + for (std::size_t i = 0u; i < 9u; ++i) { + result.emplace_back(obs->accumulate(0., i)); + } + return result; } }; diff --git a/src/core/observables/RDF.cpp b/src/core/observables/RDF.cpp index e2bf0720654..c4975ca4e66 100644 --- a/src/core/observables/RDF.cpp +++ b/src/core/observables/RDF.cpp @@ -22,47 +22,52 @@ #include "fetch_particles.hpp" #include "grid.hpp" -#include #include #include #include #include #include +#include #include -#include -#include #include namespace Observables { -std::vector RDF::operator()() const { - std::vector particles1 = fetch_particles(ids1()); - std::vector particles_ptrs1(particles1.size()); - boost::transform(particles1, particles_ptrs1.begin(), - [](auto const &p) { return std::addressof(p); }); +std::vector +RDF::operator()(boost::mpi::communicator const &comm) const { + auto const &local_particles_1 = fetch_particles(ids1()); if (ids2().empty()) { - return this->evaluate(particles_ptrs1, {}); + return this->evaluate(comm, local_particles_1, {}, {}); } - std::vector particles2 = fetch_particles(ids2()); - std::vector particles_ptrs2(particles2.size()); - boost::transform(particles2, particles_ptrs2.begin(), - [](auto const &p) { return std::addressof(p); }); - return this->evaluate(particles_ptrs1, particles_ptrs2); + auto const &local_particles_2 = fetch_particles(ids2()); + + return this->evaluate(comm, local_particles_1, local_particles_2, {}); } std::vector -RDF::evaluate(Utils::Span particles1, - Utils::Span particles2) const { +RDF::evaluate(boost::mpi::communicator const &comm, + ParticleReferenceRange const &local_particles_1, + ParticleReferenceRange const &local_particles_2, + const ParticleObservables::traits &traits) const { + auto const positions_1 = detail::get_all_particle_positions( + comm, local_particles_1, ids1(), traits, false); + auto const positions_2 = detail::get_all_particle_positions( + comm, local_particles_2, ids2(), traits, false); + + if (comm.rank() != 0) { + return {}; + } + auto const bin_width = (max_r - min_r) / static_cast(n_r_bins); auto const inv_bin_width = 1.0 / bin_width; std::vector res(n_values(), 0.0); long int cnt = 0; - auto op = [this, inv_bin_width, &cnt, &res](const Particle *const p1, - const Particle *const p2) { - auto const dist = box_geo.get_mi_vector(p1->pos(), p2->pos()).norm(); + auto op = [this, inv_bin_width, &cnt, &res](auto const &pos1, + auto const &pos2) { + auto const dist = box_geo.get_mi_vector(pos1, pos2).norm(); if (dist > min_r && dist < max_r) { auto const ind = static_cast(std::floor((dist - min_r) * inv_bin_width)); @@ -70,11 +75,27 @@ RDF::evaluate(Utils::Span particles1, } cnt++; }; - if (particles2.empty()) { - Utils::for_each_pair(particles1, op); + + if (local_particles_2.empty()) { + Utils::for_each_pair(positions_1, op); } else { - auto cmp = std::not_equal_to(); - Utils::for_each_cartesian_pair_if(particles1, particles2, op, cmp); + auto const combine_1 = boost::combine(ids1(), positions_1); + auto const combine_2 = boost::combine(ids2(), positions_2); + + auto op2 = [&op](auto const &it1, auto const &it2) { + auto const &[id1, pos1] = it1; + auto const &[id2, pos2] = it2; + + op(pos1, pos2); + }; + + auto cmp = [](auto const &it1, auto const &it2) { + auto const &[id1, pos1] = it1; + auto const &[id2, pos2] = it2; + + return id1 != id2; + }; + Utils::for_each_cartesian_pair_if(combine_1, combine_2, op2, cmp); } if (cnt == 0) return res; diff --git a/src/core/observables/RDF.hpp b/src/core/observables/RDF.hpp index ae8e7cde815..ccebb1bd37b 100644 --- a/src/core/observables/RDF.hpp +++ b/src/core/observables/RDF.hpp @@ -21,9 +21,8 @@ #define OBSERVABLES_RDF_HPP #include "Observable.hpp" -#include "Particle.hpp" -#include +#include "PidObservable.hpp" #include #include @@ -40,9 +39,11 @@ class RDF : public Observable { /** Identifiers of the distant particles */ std::vector m_ids2; - virtual std::vector - evaluate(Utils::Span particles1, - Utils::Span particles2) const; + std::vector + evaluate(boost::mpi::communicator const &comm, + ParticleReferenceRange const &local_particles_1, + ParticleReferenceRange const &local_particles_2, + const ParticleObservables::traits &traits) const; public: // Range of the profile. @@ -61,7 +62,8 @@ class RDF : public Observable { if (n_r_bins <= 0) throw std::domain_error("n_r_bins has to be >= 1"); } - std::vector operator()() const final; + std::vector + operator()(boost::mpi::communicator const &comm) const final; std::vector &ids1() { return m_ids1; } std::vector &ids2() { return m_ids2; } diff --git a/src/core/observables/TotalForce.hpp b/src/core/observables/TotalForce.hpp index e35b2577de5..7125de3f826 100644 --- a/src/core/observables/TotalForce.hpp +++ b/src/core/observables/TotalForce.hpp @@ -21,6 +21,8 @@ #include "PidObservable.hpp" +#include + #include #include @@ -31,15 +33,20 @@ class TotalForce : public PidObservable { std::vector shape() const override { return {3}; } std::vector - evaluate(ParticleReferenceRange particles, + evaluate(boost::mpi::communicator const &comm, + ParticleReferenceRange const &local_particles, const ParticleObservables::traits &) const override { - Utils::Vector3d res{}; - for (auto const &p : particles) { + Utils::Vector3d local_force{}; + for (auto const &p : local_particles) { if (p.get().is_virtual()) continue; - res += p.get().force(); + local_force += p.get().force(); } - return res.as_vector(); + + decltype(local_force) global_force; + boost::mpi::reduce(comm, local_force, global_force, std::plus<>(), 0); + + return global_force.as_vector(); } }; } // Namespace Observables diff --git a/src/core/observables/fetch_particles.hpp b/src/core/observables/fetch_particles.hpp index e6ab60a548b..af33a133df5 100644 --- a/src/core/observables/fetch_particles.hpp +++ b/src/core/observables/fetch_particles.hpp @@ -20,13 +20,12 @@ #ifndef FETCH_PARTICLES_HPP #define FETCH_PARTICLES_HPP -#include "grid.hpp" -#include "particle_node.hpp" - -#include +#include "PidObservable.hpp" +#include "cells.hpp" #include -#include +#include +#include #include /** Fetch a group of particles. @@ -34,30 +33,13 @@ * @param ids particle identifiers * @return array of particle copies, with positions in the current box. */ -inline std::vector fetch_particles(std::vector const &ids) { - std::vector particles; - particles.reserve(ids.size()); - - auto const chunk_size = fetch_cache_max_size(); - for (std::size_t offset = 0; offset < ids.size();) { - auto const this_size = std::clamp(chunk_size, std::size_t{0}, - std::size_t{ids.size() - offset}); - auto const chunk_ids = - Utils::make_const_span(ids.data() + offset, this_size); - - prefetch_particle_data(chunk_ids); - - for (auto id : chunk_ids) { - particles.push_back(get_particle_data(id)); - - auto &p = particles.back(); - p.pos() += image_shift(p.image_box(), box_geo.length()); - p.image_box() = {}; - } - - offset += this_size; - } - - return particles; +inline auto fetch_particles(std::vector const &ids) { + auto const ids_set = std::set{ids.begin(), ids.end()}; + auto const local_particles = ::cell_structure.local_particles(); + Observables::ParticleReferenceRange local_particle_refs; + std::copy_if(local_particles.begin(), local_particles.end(), + std::back_inserter(local_particle_refs), + [&ids_set](Particle &p) { return ids_set.count(p.id()) != 0; }); + return local_particle_refs; } #endif diff --git a/src/core/observables/utils_histogram.hpp b/src/core/observables/utils_histogram.hpp new file mode 100644 index 00000000000..7f4bbc003e6 --- /dev/null +++ b/src/core/observables/utils_histogram.hpp @@ -0,0 +1,61 @@ +/* + * Copyright (C) 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 . + */ + +#pragma once + +#include + +#include +#include + +namespace Observables { + +/** @brief Accumulate histogram data gathered from multiple MPI ranks. */ +template +void accumulate(Utils::Histogram &histogram, + std::vector> const &pos, + std::vector> const &val) { + for (std::size_t rank = 0u; rank < pos.size(); ++rank) { + auto const &pos_vec = pos[rank]; + auto const &val_vec = val[rank]; + for (std::size_t i = 0u; i < pos_vec.size(); ++i) { + histogram.update(pos_vec[i], val_vec[i]); + } + } +} + +struct empty_bin_exception {}; + +/** @brief Normalize histogram by the number of values in each bin. */ +template +auto normalize_by_bin_size(Utils::Histogram &histogram, + bool allow_empty_bins = true) { + auto hist_data = histogram.get_histogram(); + auto tot_count = histogram.get_tot_count(); + for (std::size_t i = 0u; i < hist_data.size(); ++i) { + if (tot_count[i] != 0u) { + hist_data[i] /= static_cast(tot_count[i]); + } else if (not allow_empty_bins) { + throw empty_bin_exception{}; + } + } + return hist_data; +} + +} // Namespace Observables diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp index da26f84200e..d8b8ab8f734 100644 --- a/src/core/pressure.cpp +++ b/src/core/pressure.cpp @@ -47,11 +47,7 @@ #include -#include -#include -#include #include -#include std::shared_ptr calculate_pressure() { @@ -125,15 +121,3 @@ std::shared_ptr calculate_pressure() { obs_pressure.mpi_reduce(); return obs_pressure_ptr; } - -REGISTER_CALLBACK_MAIN_RANK(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->accumulate(0, j); - } - return pressure_tensor; -} diff --git a/src/core/pressure.hpp b/src/core/pressure.hpp index 7fd55d37841..6ee5883575d 100644 --- a/src/core/pressure.hpp +++ b/src/core/pressure.hpp @@ -34,7 +34,4 @@ /** Parallel pressure calculation from a virial expansion. */ std::shared_ptr calculate_pressure(); -/** Helper function for @ref Observables::PressureTensor. */ -Utils::Vector9d mpi_observable_compute_pressure_tensor(); - #endif diff --git a/src/core/unit_tests/EspressoSystemStandAlone_test.cpp b/src/core/unit_tests/EspressoSystemStandAlone_test.cpp index 3c9de972eb5..abd19fb80f8 100644 --- a/src/core/unit_tests/EspressoSystemStandAlone_test.cpp +++ b/src/core/unit_tests/EspressoSystemStandAlone_test.cpp @@ -45,6 +45,7 @@ namespace utf = boost::unit_test; #include "magnetostatics/dipoles.hpp" #include "nonbonded_interactions/lj.hpp" #include "observables/ParticleVelocities.hpp" +#include "observables/PidObservable.hpp" #include "particle_node.hpp" #include "system/System.hpp" @@ -59,7 +60,9 @@ namespace utf = boost::unit_test; #include #include +#include #include +#include #include #include #include @@ -118,6 +121,46 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory) { } }; + // check observables + { + auto const pid4 = 10; + auto const pids = std::vector{pid2, pid3, pid1, pid4}; + Observables::ParticleReferenceRange particle_range{}; + for (int pid : pids) { + if (auto const p = ::cell_structure.get_local_particle(pid)) { + particle_range.emplace_back(*p); + } + } + Particle p{}; + p.id() = pid4; + p.pos() = {1., 1., 1.}; + p.image_box() = {1, -1, 0}; + if (rank == 0) { + particle_range.emplace_back(p); + } + for (auto const use_folded_positions : {true, false}) { + auto const vec = Observables::detail::get_all_particle_positions( + comm, particle_range, pids, {}, use_folded_positions); + if (rank == 0) { + BOOST_CHECK_EQUAL(vec.size(), 4ul); + for (std::size_t i = 0u; i < pids.size(); ++i) { + Utils::Vector3d dist{}; + if (pids[i] == pid4) { + dist = p.pos() - vec[i]; + if (not use_folded_positions) { + dist += p.image_box() * box_l; + } + } else { + dist = start_positions.at(pids[i]) - vec[i]; + } + BOOST_CHECK_LE(dist.norm(), tol); + } + } else { + BOOST_CHECK_EQUAL(vec.size(), 0ul); + } + } + } + // check accumulators if (n_nodes == 1) { auto const pids = std::vector{pid2}; @@ -133,12 +176,12 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory) { for (int i = 0; i < 5; ++i) { set_particle_v(pid2, {static_cast(i), 0., 0.}); - acc.update(); + acc.update(comm); auto const time_series = acc.time_series(); BOOST_REQUIRE_EQUAL(time_series.size(), i + 1); auto const acc_value = time_series.back(); - auto const obs_value = (*obs)(); + auto const obs_value = (*obs)(comm); auto const &p = get_particle_data(pid2); BOOST_TEST(obs_value == p.v(), boost::test_tools::per_element()); BOOST_TEST(acc_value == p.v(), boost::test_tools::per_element()); diff --git a/src/core/unit_tests/lb_particle_coupling_test.cpp b/src/core/unit_tests/lb_particle_coupling_test.cpp index 5c01c36b99a..25378093921 100644 --- a/src/core/unit_tests/lb_particle_coupling_test.cpp +++ b/src/core/unit_tests/lb_particle_coupling_test.cpp @@ -533,6 +533,7 @@ bool test_lb_domain_mismatch_local() { BOOST_AUTO_TEST_CASE(exceptions) { { + boost::mpi::communicator world; using std::exception; // accessing uninitialized pointers is not allowed BOOST_CHECK_THROW(lb_walberla(), std::runtime_error); @@ -541,7 +542,7 @@ BOOST_AUTO_TEST_CASE(exceptions) { BOOST_CHECK_THROW(LB::get_agrid(), exception); BOOST_CHECK_THROW(LB::get_tau(), exception); BOOST_CHECK_THROW(LB::get_kT(), exception); - BOOST_CHECK_THROW(LB::get_pressure_tensor(), exception); + BOOST_CHECK_THROW(LB::get_pressure_tensor(world), exception); BOOST_CHECK_THROW(LB::get_force_to_be_applied({-10., -10., -10.}), std::runtime_error); // coupling, interpolation, boundaries diff --git a/src/particle_observables/include/particle_observables/algorithms.hpp b/src/particle_observables/include/particle_observables/algorithms.hpp index 1c3f9dd675c..c4e9f7c04f4 100644 --- a/src/particle_observables/include/particle_observables/algorithms.hpp +++ b/src/particle_observables/include/particle_observables/algorithms.hpp @@ -45,7 +45,7 @@ template struct WeightedSum { using particle_type = typename ParticleRange::value_type; using value_op_type = decltype(ValueOp{}(std::declval())); using weight_op_type = decltype(WeightOp{}(std::declval())); - auto func = [](auto sum, auto const &p) { + auto func = [](auto const &sum, auto const &p) { auto const w = WeightOp{}(p); return std::make_pair(sum.first + ValueOp{}(p)*w, sum.second + w); }; @@ -59,14 +59,25 @@ template struct WeightedSum { template struct WeightedSum { template auto operator()(ParticleRange const &particles) const { - return detail::WeightedSum()(particles).first; + auto const ws = detail::WeightedSum()(particles); + + return std::make_pair(ws.first, ws.second); + } + + template auto operator()(T const &acc, T const &val) const { + return std::make_pair(acc.first + val.first, acc.second + val.second); } }; template struct Sum { template auto operator()(ParticleRange const &particles) const { - return detail::WeightedSum()(particles).first; + return std::make_pair( + detail::WeightedSum()(particles).first, 0.0); + } + + template auto operator()(T const &acc, T const &val) const { + return WeightedSum{}.operator()(acc, val); } }; @@ -74,7 +85,14 @@ template struct WeightedAverage { template auto operator()(ParticleRange const &particles) const { auto const ws = detail::WeightedSum()(particles); - return (ws.second) ? ws.first / ws.second : ws.first; + return std::make_pair((ws.second) ? ws.first / ws.second : ws.first, + ws.second); + } + + template auto operator()(T const &acc, T const &val) const { + auto const value = acc.first * acc.second + val.first * val.second; + auto const weight = acc.second + val.second; + return std::make_pair((weight) ? value / weight : value, weight); } }; @@ -83,6 +101,10 @@ template struct Average { auto operator()(ParticleRange const &particles) const { return WeightedAverage()(particles); } + + template auto operator()(T const &acc, T const &val) const { + return WeightedAverage{}.operator()(acc, val); + } }; template struct Map { diff --git a/src/particle_observables/include/particle_observables/observable.hpp b/src/particle_observables/include/particle_observables/observable.hpp index 0712333e846..6e66fd6532a 100644 --- a/src/particle_observables/include/particle_observables/observable.hpp +++ b/src/particle_observables/include/particle_observables/observable.hpp @@ -50,6 +50,7 @@ using Momentum = Product; using AverageMomentum = Average; using CenterOfMassPosition = WeightedAverage; using CenterOfMassVelocity = WeightedAverage; +using Identities = Map; using Forces = Map; using Positions = Map; using Velocities = Map; diff --git a/src/particle_observables/include/particle_observables/properties.hpp b/src/particle_observables/include/particle_observables/properties.hpp index 6ee1721162e..3e3510c536e 100644 --- a/src/particle_observables/include/particle_observables/properties.hpp +++ b/src/particle_observables/include/particle_observables/properties.hpp @@ -39,6 +39,14 @@ template using decay_t = typename decay::type; template using default_traits = traits>; +struct Identity { + template > + decltype(auto) operator()(Particle const &p, + Traits particle_traits = {}) const { + return particle_traits.id(p); + } +}; + struct Force { template > decltype(auto) operator()(Particle const &p, diff --git a/src/particle_observables/tests/algorithms.cpp b/src/particle_observables/tests/algorithms.cpp index 7fabb7c3ef4..bd03599e3e6 100644 --- a/src/particle_observables/tests/algorithms.cpp +++ b/src/particle_observables/tests/algorithms.cpp @@ -52,24 +52,26 @@ T average(std::vector const &value, std::size_t denominator) { BOOST_AUTO_TEST_CASE(algorithms_integer) { std::vector const values{1, 2, 3, 4}; { - auto const res = WeightedAverage()(values); + auto const res = + WeightedAverage()(values).first; BOOST_CHECK_EQUAL(res, Testing::average(values, values.size())); } { auto const res = - WeightedAverage()(values); + WeightedAverage()(values).first; BOOST_CHECK_EQUAL(res, (1 * 2 + 2 * 3 + 3 * 4 + 4 * 5) / 14); } { - auto const res = WeightedSum()(values); + auto const res = + WeightedSum()(values).first; BOOST_CHECK_EQUAL(res, (1 * 2 + 2 * 3 + 3 * 4 + 4 * 5)); } { - auto const res = Average()(values); + auto const res = Average()(values).first; BOOST_CHECK_EQUAL(res, Testing::average(values, values.size())); } { - auto const res = Sum{}(values); + auto const res = Sum{}(values).first; BOOST_CHECK_EQUAL(res, Testing::average(values, 1u)); } { @@ -82,15 +84,16 @@ BOOST_AUTO_TEST_CASE(algorithms_double) { auto constexpr tol = 8. * 100. * std::numeric_limits::epsilon(); std::vector const values{1., 2., 3., 4.}; { - auto const res = WeightedAverage()(values); + auto const res = + WeightedAverage()(values).first; BOOST_CHECK_CLOSE(res, Testing::average(values, values.size()), tol); } { - auto const res = Average()(values); + auto const res = Average()(values).first; BOOST_CHECK_EQUAL(res, Testing::average(values, values.size())); } { - auto const res = Sum{}(values); + auto const res = Sum{}(values).first; BOOST_CHECK_EQUAL(res, Testing::average(values, 1u)); } { diff --git a/src/particle_observables/tests/observables.cpp b/src/particle_observables/tests/observables.cpp index 69ce8eb179c..ec56b4a28f8 100644 --- a/src/particle_observables/tests/observables.cpp +++ b/src/particle_observables/tests/observables.cpp @@ -52,17 +52,17 @@ BOOST_AUTO_TEST_CASE(obs) { std::vector parts{p, Testing::Particle{}}; parts[1].mass = 5.; { - auto const res = AverageMomentum{}(parts); + auto const res = AverageMomentum{}(parts).first; BOOST_CHECK(res == 0.5 * (Momentum{}(parts[0]) + Momentum{}(parts[1]))); } { - auto const res = CenterOfMassPosition{}(parts); + auto const res = CenterOfMassPosition{}(parts).first; BOOST_CHECK(res == (Mass{}(parts[0]) * Position{}(parts[0]) + Mass{}(parts[1]) * Position{}(parts[1])) / (Mass{}(parts[0]) + Mass{}(parts[1]))); } { - auto const res = CenterOfMassVelocity{}(parts); + auto const res = CenterOfMassVelocity{}(parts).first; BOOST_CHECK(res == (Mass{}(parts[0]) * Velocity{}(parts[0]) + Mass{}(parts[1]) * Velocity{}(parts[1])) / (Mass{}(parts[0]) + Mass{}(parts[1]))); diff --git a/src/python/espressomd/accumulators.py b/src/python/espressomd/accumulators.py index 36c4d7542aa..e43e66a585e 100644 --- a/src/python/espressomd/accumulators.py +++ b/src/python/espressomd/accumulators.py @@ -41,7 +41,7 @@ class MeanVarianceCalculator(ScriptInterfaceHelper): "update", "shape", ) - _so_creation_policy = "LOCAL" + _so_creation_policy = "GLOBAL" def mean(self): """ @@ -90,7 +90,7 @@ class TimeSeries(ScriptInterfaceHelper): "shape", "clear" ) - _so_creation_policy = "LOCAL" + _so_creation_policy = "GLOBAL" def time_series(self): """ @@ -263,7 +263,7 @@ class Correlator(ScriptInterfaceHelper): "update", "shape", "finalize") - _so_creation_policy = "LOCAL" + _so_creation_policy = "GLOBAL" def result(self): """ @@ -307,7 +307,7 @@ class AutoUpdateAccumulators(ScriptObjectList): """ _so_name = "Accumulators::AutoUpdateAccumulators" - _so_creation_policy = "LOCAL" + _so_creation_policy = "GLOBAL" def add(self, accumulator): """ diff --git a/src/python/espressomd/observables.py b/src/python/espressomd/observables.py index 986442e8a5e..48b5614b4c6 100644 --- a/src/python/espressomd/observables.py +++ b/src/python/espressomd/observables.py @@ -32,7 +32,7 @@ class Observable(ScriptInterfaceHelper): """ _so_name = "Observables::Observable" _so_bind_methods = ("shape",) - _so_creation_policy = "LOCAL" + _so_creation_policy = "GLOBAL" def calculate(self): return np.array(self.call_method("calculate")).reshape(self.shape()) diff --git a/src/script_interface/accumulators/Correlator.hpp b/src/script_interface/accumulators/Correlator.hpp index 95d0908ae19..4e798b81dc7 100644 --- a/src/script_interface/accumulators/Correlator.hpp +++ b/src/script_interface/accumulators/Correlator.hpp @@ -67,11 +67,14 @@ class Correlator : public AccumulatorBase { auto const comp1 = get_value_or(args, "compress1", "discard2"); auto const comp2 = get_value_or(args, "compress2", comp1); - m_correlator = std::make_shared( - get_value(args, "tau_lin"), get_value(args, "tau_max"), - get_value(args, "delta_N"), comp1, comp2, - get_value(args, "corr_operation"), m_obs1->observable(), - m_obs2->observable(), get_value_or(args, "args", {})); + ObjectHandle::context()->parallel_try_catch([&]() { + m_correlator = std::make_shared( + get_value(args, "tau_lin"), get_value(args, "tau_max"), + get_value(args, "delta_N"), comp1, comp2, + get_value(args, "corr_operation"), m_obs1->observable(), + m_obs2->observable(), + get_value_or(args, "args", {})); + }); } std::shared_ptr<::Accumulators::Correlator> correlator() { @@ -80,16 +83,32 @@ class Correlator : public AccumulatorBase { Variant do_call_method(std::string const &method, VariantMap const ¶meters) override { - if (method == "update") - correlator()->update(); - if (method == "finalize") - correlator()->finalize(); - if (method == "get_correlation") - return correlator()->get_correlation(); - if (method == "get_lag_times") - return correlator()->get_lag_times(); - if (method == "get_samples_sizes") - return correlator()->get_samples_sizes(); + if (method == "update") { + ObjectHandle::context()->parallel_try_catch( + [&]() { correlator()->update(context()->get_comm()); }); + } + if (method == "finalize") { + ObjectHandle::context()->parallel_try_catch( + [&]() { correlator()->finalize(context()->get_comm()); }); + } + if (method == "get_correlation") { + if (ObjectHandle::context()->is_head_node()) { + return correlator()->get_correlation(); + } + return {}; + } + if (method == "get_lag_times") { + if (ObjectHandle::context()->is_head_node()) { + return correlator()->get_lag_times(); + } + return {}; + } + if (method == "get_samples_sizes") { + if (ObjectHandle::context()->is_head_node()) { + return correlator()->get_samples_sizes(); + } + return {}; + } return AccumulatorBase::call_method(method, parameters); } diff --git a/src/script_interface/accumulators/MeanVarianceCalculator.hpp b/src/script_interface/accumulators/MeanVarianceCalculator.hpp index e7586af9df4..81e1428b82d 100644 --- a/src/script_interface/accumulators/MeanVarianceCalculator.hpp +++ b/src/script_interface/accumulators/MeanVarianceCalculator.hpp @@ -58,8 +58,10 @@ class MeanVarianceCalculator : public AccumulatorBase { Variant do_call_method(std::string const &method, VariantMap const ¶meters) override { - if (method == "update") - mean_variance_calculator()->update(); + if (method == "update") { + ObjectHandle::context()->parallel_try_catch( + [&]() { mean_variance_calculator()->update(context()->get_comm()); }); + } if (method == "mean") return mean_variance_calculator()->mean(); if (method == "variance") diff --git a/src/script_interface/accumulators/TimeSeries.hpp b/src/script_interface/accumulators/TimeSeries.hpp index 44d9c768fd5..1c113d68856 100644 --- a/src/script_interface/accumulators/TimeSeries.hpp +++ b/src/script_interface/accumulators/TimeSeries.hpp @@ -49,7 +49,8 @@ class TimeSeries : public AccumulatorBase { Variant do_call_method(std::string const &method, VariantMap const ¶meters) override { if (method == "update") { - m_accumulator->update(); + ObjectHandle::context()->parallel_try_catch( + [&]() { m_accumulator->update(context()->get_comm()); }); } if (method == "time_series") { auto const &series = m_accumulator->time_series(); diff --git a/src/script_interface/analysis/Analysis.cpp b/src/script_interface/analysis/Analysis.cpp index 19763b20166..d7f9f20f530 100644 --- a/src/script_interface/analysis/Analysis.cpp +++ b/src/script_interface/analysis/Analysis.cpp @@ -109,6 +109,12 @@ Variant Analysis::do_call_method(std::string const &name, }); return make_unordered_map_of_variants(dict); } +#ifdef DPD + if (name == "dpd_stress") { + auto const result = dpd_stress(context()->get_comm()); + return result.as_vector(); + } +#endif // DPD if (not context()->is_head_node()) { return {}; } @@ -139,12 +145,6 @@ Variant Analysis::do_call_method(std::string const &name, auto const result = nbhood(partCfg(), pos, radius); return result; } -#ifdef DPD - if (name == "dpd_stress") { - auto const result = dpd_stress(); - return result.as_vector(); - } -#endif // DPD 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 index fef69d3f296..5c96ac36222 100644 --- a/src/script_interface/communication.hpp +++ b/src/script_interface/communication.hpp @@ -19,10 +19,8 @@ #ifndef ESPRESSO_SRC_SCRIPT_INTERFACE_COMMUNICATION_HPP #define ESPRESSO_SRC_SCRIPT_INTERFACE_COMMUNICATION_HPP -#include #include #include -#include #include #include @@ -40,30 +38,6 @@ T mpi_reduce_sum(boost::mpi::communicator const &comm, T const &result) { return out; } -/** - * @brief Reduce an optional on the head node. - * Worker nodes get a default-constructed object. - */ -template -T mpi_reduce_optional(boost::mpi::communicator const &comm, - boost::optional const &result) { - assert(1 == boost::mpi::all_reduce(comm, static_cast(!!result), - std::plus<>()) && - "Incorrect number of return values"); - if (comm.rank() == 0) { - if (result) { - return *result; - } - T value; - comm.recv(boost::mpi::any_source, 42, value); - return value; - } - if (result) { - comm.send(0, 42, *result); - } - return {}; -} - } // namespace ScriptInterface #endif diff --git a/src/script_interface/observables/CylindricalLBProfileObservable.hpp b/src/script_interface/observables/CylindricalLBProfileObservable.hpp index 3d5135ba10c..455e7f061a4 100644 --- a/src/script_interface/observables/CylindricalLBProfileObservable.hpp +++ b/src/script_interface/observables/CylindricalLBProfileObservable.hpp @@ -105,19 +105,22 @@ class CylindricalLBProfileObservable void do_construct(VariantMap const ¶ms) override { set_from_args(m_transform_params, params, "transform_params"); - if (m_transform_params) - m_observable = std::make_shared( - m_transform_params->cyl_transform_params(), - get_value_or(params, "n_r_bins", 1), - get_value_or(params, "n_phi_bins", 1), - get_value_or(params, "n_z_bins", 1), - get_value_or(params, "min_r", 0.), - get_value(params, "max_r"), - get_value_or(params, "min_phi", -Utils::pi()), - get_value_or(params, "max_phi", Utils::pi()), - get_value(params, "min_z"), - get_value(params, "max_z"), - get_value(params, "sampling_density")); + if (m_transform_params) { + ObjectHandle::context()->parallel_try_catch([&]() { + m_observable = std::make_shared( + m_transform_params->cyl_transform_params(), + get_value_or(params, "n_r_bins", 1), + get_value_or(params, "n_phi_bins", 1), + get_value_or(params, "n_z_bins", 1), + get_value_or(params, "min_r", 0.), + get_value(params, "max_r"), + get_value_or(params, "min_phi", -Utils::pi()), + get_value_or(params, "max_phi", Utils::pi()), + get_value(params, "min_z"), + get_value(params, "max_z"), + get_value(params, "sampling_density")); + }); + } } Variant do_call_method(std::string const &method, diff --git a/src/script_interface/observables/CylindricalPidProfileObservable.hpp b/src/script_interface/observables/CylindricalPidProfileObservable.hpp index b8823491e5a..629a491bf63 100644 --- a/src/script_interface/observables/CylindricalPidProfileObservable.hpp +++ b/src/script_interface/observables/CylindricalPidProfileObservable.hpp @@ -104,19 +104,22 @@ class CylindricalPidProfileObservable void do_construct(VariantMap const ¶ms) override { set_from_args(m_transform_params, params, "transform_params"); - if (m_transform_params) - m_observable = std::make_shared( - get_value>(params, "ids"), - m_transform_params->cyl_transform_params(), - get_value_or(params, "n_r_bins", 1), - get_value_or(params, "n_phi_bins", 1), - get_value_or(params, "n_z_bins", 1), - get_value_or(params, "min_r", 0.), - get_value(params, "max_r"), - get_value_or(params, "min_phi", -Utils::pi()), - get_value_or(params, "max_phi", Utils::pi()), - get_value(params, "min_z"), - get_value(params, "max_z")); + if (m_transform_params) { + ObjectHandle::context()->parallel_try_catch([&]() { + m_observable = std::make_shared( + get_value>(params, "ids"), + m_transform_params->cyl_transform_params(), + get_value_or(params, "n_r_bins", 1), + get_value_or(params, "n_phi_bins", 1), + get_value_or(params, "n_z_bins", 1), + get_value_or(params, "min_r", 0.), + get_value(params, "max_r"), + get_value_or(params, "min_phi", -Utils::pi()), + get_value_or(params, "max_phi", Utils::pi()), + get_value(params, "min_z"), + get_value(params, "max_z")); + }); + } } Variant do_call_method(std::string const &method, diff --git a/src/script_interface/observables/LBProfileObservable.hpp b/src/script_interface/observables/LBProfileObservable.hpp index eb648da959a..3ebe9f02780 100644 --- a/src/script_interface/observables/LBProfileObservable.hpp +++ b/src/script_interface/observables/LBProfileObservable.hpp @@ -91,14 +91,16 @@ class LBProfileObservable } void do_construct(VariantMap const ¶ms) override { - m_observable = - make_shared_from_args( - params, "sampling_delta_x", "sampling_delta_y", "sampling_delta_z", - "sampling_offset_x", "sampling_offset_y", "sampling_offset_z", - "n_x_bins", "n_y_bins", "n_z_bins", "min_x", "max_x", "min_y", - "max_y", "min_z", "max_z", "allow_empty_bins"); + ObjectHandle::context()->parallel_try_catch([&]() { + m_observable = + make_shared_from_args( + params, "sampling_delta_x", "sampling_delta_y", + "sampling_delta_z", "sampling_offset_x", "sampling_offset_y", + "sampling_offset_z", "n_x_bins", "n_y_bins", "n_z_bins", "min_x", + "max_x", "min_y", "max_y", "min_z", "max_z", "allow_empty_bins"); + }); } Variant do_call_method(std::string const &method, diff --git a/src/script_interface/observables/Observable.hpp b/src/script_interface/observables/Observable.hpp index fa27e60b148..238b5dcc258 100644 --- a/src/script_interface/observables/Observable.hpp +++ b/src/script_interface/observables/Observable.hpp @@ -41,7 +41,11 @@ class Observable : public ObjectHandle { Variant do_call_method(std::string const &method, VariantMap const ¶meters) override { if (method == "calculate") { - return observable()->operator()(); + std::vector out{}; + context()->parallel_try_catch([this, &out]() { + out = observable()->operator()(context()->get_comm()); + }); + return out; } if (method == "shape") { auto const shape = observable()->shape(); diff --git a/src/script_interface/observables/PidObservable.hpp b/src/script_interface/observables/PidObservable.hpp index 57b0cc7a4d3..7431d171602 100644 --- a/src/script_interface/observables/PidObservable.hpp +++ b/src/script_interface/observables/PidObservable.hpp @@ -52,8 +52,10 @@ class PidObservable } void do_construct(VariantMap const ¶ms) override { - m_observable = - make_shared_from_args>(params, "ids"); + ObjectHandle::context()->parallel_try_catch([&]() { + m_observable = + make_shared_from_args>(params, "ids"); + }); } std::shared_ptr<::Observables::Observable> observable() const override { diff --git a/src/script_interface/observables/PidProfileObservable.hpp b/src/script_interface/observables/PidProfileObservable.hpp index 8794d253479..4e6c1ebe947 100644 --- a/src/script_interface/observables/PidProfileObservable.hpp +++ b/src/script_interface/observables/PidProfileObservable.hpp @@ -80,11 +80,13 @@ class PidProfileObservable } void do_construct(VariantMap const ¶ms) override { - m_observable = - make_shared_from_args, int, int, int, double, - double, double, double, double, double>( - params, "ids", "n_x_bins", "n_y_bins", "n_z_bins", "min_x", "max_x", - "min_y", "max_y", "min_z", "max_z"); + ObjectHandle::context()->parallel_try_catch([&]() { + m_observable = + make_shared_from_args, int, int, int, + double, double, double, double, double, double>( + params, "ids", "n_x_bins", "n_y_bins", "n_z_bins", "min_x", + "max_x", "min_y", "max_y", "min_z", "max_z"); + }); } Variant do_call_method(std::string const &method, diff --git a/src/script_interface/observables/RDF.hpp b/src/script_interface/observables/RDF.hpp index 31c1e8b17c2..8d912c3ad7c 100644 --- a/src/script_interface/observables/RDF.hpp +++ b/src/script_interface/observables/RDF.hpp @@ -55,9 +55,12 @@ class RDF : public AutoParameters { } void do_construct(VariantMap const ¶ms) override { - m_observable = make_shared_from_args<::Observables::RDF, std::vector, - std::vector, int, double, double>( - params, "ids1", "ids2", "n_r_bins", "min_r", "max_r"); + context()->parallel_try_catch([&]() { + m_observable = + make_shared_from_args<::Observables::RDF, std::vector, + std::vector, int, double, double>( + params, "ids1", "ids2", "n_r_bins", "min_r", "max_r"); + }); } std::shared_ptr<::Observables::RDF> rdf_observable() const { diff --git a/src/script_interface/particle_data/ParticleHandle.cpp b/src/script_interface/particle_data/ParticleHandle.cpp index 3e7ad5a0df4..2a89b544c24 100644 --- a/src/script_interface/particle_data/ParticleHandle.cpp +++ b/src/script_interface/particle_data/ParticleHandle.cpp @@ -22,7 +22,6 @@ #include "ParticleHandle.hpp" #include "script_interface/Variant.hpp" -#include "script_interface/communication.hpp" #include "script_interface/get_value.hpp" #include "core/bonded_interactions/bonded_interaction_data.hpp" @@ -36,6 +35,7 @@ #include "core/virtual_sites.hpp" #include +#include #include #include @@ -140,7 +140,7 @@ T ParticleHandle::get_particle_property(F const &fun) const { } else { ret = {fun(*ptr)}; } - return mpi_reduce_optional(comm, ret); + return Utils::Mpi::reduce_optional(comm, ret); } template diff --git a/src/script_interface/tests/Accumulators_test.cpp b/src/script_interface/tests/Accumulators_test.cpp index 0a35c371c8f..f5dafa478ab 100644 --- a/src/script_interface/tests/Accumulators_test.cpp +++ b/src/script_interface/tests/Accumulators_test.cpp @@ -17,12 +17,15 @@ * along with this program. If not, see . */ +#define BOOST_TEST_NO_MAIN #define BOOST_TEST_MODULE Accumulators test #define BOOST_TEST_DYN_LINK #include #include +#include "script_interface/GlobalContext.hpp" + #include "script_interface/accumulators/Correlator.hpp" #include "script_interface/accumulators/MeanVarianceCalculator.hpp" #include "script_interface/accumulators/TimeSeries.hpp" @@ -31,6 +34,12 @@ #include "core/observables/Observable.hpp" +#include "core/MpiCallbacks.hpp" +#include "core/errorhandling.hpp" + +#include +#include + #include #include #include @@ -40,7 +49,10 @@ namespace Observables { class MockObservable : public Observable { public: - std::vector operator()() const override { return {1., 2., 3., 4.}; } + std::vector + operator()(boost::mpi::communicator const &comm) const override { + return {1., 2., 3., 4.}; + } std::vector shape() const override { return {2u, 2u}; } }; } // namespace Observables @@ -53,130 +65,202 @@ NEW_PARAMLESS_OBSERVABLE(MockObservable) using TestObs = ScriptInterface::Observables::MockObservable; using TestObsPtr = std::shared_ptr; +using TimeSeries = ScriptInterface::Accumulators::TimeSeries; +using Correlator = ScriptInterface::Accumulators::Correlator; +using MeanVarianceCalculator = + ScriptInterface::Accumulators::MeanVarianceCalculator; + using namespace ScriptInterface; +auto make_global_context(boost::mpi::communicator const &comm, + Communication::MpiCallbacks &cb) { + Utils::Factory factory; + factory.register_new("TestObs"); + factory.register_new("TimeSeries"); + factory.register_new("Correlator"); + factory.register_new("MeanVarianceCalculator"); + + return std::make_shared( + cb, std::make_shared(factory, comm)); +} + BOOST_AUTO_TEST_CASE(time_series) { - auto const obs = std::make_shared(); - ScriptInterface::Accumulators::TimeSeries acc; - acc.do_construct({{"obs", obs}, {"delta_N", 2}}); - acc.do_call_method("update", VariantMap{}); - { - BOOST_CHECK_EQUAL(get_value(acc.get_parameter("delta_N")), 2); - BOOST_CHECK_EQUAL(get_value(acc.get_parameter("obs")), obs); - } - { - auto const shape_ref = std::vector{{1, 2, 2}}; - // check non-const access - auto const variant = acc.do_call_method("shape", VariantMap{}); - auto const shape = get_value>(variant); - BOOST_TEST(shape == shape_ref, boost::test_tools::per_element()); - // check const access - auto const shape_const = std::as_const(acc).accumulator()->shape(); - BOOST_TEST(shape_const == shape_ref, boost::test_tools::per_element()); - } - { - auto const variant = acc.do_call_method("time_series", VariantMap{}); - auto const time_series = get_value>(variant); - BOOST_REQUIRE_EQUAL(time_series.size(), 1u); - auto const series = get_value>(time_series[0]); - auto const series_ref = std::vector{1., 2., 3., 4.}; - BOOST_TEST(series == series_ref, boost::test_tools::per_element()); + boost::mpi::communicator world; + Communication::MpiCallbacks cb(world); + auto ctx = make_global_context(world, cb); + + if (world.rank() != 0) { + cb.loop(); + } else { + auto const obs = ctx->make_shared("TestObs", {}); + BOOST_REQUIRE(obs != nullptr); + BOOST_CHECK_EQUAL(obs->context(), ctx.get()); + + auto const acc_si = + ctx->make_shared("TimeSeries", {{"obs", obs}, {"delta_N", 2}}); + BOOST_REQUIRE(acc_si != nullptr); + auto const acc = std::dynamic_pointer_cast(acc_si); + BOOST_REQUIRE(acc != nullptr); + + acc_si->call_method("update", VariantMap{}); + { + BOOST_CHECK_EQUAL(get_value(acc_si->get_parameter("delta_N")), 2); + BOOST_CHECK_EQUAL(get_value(acc_si->get_parameter("obs")), + obs); + } + { + auto const shape_ref = std::vector{{1, 2, 2}}; + // check non-const access + auto const variant = acc_si->call_method("shape", VariantMap{}); + auto const shape = get_value>(variant); + BOOST_TEST(shape == shape_ref, boost::test_tools::per_element()); + // check const access + auto const shape_const = std::as_const(*acc).accumulator()->shape(); + BOOST_TEST(shape_const == shape_ref, boost::test_tools::per_element()); + } + { + auto const variant = acc_si->call_method("time_series", VariantMap{}); + auto const time_series = get_value>(variant); + BOOST_REQUIRE_EQUAL(time_series.size(), 1u); + auto const series = get_value>(time_series[0]); + auto const series_ref = std::vector{1., 2., 3., 4.}; + BOOST_TEST(series == series_ref, boost::test_tools::per_element()); + } } } BOOST_AUTO_TEST_CASE(correlator) { - auto const obs = std::make_shared(); - ScriptInterface::Accumulators::Correlator acc; - acc.do_construct({{"obs1", obs}, - {"delta_N", 2}, - {"tau_lin", 4}, - {"tau_max", 2.}, - {"corr_operation", std::string("componentwise_product")}}); - acc.do_call_method("update", VariantMap{}); - acc.do_call_method("finalize", VariantMap{}); - { - BOOST_CHECK_EQUAL(get_value(acc.get_parameter("delta_N")), 2); - BOOST_CHECK_EQUAL(get_value(acc.get_parameter("tau_lin")), 4); - BOOST_CHECK_EQUAL(get_value(acc.get_parameter("tau_max")), 2.); - BOOST_CHECK_EQUAL( - get_value(acc.get_parameter("corr_operation")), - std::string("componentwise_product")); - BOOST_CHECK_EQUAL(get_value(acc.get_parameter("obs1")), obs); - BOOST_CHECK_EQUAL(get_value(acc.get_parameter("obs2")), obs); - } - { - auto const shape_ref = std::vector{{5, 2, 2}}; - // check non-const access - auto const variant = acc.do_call_method("shape", VariantMap{}); - auto const shape = get_value>(variant); - BOOST_TEST(shape == shape_ref, boost::test_tools::per_element()); - // check const access - auto const shape_const = std::as_const(acc).accumulator()->shape(); - BOOST_TEST(shape_const == shape_ref, boost::test_tools::per_element()); - } - { - auto const variant = acc.do_call_method("get_correlation", VariantMap{}); - auto const correlation = get_value>(variant); - BOOST_REQUIRE_EQUAL(correlation.size(), 20u); - auto corr_ref = std::vector{1., 4., 9., 16.}; - corr_ref.resize(correlation.size()); - BOOST_TEST(correlation == corr_ref, boost::test_tools::per_element()); - } - { - auto const variant = acc.do_call_method("get_lag_times", VariantMap{}); - auto const lag_times = get_value>(variant); - BOOST_REQUIRE_EQUAL(lag_times.size(), 5u); - auto const lag_times_ref = std::vector{0., -2., -4., -6., -8.}; - BOOST_TEST(lag_times == lag_times_ref, boost::test_tools::per_element()); - } - { - auto const variant = acc.do_call_method("get_samples_sizes", VariantMap{}); - auto const samples_n = get_value>(variant); - BOOST_REQUIRE_EQUAL(samples_n.size(), 5u); - auto const samples_n_ref = std::vector{1, 0, 0, 0, 0}; - BOOST_TEST(samples_n == samples_n_ref, boost::test_tools::per_element()); + boost::mpi::communicator world; + Communication::MpiCallbacks cb(world); + auto ctx = make_global_context(world, cb); + + if (world.rank() != 0) { + cb.loop(); + } else { + auto const obs = ctx->make_shared("TestObs", {}); + BOOST_REQUIRE(obs != nullptr); + + auto const acc_si = ctx->make_shared( + "Correlator", + {{"obs1", obs}, + {"delta_N", 2}, + {"tau_lin", 4}, + {"tau_max", 2.}, + {"corr_operation", std::string("componentwise_product")}}); + BOOST_REQUIRE(acc_si != nullptr); + auto acc = std::dynamic_pointer_cast(acc_si); + BOOST_REQUIRE(acc != nullptr); + + acc_si->call_method("update", VariantMap{}); + acc_si->call_method("finalize", VariantMap{}); + { + BOOST_CHECK_EQUAL(get_value(acc_si->get_parameter("delta_N")), 2); + BOOST_CHECK_EQUAL(get_value(acc_si->get_parameter("tau_lin")), 4); + BOOST_CHECK_EQUAL(get_value(acc_si->get_parameter("tau_max")), + 2.); + BOOST_CHECK_EQUAL( + get_value(acc_si->get_parameter("corr_operation")), + std::string("componentwise_product")); + BOOST_CHECK_EQUAL(get_value(acc_si->get_parameter("obs1")), + obs); + BOOST_CHECK_EQUAL(get_value(acc_si->get_parameter("obs2")), + obs); + } + { + auto const shape_ref = std::vector{{5, 2, 2}}; + // check non-const access + auto const variant = acc_si->call_method("shape", VariantMap{}); + auto const shape = get_value>(variant); + BOOST_TEST(shape == shape_ref, boost::test_tools::per_element()); + // check const access + auto const shape_const = std::as_const(*acc).accumulator()->shape(); + BOOST_TEST(shape_const == shape_ref, boost::test_tools::per_element()); + } + { + auto const variant = acc_si->call_method("get_correlation", VariantMap{}); + auto const correlation = get_value>(variant); + BOOST_REQUIRE_EQUAL(correlation.size(), 20u); + auto corr_ref = std::vector{1., 4., 9., 16.}; + corr_ref.resize(correlation.size()); + BOOST_TEST(correlation == corr_ref, boost::test_tools::per_element()); + } + { + auto const variant = acc_si->call_method("get_lag_times", VariantMap{}); + auto const lag_times = get_value>(variant); + BOOST_REQUIRE_EQUAL(lag_times.size(), 5u); + auto const lag_times_ref = std::vector{0., -2., -4., -6., -8.}; + BOOST_TEST(lag_times == lag_times_ref, boost::test_tools::per_element()); + } + { + auto const variant = + acc_si->call_method("get_samples_sizes", VariantMap{}); + auto const samples_n = get_value>(variant); + BOOST_REQUIRE_EQUAL(samples_n.size(), 5u); + auto const samples_n_ref = std::vector{1, 0, 0, 0, 0}; + BOOST_TEST(samples_n == samples_n_ref, boost::test_tools::per_element()); + } } } BOOST_AUTO_TEST_CASE(mean_variance) { - auto const obs = std::make_shared(); - ScriptInterface::Accumulators::MeanVarianceCalculator acc; - acc.do_construct({{"obs", obs}, {"delta_N", 2}}); - acc.do_call_method("update", VariantMap{}); - acc.do_call_method("update", VariantMap{}); - { - BOOST_CHECK_EQUAL(get_value(acc.get_parameter("delta_N")), 2); - BOOST_CHECK_EQUAL(get_value(acc.get_parameter("obs")), obs); - } - { - auto const shape_ref = std::vector{{2, 2}}; - // check non-const access - auto const variant = acc.do_call_method("shape", VariantMap{}); - auto const shape = get_value>(variant); - BOOST_TEST(shape == shape_ref, boost::test_tools::per_element()); - // check const access - auto const shape_const = std::as_const(acc).accumulator()->shape(); - BOOST_TEST(shape_const == shape_ref, boost::test_tools::per_element()); - } - { - auto const variant = acc.do_call_method("mean", VariantMap{}); - auto const mean = get_value>(variant); - BOOST_REQUIRE_EQUAL(mean.size(), 4u); - auto const mean_ref = std::vector{1., 2., 3., 4.}; - BOOST_TEST(mean == mean_ref, boost::test_tools::per_element()); - } - { - auto const variant = acc.do_call_method("variance", VariantMap{}); - auto const variance = get_value>(variant); - BOOST_REQUIRE_EQUAL(variance.size(), 4u); - auto const variance_ref = std::vector{0., 0., 0., 0.}; - BOOST_TEST(variance == variance_ref, boost::test_tools::per_element()); - } - { - auto const variant = acc.do_call_method("std_error", VariantMap{}); - auto const stderror = get_value>(variant); - BOOST_REQUIRE_EQUAL(stderror.size(), 4u); - auto const stderror_ref = std::vector{0., 0., 0., 0.}; - BOOST_TEST(stderror == stderror_ref, boost::test_tools::per_element()); + boost::mpi::communicator world; + Communication::MpiCallbacks cb(world); + auto ctx = make_global_context(world, cb); + + if (world.rank() != 0) { + cb.loop(); + } else { + auto const obs = ctx->make_shared("TestObs", {}); + + auto const acc_si = ctx->make_shared("MeanVarianceCalculator", + {{"obs", obs}, {"delta_N", 2}}); + BOOST_REQUIRE(acc_si != nullptr); + auto acc = std::dynamic_pointer_cast(acc_si); + BOOST_REQUIRE(acc != nullptr); + + acc_si->call_method("update", VariantMap{}); + acc_si->call_method("update", VariantMap{}); + { + BOOST_CHECK_EQUAL(get_value(acc_si->get_parameter("delta_N")), 2); + BOOST_CHECK_EQUAL(get_value(acc_si->get_parameter("obs")), + obs); + } + { + auto const shape_ref = std::vector{{2, 2}}; + // check non-const access + auto const variant = acc->call_method("shape", VariantMap{}); + auto const shape = get_value>(variant); + BOOST_TEST(shape == shape_ref, boost::test_tools::per_element()); + // check const access + auto const shape_const = std::as_const(*acc).accumulator()->shape(); + BOOST_TEST(shape_const == shape_ref, boost::test_tools::per_element()); + } + { + auto const variant = acc_si->call_method("mean", VariantMap{}); + auto const mean = get_value>(variant); + BOOST_REQUIRE_EQUAL(mean.size(), 4u); + auto const mean_ref = std::vector{1., 2., 3., 4.}; + BOOST_TEST(mean == mean_ref, boost::test_tools::per_element()); + } + { + auto const variant = acc_si->call_method("variance", VariantMap{}); + auto const variance = get_value>(variant); + BOOST_REQUIRE_EQUAL(variance.size(), 4u); + auto const variance_ref = std::vector{0., 0., 0., 0.}; + BOOST_TEST(variance == variance_ref, boost::test_tools::per_element()); + } + { + auto const variant = acc_si->call_method("std_error", VariantMap{}); + auto const stderror = get_value>(variant); + BOOST_REQUIRE_EQUAL(stderror.size(), 4u); + auto const stderror_ref = std::vector{0., 0., 0., 0.}; + BOOST_TEST(stderror == stderror_ref, boost::test_tools::per_element()); + } } } + +int main(int argc, char **argv) { + boost::mpi::environment mpi_env(argc, argv); + + return boost::unit_test::unit_test_main(init_unit_test, argc, argv); +} diff --git a/src/script_interface/tests/CMakeLists.txt b/src/script_interface/tests/CMakeLists.txt index 0fc25729ec1..399e44255a7 100644 --- a/src/script_interface/tests/CMakeLists.txt +++ b/src/script_interface/tests/CMakeLists.txt @@ -48,7 +48,8 @@ unit_test(NAME ObjectList_test SRC ObjectList_test.cpp DEPENDS unit_test(NAME ObjectMap_test SRC ObjectMap_test.cpp DEPENDS espresso::script_interface espresso::core Boost::mpi) unit_test(NAME Accumulators_test SRC Accumulators_test.cpp DEPENDS - espresso::script_interface espresso::core) + espresso::script_interface espresso::core Boost::mpi MPI::MPI_CXX + NUM_PROC 2) unit_test(NAME Constraints_test SRC Constraints_test.cpp DEPENDS espresso::script_interface espresso::core) unit_test(NAME Actors_test SRC Actors_test.cpp DEPENDS diff --git a/src/script_interface/tests/GlobalContext_test.cpp b/src/script_interface/tests/GlobalContext_test.cpp index b5556b2d18e..a1ef31d0a8d 100644 --- a/src/script_interface/tests/GlobalContext_test.cpp +++ b/src/script_interface/tests/GlobalContext_test.cpp @@ -54,10 +54,10 @@ struct Dummy : si::ObjectHandle { } }; -auto make_global_context(Communication::MpiCallbacks &cb) { +auto make_global_context(boost::mpi::communicator const &comm, + Communication::MpiCallbacks &cb) { Utils::Factory factory; factory.register_new("Dummy"); - boost::mpi::communicator comm; return std::make_shared( cb, std::make_shared(factory, comm)); @@ -66,7 +66,7 @@ auto make_global_context(Communication::MpiCallbacks &cb) { BOOST_AUTO_TEST_CASE(GlobalContext_make_shared) { boost::mpi::communicator world; Communication::MpiCallbacks cb{world}; - auto ctx = make_global_context(cb); + auto ctx = make_global_context(world, cb); if (world.rank() == 0) { auto res = ctx->make_shared("Dummy", {}); @@ -81,7 +81,7 @@ BOOST_AUTO_TEST_CASE(GlobalContext_make_shared) { BOOST_AUTO_TEST_CASE(GlobalContext_serialization) { boost::mpi::communicator world; Communication::MpiCallbacks cb{world}; - auto ctx = make_global_context(cb); + auto ctx = make_global_context(world, cb); if (world.rank() == 0) { auto const serialized = [&]() { diff --git a/src/script_interface/tests/reduction_test.cpp b/src/script_interface/tests/reduction_test.cpp index 24b8354c4ff..d5613e440f4 100644 --- a/src/script_interface/tests/reduction_test.cpp +++ b/src/script_interface/tests/reduction_test.cpp @@ -25,7 +25,6 @@ #include "script_interface/communication.hpp" #include -#include BOOST_AUTO_TEST_CASE(reduce_sum) { boost::mpi::communicator comm; @@ -39,24 +38,6 @@ BOOST_AUTO_TEST_CASE(reduce_sum) { } } -BOOST_AUTO_TEST_CASE(reduce_optional) { - boost::mpi::communicator comm; - - for (int rank = 0; rank < comm.size(); ++rank) { - boost::optional maybe; - if (comm.rank() == rank) { - maybe = 42; - } - auto const sum = ScriptInterface::mpi_reduce_optional(comm, maybe); - - if (comm.rank() == 0) { - BOOST_CHECK_EQUAL(sum, 42); - } else { - BOOST_CHECK_EQUAL(sum, 0); - } - } -} - int main(int argc, char **argv) { boost::mpi::environment mpi_env(argc, argv); diff --git a/src/script_interface/walberla/EKReaction.hpp b/src/script_interface/walberla/EKReaction.hpp index 72ba11569ba..3adc5f88479 100644 --- a/src/script_interface/walberla/EKReaction.hpp +++ b/src/script_interface/walberla/EKReaction.hpp @@ -33,7 +33,8 @@ #include #include -#include + +#include #include #include @@ -163,7 +164,7 @@ class EKIndexedReaction : public EKReaction { get_value(parameters, "node"), get_instance()->get_lattice()->get_grid_dimensions()); auto const result = m_ekreaction_impl->get_node_is_boundary(index); - return mpi_reduce_optional(context()->get_comm(), result); + return Utils::Mpi::reduce_optional(context()->get_comm(), result); } return {}; } diff --git a/src/script_interface/walberla/EKSpeciesNode.cpp b/src/script_interface/walberla/EKSpeciesNode.cpp index 594f1a73cd0..b367deda65a 100644 --- a/src/script_interface/walberla/EKSpeciesNode.cpp +++ b/src/script_interface/walberla/EKSpeciesNode.cpp @@ -25,12 +25,11 @@ #include "LatticeIndices.hpp" -#include - #include #include #include +#include #include #include @@ -69,18 +68,20 @@ Variant EKSpeciesNode::do_call_method(std::string const &name, } if (name == "get_density") { auto const result = m_ek_species->get_node_density(m_index); - return mpi_reduce_optional(context()->get_comm(), result) / m_conv_dens; + return Utils::Mpi::reduce_optional(context()->get_comm(), result) / + m_conv_dens; } if (name == "get_is_boundary") { auto const result = m_ek_species->get_node_is_boundary(m_index); - return mpi_reduce_optional(context()->get_comm(), result); + return Utils::Mpi::reduce_optional(context()->get_comm(), result); } if (name == "get_node_density_at_boundary") { auto const boundary_opt = m_ek_species->get_node_is_density_boundary(m_index); if (is_boundary_all_reduce(context()->get_comm(), boundary_opt)) { auto const result = m_ek_species->get_node_density_at_boundary(m_index); - return mpi_reduce_optional(context()->get_comm(), result) / m_conv_dens; + return Utils::Mpi::reduce_optional(context()->get_comm(), result) / + m_conv_dens; } return Variant{None{}}; } @@ -97,7 +98,8 @@ Variant EKSpeciesNode::do_call_method(std::string const &name, auto const boundary_opt = m_ek_species->get_node_is_flux_boundary(m_index); if (is_boundary_all_reduce(context()->get_comm(), boundary_opt)) { auto const result = m_ek_species->get_node_flux_at_boundary(m_index); - return mpi_reduce_optional(context()->get_comm(), result) / m_conv_flux; + return Utils::Mpi::reduce_optional(context()->get_comm(), result) / + m_conv_flux; } return Variant{None{}}; } diff --git a/src/script_interface/walberla/LBFluid.cpp b/src/script_interface/walberla/LBFluid.cpp index 42e6ecbb4c5..c70c3f82a39 100644 --- a/src/script_interface/walberla/LBFluid.cpp +++ b/src/script_interface/walberla/LBFluid.cpp @@ -40,6 +40,7 @@ #include #include +#include #include #include @@ -208,7 +209,8 @@ std::vector LBFluid::get_average_pressure_tensor() const { Variant LBFluid::get_interpolated_velocity(Utils::Vector3d const &pos) const { auto const lb_pos = folded_position(pos, box_geo) * m_conv_dist; auto const result = m_instance->get_velocity_at_pos(lb_pos); - return mpi_reduce_optional(context()->get_comm(), result) / m_conv_speed; + return Utils::Mpi::reduce_optional(context()->get_comm(), result) / + m_conv_speed; } void LBFluid::load_checkpoint(std::string const &filename, int mode) { diff --git a/src/script_interface/walberla/LBFluidNode.cpp b/src/script_interface/walberla/LBFluidNode.cpp index 600fa9a74b1..aeb41837e6f 100644 --- a/src/script_interface/walberla/LBFluidNode.cpp +++ b/src/script_interface/walberla/LBFluidNode.cpp @@ -23,10 +23,9 @@ #include "LBFluidNode.hpp" -#include - #include #include +#include #include #include @@ -61,14 +60,15 @@ Variant LBFluidNode::do_call_method(std::string const &name, auto const boundary_opt = m_lb_fluid->get_node_is_boundary(m_index); if (is_boundary_all_reduce(context()->get_comm(), boundary_opt)) { auto const result = m_lb_fluid->get_node_velocity_at_boundary(m_index); - return mpi_reduce_optional(context()->get_comm(), result) / + return Utils::Mpi::reduce_optional(context()->get_comm(), result) / m_conv_velocity; } return Variant{None{}}; } if (name == "get_density") { auto const result = m_lb_fluid->get_node_density(m_index); - return mpi_reduce_optional(context()->get_comm(), result) / m_conv_dens; + return Utils::Mpi::reduce_optional(context()->get_comm(), result) / + m_conv_dens; } if (name == "set_density") { auto const dens = get_value(params, "value"); @@ -78,7 +78,7 @@ Variant LBFluidNode::do_call_method(std::string const &name, } if (name == "get_population") { auto const result = m_lb_fluid->get_node_population(m_index); - return mpi_reduce_optional(context()->get_comm(), result); + return Utils::Mpi::reduce_optional(context()->get_comm(), result); } if (name == "set_population") { auto const pop = get_value>(params, "value"); @@ -88,7 +88,8 @@ Variant LBFluidNode::do_call_method(std::string const &name, } if (name == "get_velocity") { auto const result = m_lb_fluid->get_node_velocity(m_index); - return mpi_reduce_optional(context()->get_comm(), result) / m_conv_velocity; + return Utils::Mpi::reduce_optional(context()->get_comm(), result) / + m_conv_velocity; } if (name == "set_velocity") { auto const u = @@ -99,13 +100,14 @@ Variant LBFluidNode::do_call_method(std::string const &name, } if (name == "get_is_boundary") { auto const result = m_lb_fluid->get_node_is_boundary(m_index); - return mpi_reduce_optional(context()->get_comm(), result); + return Utils::Mpi::reduce_optional(context()->get_comm(), result); } if (name == "get_boundary_force") { auto const boundary_opt = m_lb_fluid->get_node_is_boundary(m_index); if (is_boundary_all_reduce(context()->get_comm(), boundary_opt)) { auto result = m_lb_fluid->get_node_boundary_force(m_index); - return mpi_reduce_optional(context()->get_comm(), result) / m_conv_force; + return Utils::Mpi::reduce_optional(context()->get_comm(), result) / + m_conv_force; } return Variant{None{}}; } @@ -115,7 +117,7 @@ Variant LBFluidNode::do_call_method(std::string const &name, if (result) { value = (*result / m_conv_press).as_vector(); } - auto vec = mpi_reduce_optional(context()->get_comm(), value); + auto vec = Utils::Mpi::reduce_optional(context()->get_comm(), value); if (context()->is_head_node()) { if (name == "get_pressure_tensor_neq") { auto constexpr c_sound_sq = 1. / 3.; @@ -135,7 +137,8 @@ Variant LBFluidNode::do_call_method(std::string const &name, } if (name == "get_last_applied_force") { auto const result = m_lb_fluid->get_node_last_applied_force(m_index); - return mpi_reduce_optional(context()->get_comm(), result) / m_conv_force; + return Utils::Mpi::reduce_optional(context()->get_comm(), result) / + m_conv_force; } if (name == "set_last_applied_force") { auto const f = get_value(params, "value"); diff --git a/src/utils/include/utils/mpi/reduce_optional.hpp b/src/utils/include/utils/mpi/reduce_optional.hpp new file mode 100644 index 00000000000..216a234c364 --- /dev/null +++ b/src/utils/include/utils/mpi/reduce_optional.hpp @@ -0,0 +1,57 @@ +/* + * 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 UTILS_MPI_REDUCE_OPTIONAL_HPP +#define UTILS_MPI_REDUCE_OPTIONAL_HPP + +#include +#include +#include + +#include +#include + +namespace Utils::Mpi { + +/** + * @brief Reduce an optional on the head node. + * Worker nodes get a default-constructed object. + */ +template +T reduce_optional(boost::mpi::communicator const &comm, + boost::optional const &result) { + assert(1 == boost::mpi::all_reduce(comm, static_cast(!!result), + std::plus<>()) && + "Incorrect number of return values"); + if (comm.rank() == 0) { + if (result) { + return *result; + } + T value; + comm.recv(boost::mpi::any_source, 42, value); + return value; + } + if (result) { + comm.send(0, 42, *result); + } + return {}; +} + +} // namespace Utils::Mpi + +#endif diff --git a/src/utils/tests/CMakeLists.txt b/src/utils/tests/CMakeLists.txt index 092d891cfb6..d7935217f6e 100644 --- a/src/utils/tests/CMakeLists.txt +++ b/src/utils/tests/CMakeLists.txt @@ -109,3 +109,5 @@ unit_test(NAME matrix_test SRC matrix_test.cpp DEPENDS espresso::utils Boost::serialization NUM_PROC 1) unit_test(NAME orthonormal_vec_test SRC orthonormal_vec_test.cpp DEPENDS espresso::utils Boost::serialization NUM_PROC 1) +unit_test(NAME reduce_optional_test SRC reduce_optional_test.cpp DEPENDS + espresso::utils::mpi Boost::mpi MPI::MPI_CXX NUM_PROC 4) diff --git a/src/utils/tests/reduce_optional_test.cpp b/src/utils/tests/reduce_optional_test.cpp new file mode 100644 index 00000000000..1a95151d8d8 --- /dev/null +++ b/src/utils/tests/reduce_optional_test.cpp @@ -0,0 +1,52 @@ +/* + * 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 . + */ + +#define BOOST_TEST_NO_MAIN +#define BOOST_TEST_MODULE MPI reduction algorithms +#define BOOST_TEST_DYN_LINK +#include + +#include + +#include +#include + +BOOST_AUTO_TEST_CASE(reduce_optional) { + boost::mpi::communicator comm; + + for (int rank = 0; rank < comm.size(); ++rank) { + boost::optional maybe; + if (comm.rank() == rank) { + maybe = 42; + } + auto const sum = Utils::Mpi::reduce_optional(comm, maybe); + + if (comm.rank() == 0) { + BOOST_CHECK_EQUAL(sum, 42); + } else { + BOOST_CHECK_EQUAL(sum, 0); + } + } +} + +int main(int argc, char **argv) { + boost::mpi::environment mpi_env(argc, argv); + + return boost::unit_test::unit_test_main(init_unit_test, argc, argv); +} diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index 3730c3190a6..f504423711d 100644 --- a/testsuite/python/CMakeLists.txt +++ b/testsuite/python/CMakeLists.txt @@ -300,7 +300,7 @@ python_test(FILE dpd_stats.py MAX_NUM_PROC 4 LABELS long) python_test(FILE hat.py MAX_NUM_PROC 4) python_test(FILE analyze_energy.py MAX_NUM_PROC 2 GPU_SLOTS 1) python_test(FILE analyze_mass_related.py MAX_NUM_PROC 4) -python_test(FILE rdf.py MAX_NUM_PROC 1) +python_test(FILE rdf.py MAX_NUM_PROC 2) python_test(FILE sf_simple_lattice.py MAX_NUM_PROC 1) python_test(FILE coulomb_mixed_periodicity.py MAX_NUM_PROC 4) python_test(FILE coulomb_cloud_wall_duplicated.py MAX_NUM_PROC 4 GPU_SLOTS 3) diff --git a/testsuite/python/observable_profile.py b/testsuite/python/observable_profile.py index 0a7b7eafff5..efa6f8627cb 100644 --- a/testsuite/python/observable_profile.py +++ b/testsuite/python/observable_profile.py @@ -31,6 +31,7 @@ class ProfileObservablesTest(ut.TestCase): system.time_step = 0.01 system.part.add(pos=[4.0, 4.0, 6.0], v=[0.0, 0.0, 1.0]) system.part.add(pos=[4.0, 4.0, 6.0], v=[0.0, 0.0, 1.0]) + system.part.add(pos=[6.0, 8.0, 16.0], v=[0.0, 1.0, 0.0]) bin_volume = 5.0**3 kwargs = {'ids': list(system.part.all().id), 'n_x_bins': 2, @@ -79,7 +80,8 @@ def test_force_density_profile(self): obs_data.shape, [self.kwargs['n_x_bins'], self.kwargs['n_y_bins'], self.kwargs['n_z_bins'], 3]) self.assertEqual(obs_data[0, 0, 1, 2], 2.0 / self.bin_volume) - self.assertEqual(np.sum(np.abs(obs_data)), 2.0 / self.bin_volume) + self.assertEqual(obs_data[1, 1, 3, 2], 1.0 / self.bin_volume) + self.assertEqual(np.sum(np.abs(obs_data)), 3.0 / self.bin_volume) def test_flux_density_profile(self): density_profile = espressomd.observables.FluxDensityProfile( @@ -90,7 +92,8 @@ def test_flux_density_profile(self): obs_data.shape, [self.kwargs['n_x_bins'], self.kwargs['n_y_bins'], self.kwargs['n_z_bins'], 3]) self.assertEqual(obs_data[0, 0, 1, 2], 2.0 / self.bin_volume) - self.assertEqual(np.sum(np.abs(obs_data)), 2.0 / self.bin_volume) + self.assertEqual(obs_data[1, 1, 3, 1], 1.0 / self.bin_volume) + self.assertEqual(np.sum(np.abs(obs_data)), 3.0 / self.bin_volume) def test_pid_profile_interface(self): # test setters and getters