From 27fbfcd1e161ebe5d1330e31c49e9fb31cc70e5e Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Thu, 9 Apr 2020 15:58:23 +0200 Subject: [PATCH 01/27] core: short_range_loop: Split particle kernels Split single particle kernel application from pair kernel application. --- src/core/algorithm/for_each_pair.hpp | 38 +++++++++++++------------- src/core/algorithm/link_cell.hpp | 8 ++---- src/core/algorithm/verlet_ia.hpp | 29 ++++++-------------- src/core/short_range_loop.hpp | 32 ++++++++++------------ src/core/unit_tests/link_cell_test.cpp | 7 ----- src/core/unit_tests/verlet_ia_test.cpp | 21 -------------- 6 files changed, 44 insertions(+), 91 deletions(-) diff --git a/src/core/algorithm/for_each_pair.hpp b/src/core/algorithm/for_each_pair.hpp index 20d35f64b12..f64bf699252 100644 --- a/src/core/algorithm/for_each_pair.hpp +++ b/src/core/algorithm/for_each_pair.hpp @@ -26,15 +26,14 @@ namespace Algorithm { /** - * @brief Run single and pair kernel for each particle (pair) from cell range. + * @brief Run pair kernel for each particle pair from a cell range. * - * Iterates over all cells in [first, last), and calls @p particle_kernel for - * each particle in the cells. Then, for every particle pair within the - * cell and for each pair with the cells neighbors, @p distance_function is + * Iterates over all cells in [first, last), and for every particle pair within + * the cell and for each pair with the cells neighbors, @p distance_function is * evaluated and @p verlet_criterion is evaluated with the calculated distance. - * Iff true, the pair_kernel is called. + * Iff true, the @p pair_kernel is called. * - * For details see verlet_ia and link_cell. + * For details see @ref verlet_ia and @ref link_cell. * * Requirements on the types: * The Cell type has to provide a function %neighbors() that returns @@ -43,29 +42,30 @@ namespace Algorithm { * container that can be used to store particle pairs. It can be empty and is * not touched if @p use_verlet_list is false. * - * verlet_criterion(p1, p2, distance_function(p1, p2)) has to be valid and - * convertible to bool. + * verlet_criterion(p1, p2, distance_function(p1, p2)) has to be valid + * and convertible to bool. * - * ParticleKernel has to provide an %operator() member that can be called - * with a particle reference. - * PairKernel has to provide an %operator() member that can be called - * with two particle references and a distance. + * @tparam CellIterator cell structure iterator + * @tparam PairKernel has to provide an %operator() member that can be + * called with two particle references and a distance + * @tparam DistanceFunction calculates distances between particles + * @tparam VerletCriterion determines whether a particle pair is interacting, + * e.g. by comparing the particle pair distance to + * the interaction range */ -template +template void for_each_pair(CellIterator first, CellIterator last, - ParticleKernel &&particle_kernel, PairKernel &&pair_kernel, + PairKernel &&pair_kernel, DistanceFunction &&distance_function, VerletCriterion &&verlet_criterion, bool use_verlet_list, bool rebuild) { if (use_verlet_list) { - verlet_ia(first, last, std::forward(particle_kernel), - std::forward(pair_kernel), + verlet_ia(first, last, std::forward(pair_kernel), std::forward(distance_function), std::forward(verlet_criterion), rebuild); } else { - link_cell(first, last, std::forward(particle_kernel), - std::forward(pair_kernel), + link_cell(first, last, std::forward(pair_kernel), std::forward(distance_function)); } } diff --git a/src/core/algorithm/link_cell.hpp b/src/core/algorithm/link_cell.hpp index ef4c6bc6e13..210037f165d 100644 --- a/src/core/algorithm/link_cell.hpp +++ b/src/core/algorithm/link_cell.hpp @@ -26,18 +26,14 @@ namespace Algorithm { * and over all pairs within the cells and with * their neighbors. */ -template -void link_cell(CellIterator first, CellIterator last, - ParticleKernel &&particle_kernel, PairKernel &&pair_kernel, +template +void link_cell(CellIterator first, CellIterator last, PairKernel &&pair_kernel, DistanceFunction &&distance_function) { for (; first != last; ++first) { for (auto it = first->particles().begin(); it != first->particles().end(); ++it) { auto &p1 = *it; - particle_kernel(p1); - /* Pairs in this cell */ for (auto jt = std::next(it); jt != first->particles().end(); ++jt) { auto const dist = distance_function(p1, *jt); diff --git a/src/core/algorithm/verlet_ia.hpp b/src/core/algorithm/verlet_ia.hpp index 55c08639cee..0941c0e98d8 100644 --- a/src/core/algorithm/verlet_ia.hpp +++ b/src/core/algorithm/verlet_ia.hpp @@ -23,11 +23,9 @@ namespace Algorithm { namespace detail { - -template +template void update_and_kernel(CellIterator first, CellIterator last, - ParticleKernel &&particle_kernel, PairKernel &&pair_kernel, DistanceFunction &&distance_function, VerletCriterion &&verlet_criterion) { @@ -39,8 +37,6 @@ void update_and_kernel(CellIterator first, CellIterator last, ++it) { auto &p1 = *it; - particle_kernel(p1); - /* Pairs in this cell */ for (auto jt = std::next(it); jt != first->particles().end(); ++jt) { auto const dist = distance_function(p1, *jt); @@ -64,16 +60,10 @@ void update_and_kernel(CellIterator first, CellIterator last, } } -template -void kernel(CellIterator first, CellIterator last, - ParticleKernel &&particle_kernel, PairKernel &&pair_kernel, +template +void kernel(CellIterator first, CellIterator last, PairKernel &&pair_kernel, DistanceFunction &&distance_function) { for (; first != last; ++first) { - for (auto &p : first->particles()) { - particle_kernel(p); - } - for (auto &pair : first->m_verlet_list) { auto const dist = distance_function(*pair.first, *pair.second); pair_kernel(*pair.first, *pair.second, dist); @@ -88,21 +78,18 @@ void kernel(CellIterator first, CellIterator last, * If rebuild is true, all neighbor cells are iterated * and the Verlet lists are updated with the so found pairs. */ -template -void verlet_ia(CellIterator first, CellIterator last, - ParticleKernel &&particle_kernel, PairKernel &&pair_kernel, +template +void verlet_ia(CellIterator first, CellIterator last, PairKernel &&pair_kernel, DistanceFunction &&distance_function, VerletCriterion &&verlet_criterion, bool rebuild) { if (rebuild) { detail::update_and_kernel(first, last, - std::forward(particle_kernel), std::forward(pair_kernel), std::forward(distance_function), std::forward(verlet_criterion)); } else { - detail::kernel(first, last, std::forward(particle_kernel), - std::forward(pair_kernel), + detail::kernel(first, last, std::forward(pair_kernel), std::forward(distance_function)); } } diff --git a/src/core/short_range_loop.hpp b/src/core/short_range_loop.hpp index f4d5c569b74..05c37a3582e 100644 --- a/src/core/short_range_loop.hpp +++ b/src/core/short_range_loop.hpp @@ -59,21 +59,19 @@ struct EuclidianDistance { * @brief Decide which distance function to use depending on the * cell system, and call the pair code. */ -template +template void decide_distance(CellIterator first, CellIterator last, - ParticleKernel &&particle_kernel, PairKernel &&pair_kernel, + PairKernel &&pair_kernel, VerletCriterion &&verlet_criterion) { if (cell_structure.minimum_image_distance()) { - Algorithm::for_each_pair( - first, last, std::forward(particle_kernel), - std::forward(pair_kernel), MinimalImageDistance{box_geo}, - std::forward(verlet_criterion), - cell_structure.use_verlet_list, rebuild_verletlist); + Algorithm::for_each_pair(first, last, std::forward(pair_kernel), + MinimalImageDistance{box_geo}, + std::forward(verlet_criterion), + cell_structure.use_verlet_list, + rebuild_verletlist); } else { Algorithm::for_each_pair( - first, last, std::forward(particle_kernel), - std::forward(pair_kernel), EuclidianDistance{}, + first, last, std::forward(pair_kernel), EuclidianDistance{}, std::forward(verlet_criterion), cell_structure.use_verlet_list, rebuild_verletlist); } @@ -96,21 +94,21 @@ void short_range_loop(ParticleKernel &&particle_kernel, assert(cell_structure.get_resort_particles() == Cells::RESORT_NONE); + auto local_particles = cell_structure.local_particles(); + for (auto &p : local_particles) { + particle_kernel(p); + } + if (interaction_range() != INACTIVE_CUTOFF) { auto first = boost::make_indirect_iterator(cell_structure.local_cells().begin()); auto last = boost::make_indirect_iterator(cell_structure.local_cells().end()); - detail::decide_distance( - first, last, std::forward(particle_kernel), - std::forward(pair_kernel), verlet_criterion); + detail::decide_distance(first, last, std::forward(pair_kernel), + verlet_criterion); rebuild_verletlist = false; - } else { - for (auto &p : cell_structure.local_particles()) { - particle_kernel(p); - } } } diff --git a/src/core/unit_tests/link_cell_test.cpp b/src/core/unit_tests/link_cell_test.cpp index 55456bcf32b..aac98d9580d 100644 --- a/src/core/unit_tests/link_cell_test.cpp +++ b/src/core/unit_tests/link_cell_test.cpp @@ -54,11 +54,9 @@ BOOST_AUTO_TEST_CASE(link_cell) { std::vector> lc_pairs; lc_pairs.reserve((n_part * (n_part - 1)) / 2); - std::vector id_counts(n_part, 0u); Algorithm::link_cell( cells.begin(), cells.end(), - [&id_counts](Particle const &p) { id_counts[p.p.identity]++; }, [&lc_pairs](Particle const &p1, Particle const &p2, std::pair d) { /* Check that the "distance function" has been called with the correct @@ -71,11 +69,6 @@ BOOST_AUTO_TEST_CASE(link_cell) { return std::make_pair(p1.p.identity, p2.p.identity); }); - /* Check that the particle kernel has been executed exactly once for every - * particle. */ - BOOST_CHECK(std::all_of(id_counts.begin(), id_counts.end(), - [](int count) { return count == 1; })); - BOOST_CHECK(lc_pairs.size() == (n_part * (n_part - 1) / 2)); auto it = lc_pairs.begin(); diff --git a/src/core/unit_tests/verlet_ia_test.cpp b/src/core/unit_tests/verlet_ia_test.cpp index c1360eec77f..dc2c2e1547c 100644 --- a/src/core/unit_tests/verlet_ia_test.cpp +++ b/src/core/unit_tests/verlet_ia_test.cpp @@ -78,12 +78,10 @@ BOOST_AUTO_TEST_CASE(verlet_ia) { std::vector> pairs; pairs.reserve((n_part * (n_part - 1)) / 2); - std::vector id_counts(n_part, 0u); /* Build VL */ Algorithm::verlet_ia( cells.begin(), cells.end(), - [&id_counts](Particle const &p) { id_counts[p.p.identity]++; }, [&pairs](Particle const &p1, Particle const &p2, Distance const &) { pairs.emplace_back(p1.p.identity, p2.p.identity); }, @@ -92,21 +90,14 @@ BOOST_AUTO_TEST_CASE(verlet_ia) { }, VerletCriterion{}, /* rebuild */ true); - /* Check that the particle kernel has been executed exactly once for every - * particle. */ - BOOST_CHECK(std::all_of(id_counts.begin(), id_counts.end(), - [](int count) { return count == 1; })); - check_pairs(n_part, pairs); /* Reset everything */ pairs.clear(); - std::fill(id_counts.begin(), id_counts.end(), 0); /* Now check the Verlet lists */ Algorithm::verlet_ia( cells.begin(), cells.end(), - [&id_counts](Particle const &p) { id_counts[p.p.identity]++; }, [&pairs](Particle const &p1, Particle const &p2, Distance const &) { pairs.emplace_back(p1.p.identity, p2.p.identity); }, @@ -115,21 +106,14 @@ BOOST_AUTO_TEST_CASE(verlet_ia) { }, VerletCriterion{}, /* rebuild */ false); - /* Check that the particle kernel has been executed exactly once for every - * particle. */ - BOOST_CHECK(std::all_of(id_counts.begin(), id_counts.end(), - [](int count) { return count == 1; })); - check_pairs(n_part, pairs); /* Reset everything */ pairs.clear(); - std::fill(id_counts.begin(), id_counts.end(), 0); /* Rebuild again */ Algorithm::verlet_ia( cells.begin(), cells.end(), - [&id_counts](Particle const &p) { id_counts[p.p.identity]++; }, [&pairs](Particle const &p1, Particle const &p2, Distance const &) { pairs.emplace_back(p1.p.identity, p2.p.identity); }, @@ -138,10 +122,5 @@ BOOST_AUTO_TEST_CASE(verlet_ia) { }, VerletCriterion{}, /* rebuild */ true); - /* Check that the particle kernel has been executed exactly once for every - * particle. */ - BOOST_CHECK(std::all_of(id_counts.begin(), id_counts.end(), - [](int count) { return count == 1; })); - check_pairs(n_part, pairs); } From 4bab12f4801a690d9a0306ba873754a0f0566a58 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Thu, 16 Apr 2020 22:57:30 +0200 Subject: [PATCH 02/27] core: short_range_loop: Refactor code Create a new Verlet list at the CellStructure level. Remove global variable `rebuild_verletlist`. Properly detect cutoff of coldet. --- src/core/CellStructure.cpp | 2 + src/core/CellStructure.hpp | 4 +- src/core/algorithm/verlet_ia.hpp | 2 +- src/core/cells.cpp | 4 - src/core/cells.hpp | 5 -- src/core/collision.cpp | 9 +- src/core/collision.hpp | 2 +- .../nonbonded_interaction_data.cpp | 2 + src/core/short_range_loop.hpp | 84 ++++++++++++------- src/core/unit_tests/link_cell_test.cpp | 26 +++--- 10 files changed, 76 insertions(+), 64 deletions(-) diff --git a/src/core/CellStructure.cpp b/src/core/CellStructure.cpp index 630e11f2946..5a5b2a23892 100644 --- a/src/core/CellStructure.cpp +++ b/src/core/CellStructure.cpp @@ -166,6 +166,8 @@ void CellStructure::resort_particles(int global_flag) { for (auto d : diff) { boost::apply_visitor(UpdateParticleIndexVisitor{this}, d); } + + m_rebuild_verlet_list = true; } void CellStructure::set_atom_decomposition(boost::mpi::communicator const &comm, diff --git a/src/core/CellStructure.hpp b/src/core/CellStructure.hpp index 5d12e08e9a3..74cd0cb7cb4 100644 --- a/src/core/CellStructure.hpp +++ b/src/core/CellStructure.hpp @@ -103,6 +103,9 @@ struct CellStructure { public: bool use_verlet_list = true; + bool m_rebuild_verlet_list = true; + std::vector> m_verlet_list; + /** * @brief Update local particle index. * @@ -293,7 +296,6 @@ struct CellStructure { return assert(m_decomposition), *m_decomposition; } -private: ParticleDecomposition &decomposition() { return assert(m_decomposition), *m_decomposition; } diff --git a/src/core/algorithm/verlet_ia.hpp b/src/core/algorithm/verlet_ia.hpp index 0941c0e98d8..1e88d813722 100644 --- a/src/core/algorithm/verlet_ia.hpp +++ b/src/core/algorithm/verlet_ia.hpp @@ -49,7 +49,7 @@ void update_and_kernel(CellIterator first, CellIterator last, /* Pairs with neighbors */ for (auto &neighbor : first->neighbors().red()) { for (auto &p2 : neighbor->particles()) { - auto dist = distance_function(p1, p2); + auto const dist = distance_function(p1, p2); if (verlet_criterion(p1, p2, dist)) { pair_kernel(p1, p2, dist); first->m_verlet_list.emplace_back(&p1, &p2); diff --git a/src/core/cells.cpp b/src/core/cells.cpp index a5efb55e2df..2a71e1f9cd0 100644 --- a/src/core/cells.cpp +++ b/src/core/cells.cpp @@ -49,8 +49,6 @@ /** Type of cell structure in use */ CellStructure cell_structure; -bool rebuild_verletlist = true; - /** * @brief Get pairs closer than distance from the cells. * @@ -144,8 +142,6 @@ void cells_resort_particles(int global_flag) { check_particle_consistency(); check_particle_sorting(); #endif - - rebuild_verletlist = true; } /*************************************************/ diff --git a/src/core/cells.hpp b/src/core/cells.hpp index d5897a580a2..f123d5e9cdf 100644 --- a/src/core/cells.hpp +++ b/src/core/cells.hpp @@ -67,11 +67,6 @@ enum { /** Type of cell structure in use. */ extern CellStructure cell_structure; -/** If non-zero, cell systems should reset the position for checking - * the Verlet criterion. Moreover, the Verlet list has to be rebuilt. - */ -extern bool rebuild_verletlist; - /*@}*/ /************************************************************/ diff --git a/src/core/collision.cpp b/src/core/collision.cpp index d31f443d24f..526bad17641 100644 --- a/src/core/collision.cpp +++ b/src/core/collision.cpp @@ -102,12 +102,6 @@ bool validate_collision_parameters() { // Cache square of cutoff collision_params.distance2 = Utils::sqr(collision_params.distance); - - if (collision_params.distance > min_global_cut) { - runtimeErrorMsg() << "The minimum global cutoff (System.min_global_cut) " - "must be larger or equal the collision detection " - "distance."; - } } #ifndef VIRTUAL_SITES_RELATIVE @@ -241,8 +235,7 @@ bool validate_collision_parameters() { make_particle_type_exist(collision_params.part_type_after_glueing); } - recalc_forces = true; - rebuild_verletlist = true; + on_short_range_ia_change(); return true; } diff --git a/src/core/collision.hpp b/src/core/collision.hpp index 829caa9a24e..77847d65d64 100644 --- a/src/core/collision.hpp +++ b/src/core/collision.hpp @@ -160,7 +160,7 @@ inline double collision_detection_cutoff() { if (collision_params.mode) return collision_params.distance; #endif - return 0.; + return -1.; } #endif diff --git a/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp b/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp index b774c3aca31..2a3dad1e86b 100644 --- a/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp +++ b/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp @@ -23,6 +23,7 @@ */ #include "nonbonded_interactions/nonbonded_interaction_data.hpp" #include "bonded_interactions/bonded_interaction_data.hpp" +#include "collision.hpp" #include "communication.hpp" #include "errorhandling.hpp" #include "grid.hpp" @@ -198,6 +199,7 @@ double maximal_cutoff() { max_cut = std::max(max_cut, max_cut_long_range); max_cut = std::max(max_cut, max_cut_bonded); max_cut = std::max(max_cut, max_cut_nonbonded); + max_cut = std::max(max_cut, collision_detection_cutoff()); return max_cut; } diff --git a/src/core/short_range_loop.hpp b/src/core/short_range_loop.hpp index 05c37a3582e..8c9dfc50d60 100644 --- a/src/core/short_range_loop.hpp +++ b/src/core/short_range_loop.hpp @@ -56,33 +56,56 @@ struct EuclidianDistance { }; /** - * @brief Decide which distance function to use depending on the - * cell system, and call the pair code. + * @brief Functor that returns true for + * any arguments. */ -template -void decide_distance(CellIterator first, CellIterator last, - PairKernel &&pair_kernel, - VerletCriterion &&verlet_criterion) { - if (cell_structure.minimum_image_distance()) { - Algorithm::for_each_pair(first, last, std::forward(pair_kernel), - MinimalImageDistance{box_geo}, - std::forward(verlet_criterion), - cell_structure.use_verlet_list, - rebuild_verletlist); +struct True { + template bool operator()(T...) const { return true; } +}; + +template +void pair_loop(PairKernel &&pair_kernel, DistanceFunction df, + const VerletCriterion &verlet_criterion) { + auto first = + boost::make_indirect_iterator(cell_structure.local_cells().begin()); + auto last = boost::make_indirect_iterator(cell_structure.local_cells().end()); + + if (cell_structure.use_verlet_list && cell_structure.m_rebuild_verlet_list) { + cell_structure.m_verlet_list.clear(); + + Algorithm::link_cell( + first, last, + [&pair_kernel, &df, &verlet_criterion](Particle &p1, Particle &p2, + Distance const &) { + auto const d = df(p1, p2); + if (verlet_criterion(p1, p2, d)) { + cell_structure.m_verlet_list.emplace_back(&p1, &p2); + pair_kernel(p1, p2, d); + } + }, + df); + + cell_structure.m_rebuild_verlet_list = false; + } else if (cell_structure.use_verlet_list && + not cell_structure.m_rebuild_verlet_list) { + for (auto &pair : cell_structure.m_verlet_list) { + pair_kernel(*pair.first, *pair.second, df(*pair.first, *pair.second)); + } } else { - Algorithm::for_each_pair( - first, last, std::forward(pair_kernel), EuclidianDistance{}, - std::forward(verlet_criterion), - cell_structure.use_verlet_list, rebuild_verletlist); + Algorithm::link_cell( + first, last, + [&pair_kernel, &df, &verlet_criterion](Particle &p1, Particle &p2, + Distance const &) { + auto const d = df(p1, p2); + if (verlet_criterion(p1, p2, d)) { + cell_structure.m_verlet_list.emplace_back(&p1, &p2); + pair_kernel(p1, p2, d); + } + }, + df); } } -/** - * @brief Functor that returns true for any argument. - */ -struct True { - template bool operator()(T...) const { return true; } -}; } // namespace detail template (pair_kernel), - verlet_criterion); - - rebuild_verletlist = false; + if (cell_structure.decomposition().minimum_image_distance()) { + detail::pair_loop(std::forward(pair_kernel), + detail::MinimalImageDistance{box_geo}, + verlet_criterion); + } else { + detail::pair_loop(std::forward(pair_kernel), + detail::EuclidianDistance{}, verlet_criterion); + } } } diff --git a/src/core/unit_tests/link_cell_test.cpp b/src/core/unit_tests/link_cell_test.cpp index aac98d9580d..4c9ff876a88 100644 --- a/src/core/unit_tests/link_cell_test.cpp +++ b/src/core/unit_tests/link_cell_test.cpp @@ -55,19 +55,19 @@ BOOST_AUTO_TEST_CASE(link_cell) { std::vector> lc_pairs; lc_pairs.reserve((n_part * (n_part - 1)) / 2); - Algorithm::link_cell( - cells.begin(), cells.end(), - [&lc_pairs](Particle const &p1, Particle const &p2, - std::pair d) { - /* Check that the "distance function" has been called with the correct - * arguments */ - BOOST_CHECK((d.first == p1.p.identity) && (d.second == p2.p.identity)); - if (p1.p.identity <= p2.p.identity) - lc_pairs.emplace_back(p1.p.identity, p2.p.identity); - }, - [](Particle const &p1, Particle const &p2) { - return std::make_pair(p1.p.identity, p2.p.identity); - }); + Algorithm::link_cell(cells.begin(), cells.end(), + [&lc_pairs](Particle const &p1, Particle const &p2, + std::pair d) { + /* Check that the "distance function" has been called + * with the correct arguments */ + BOOST_CHECK((d.first == p1.p.identity) && + (d.second == p2.p.identity)); + if (p1.p.identity <= p2.p.identity) + lc_pairs.emplace_back(p1.p.identity, p2.p.identity); + }, + [](Particle const &p1, Particle const &p2) { + return std::make_pair(p1.p.identity, p2.p.identity); + }); BOOST_CHECK(lc_pairs.size() == (n_part * (n_part - 1) / 2)); From 9832504542cb999f5e8b20f03c55703a7186d6bd Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Thu, 16 Apr 2020 23:40:57 +0200 Subject: [PATCH 03/27] core: short_range_loop: Pull up DistanceFunction --- src/core/algorithm/for_each_pair.hpp | 74 --------------- src/core/algorithm/link_cell.hpp | 12 +-- src/core/algorithm/verlet_ia.hpp | 98 ------------------- src/core/short_range_loop.hpp | 15 +-- src/core/unit_tests/CMakeLists.txt | 1 - src/core/unit_tests/link_cell_test.cpp | 10 +- src/core/unit_tests/verlet_ia_test.cpp | 126 ------------------------- 7 files changed, 11 insertions(+), 325 deletions(-) delete mode 100644 src/core/algorithm/for_each_pair.hpp delete mode 100644 src/core/algorithm/verlet_ia.hpp delete mode 100644 src/core/unit_tests/verlet_ia_test.cpp diff --git a/src/core/algorithm/for_each_pair.hpp b/src/core/algorithm/for_each_pair.hpp deleted file mode 100644 index f64bf699252..00000000000 --- a/src/core/algorithm/for_each_pair.hpp +++ /dev/null @@ -1,74 +0,0 @@ -/* - * Copyright (C) 2010-2019 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 CORE_ALGORITHM_PAIR_LOOP_HPP -#define CORE_ALGORITHM_PAIR_LOOP_HPP - -#include - -#include "link_cell.hpp" -#include "verlet_ia.hpp" - -namespace Algorithm { -/** - * @brief Run pair kernel for each particle pair from a cell range. - * - * Iterates over all cells in [first, last), and for every particle pair within - * the cell and for each pair with the cells neighbors, @p distance_function is - * evaluated and @p verlet_criterion is evaluated with the calculated distance. - * Iff true, the @p pair_kernel is called. - * - * For details see @ref verlet_ia and @ref link_cell. - * - * Requirements on the types: - * The Cell type has to provide a function %neighbors() that returns - * a cell range comprised of the topological neighbors of the cell, - * excluding the cell itself. The cells have to provide a %m_verlet_list - * container that can be used to store particle pairs. It can be empty and is - * not touched if @p use_verlet_list is false. - * - * verlet_criterion(p1, p2, distance_function(p1, p2)) has to be valid - * and convertible to bool. - * - * @tparam CellIterator cell structure iterator - * @tparam PairKernel has to provide an %operator() member that can be - * called with two particle references and a distance - * @tparam DistanceFunction calculates distances between particles - * @tparam VerletCriterion determines whether a particle pair is interacting, - * e.g. by comparing the particle pair distance to - * the interaction range - */ -template -void for_each_pair(CellIterator first, CellIterator last, - PairKernel &&pair_kernel, - DistanceFunction &&distance_function, - VerletCriterion &&verlet_criterion, bool use_verlet_list, - bool rebuild) { - if (use_verlet_list) { - verlet_ia(first, last, std::forward(pair_kernel), - std::forward(distance_function), - std::forward(verlet_criterion), rebuild); - } else { - link_cell(first, last, std::forward(pair_kernel), - std::forward(distance_function)); - } -} -} // namespace Algorithm - -#endif diff --git a/src/core/algorithm/link_cell.hpp b/src/core/algorithm/link_cell.hpp index 210037f165d..28a30f19370 100644 --- a/src/core/algorithm/link_cell.hpp +++ b/src/core/algorithm/link_cell.hpp @@ -26,9 +26,9 @@ namespace Algorithm { * and over all pairs within the cells and with * their neighbors. */ -template -void link_cell(CellIterator first, CellIterator last, PairKernel &&pair_kernel, - DistanceFunction &&distance_function) { +template +void link_cell(CellIterator first, CellIterator last, + PairKernel &&pair_kernel) { for (; first != last; ++first) { for (auto it = first->particles().begin(); it != first->particles().end(); ++it) { @@ -36,15 +36,13 @@ void link_cell(CellIterator first, CellIterator last, PairKernel &&pair_kernel, /* Pairs in this cell */ for (auto jt = std::next(it); jt != first->particles().end(); ++jt) { - auto const dist = distance_function(p1, *jt); - pair_kernel(p1, *jt, dist); + pair_kernel(p1, *jt); } /* Pairs with neighbors */ for (auto &neighbor : first->neighbors().red()) { for (auto &p2 : neighbor->particles()) { - auto const dist = distance_function(p1, p2); - pair_kernel(p1, p2, dist); + pair_kernel(p1, p2); } } } diff --git a/src/core/algorithm/verlet_ia.hpp b/src/core/algorithm/verlet_ia.hpp deleted file mode 100644 index 1e88d813722..00000000000 --- a/src/core/algorithm/verlet_ia.hpp +++ /dev/null @@ -1,98 +0,0 @@ -/* - * Copyright (C) 2010-2019 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 CORE_ALGORITHM_VERLET_IA_HPP -#define CORE_ALGORITHM_VERLET_IA_HPP - -#include - -namespace Algorithm { -namespace detail { -template -void update_and_kernel(CellIterator first, CellIterator last, - PairKernel &&pair_kernel, - DistanceFunction &&distance_function, - VerletCriterion &&verlet_criterion) { - for (; first != last; ++first) { - /* Clear the VL */ - first->m_verlet_list.clear(); - - for (auto it = first->particles().begin(); it != first->particles().end(); - ++it) { - auto &p1 = *it; - - /* Pairs in this cell */ - for (auto jt = std::next(it); jt != first->particles().end(); ++jt) { - auto const dist = distance_function(p1, *jt); - if (verlet_criterion(p1, *jt, dist)) { - pair_kernel(p1, *jt, dist); - first->m_verlet_list.emplace_back(&p1, &(*jt)); - } - } - - /* Pairs with neighbors */ - for (auto &neighbor : first->neighbors().red()) { - for (auto &p2 : neighbor->particles()) { - auto const dist = distance_function(p1, p2); - if (verlet_criterion(p1, p2, dist)) { - pair_kernel(p1, p2, dist); - first->m_verlet_list.emplace_back(&p1, &p2); - } - } - } - } - } -} - -template -void kernel(CellIterator first, CellIterator last, PairKernel &&pair_kernel, - DistanceFunction &&distance_function) { - for (; first != last; ++first) { - for (auto &pair : first->m_verlet_list) { - auto const dist = distance_function(*pair.first, *pair.second); - pair_kernel(*pair.first, *pair.second, dist); - } - } -} -} // namespace detail - -/** - * @brief Iterates over all particles in the cell range - * and all pairs in the Verlet list of the cells. - * If rebuild is true, all neighbor cells are iterated - * and the Verlet lists are updated with the so found pairs. - */ -template -void verlet_ia(CellIterator first, CellIterator last, PairKernel &&pair_kernel, - DistanceFunction &&distance_function, - VerletCriterion &&verlet_criterion, bool rebuild) { - if (rebuild) { - detail::update_and_kernel(first, last, - std::forward(pair_kernel), - std::forward(distance_function), - std::forward(verlet_criterion)); - } else { - detail::kernel(first, last, std::forward(pair_kernel), - std::forward(distance_function)); - } -} -} // namespace Algorithm - -#endif diff --git a/src/core/short_range_loop.hpp b/src/core/short_range_loop.hpp index 8c9dfc50d60..99745b5c521 100644 --- a/src/core/short_range_loop.hpp +++ b/src/core/short_range_loop.hpp @@ -19,7 +19,7 @@ #ifndef CORE_SHORT_RANGE_HPP #define CORE_SHORT_RANGE_HPP -#include "algorithm/for_each_pair.hpp" +#include "algorithm/link_cell.hpp" #include "cells.hpp" #include "grid.hpp" #include "integrate.hpp" @@ -75,16 +75,13 @@ void pair_loop(PairKernel &&pair_kernel, DistanceFunction df, Algorithm::link_cell( first, last, - [&pair_kernel, &df, &verlet_criterion](Particle &p1, Particle &p2, - Distance const &) { + [&pair_kernel, &df, &verlet_criterion](Particle &p1, Particle &p2) { auto const d = df(p1, p2); if (verlet_criterion(p1, p2, d)) { cell_structure.m_verlet_list.emplace_back(&p1, &p2); pair_kernel(p1, p2, d); } - }, - df); - + }); cell_structure.m_rebuild_verlet_list = false; } else if (cell_structure.use_verlet_list && not cell_structure.m_rebuild_verlet_list) { @@ -94,15 +91,13 @@ void pair_loop(PairKernel &&pair_kernel, DistanceFunction df, } else { Algorithm::link_cell( first, last, - [&pair_kernel, &df, &verlet_criterion](Particle &p1, Particle &p2, - Distance const &) { + [&pair_kernel, &df, &verlet_criterion](Particle &p1, Particle &p2) { auto const d = df(p1, p2); if (verlet_criterion(p1, p2, d)) { cell_structure.m_verlet_list.emplace_back(&p1, &p2); pair_kernel(p1, p2, d); } - }, - df); + }); } } diff --git a/src/core/unit_tests/CMakeLists.txt b/src/core/unit_tests/CMakeLists.txt index 735f3145ccf..abcc0a53e8b 100644 --- a/src/core/unit_tests/CMakeLists.txt +++ b/src/core/unit_tests/CMakeLists.txt @@ -41,7 +41,6 @@ unit_test(NAME Variant_test SRC Variant_test.cpp DEPENDS ScriptInterface) unit_test(NAME ParticleIterator_test SRC ParticleIterator_test.cpp DEPENDS EspressoUtils) unit_test(NAME link_cell_test SRC link_cell_test.cpp DEPENDS EspressoUtils) -unit_test(NAME verlet_ia_test SRC verlet_ia_test.cpp DEPENDS EspressoUtils) unit_test(NAME Particle_test SRC Particle_test.cpp DEPENDS EspressoUtils Boost::serialization) unit_test(NAME get_value SRC get_value_test.cpp DEPENDS ScriptInterface) diff --git a/src/core/unit_tests/link_cell_test.cpp b/src/core/unit_tests/link_cell_test.cpp index 4c9ff876a88..8f711f14d98 100644 --- a/src/core/unit_tests/link_cell_test.cpp +++ b/src/core/unit_tests/link_cell_test.cpp @@ -56,17 +56,9 @@ BOOST_AUTO_TEST_CASE(link_cell) { lc_pairs.reserve((n_part * (n_part - 1)) / 2); Algorithm::link_cell(cells.begin(), cells.end(), - [&lc_pairs](Particle const &p1, Particle const &p2, - std::pair d) { - /* Check that the "distance function" has been called - * with the correct arguments */ - BOOST_CHECK((d.first == p1.p.identity) && - (d.second == p2.p.identity)); + [&lc_pairs](Particle const &p1, Particle const &p2) { if (p1.p.identity <= p2.p.identity) lc_pairs.emplace_back(p1.p.identity, p2.p.identity); - }, - [](Particle const &p1, Particle const &p2) { - return std::make_pair(p1.p.identity, p2.p.identity); }); BOOST_CHECK(lc_pairs.size() == (n_part * (n_part - 1) / 2)); diff --git a/src/core/unit_tests/verlet_ia_test.cpp b/src/core/unit_tests/verlet_ia_test.cpp deleted file mode 100644 index dc2c2e1547c..00000000000 --- a/src/core/unit_tests/verlet_ia_test.cpp +++ /dev/null @@ -1,126 +0,0 @@ -/* - * Copyright (C) 2010-2019 The ESPResSo project - * - * This file is part of ESPResSo. - * - * ESPResSo is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * ESPResSo is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see . - */ -#include -#include -#include - -#define BOOST_TEST_MODULE link_cell test -#define BOOST_TEST_DYN_LINK -#include - -#include "Cell.hpp" -#include "algorithm/verlet_ia.hpp" - -void check_pairs(int n_part, std::vector> const &pairs) { - BOOST_CHECK(pairs.size() == (n_part * (n_part - 1) / 2)); - - auto it = pairs.begin(); - for (int i = 0; i < n_part; i++) - for (int j = i + 1; j < n_part; j++) { - BOOST_CHECK((it->first == i) && (it->second == j)); - ++it; - } -} - -/* Dummy distance */ -struct Distance { - bool interact; -}; - -/* Dummy interaction criterion */ -struct VerletCriterion { - bool operator()(Particle const &p1, Particle const &p2, - Distance const &d) const { - return d.interact; - } -}; - -BOOST_AUTO_TEST_CASE(verlet_ia) { - const unsigned n_cells = 100; - const auto n_part_per_cell = 10; - const auto n_part = n_cells * n_part_per_cell; - - std::vector cells(n_cells); - - auto id = 0; - for (auto &c : cells) { - std::vector neighbors; - - for (auto &n : cells) { - if (&c != &n) - neighbors.push_back(&n); - } - - c.m_neighbors = Neighbors(neighbors, {}); - - c.particles().resize(n_part_per_cell); - - for (auto &p : c.particles()) { - p.p.identity = id++; - } - } - - std::vector> pairs; - pairs.reserve((n_part * (n_part - 1)) / 2); - - /* Build VL */ - Algorithm::verlet_ia( - cells.begin(), cells.end(), - [&pairs](Particle const &p1, Particle const &p2, Distance const &) { - pairs.emplace_back(p1.p.identity, p2.p.identity); - }, - [](Particle const &p1, Particle const &p2) { - return Distance{p1.p.identity <= p2.p.identity}; - }, - VerletCriterion{}, /* rebuild */ true); - - check_pairs(n_part, pairs); - - /* Reset everything */ - pairs.clear(); - - /* Now check the Verlet lists */ - Algorithm::verlet_ia( - cells.begin(), cells.end(), - [&pairs](Particle const &p1, Particle const &p2, Distance const &) { - pairs.emplace_back(p1.p.identity, p2.p.identity); - }, - [](Particle const &p1, Particle const &p2) { - return Distance{p1.p.identity <= p2.p.identity}; - }, - VerletCriterion{}, /* rebuild */ false); - - check_pairs(n_part, pairs); - - /* Reset everything */ - pairs.clear(); - - /* Rebuild again */ - Algorithm::verlet_ia( - cells.begin(), cells.end(), - [&pairs](Particle const &p1, Particle const &p2, Distance const &) { - pairs.emplace_back(p1.p.identity, p2.p.identity); - }, - [](Particle const &p1, Particle const &p2) { - return Distance{p1.p.identity <= p2.p.identity}; - }, - VerletCriterion{}, /* rebuild */ true); - - check_pairs(n_part, pairs); -} From 6396d302d96b7c71981eb42aee2e337bb1822c76 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Fri, 17 Apr 2020 00:12:58 +0200 Subject: [PATCH 04/27] core: short_range_loop: Specialization without particle kernel --- src/core/dpd.cpp | 1 - .../electrostatics_magnetostatics/icc.cpp | 5 ++--- src/core/short_range_loop.hpp | 19 +++++++++++++++++++ 3 files changed, 21 insertions(+), 4 deletions(-) diff --git a/src/core/dpd.cpp b/src/core/dpd.cpp index 87846abacfd..e897c2e2475 100644 --- a/src/core/dpd.cpp +++ b/src/core/dpd.cpp @@ -143,7 +143,6 @@ static auto dpd_viscous_stress_local() { Utils::Vector stress{}; short_range_loop( - Utils::NoOp{}, [&stress](const Particle &p1, const Particle &p2, Distance const &d) { auto const v21 = p1.m.v - p2.m.v; diff --git a/src/core/electrostatics_magnetostatics/icc.cpp b/src/core/electrostatics_magnetostatics/icc.cpp index ec77dbbd577..66f227270ed 100644 --- a/src/core/electrostatics_magnetostatics/icc.cpp +++ b/src/core/electrostatics_magnetostatics/icc.cpp @@ -197,9 +197,8 @@ void force_calc_iccp3m(const ParticleRange &particles, const ParticleRange &ghost_particles) { init_forces_iccp3m(particles, ghost_particles); - short_range_loop(Utils::NoOp{}, [](Particle &p1, Particle &p2, - Distance const &d) { - /* calc non-bonded interactions */ + short_range_loop([](Particle &p1, Particle &p2, Distance const &d) { + /* calc non bonded interactions */ add_non_bonded_pair_force_iccp3m(p1, p2, d.vec21, sqrt(d.dist2), d.dist2); }); diff --git a/src/core/short_range_loop.hpp b/src/core/short_range_loop.hpp index 99745b5c521..4e69f265ff6 100644 --- a/src/core/short_range_loop.hpp +++ b/src/core/short_range_loop.hpp @@ -103,6 +103,25 @@ void pair_loop(PairKernel &&pair_kernel, DistanceFunction df, } // namespace detail +template +void short_range_loop(PairKernel &&pair_kernel, + const VerletCriterion &verlet_criterion = {}) { + ESPRESSO_PROFILER_CXX_MARK_FUNCTION; + + assert(cell_structure.get_resort_particles() == Cells::RESORT_NONE); + + if (interaction_range() != INACTIVE_CUTOFF) { + if (cell_structure.decomposition().minimum_image_distance()) { + detail::pair_loop(std::forward(pair_kernel), + detail::MinimalImageDistance{box_geo}, + verlet_criterion); + } else { + detail::pair_loop(std::forward(pair_kernel), + detail::EuclidianDistance{}, verlet_criterion); + } + } +} + template void short_range_loop(ParticleKernel &&particle_kernel, From ef5999dbbd0f77720134041d499eb5c80799d68a Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Fri, 17 Apr 2020 01:01:40 +0200 Subject: [PATCH 05/27] core: Reduce use of CellStructure.local_particles Replace calls to CellStructure::execute_bond_handler() by the convenience method CellStructure::bond_loop(). --- src/core/CellStructure.hpp | 6 + src/core/event.cpp | 10 +- .../immersed_boundary/ImmersedBoundaries.cpp | 118 ++++++++--------- .../immersed_boundary/ImmersedBoundaries.hpp | 3 +- .../object-in-fluid/oif_global_forces.cpp | 125 +++++++++--------- src/core/rattle.cpp | 56 ++++---- 6 files changed, 157 insertions(+), 161 deletions(-) diff --git a/src/core/CellStructure.hpp b/src/core/CellStructure.hpp index 74cd0cb7cb4..82d70910cd6 100644 --- a/src/core/CellStructure.hpp +++ b/src/core/CellStructure.hpp @@ -451,6 +451,12 @@ struct CellStructure { bool minimum_image_distance() const { return m_decomposition->minimum_image_distance(); } + + template void bond_loop(BondKernel const &bond_kernel) { + for (auto &p : local_particles()) { + execute_bond_handler(p, bond_kernel); + } + } }; #endif // ESPRESSO_CELLSTRUCTURE_HPP diff --git a/src/core/event.cpp b/src/core/event.cpp index 15395fd6418..5bc622c8904 100644 --- a/src/core/event.cpp +++ b/src/core/event.cpp @@ -113,11 +113,6 @@ void on_integration_start() { sizeof(CUDA_global_part_vars), MPI_BYTE, 0, comm_cart); #endif - // Here we initialize volume conservation - // This function checks if the reference volumes have been set and if - // necessary calculates them - immersed_boundaries.init_volume_conservation(cell_structure); - /* Prepare the thermostat */ if (reinit_thermo) { thermo_init(); @@ -373,4 +368,9 @@ void update_dependent_particles() { #ifdef ELECTROSTATICS Coulomb::update_dependent_particles(); #endif + + // Here we initialize volume conservation + // This function checks if the reference volumes have been set and if + // necessary calculates them + immersed_boundaries.init_volume_conservation(cell_structure); } diff --git a/src/core/immersed_boundary/ImmersedBoundaries.cpp b/src/core/immersed_boundary/ImmersedBoundaries.cpp index 6059b54cab0..943f7cbf817 100644 --- a/src/core/immersed_boundary/ImmersedBoundaries.cpp +++ b/src/core/immersed_boundary/ImmersedBoundaries.cpp @@ -45,16 +45,19 @@ void ImmersedBoundaries::volume_conservation(CellStructure &cs) { /** Initialize volume conservation */ void ImmersedBoundaries::init_volume_conservation(CellStructure &cs) { - // Check since this function is called at the start of every integrate loop // Also check if volume has been set due to reading of a checkpoint - if (!VolumeInitDone) { + if (not BoundariesFound) { + BoundariesFound = boost::algorithm::any_of( + bonded_ia_params, [](auto const &bonded_ia_param) { + return bonded_ia_param.type == BONDED_IA_IBM_VOLUME_CONSERVATION; + }); + } + if (!VolumeInitDone && BoundariesFound) { // Calculate volumes calc_volumes(cs); - // numWriteCOM = 0; - // Loop through all bonded interactions and check if we need to set the // reference volume for (auto &bonded_ia_param : bonded_ia_params) { @@ -70,9 +73,9 @@ void ImmersedBoundaries::init_volume_conservation(CellStructure &cs) { } } } - } - VolumeInitDone = true; + VolumeInitDone = true; + } } /** Set parameters of volume conservation */ @@ -188,58 +191,53 @@ void ImmersedBoundaries::calc_volumes(CellStructure &cs) { /** Calculate and add the volume force to each node */ void ImmersedBoundaries::calc_volume_force(CellStructure &cs) { - // Loop over all particles on local node - for (auto &p1 : cs.local_particles()) { - // Check if particle has a BONDED_IA_IBM_TRIEL and a - // BONDED_IA_IBM_VOLUME_CONSERVATION. Basically this loops over all - // triangles, not all particles. First round to check for volume - // conservation. - const IBM_VolCons_Parameters *ibmVolConsParameters = - vol_cons_parameters(p1); - if (not ibmVolConsParameters) - return; - - auto current_volume = VolumesCurrent[ibmVolConsParameters->softID]; - - // Second round for triel - cs.execute_bond_handler(p1, [ibmVolConsParameters, current_volume]( - Particle &p1, int bond_id, - Utils::Span partners) { - auto const &iaparams = bonded_ia_params[bond_id]; - - if (iaparams.type == BONDED_IA_IBM_TRIEL) { - // Our particle is the leading particle of a triel - // Get second and third particle of the triangle - Particle &p2 = *partners[0]; - Particle &p3 = *partners[1]; - - // Unfold position of first node. - // This is to get a continuous trajectory with no jumps when box - // boundaries are crossed. - auto const x1 = unfolded_position(p1.r.p, p1.l.i, box_geo.length()); - - // Unfolding seems to work only for the first particle of a triel - // so get the others from relative vectors considering PBC - auto const a12 = get_mi_vector(p2.r.p, x1, box_geo); - auto const a13 = get_mi_vector(p3.r.p, x1, box_geo); - - // Now we have the true and good coordinates - // This is eq. (9) in @cite dupin08a. - auto const n = vector_product(a12, a13); - const double ln = n.norm(); - const double A = 0.5 * ln; - const double fact = ibmVolConsParameters->kappaV * - (current_volume - ibmVolConsParameters->volRef) / - current_volume; - - auto const nHat = n / ln; - auto const force = -fact * A * nHat; - - p1.f.f += force; - p2.f.f += force; - p3.f.f += force; - } - return false; - }); - } + cs.bond_loop( + [this](Particle &p1, int bond_id, Utils::Span partners) { + auto const &iaparams = bonded_ia_params[bond_id]; + + if (iaparams.type == BONDED_IA_IBM_TRIEL) { + // Check if particle has a BONDED_IA_IBM_TRIEL and a + // BONDED_IA_IBM_VOLUME_CONSERVATION. Basically this loops over all + // triangles, not all particles. First round to check for volume + // conservation. + const IBM_VolCons_Parameters *ibmVolConsParameters = + vol_cons_parameters(p1); + if (not ibmVolConsParameters) + return false; + + auto current_volume = VolumesCurrent[ibmVolConsParameters->softID]; + + // Our particle is the leading particle of a triel + // Get second and third particle of the triangle + Particle &p2 = *partners[0]; + Particle &p3 = *partners[1]; + + // Unfold position of first node. + // This is to get a continuous trajectory with no jumps when box + // boundaries are crossed. + auto const x1 = unfolded_position(p1.r.p, p1.l.i, box_geo.length()); + + // Unfolding seems to work only for the first particle of a triel + // so get the others from relative vectors considering PBC + auto const a12 = get_mi_vector(p2.r.p, x1, box_geo); + auto const a13 = get_mi_vector(p3.r.p, x1, box_geo); + + // Now we have the true and good coordinates + // This is eq. (9) in @cite dupin08a. + auto const n = vector_product(a12, a13); + const double ln = n.norm(); + const double A = 0.5 * ln; + const double fact = ibmVolConsParameters->kappaV * + (current_volume - ibmVolConsParameters->volRef) / + current_volume; + + auto const nHat = n / ln; + auto const force = -fact * A * nHat; + + p1.f.f += force; + p2.f.f += force; + p3.f.f += force; + } + return false; + }); } diff --git a/src/core/immersed_boundary/ImmersedBoundaries.hpp b/src/core/immersed_boundary/ImmersedBoundaries.hpp index 839927a9154..a89a9ad7873 100644 --- a/src/core/immersed_boundary/ImmersedBoundaries.hpp +++ b/src/core/immersed_boundary/ImmersedBoundaries.hpp @@ -32,10 +32,11 @@ class ImmersedBoundaries { void init_volume_conservation(CellStructure &cs); void volume_conservation(CellStructure &cs); int volume_conservation_set_params(int bond_type, int softID, double kappaV); + +private: void calc_volumes(CellStructure &cs); void calc_volume_force(CellStructure &cs); -private: const int MaxNumIBM; std::vector VolumesCurrent; bool VolumeInitDone = false; diff --git a/src/core/object-in-fluid/oif_global_forces.cpp b/src/core/object-in-fluid/oif_global_forces.cpp index 7708a246004..ff50fda95a5 100644 --- a/src/core/object-in-fluid/oif_global_forces.cpp +++ b/src/core/object-in-fluid/oif_global_forces.cpp @@ -63,34 +63,32 @@ void calc_oif_global(Utils::Vector2d &area_volume, int molType, Utils::Vector2d part_area_volume; // added - for (auto &p : cs.local_particles()) { - if (p.p.mol_id != molType) - continue; - - cs.execute_bond_handler(p, [&partArea, &VOL_partVol]( - Particle &p1, int bond_id, - Utils::Span partners) { - auto const &iaparams = bonded_ia_params[bond_id]; - - if (iaparams.type == BONDED_IA_OIF_GLOBAL_FORCES) { - // remaining neighbors fetched - auto const p11 = unfolded_position(p1.r.p, p1.l.i, box_geo.length()); - auto const p22 = p11 + get_mi_vector(partners[0]->r.p, p11, box_geo); - auto const p33 = p11 + get_mi_vector(partners[1]->r.p, p11, box_geo); - - // unfolded positions correct - auto const VOL_A = area_triangle(p11, p22, p33); - partArea += VOL_A; - - auto const VOL_norm = get_n_triangle(p11, p22, p33); - auto const VOL_dn = VOL_norm.norm(); - auto const VOL_hz = 1.0 / 3.0 * (p11[2] + p22[2] + p33[2]); - VOL_partVol += VOL_A * -1 * VOL_norm[2] / VOL_dn * VOL_hz; - } + cs.bond_loop( + [&partArea, &VOL_partVol, molType](Particle &p1, int bond_id, + Utils::Span partners) { + if (p1.p.mol_id != molType) + return false; - return false; - }); - } + auto const &iaparams = bonded_ia_params[bond_id]; + + if (iaparams.type == BONDED_IA_OIF_GLOBAL_FORCES) { + // remaining neighbors fetched + auto const p11 = unfolded_position(p1.r.p, p1.l.i, box_geo.length()); + auto const p22 = p11 + get_mi_vector(partners[0]->r.p, p11, box_geo); + auto const p33 = p11 + get_mi_vector(partners[1]->r.p, p11, box_geo); + + // unfolded positions correct + auto const VOL_A = area_triangle(p11, p22, p33); + partArea += VOL_A; + + auto const VOL_norm = get_n_triangle(p11, p22, p33); + auto const VOL_dn = VOL_norm.norm(); + auto const VOL_hz = 1.0 / 3.0 * (p11[2] + p22[2] + p33[2]); + VOL_partVol += VOL_A * -1 * VOL_norm[2] / VOL_dn * VOL_hz; + } + + return false; + }); part_area_volume[0] = partArea; part_area_volume[1] = VOL_partVol; @@ -105,55 +103,52 @@ void add_oif_global_forces(Utils::Vector2d const &area_volume, int molType, double area = area_volume[0]; double VOL_volume = area_volume[1]; - for (auto &p : cs.local_particles()) { - if (p.p.mol_id != molType) - continue; + cs.bond_loop([area, VOL_volume, molType](Particle &p1, int bond_id, + Utils::Span partners) { + if (p1.p.mol_id != molType) + return false; - cs.execute_bond_handler(p, [area, - VOL_volume](Particle &p1, int bond_id, - Utils::Span partners) { - auto const &iaparams = bonded_ia_params[bond_id]; + auto const &iaparams = bonded_ia_params[bond_id]; - if (iaparams.type == BONDED_IA_OIF_GLOBAL_FORCES) { - auto const p11 = unfolded_position(p1.r.p, p1.l.i, box_geo.length()); - auto const p22 = p11 + get_mi_vector(partners[0]->r.p, p11, box_geo); - auto const p33 = p11 + get_mi_vector(partners[1]->r.p, p11, box_geo); + if (iaparams.type == BONDED_IA_OIF_GLOBAL_FORCES) { + auto const p11 = unfolded_position(p1.r.p, p1.l.i, box_geo.length()); + auto const p22 = p11 + get_mi_vector(partners[0]->r.p, p11, box_geo); + auto const p33 = p11 + get_mi_vector(partners[1]->r.p, p11, box_geo); - // unfolded positions correct - // starting code from volume force - auto const VOL_norm = get_n_triangle(p11, p22, p33).normalize(); - auto const VOL_A = area_triangle(p11, p22, p33); - auto const VOL_vv = (VOL_volume - iaparams.p.oif_global_forces.V0) / - iaparams.p.oif_global_forces.V0; + // unfolded positions correct + // starting code from volume force + auto const VOL_norm = get_n_triangle(p11, p22, p33).normalize(); + auto const VOL_A = area_triangle(p11, p22, p33); + auto const VOL_vv = (VOL_volume - iaparams.p.oif_global_forces.V0) / + iaparams.p.oif_global_forces.V0; - auto const VOL_force = (1.0 / 3.0) * iaparams.p.oif_global_forces.kv * - VOL_vv * VOL_A * VOL_norm; + auto const VOL_force = (1.0 / 3.0) * iaparams.p.oif_global_forces.kv * + VOL_vv * VOL_A * VOL_norm; - auto const h = (1. / 3.) * (p11 + p22 + p33); + auto const h = (1. / 3.) * (p11 + p22 + p33); - auto const deltaA = (area - iaparams.p.oif_global_forces.A0_g) / - iaparams.p.oif_global_forces.A0_g; + auto const deltaA = (area - iaparams.p.oif_global_forces.A0_g) / + iaparams.p.oif_global_forces.A0_g; - auto const m1 = h - p11; - auto const m2 = h - p22; - auto const m3 = h - p33; + auto const m1 = h - p11; + auto const m2 = h - p22; + auto const m3 = h - p33; - auto const m1_length = m1.norm(); - auto const m2_length = m2.norm(); - auto const m3_length = m3.norm(); + auto const m1_length = m1.norm(); + auto const m2_length = m2.norm(); + auto const m3_length = m3.norm(); - auto const fac = iaparams.p.oif_global_forces.ka_g * VOL_A * deltaA / - (m1_length * m1_length + m2_length * m2_length + - m3_length * m3_length); + auto const fac = iaparams.p.oif_global_forces.ka_g * VOL_A * deltaA / + (m1_length * m1_length + m2_length * m2_length + + m3_length * m3_length); - p1.f.f += fac * m1 + VOL_force; - partners[0]->f.f += fac * m2 + VOL_force; - partners[1]->f.f += fac * m3 + VOL_force; - } + p1.f.f += fac * m1 + VOL_force; + partners[0]->f.f += fac * m2 + VOL_force; + partners[1]->f.f += fac * m3 + VOL_force; + } - return false; - }); - } + return false; + }); } int max_oif_objects = 0; diff --git a/src/core/rattle.cpp b/src/core/rattle.cpp index 3e940a7fe82..cfa9e24c15c 100644 --- a/src/core/rattle.cpp +++ b/src/core/rattle.cpp @@ -136,21 +136,19 @@ static bool add_pos_corr_vec(Bonded_ia_parameters const &ia_params, } /**Compute positional corrections*/ static void compute_pos_corr_vec(int *repeat_, CellStructure &cs) { - for (auto &p1 : cs.local_particles()) { - cs.execute_bond_handler(p1, [repeat_](Particle &p1, int bond_id, - Utils::Span partners) { - auto const &iaparams = bonded_ia_params[bond_id]; - - if (iaparams.type == BONDED_IA_RIGID_BOND) { - auto const corrected = add_pos_corr_vec(iaparams, p1, *partners[0]); - if (corrected) - *repeat_ += 1; - } - - /* Rigid bonds cannot break */ - return false; - }); - } + cs.bond_loop( + [repeat_](Particle &p1, int bond_id, Utils::Span partners) { + auto const &iaparams = bonded_ia_params[bond_id]; + + if (iaparams.type == BONDED_IA_RIGID_BOND) { + auto const corrected = add_pos_corr_vec(iaparams, p1, *partners[0]); + if (corrected) + *repeat_ += 1; + } + + /* Rigid bonds cannot break */ + return false; + }); } /**Apply corrections to each particle**/ @@ -253,21 +251,19 @@ static bool add_vel_corr_vec(Bonded_ia_parameters const &ia_params, /** Velocity correction vectors are computed*/ static void compute_vel_corr_vec(int *repeat_, CellStructure &cs) { - for (auto &p1 : cs.local_particles()) { - cs.execute_bond_handler(p1, [repeat_](Particle &p1, int bond_id, - Utils::Span partners) { - auto const &iaparams = bonded_ia_params[bond_id]; - - if (iaparams.type == BONDED_IA_RIGID_BOND) { - auto const corrected = add_vel_corr_vec(iaparams, p1, *partners[0]); - if (corrected) - *repeat_ += 1; - } - - /* Rigid bonds can not break */ - return false; - }); - } + cs.bond_loop( + [repeat_](Particle &p1, int bond_id, Utils::Span partners) { + auto const &iaparams = bonded_ia_params[bond_id]; + + if (iaparams.type == BONDED_IA_RIGID_BOND) { + auto const corrected = add_vel_corr_vec(iaparams, p1, *partners[0]); + if (corrected) + *repeat_ += 1; + } + + /* Rigid bonds can not break */ + return false; + }); } /**Apply velocity corrections*/ From c490e5d0f6d43c77a88641b1fbfd72b4af3501c9 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Thu, 23 Apr 2020 10:57:45 +0200 Subject: [PATCH 06/27] core: cells: Move pair_loop() to CellStructure Pass BondKernel as argument to short_range_loop and delegate pair loop to the CellStructure. Rewrite pressure/force/energy kernels from free functions to lambdas. --- src/core/CellStructure.hpp | 33 ++++++++ src/core/cells.cpp | 5 +- src/core/dpd.cpp | 3 +- .../electrostatics_magnetostatics/icc.cpp | 11 ++- src/core/energy.cpp | 13 +++- src/core/energy_inline.hpp | 18 ----- src/core/forces.cpp | 2 +- src/core/forces_inline.hpp | 6 -- src/core/pressure.cpp | 26 +++++-- src/core/pressure_inline.hpp | 23 ------ src/core/short_range_loop.hpp | 77 +++---------------- 11 files changed, 88 insertions(+), 129 deletions(-) diff --git a/src/core/CellStructure.hpp b/src/core/CellStructure.hpp index 82d70910cd6..8e41a7aebba 100644 --- a/src/core/CellStructure.hpp +++ b/src/core/CellStructure.hpp @@ -30,11 +30,13 @@ #include "ParticleDecomposition.hpp" #include "ParticleList.hpp" #include "ParticleRange.hpp" +#include "algorithm/link_cell.hpp" #include "bond_error.hpp" #include "ghosts.hpp" #include #include +#include #include #include @@ -457,6 +459,37 @@ struct CellStructure { execute_bond_handler(p, bond_kernel); } } + + template + void pair_loop(PairKernel &&pair_kernel, DistanceFunction df, + const VerletCriterion &verlet_criterion) { + auto first = boost::make_indirect_iterator(local_cells().begin()); + auto last = boost::make_indirect_iterator(local_cells().end()); + + if (use_verlet_list && m_rebuild_verlet_list) { + m_verlet_list.clear(); + + Algorithm::link_cell(first, last, + [&pair_kernel, &df, &verlet_criterion, + this](Particle &p1, Particle &p2) { + auto const d = df(p1, p2); + if (verlet_criterion(p1, p2, d)) { + m_verlet_list.emplace_back(&p1, &p2); + pair_kernel(p1, p2, d); + } + }); + m_rebuild_verlet_list = false; + } else if (use_verlet_list && not m_rebuild_verlet_list) { + for (auto &pair : m_verlet_list) { + pair_kernel(*pair.first, *pair.second, df(*pair.first, *pair.second)); + } + } else { + Algorithm::link_cell(first, last, + [&pair_kernel, &df](Particle &p1, Particle &p2) { + pair_kernel(p1, p2, df(p1, p2)); + }); + } + } }; #endif // ESPRESSO_CELLSTRUCTURE_HPP diff --git a/src/core/cells.cpp b/src/core/cells.cpp index 2a71e1f9cd0..d72e6b32eb8 100644 --- a/src/core/cells.cpp +++ b/src/core/cells.cpp @@ -38,7 +38,6 @@ #include "AtomDecomposition.hpp" #include "DomainDecomposition.hpp" -#include #include #include @@ -70,7 +69,7 @@ std::vector> get_pairs(double distance) { ret.emplace_back(p1.p.identity, p2.p.identity); }; - short_range_loop(Utils::NoOp{}, pair_kernel); + short_range_loop(pair_kernel, detail::True{}); /* Sort pairs */ for (auto &pair : ret) { @@ -213,4 +212,4 @@ const DomainDecomposition *get_domain_decomposition() { void cells_set_use_verlet_lists(bool use_verlet_lists) { cell_structure.use_verlet_list = use_verlet_lists; -} \ No newline at end of file +} diff --git a/src/core/dpd.cpp b/src/core/dpd.cpp index e897c2e2475..1f064e48c5a 100644 --- a/src/core/dpd.cpp +++ b/src/core/dpd.cpp @@ -159,7 +159,8 @@ static auto dpd_viscous_stress_local() { auto const f = P * (f_r - f_t) + f_t; stress += tensor_product(d.vec21, f); - }); + }, + detail::True{}); return stress; } diff --git a/src/core/electrostatics_magnetostatics/icc.cpp b/src/core/electrostatics_magnetostatics/icc.cpp index 66f227270ed..27f7735a332 100644 --- a/src/core/electrostatics_magnetostatics/icc.cpp +++ b/src/core/electrostatics_magnetostatics/icc.cpp @@ -197,10 +197,13 @@ void force_calc_iccp3m(const ParticleRange &particles, const ParticleRange &ghost_particles) { init_forces_iccp3m(particles, ghost_particles); - short_range_loop([](Particle &p1, Particle &p2, Distance const &d) { - /* calc non bonded interactions */ - add_non_bonded_pair_force_iccp3m(p1, p2, d.vec21, sqrt(d.dist2), d.dist2); - }); + short_range_loop( + [](Particle &p1, Particle &p2, Distance const &d) { + /* calc non bonded interactions */ + add_non_bonded_pair_force_iccp3m(p1, p2, d.vec21, sqrt(d.dist2), + d.dist2); + }, + detail::True{}); Coulomb::calc_long_range_force(particles); } diff --git a/src/core/energy.cpp b/src/core/energy.cpp index af3e3485e15..99a14171994 100644 --- a/src/core/energy.cpp +++ b/src/core/energy.cpp @@ -71,11 +71,20 @@ void energy_calc(const double time) { } short_range_loop( - [](Particle &p) { add_bonded_energy(p, obs_energy); }, + [](Particle &p1, int bond_id, Utils::Span partners) { + auto const &iaparams = bonded_ia_params[bond_id]; + auto const result = calc_bonded_energy(iaparams, p1, partners); + if (result) { + obs_energy.bonded_contribution(bond_id)[0] += result.get(); + return false; + } + return true; + }, [](Particle const &p1, Particle const &p2, Distance const &d) { add_non_bonded_pair_energy(p1, p2, d.vec21, sqrt(d.dist2), d.dist2, obs_energy); - }); + }, + detail::True{}); calc_long_range_energies(cell_structure.local_particles()); diff --git a/src/core/energy_inline.hpp b/src/core/energy_inline.hpp index 26242f87633..942556e3c14 100644 --- a/src/core/energy_inline.hpp +++ b/src/core/energy_inline.hpp @@ -301,22 +301,4 @@ inline double calc_kinetic_energy(Particle const &p1) { return res; } -/** Add bonded energies for one particle to the energy observable. - * @param[in] p particle for which to calculate energies - * @param[out] obs_energy energy observable - */ -inline void add_bonded_energy(Particle &p, Observable_stat &obs_energy) { - cell_structure.execute_bond_handler( - p, [&obs_energy](Particle &p1, int bond_id, - Utils::Span partners) { - auto const &iaparams = bonded_ia_params[bond_id]; - auto const result = calc_bonded_energy(iaparams, p1, partners); - if (result) { - obs_energy.bonded_contribution(bond_id)[0] += result.get(); - return false; - } - return true; - }); -} - #endif // ENERGY_INLINE_HPP diff --git a/src/core/forces.cpp b/src/core/forces.cpp index eaba1e80128..b616bdb81b6 100644 --- a/src/core/forces.cpp +++ b/src/core/forces.cpp @@ -117,7 +117,7 @@ void force_calc(CellStructure &cell_structure) { #endif short_range_loop( - [](Particle &p) { add_single_particle_force(p); }, + add_bonded_force, [](Particle &p1, Particle &p2, Distance const &d) { add_non_bonded_pair_force(p1, p2, d.vec21, sqrt(d.dist2), d.dist2); #ifdef COLLISION_DETECTION diff --git a/src/core/forces_inline.hpp b/src/core/forces_inline.hpp index 4f8dcb2a29c..04571b60941 100644 --- a/src/core/forces_inline.hpp +++ b/src/core/forces_inline.hpp @@ -512,10 +512,4 @@ inline bool add_bonded_force(Particle &p1, int bond_id, } } -/** Calculate bonded forces for one particle. - * @param p particle for which to calculate forces - */ -inline void add_single_particle_force(Particle &p) { - cell_structure.execute_bond_handler(p, add_bonded_force); -} #endif diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp index db90822ff55..50aa164364d 100644 --- a/src/core/pressure.cpp +++ b/src/core/pressure.cpp @@ -83,11 +83,27 @@ void pressure_calc() { add_kinetic_virials(p, obs_pressure); } - short_range_loop([](Particle &p) { add_bonded_virials(p, obs_pressure); }, - [](Particle &p1, Particle &p2, Distance const &d) { - add_non_bonded_pair_virials(p1, p2, d.vec21, sqrt(d.dist2), - obs_pressure); - }); + short_range_loop( + [](Particle &p1, int bond_id, + Utils::Span partners) { + auto const &iaparams = bonded_ia_params[bond_id]; + auto const result = calc_bonded_pressure_tensor(iaparams, p1, partners); + if (result) { + auto const &tensor = result.get(); + /* pressure tensor part */ + for (int k = 0; k < 3; k++) + for (int l = 0; l < 3; l++) + obs_pressure.bonded_contribution(bond_id)[k * 3 + l] += + tensor[k][l]; + + return false; + } + return true; + }, + [](Particle &p1, Particle &p2, Distance const &d) { + add_non_bonded_pair_virials(p1, p2, d.vec21, sqrt(d.dist2), + obs_pressure); + }, detail::True{}); calc_long_range_virials(cell_structure.local_particles()); diff --git a/src/core/pressure_inline.hpp b/src/core/pressure_inline.hpp index 3b4f7050c01..b77faa107c7 100644 --- a/src/core/pressure_inline.hpp +++ b/src/core/pressure_inline.hpp @@ -151,27 +151,4 @@ inline void add_kinetic_virials(Particle const &p1, obs_pressure.kinetic[k * 3 + l] += p1.m.v[k] * p1.m.v[l] * p1.p.mass; } -/** Add bonded energies for one particle to the energy observable. - * @param[in] p particle for which to calculate pressure - * @param[out] obs_pressure pressure observable - */ -inline void add_bonded_virials(Particle &p, Observable_stat &obs_pressure) { - cell_structure.execute_bond_handler(p, [&obs_pressure]( - Particle &p1, int bond_id, - Utils::Span partners) { - auto const &iaparams = bonded_ia_params[bond_id]; - auto const result = calc_bonded_pressure_tensor(iaparams, p1, partners); - if (result) { - auto const &tensor = result.get(); - /* pressure tensor part */ - for (int k = 0; k < 3; k++) - for (int l = 0; l < 3; l++) - obs_pressure.bonded_contribution(bond_id)[k * 3 + l] += tensor[k][l]; - - return false; - } - return true; - }); -} - #endif diff --git a/src/core/short_range_loop.hpp b/src/core/short_range_loop.hpp index 4e69f265ff6..71b428532e9 100644 --- a/src/core/short_range_loop.hpp +++ b/src/core/short_range_loop.hpp @@ -19,12 +19,10 @@ #ifndef CORE_SHORT_RANGE_HPP #define CORE_SHORT_RANGE_HPP -#include "algorithm/link_cell.hpp" #include "cells.hpp" #include "grid.hpp" #include "integrate.hpp" -#include #include #include @@ -62,90 +60,37 @@ struct EuclidianDistance { struct True { template bool operator()(T...) const { return true; } }; - -template -void pair_loop(PairKernel &&pair_kernel, DistanceFunction df, - const VerletCriterion &verlet_criterion) { - auto first = - boost::make_indirect_iterator(cell_structure.local_cells().begin()); - auto last = boost::make_indirect_iterator(cell_structure.local_cells().end()); - - if (cell_structure.use_verlet_list && cell_structure.m_rebuild_verlet_list) { - cell_structure.m_verlet_list.clear(); - - Algorithm::link_cell( - first, last, - [&pair_kernel, &df, &verlet_criterion](Particle &p1, Particle &p2) { - auto const d = df(p1, p2); - if (verlet_criterion(p1, p2, d)) { - cell_structure.m_verlet_list.emplace_back(&p1, &p2); - pair_kernel(p1, p2, d); - } - }); - cell_structure.m_rebuild_verlet_list = false; - } else if (cell_structure.use_verlet_list && - not cell_structure.m_rebuild_verlet_list) { - for (auto &pair : cell_structure.m_verlet_list) { - pair_kernel(*pair.first, *pair.second, df(*pair.first, *pair.second)); - } - } else { - Algorithm::link_cell( - first, last, - [&pair_kernel, &df, &verlet_criterion](Particle &p1, Particle &p2) { - auto const d = df(p1, p2); - if (verlet_criterion(p1, p2, d)) { - cell_structure.m_verlet_list.emplace_back(&p1, &p2); - pair_kernel(p1, p2, d); - } - }); - } -} - } // namespace detail template void short_range_loop(PairKernel &&pair_kernel, - const VerletCriterion &verlet_criterion = {}) { + const VerletCriterion &verlet_criterion) { ESPRESSO_PROFILER_CXX_MARK_FUNCTION; assert(cell_structure.get_resort_particles() == Cells::RESORT_NONE); if (interaction_range() != INACTIVE_CUTOFF) { if (cell_structure.decomposition().minimum_image_distance()) { - detail::pair_loop(std::forward(pair_kernel), - detail::MinimalImageDistance{box_geo}, - verlet_criterion); + cell_structure.pair_loop(std::forward(pair_kernel), + detail::MinimalImageDistance{box_geo}, + verlet_criterion); } else { - detail::pair_loop(std::forward(pair_kernel), - detail::EuclidianDistance{}, verlet_criterion); + cell_structure.pair_loop(std::forward(pair_kernel), + detail::EuclidianDistance{}, verlet_criterion); } } } -template -void short_range_loop(ParticleKernel &&particle_kernel, - PairKernel &&pair_kernel, - const VerletCriterion &verlet_criterion = {}) { +void short_range_loop(BondKernel bond_kernel, PairKernel pair_kernel, + const VerletCriterion &verlet_criterion) { ESPRESSO_PROFILER_CXX_MARK_FUNCTION; assert(cell_structure.get_resort_particles() == Cells::RESORT_NONE); - auto local_particles = cell_structure.local_particles(); - for (auto &p : local_particles) { - particle_kernel(p); - } + cell_structure.bond_loop(bond_kernel); - if (interaction_range() != INACTIVE_CUTOFF) { - if (cell_structure.decomposition().minimum_image_distance()) { - detail::pair_loop(std::forward(pair_kernel), - detail::MinimalImageDistance{box_geo}, - verlet_criterion); - } else { - detail::pair_loop(std::forward(pair_kernel), - detail::EuclidianDistance{}, verlet_criterion); - } - } + short_range_loop(pair_kernel, verlet_criterion); } - #endif From 1f8c6534b3883f15ff857ff2c49bec0121c8265e Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Thu, 23 Apr 2020 22:49:48 +0200 Subject: [PATCH 07/27] core: Clear particle node on CellStructure change --- src/core/event.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/core/event.cpp b/src/core/event.cpp index 5bc622c8904..48dd925d567 100644 --- a/src/core/event.cpp +++ b/src/core/event.cpp @@ -251,6 +251,8 @@ void on_boxl_change() { } void on_cell_structure_change() { + clear_particle_node(); + /* Now give methods a chance to react to the change in cell structure. Most ES methods need to reinitialize, as they depend on skin, node grid and so on. */ From 99990f869311ab3d0622d915d67af3170d3e0c9b Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Mon, 27 Apr 2020 14:51:18 +0200 Subject: [PATCH 08/27] core: cells: Hide implementation details Removed final use of CellStructure::execute_bond_handler(), made it private. --- src/core/CellStructure.hpp | 4 ++-- src/core/dpd.cpp | 2 +- src/core/electrostatics_magnetostatics/icc.cpp | 2 +- .../immersed_boundary/ImmersedBoundaries.cpp | 16 +++++----------- 4 files changed, 9 insertions(+), 15 deletions(-) diff --git a/src/core/CellStructure.hpp b/src/core/CellStructure.hpp index 8e41a7aebba..707790b9914 100644 --- a/src/core/CellStructure.hpp +++ b/src/core/CellStructure.hpp @@ -331,6 +331,7 @@ struct CellStructure { * Cells::DataPart */ void ghosts_update(unsigned data_parts); + /** * @brief Add forces from ghost particles to real particles. */ @@ -359,7 +360,6 @@ struct CellStructure { return partners; } -public: /** * @brief Execute kernel for every bond on particle. * @tparam Handler Callable, which can be invoked with @@ -391,7 +391,6 @@ struct CellStructure { } } -private: /** Go through ghost cells and remove the ghost entries from the local particle index. */ void invalidate_ghosts() { @@ -454,6 +453,7 @@ struct CellStructure { return m_decomposition->minimum_image_distance(); } +public: template void bond_loop(BondKernel const &bond_kernel) { for (auto &p : local_particles()) { execute_bond_handler(p, bond_kernel); diff --git a/src/core/dpd.cpp b/src/core/dpd.cpp index 1f064e48c5a..e288879400a 100644 --- a/src/core/dpd.cpp +++ b/src/core/dpd.cpp @@ -160,7 +160,7 @@ static auto dpd_viscous_stress_local() { stress += tensor_product(d.vec21, f); }, - detail::True{}); + {}); return stress; } diff --git a/src/core/electrostatics_magnetostatics/icc.cpp b/src/core/electrostatics_magnetostatics/icc.cpp index 27f7735a332..40be1caf448 100644 --- a/src/core/electrostatics_magnetostatics/icc.cpp +++ b/src/core/electrostatics_magnetostatics/icc.cpp @@ -203,7 +203,7 @@ void force_calc_iccp3m(const ParticleRange &particles, add_non_bonded_pair_force_iccp3m(p1, p2, d.vec21, sqrt(d.dist2), d.dist2); }, - detail::True{}); + {}); Coulomb::calc_long_range_force(particles); } diff --git a/src/core/immersed_boundary/ImmersedBoundaries.cpp b/src/core/immersed_boundary/ImmersedBoundaries.cpp index 943f7cbf817..5dc64c7003a 100644 --- a/src/core/immersed_boundary/ImmersedBoundaries.cpp +++ b/src/core/immersed_boundary/ImmersedBoundaries.cpp @@ -134,16 +134,12 @@ void ImmersedBoundaries::calc_volumes(CellStructure &cs) { std::vector tempVol(MaxNumIBM); // Loop over all particles on local node - for (auto &p1 : cs.local_particles()) { - auto vol_cons_params = vol_cons_parameters(p1); - - if (vol_cons_params) { - cs.execute_bond_handler(p1, [softID = vol_cons_params->softID, - &tempVol](Particle &p1, int bond_id, - Utils::Span partners) { + cs.bond_loop( + [&tempVol](Particle &p1, int bond_id, Utils::Span partners) { auto const &iaparams = bonded_ia_params[bond_id]; + auto vol_cons_params = vol_cons_parameters(p1); - if (iaparams.type == BONDED_IA_IBM_TRIEL) { + if (vol_cons_params && iaparams.type == BONDED_IA_IBM_TRIEL) { // Our particle is the leading particle of a triel // Get second and third particle of the triangle Particle &p2 = *partners[0]; @@ -173,13 +169,11 @@ void ImmersedBoundaries::calc_volumes(CellStructure &cs) { const double v213 = x2[0] * x1[1] * x3[2]; const double v123 = x1[0] * x2[1] * x3[2]; - tempVol[softID] += + tempVol[vol_cons_params->softID] += 1.0 / 6.0 * (-v321 + v231 + v312 - v132 - v213 + v123); } return false; }); - } - } for (int i = 0; i < MaxNumIBM; i++) VolumesCurrent[i] = 0; From 1a4cfec9da930c0e6ed1edde8747be42d4afb062 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Mon, 27 Apr 2020 15:58:33 +0200 Subject: [PATCH 09/27] core: Style and documentation --- src/core/BoxGeometry.hpp | 4 ++++ src/core/CellStructure.cpp | 2 +- src/core/CellStructure.hpp | 6 ++++-- src/core/DomainDecomposition.cpp | 1 + src/core/DomainDecomposition.hpp | 1 - 5 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/core/BoxGeometry.hpp b/src/core/BoxGeometry.hpp index 5083447a322..ed8487cd5df 100644 --- a/src/core/BoxGeometry.hpp +++ b/src/core/BoxGeometry.hpp @@ -25,6 +25,7 @@ #include #include +#include class BoxGeometry { public: @@ -92,6 +93,9 @@ template T get_mi_coord(T a, T b, T box_length, bool periodic) { /** * @brief Get the minimum-image vector between two coordinates. + * + * @tparam T Floating point type. + * * @param a Coordinate of the terminal point. * @param b Coordinate of the initial point. * @param box Box parameters (side lengths, periodicity). diff --git a/src/core/CellStructure.cpp b/src/core/CellStructure.cpp index 5a5b2a23892..d94817e8309 100644 --- a/src/core/CellStructure.cpp +++ b/src/core/CellStructure.cpp @@ -182,4 +182,4 @@ void CellStructure::set_domain_decomposition( set_particle_decomposition( std::make_unique(comm, range, box, local_geo)); m_type = CELL_STRUCTURE_DOMDEC; -} \ No newline at end of file +} diff --git a/src/core/CellStructure.hpp b/src/core/CellStructure.hpp index 707790b9914..1afbcfef02b 100644 --- a/src/core/CellStructure.hpp +++ b/src/core/CellStructure.hpp @@ -391,8 +391,10 @@ struct CellStructure { } } - /** Go through ghost cells and remove the ghost entries from the - local particle index. */ + /** + * @brief Go through ghost cells and remove the ghost entries from the + * local particle index. + */ void invalidate_ghosts() { for (auto const &p : ghost_particles()) { if (get_local_particle(p.identity()) == &p) { diff --git a/src/core/DomainDecomposition.cpp b/src/core/DomainDecomposition.cpp index 748c2fa7b41..7dafdace213 100644 --- a/src/core/DomainDecomposition.cpp +++ b/src/core/DomainDecomposition.cpp @@ -224,6 +224,7 @@ void DomainDecomposition::resort(bool global, } } } + void DomainDecomposition::mark_cells() { int cnt_c = 0; diff --git a/src/core/DomainDecomposition.hpp b/src/core/DomainDecomposition.hpp index 617beaf5582..990091cbcdd 100644 --- a/src/core/DomainDecomposition.hpp +++ b/src/core/DomainDecomposition.hpp @@ -81,7 +81,6 @@ struct DomainDecomposition : public ParticleDecomposition { const BoxGeometry &box_geo, const LocalBox &local_geo); -public: GhostCommunicator const &exchange_ghosts_comm() const override { return m_exchange_ghosts_comm; } From 104a722d8593f3825dafdef54087b2c5e0c575d5 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Sun, 3 May 2020 16:01:00 +0200 Subject: [PATCH 10/27] core: cells: Move distance functions --- src/core/CellStructure.hpp | 26 ++++++++++++++++++++++++++ src/core/short_range_loop.hpp | 25 ------------------------- 2 files changed, 26 insertions(+), 25 deletions(-) diff --git a/src/core/CellStructure.hpp b/src/core/CellStructure.hpp index 1afbcfef02b..fcd4f61cd15 100644 --- a/src/core/CellStructure.hpp +++ b/src/core/CellStructure.hpp @@ -82,6 +82,32 @@ inline ParticleRange particles(Utils::Span cells) { } } // namespace Cells +/** + * @brief Distance vector and length handed to pair kernels. + */ +struct Distance { + explicit Distance(Utils::Vector3d const &vec21) + : vec21(vec21), dist2(vec21.norm2()) {} + + Utils::Vector3d vec21; + double dist2; +}; +namespace detail { +struct MinimalImageDistance { + const BoxGeometry box; + + Distance operator()(Particle const &p1, Particle const &p2) const { + return Distance(get_mi_vector(p1.r.p, p2.r.p, box)); + } +}; + +struct EuclidianDistance { + Distance operator()(Particle const &p1, Particle const &p2) const { + return Distance(p1.r.p - p2.r.p); + } +}; +} // namespace detail + /** Describes a cell structure / cell system. Contains information * about the communication of cell contents (particles, ghosts, ...) * between different nodes and the relation between particle diff --git a/src/core/short_range_loop.hpp b/src/core/short_range_loop.hpp index 71b428532e9..cde62f526ef 100644 --- a/src/core/short_range_loop.hpp +++ b/src/core/short_range_loop.hpp @@ -27,32 +27,7 @@ #include -/** - * @brief Distance vector and length handed to pair kernels. - */ -struct Distance { - explicit Distance(Utils::Vector3d const &vec21) - : vec21(vec21), dist2(vec21.norm2()) {} - - Utils::Vector3d vec21; - double dist2; -}; - namespace detail { -struct MinimalImageDistance { - const BoxGeometry box; - - Distance operator()(Particle const &p1, Particle const &p2) const { - return Distance(get_mi_vector(p1.r.p, p2.r.p, box)); - } -}; - -struct EuclidianDistance { - Distance operator()(Particle const &p1, Particle const &p2) const { - return Distance(p1.r.p - p2.r.p); - } -}; - /** * @brief Functor that returns true for * any arguments. From 9446ffb667510ecb806ff0b488f45e955dcde420 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Sun, 3 May 2020 16:21:21 +0200 Subject: [PATCH 11/27] core: cells: Return box object --- src/core/AtomDecomposition.hpp | 4 +++- src/core/CellStructure.hpp | 7 ------- src/core/DomainDecomposition.hpp | 5 ++++- src/core/ParticleDecomposition.hpp | 7 ++++--- 4 files changed, 11 insertions(+), 12 deletions(-) diff --git a/src/core/AtomDecomposition.hpp b/src/core/AtomDecomposition.hpp index 8961fdbfd45..ce1c89c10d0 100644 --- a/src/core/AtomDecomposition.hpp +++ b/src/core/AtomDecomposition.hpp @@ -90,7 +90,9 @@ class AtomDecomposition : public ParticleDecomposition { Utils::Vector3d max_range() const override; /* Return true if minimum image convention is * needed for distance calculation. */ - bool minimum_image_distance() const override { return true; } + boost::optional minimum_image_distance() const override { + return m_box; + } private: /** diff --git a/src/core/CellStructure.hpp b/src/core/CellStructure.hpp index fcd4f61cd15..43c106e2d9f 100644 --- a/src/core/CellStructure.hpp +++ b/src/core/CellStructure.hpp @@ -474,13 +474,6 @@ struct CellStructure { double range, BoxGeometry const &box, LocalBox const &local_geo); - /** - * @brief Return true if minimum image convention is - * needed for distance calculation. */ - bool minimum_image_distance() const { - return m_decomposition->minimum_image_distance(); - } - public: template void bond_loop(BondKernel const &bond_kernel) { for (auto &p : local_particles()) { diff --git a/src/core/DomainDecomposition.hpp b/src/core/DomainDecomposition.hpp index 990091cbcdd..6577da19080 100644 --- a/src/core/DomainDecomposition.hpp +++ b/src/core/DomainDecomposition.hpp @@ -99,10 +99,13 @@ struct DomainDecomposition : public ParticleDecomposition { return position_to_cell(p.r.p); } - bool minimum_image_distance() const override { return false; } void resort(bool global, std::vector &diff) override; Utils::Vector3d max_range() const override; + boost::optional minimum_image_distance() const override { + return {}; + } + private: /** Fill local_cells list and ghost_cells list for use with domain * decomposition. diff --git a/src/core/ParticleDecomposition.hpp b/src/core/ParticleDecomposition.hpp index 7fd0a282b39..f06ff2ee168 100644 --- a/src/core/ParticleDecomposition.hpp +++ b/src/core/ParticleDecomposition.hpp @@ -114,10 +114,11 @@ class ParticleDecomposition { */ virtual Utils::Vector3d max_range() const = 0; /** - * @brief Return true if minimum image convention is - * needed for distance calculation. + * @brief Return the box geometry needed for distance calculation + * if minimum image convention should be used needed for + * distance calculation. */ - virtual bool minimum_image_distance() const = 0; + virtual boost::optional minimum_image_distance() const = 0; virtual ~ParticleDecomposition() = default; }; From da9068d71a47bd13cd54b3cb13cf02ed08d4e295 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Sun, 3 May 2020 16:49:15 +0200 Subject: [PATCH 12/27] core: short_range_loop: Pull distance function out Delegate the DistanceFunction logic to the CellStructure. --- src/core/CellStructure.hpp | 75 +++++++++++++++++++++++++---------- src/core/short_range_loop.hpp | 11 +---- 2 files changed, 57 insertions(+), 29 deletions(-) diff --git a/src/core/CellStructure.hpp b/src/core/CellStructure.hpp index 43c106e2d9f..420e86ca9c6 100644 --- a/src/core/CellStructure.hpp +++ b/src/core/CellStructure.hpp @@ -481,34 +481,69 @@ struct CellStructure { } } - template - void pair_loop(PairKernel &&pair_kernel, DistanceFunction df, - const VerletCriterion &verlet_criterion) { - auto first = boost::make_indirect_iterator(local_cells().begin()); - auto last = boost::make_indirect_iterator(local_cells().end()); +private: + /** + * @brief Run link_cell algorithm for local cells. + * + * @tparam Kernel Needs to be callable with (Particle, Particle, Distance). + * @param kernel Pair kernel functor. + */ + template void link_cell(Kernel kernel) { + auto const maybe_box = decomposition().minimum_image_distance(); + auto const first = boost::make_indirect_iterator(local_cells().begin()); + auto const last = boost::make_indirect_iterator(local_cells().end()); + + if (maybe_box) { + Algorithm::link_cell( + first, last, + [&kernel, df = detail::MinimalImageDistance{*maybe_box}]( + Particle &p1, Particle &p2) { kernel(p1, p2, df(p1, p2)); }); + } else { + Algorithm::link_cell( + first, last, + [&kernel, df = detail::EuclidianDistance{}]( + Particle &p1, Particle &p2) { kernel(p1, p2, df(p1, p2)); }); + } + } +public: + template + void pair_loop(PairKernel &&pair_kernel, + const VerletCriterion &verlet_criterion) { + /* In this case the verlet list update is attached to + * the pair kernel, and the verlet list is rebuilt as + * we go. */ if (use_verlet_list && m_rebuild_verlet_list) { m_verlet_list.clear(); - Algorithm::link_cell(first, last, - [&pair_kernel, &df, &verlet_criterion, - this](Particle &p1, Particle &p2) { - auto const d = df(p1, p2); - if (verlet_criterion(p1, p2, d)) { - m_verlet_list.emplace_back(&p1, &p2); - pair_kernel(p1, p2, d); - } - }); + link_cell([&pair_kernel, &verlet_criterion, + this](Particle &p1, Particle &p2, Distance const &d) { + if (verlet_criterion(p1, p2, d)) { + m_verlet_list.emplace_back(&p1, &p2); + pair_kernel(p1, p2, d); + } + }); + m_rebuild_verlet_list = false; } else if (use_verlet_list && not m_rebuild_verlet_list) { - for (auto &pair : m_verlet_list) { - pair_kernel(*pair.first, *pair.second, df(*pair.first, *pair.second)); + auto const maybe_box = decomposition().minimum_image_distance(); + /* In this case the pair kernel is just run over the verlet list. */ + if (maybe_box) { + auto const distance_function = detail::MinimalImageDistance{*maybe_box}; + for (auto &pair : m_verlet_list) { + pair_kernel(*pair.first, *pair.second, + distance_function(*pair.first, *pair.second)); + } + } else { + auto const distance_function = detail::EuclidianDistance{}; + for (auto &pair : m_verlet_list) { + pair_kernel(*pair.first, *pair.second, + distance_function(*pair.first, *pair.second)); + } } } else { - Algorithm::link_cell(first, last, - [&pair_kernel, &df](Particle &p1, Particle &p2) { - pair_kernel(p1, p2, df(p1, p2)); - }); + /* No verlet lists, just run the kernel with pairs from the cells. */ + link_cell(pair_kernel); } } }; diff --git a/src/core/short_range_loop.hpp b/src/core/short_range_loop.hpp index cde62f526ef..4d4fc860937 100644 --- a/src/core/short_range_loop.hpp +++ b/src/core/short_range_loop.hpp @@ -38,21 +38,14 @@ struct True { } // namespace detail template -void short_range_loop(PairKernel &&pair_kernel, +void short_range_loop(PairKernel pair_kernel, const VerletCriterion &verlet_criterion) { ESPRESSO_PROFILER_CXX_MARK_FUNCTION; assert(cell_structure.get_resort_particles() == Cells::RESORT_NONE); if (interaction_range() != INACTIVE_CUTOFF) { - if (cell_structure.decomposition().minimum_image_distance()) { - cell_structure.pair_loop(std::forward(pair_kernel), - detail::MinimalImageDistance{box_geo}, - verlet_criterion); - } else { - cell_structure.pair_loop(std::forward(pair_kernel), - detail::EuclidianDistance{}, verlet_criterion); - } + cell_structure.pair_loop(pair_kernel, verlet_criterion); } } From e0ab386c2422c5475f5e490df34a5a7ae61dd969 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Sun, 3 May 2020 21:58:53 +0200 Subject: [PATCH 13/27] core: cells: Split CellStructure::pair_loop() Rename it to CellStructure::non_bonded_loop() and create a specialization without Verlet list. Pass CellStructure refs as function arguments instead of relying on a global variable. --- src/core/CellStructure.hpp | 17 +++++++++++++++-- src/core/cells.cpp | 4 ++-- src/core/debug.cpp | 22 ++++++++++------------ src/core/debug.hpp | 5 +++-- src/core/pressure.cpp | 6 +++--- src/core/short_range_loop.hpp | 2 +- 6 files changed, 34 insertions(+), 22 deletions(-) diff --git a/src/core/CellStructure.hpp b/src/core/CellStructure.hpp index 420e86ca9c6..7b4a2df08d5 100644 --- a/src/core/CellStructure.hpp +++ b/src/core/CellStructure.hpp @@ -507,9 +507,22 @@ struct CellStructure { } public: + /** Non-bonded pair loop with potential use + * of verlet lists. + * @param pair_kernel Kernel to apply + */ + template void non_bonded_loop(PairKernel pair_kernel) { + link_cell(pair_kernel); + } + + /** Non-bonded pair loop with potential use + * of verlet lists. + * @param pair_kernel Kernel to apply + * @param verlet_criterion Filter for verlet lists. + */ template - void pair_loop(PairKernel &&pair_kernel, - const VerletCriterion &verlet_criterion) { + void non_bonded_loop(PairKernel &&pair_kernel, + const VerletCriterion &verlet_criterion) { /* In this case the verlet list update is attached to * the pair kernel, and the verlet list is rebuilt as * we go. */ diff --git a/src/core/cells.cpp b/src/core/cells.cpp index d72e6b32eb8..d049f3c97c0 100644 --- a/src/core/cells.cpp +++ b/src/core/cells.cpp @@ -138,8 +138,8 @@ void cells_resort_particles(int global_flag) { #ifdef ADDITIONAL_CHECKS /* at the end of the day, everything should be consistent again */ - check_particle_consistency(); - check_particle_sorting(); + check_particle_consistency(cell_structure); + check_particle_sorting(cell_structure); #endif } diff --git a/src/core/debug.cpp b/src/core/debug.cpp index fbbb54b6f71..0174883d3e4 100644 --- a/src/core/debug.cpp +++ b/src/core/debug.cpp @@ -24,15 +24,13 @@ #include "debug.hpp" -#include "cells.hpp" - #include -void check_particle_consistency() { - auto local_particles = cell_structure.local_particles(); +void check_particle_consistency(CellStructure &cs) { + auto local_particles = cs.local_particles(); auto const cell_part_cnt = local_particles.size(); - auto const max_id = cell_structure.get_max_local_particle_id(); + auto const max_id = cs.get_max_local_particle_id(); for (auto const &p : local_particles) { auto const id = p.identity(); @@ -41,17 +39,17 @@ void check_particle_consistency() { throw std::runtime_error("Particle id out of bounds."); } - if (cell_structure.get_local_particle(id) != &p) { + if (cs.get_local_particle(id) != &p) { throw std::runtime_error("Invalid local particle index entry."); } } /* checks: local particle id */ int local_part_cnt = 0; - for (int n = 0; n < cell_structure.get_max_local_particle_id() + 1; n++) { - if (cell_structure.get_local_particle(n) != nullptr) { + for (int n = 0; n < cs.get_max_local_particle_id() + 1; n++) { + if (cs.get_local_particle(n) != nullptr) { local_part_cnt++; - if (cell_structure.get_local_particle(n)->p.identity != n) { + if (cs.get_local_particle(n)->p.identity != n) { throw std::runtime_error("local_particles part has corrupted id."); } } @@ -64,10 +62,10 @@ void check_particle_consistency() { } } -void check_particle_sorting() { - for (auto cell : cell_structure.local_cells()) { +void check_particle_sorting(CellStructure &cs) { + for (auto cell : cs.local_cells()) { for (auto const &p : cell->particles()) { - if (cell_structure.particle_to_cell(p) != cell) { + if (cs.particle_to_cell(p) != cell) { throw std::runtime_error("misplaced particle with id " + std::to_string(p.identity())); } diff --git a/src/core/debug.hpp b/src/core/debug.hpp index 2d48d12852c..bba556d69bc 100644 --- a/src/core/debug.hpp +++ b/src/core/debug.hpp @@ -26,10 +26,11 @@ #ifndef DEBUG_HPP #define DEBUG_HPP +#include "CellStructure.hpp" /** this performs a lot of tests which will very likely detect corruptions of * the cell structure. */ -void check_particle_consistency(); +void check_particle_consistency(CellStructure &cs); /** * @brief Check if particles are in correct cells. @@ -37,6 +38,6 @@ void check_particle_consistency(); * Check that particles are in the cells the cellsystem says * they should be. */ -void check_particle_sorting(); +void check_particle_sorting(CellStructure &cs); #endif diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp index 50aa164364d..b7ce0d07688 100644 --- a/src/core/pressure.cpp +++ b/src/core/pressure.cpp @@ -84,8 +84,7 @@ void pressure_calc() { } short_range_loop( - [](Particle &p1, int bond_id, - Utils::Span partners) { + [](Particle &p1, int bond_id, Utils::Span partners) { auto const &iaparams = bonded_ia_params[bond_id]; auto const result = calc_bonded_pressure_tensor(iaparams, p1, partners); if (result) { @@ -103,7 +102,8 @@ void pressure_calc() { [](Particle &p1, Particle &p2, Distance const &d) { add_non_bonded_pair_virials(p1, p2, d.vec21, sqrt(d.dist2), obs_pressure); - }, detail::True{}); + }, + detail::True{}); calc_long_range_virials(cell_structure.local_particles()); diff --git a/src/core/short_range_loop.hpp b/src/core/short_range_loop.hpp index 4d4fc860937..3e88f66efc8 100644 --- a/src/core/short_range_loop.hpp +++ b/src/core/short_range_loop.hpp @@ -45,7 +45,7 @@ void short_range_loop(PairKernel pair_kernel, assert(cell_structure.get_resort_particles() == Cells::RESORT_NONE); if (interaction_range() != INACTIVE_CUTOFF) { - cell_structure.pair_loop(pair_kernel, verlet_criterion); + cell_structure.non_bonded_loop(pair_kernel, verlet_criterion); } } From bf00e4912e850a8a1c80275d459b675b301e49f8 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Sun, 3 May 2020 22:14:32 +0200 Subject: [PATCH 14/27] core: short_range_loop: Factor out non-bonded loop Delegate non-bonded kernels to CellStructure::non_bonded_loop(). --- src/core/EspressoSystemInterface.hpp | 1 - src/core/cells.cpp | 2 +- src/core/cuda_common_cuda.cu | 2 -- src/core/cuda_init_cuda.cu | 1 - src/core/dpd.cpp | 5 ++--- src/core/electrostatics_magnetostatics/icc.cpp | 12 +++++------- src/core/grid_based_algorithms/lbgpu_cuda.cu | 1 - src/core/particle_data.cpp | 1 - src/core/short_range_loop.hpp | 15 +-------------- 9 files changed, 9 insertions(+), 31 deletions(-) diff --git a/src/core/EspressoSystemInterface.hpp b/src/core/EspressoSystemInterface.hpp index 2a23fe80d67..b6184fd2f35 100644 --- a/src/core/EspressoSystemInterface.hpp +++ b/src/core/EspressoSystemInterface.hpp @@ -21,7 +21,6 @@ #include "SystemInterface.hpp" #include "cuda_interface.hpp" -#include "debug.hpp" /** Syntactic sugar */ #define espressoSystemInterface EspressoSystemInterface::Instance() diff --git a/src/core/cells.cpp b/src/core/cells.cpp index d049f3c97c0..ad3c514c964 100644 --- a/src/core/cells.cpp +++ b/src/core/cells.cpp @@ -69,7 +69,7 @@ std::vector> get_pairs(double distance) { ret.emplace_back(p1.p.identity, p2.p.identity); }; - short_range_loop(pair_kernel, detail::True{}); + cell_structure.non_bonded_loop(pair_kernel); /* Sort pairs */ for (auto &pair : ret) { diff --git a/src/core/cuda_common_cuda.cu b/src/core/cuda_common_cuda.cu index 02b36554453..7aa9b8c8e18 100644 --- a/src/core/cuda_common_cuda.cu +++ b/src/core/cuda_common_cuda.cu @@ -20,8 +20,6 @@ #include "config.hpp" -#include "debug.hpp" - #include "ParticleRange.hpp" #include "cuda_init.hpp" #include "cuda_interface.hpp" diff --git a/src/core/cuda_init_cuda.cu b/src/core/cuda_init_cuda.cu index 43bdf2cc8eb..a71ed74b3c0 100644 --- a/src/core/cuda_init_cuda.cu +++ b/src/core/cuda_init_cuda.cu @@ -21,7 +21,6 @@ #include "cuda_init.hpp" #include "cuda_utils.hpp" -#include "debug.hpp" #include diff --git a/src/core/dpd.cpp b/src/core/dpd.cpp index e288879400a..9df44a1e628 100644 --- a/src/core/dpd.cpp +++ b/src/core/dpd.cpp @@ -142,7 +142,7 @@ static auto dpd_viscous_stress_local() { on_observable_calc(); Utils::Vector stress{}; - short_range_loop( + cell_structure.non_bonded_loop( [&stress](const Particle &p1, const Particle &p2, Distance const &d) { auto const v21 = p1.m.v - p2.m.v; @@ -159,8 +159,7 @@ static auto dpd_viscous_stress_local() { auto const f = P * (f_r - f_t) + f_t; stress += tensor_product(d.vec21, f); - }, - {}); + }); return stress; } diff --git a/src/core/electrostatics_magnetostatics/icc.cpp b/src/core/electrostatics_magnetostatics/icc.cpp index 40be1caf448..73128a7bf07 100644 --- a/src/core/electrostatics_magnetostatics/icc.cpp +++ b/src/core/electrostatics_magnetostatics/icc.cpp @@ -197,13 +197,11 @@ void force_calc_iccp3m(const ParticleRange &particles, const ParticleRange &ghost_particles) { init_forces_iccp3m(particles, ghost_particles); - short_range_loop( - [](Particle &p1, Particle &p2, Distance const &d) { - /* calc non bonded interactions */ - add_non_bonded_pair_force_iccp3m(p1, p2, d.vec21, sqrt(d.dist2), - d.dist2); - }, - {}); + cell_structure.non_bonded_loop([](Particle &p1, Particle &p2, + Distance const &d) { + /* calc non bonded interactions */ + add_non_bonded_pair_force_iccp3m(p1, p2, d.vec21, sqrt(d.dist2), d.dist2); + }); Coulomb::calc_long_range_force(particles); } diff --git a/src/core/grid_based_algorithms/lbgpu_cuda.cu b/src/core/grid_based_algorithms/lbgpu_cuda.cu index 838590e9f96..c133bf97b3b 100644 --- a/src/core/grid_based_algorithms/lbgpu_cuda.cu +++ b/src/core/grid_based_algorithms/lbgpu_cuda.cu @@ -37,7 +37,6 @@ extern int this_node; #include "cuda_interface.hpp" #include "cuda_utils.hpp" -#include "debug.hpp" #include "errorhandling.hpp" #include "grid_based_algorithms/OptionalCounter.hpp" diff --git a/src/core/particle_data.cpp b/src/core/particle_data.cpp index b5a63663f49..8f9b875f2a2 100644 --- a/src/core/particle_data.cpp +++ b/src/core/particle_data.cpp @@ -28,7 +28,6 @@ #include "Particle.hpp" #include "cells.hpp" #include "communication.hpp" -#include "debug.hpp" #include "event.hpp" #include "grid.hpp" #include "nonbonded_interactions/nonbonded_interaction_data.hpp" diff --git a/src/core/short_range_loop.hpp b/src/core/short_range_loop.hpp index 3e88f66efc8..8fd46910b25 100644 --- a/src/core/short_range_loop.hpp +++ b/src/core/short_range_loop.hpp @@ -37,18 +37,6 @@ struct True { }; } // namespace detail -template -void short_range_loop(PairKernel pair_kernel, - const VerletCriterion &verlet_criterion) { - ESPRESSO_PROFILER_CXX_MARK_FUNCTION; - - assert(cell_structure.get_resort_particles() == Cells::RESORT_NONE); - - if (interaction_range() != INACTIVE_CUTOFF) { - cell_structure.non_bonded_loop(pair_kernel, verlet_criterion); - } -} - template void short_range_loop(BondKernel bond_kernel, PairKernel pair_kernel, @@ -58,7 +46,6 @@ void short_range_loop(BondKernel bond_kernel, PairKernel pair_kernel, assert(cell_structure.get_resort_particles() == Cells::RESORT_NONE); cell_structure.bond_loop(bond_kernel); - - short_range_loop(pair_kernel, verlet_criterion); + cell_structure.non_bonded_loop(pair_kernel, verlet_criterion); } #endif From e16556eadbca6413ca724d5ed18d3b0d6430d9cc Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Thu, 2 Jul 2020 22:42:54 +0200 Subject: [PATCH 15/27] core: Header cleanup --- src/core/dpd.cpp | 1 + src/core/short_range_loop.hpp | 4 +--- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/src/core/dpd.cpp b/src/core/dpd.cpp index 9df44a1e628..a5c7b0bb93d 100644 --- a/src/core/dpd.cpp +++ b/src/core/dpd.cpp @@ -26,6 +26,7 @@ #ifdef DPD #include "communication.hpp" #include "event.hpp" +#include "grid.hpp" #include "integrate.hpp" #include "nonbonded_interactions/nonbonded_interaction_data.hpp" #include "random.hpp" diff --git a/src/core/short_range_loop.hpp b/src/core/short_range_loop.hpp index 8fd46910b25..ef2d8d460d0 100644 --- a/src/core/short_range_loop.hpp +++ b/src/core/short_range_loop.hpp @@ -20,12 +20,10 @@ #define CORE_SHORT_RANGE_HPP #include "cells.hpp" -#include "grid.hpp" -#include "integrate.hpp" #include -#include +#include namespace detail { /** From b535757793e568464c32f5e6d929a55541fa251e Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Thu, 2 Jul 2020 23:03:22 +0200 Subject: [PATCH 16/27] core: short_range_loop: Use default argument --- src/core/energy.cpp | 3 +-- src/core/pressure.cpp | 3 +-- src/core/short_range_loop.hpp | 2 +- 3 files changed, 3 insertions(+), 5 deletions(-) diff --git a/src/core/energy.cpp b/src/core/energy.cpp index 99a14171994..8df39a5b9a3 100644 --- a/src/core/energy.cpp +++ b/src/core/energy.cpp @@ -83,8 +83,7 @@ void energy_calc(const double time) { [](Particle const &p1, Particle const &p2, Distance const &d) { add_non_bonded_pair_energy(p1, p2, d.vec21, sqrt(d.dist2), d.dist2, obs_energy); - }, - detail::True{}); + }); calc_long_range_energies(cell_structure.local_particles()); diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp index b7ce0d07688..d3b1b3027ef 100644 --- a/src/core/pressure.cpp +++ b/src/core/pressure.cpp @@ -102,8 +102,7 @@ void pressure_calc() { [](Particle &p1, Particle &p2, Distance const &d) { add_non_bonded_pair_virials(p1, p2, d.vec21, sqrt(d.dist2), obs_pressure); - }, - detail::True{}); + }); calc_long_range_virials(cell_structure.local_particles()); diff --git a/src/core/short_range_loop.hpp b/src/core/short_range_loop.hpp index ef2d8d460d0..937bcc57331 100644 --- a/src/core/short_range_loop.hpp +++ b/src/core/short_range_loop.hpp @@ -38,7 +38,7 @@ struct True { template void short_range_loop(BondKernel bond_kernel, PairKernel pair_kernel, - const VerletCriterion &verlet_criterion) { + const VerletCriterion &verlet_criterion = {}) { ESPRESSO_PROFILER_CXX_MARK_FUNCTION; assert(cell_structure.get_resort_particles() == Cells::RESORT_NONE); From 605bab2a3c6644db580db0e18c8996dc670a0828 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Mon, 13 Jul 2020 11:08:08 +0200 Subject: [PATCH 17/27] core: cells: Moved checking functions to CellStructure --- src/core/CMakeLists.txt | 1 - src/core/CellStructure.cpp | 51 ++++++++++++++++++++++++++ src/core/CellStructure.hpp | 4 +++ src/core/cells.cpp | 7 ---- src/core/debug.cpp | 74 -------------------------------------- src/core/debug.hpp | 43 ---------------------- 6 files changed, 55 insertions(+), 125 deletions(-) delete mode 100644 src/core/debug.cpp delete mode 100644 src/core/debug.hpp diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index 034e3f4468f..c3b992ccf0f 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -6,7 +6,6 @@ set(EspressoCore_SRC comfixed_global.cpp communication.cpp constraints.cpp - debug.cpp dpd.cpp energy.cpp errorhandling.cpp diff --git a/src/core/CellStructure.cpp b/src/core/CellStructure.cpp index d94817e8309..515bd6d9aa3 100644 --- a/src/core/CellStructure.cpp +++ b/src/core/CellStructure.cpp @@ -26,6 +26,52 @@ #include +#include + +void CellStructure::check_particle_consistency() { + auto const max_id = get_max_local_particle_id(); + + for (auto const &p : local_particles()) { + auto const id = p.identity(); + + if (id < 0 || id > max_id) { + throw std::runtime_error("Particle id out of bounds."); + } + + if (get_local_particle(id) != &p) { + throw std::runtime_error("Invalid local particle index entry."); + } + } + + /* checks: local particle id */ + int local_part_cnt = 0; + for (int n = 0; n < get_max_local_particle_id() + 1; n++) { + if (get_local_particle(n) != nullptr) { + local_part_cnt++; + if (get_local_particle(n)->p.identity != n) { + throw std::runtime_error("local_particles part has corrupted id."); + } + } + } + + if (local_part_cnt != local_particles().size()) { + throw std::runtime_error( + std::to_string(local_particles().size()) + " parts in cells but " + + std::to_string(local_part_cnt) + " parts in local_particles"); + } +} + +void CellStructure::check_particle_sorting() { + for (auto cell : local_cells()) { + for (auto const &p : cell->particles()) { + if (particle_to_cell(p) != cell) { + throw std::runtime_error("misplaced particle with id " + + std::to_string(p.identity())); + } + } + } +} + Cell *CellStructure::particle_to_cell(const Particle &p) { return decomposition().particle_to_cell(p); } @@ -168,6 +214,11 @@ void CellStructure::resort_particles(int global_flag) { } m_rebuild_verlet_list = true; + +#ifdef ADDITIONAL_CHECKS + check_particle_consistency(); + check_particle_sorting(); +#endif } void CellStructure::set_atom_decomposition(boost::mpi::communicator const &comm, diff --git a/src/core/CellStructure.hpp b/src/core/CellStructure.hpp index 7b4a2df08d5..9d66010398d 100644 --- a/src/core/CellStructure.hpp +++ b/src/core/CellStructure.hpp @@ -559,6 +559,10 @@ struct CellStructure { link_cell(pair_kernel); } } + +private: + void check_particle_consistency(); + void check_particle_sorting(); }; #endif // ESPRESSO_CELLSTRUCTURE_HPP diff --git a/src/core/cells.cpp b/src/core/cells.cpp index ad3c514c964..4f4fbb96e14 100644 --- a/src/core/cells.cpp +++ b/src/core/cells.cpp @@ -27,7 +27,6 @@ #include "cells.hpp" #include "Particle.hpp" #include "communication.hpp" -#include "debug.hpp" #include "errorhandling.hpp" #include "event.hpp" #include "grid.hpp" @@ -135,12 +134,6 @@ void cells_resort_particles(int global_flag) { n_verlet_updates++; cell_structure.resort_particles(global_flag); - -#ifdef ADDITIONAL_CHECKS - /* at the end of the day, everything should be consistent again */ - check_particle_consistency(cell_structure); - check_particle_sorting(cell_structure); -#endif } /*************************************************/ diff --git a/src/core/debug.cpp b/src/core/debug.cpp deleted file mode 100644 index 0174883d3e4..00000000000 --- a/src/core/debug.cpp +++ /dev/null @@ -1,74 +0,0 @@ -/* - * Copyright (C) 2010-2019 The ESPResSo project - * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 - * Max-Planck-Institute for Polymer Research, Theory Group - * - * 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 . - */ -/** \file - * Implements the malloc replacements as described in debug.hpp. - */ - -#include "debug.hpp" - -#include - -void check_particle_consistency(CellStructure &cs) { - auto local_particles = cs.local_particles(); - auto const cell_part_cnt = local_particles.size(); - - auto const max_id = cs.get_max_local_particle_id(); - - for (auto const &p : local_particles) { - auto const id = p.identity(); - - if (id < 0 || id > max_id) { - throw std::runtime_error("Particle id out of bounds."); - } - - if (cs.get_local_particle(id) != &p) { - throw std::runtime_error("Invalid local particle index entry."); - } - } - - /* checks: local particle id */ - int local_part_cnt = 0; - for (int n = 0; n < cs.get_max_local_particle_id() + 1; n++) { - if (cs.get_local_particle(n) != nullptr) { - local_part_cnt++; - if (cs.get_local_particle(n)->p.identity != n) { - throw std::runtime_error("local_particles part has corrupted id."); - } - } - } - - if (local_part_cnt != cell_part_cnt) { - throw std::runtime_error( - std::to_string(cell_part_cnt) + " parts in cells but " + - std::to_string(local_part_cnt) + " parts in local_particles"); - } -} - -void check_particle_sorting(CellStructure &cs) { - for (auto cell : cs.local_cells()) { - for (auto const &p : cell->particles()) { - if (cs.particle_to_cell(p) != cell) { - throw std::runtime_error("misplaced particle with id " + - std::to_string(p.identity())); - } - } - } -} diff --git a/src/core/debug.hpp b/src/core/debug.hpp deleted file mode 100644 index bba556d69bc..00000000000 --- a/src/core/debug.hpp +++ /dev/null @@ -1,43 +0,0 @@ -/* - * Copyright (C) 2010-2019 The ESPResSo project - * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 - * Max-Planck-Institute for Polymer Research, Theory Group - * - * 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 . - */ -/** \file - * This file controls debug facilities. - * - * Implementation in debug.cpp. - */ -#ifndef DEBUG_HPP -#define DEBUG_HPP - -#include "CellStructure.hpp" -/** this performs a lot of tests which will very likely detect corruptions of - * the cell structure. - */ -void check_particle_consistency(CellStructure &cs); - -/** - * @brief Check if particles are in correct cells. - * - * Check that particles are in the cells the cellsystem says - * they should be. - */ -void check_particle_sorting(CellStructure &cs); - -#endif From 01aa6772b76b987197b736a4d973899b7deb9c53 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Mon, 13 Jul 2020 11:13:25 +0200 Subject: [PATCH 18/27] core: cells: Docstrings and better name --- src/core/CellStructure.cpp | 4 ++-- src/core/CellStructure.hpp | 18 +++++++++++++++++- 2 files changed, 19 insertions(+), 3 deletions(-) diff --git a/src/core/CellStructure.cpp b/src/core/CellStructure.cpp index 515bd6d9aa3..a96c1c43f66 100644 --- a/src/core/CellStructure.cpp +++ b/src/core/CellStructure.cpp @@ -28,7 +28,7 @@ #include -void CellStructure::check_particle_consistency() { +void CellStructure::check_particle_index() { auto const max_id = get_max_local_particle_id(); for (auto const &p : local_particles()) { @@ -216,7 +216,7 @@ void CellStructure::resort_particles(int global_flag) { m_rebuild_verlet_list = true; #ifdef ADDITIONAL_CHECKS - check_particle_consistency(); + check_particle_index(); check_particle_sorting(); #endif } diff --git a/src/core/CellStructure.hpp b/src/core/CellStructure.hpp index 9d66010398d..1b1decdd006 100644 --- a/src/core/CellStructure.hpp +++ b/src/core/CellStructure.hpp @@ -561,7 +561,23 @@ struct CellStructure { } private: - void check_particle_consistency(); + /** + * @brief Check that particle index is commensurate with particles. + * + * For each local particles is checked that has a correct entry + * in the particles index, and that there are no excess (non-existing) + * particles in the index. + */ + void check_particle_index(); + + /** + * @brief Check that particles are in the correct cell. + * + * This checks for all local particles that the result + * of particles_to_cell is the cell the particles is + * actually in, e.g. that the particles are sorted according + * to particles_to_cell. + */ void check_particle_sorting(); }; From 6a3e36457f3be52fc82fc09964f761c570b0ce25 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Mon, 13 Jul 2020 11:30:01 +0200 Subject: [PATCH 19/27] core: Remove n_verlet_updates global --- src/core/cells.cpp | 5 ----- src/core/integrate.cpp | 8 +++++--- src/core/integrate.hpp | 3 --- 3 files changed, 5 insertions(+), 11 deletions(-) diff --git a/src/core/cells.cpp b/src/core/cells.cpp index 4f4fbb96e14..a0e695619ef 100644 --- a/src/core/cells.cpp +++ b/src/core/cells.cpp @@ -32,7 +32,6 @@ #include "grid.hpp" #include "integrate.hpp" #include "nonbonded_interactions/nonbonded_interaction_data.hpp" -#include "short_range_loop.hpp" #include "AtomDecomposition.hpp" #include "DomainDecomposition.hpp" @@ -42,8 +41,6 @@ #include #include -#include - /** Type of cell structure in use */ CellStructure cell_structure; @@ -131,8 +128,6 @@ void cells_re_init(int new_cs) { /*************************************************/ void cells_resort_particles(int global_flag) { - n_verlet_updates++; - cell_structure.resort_particles(global_flag); } diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index 3f0e921d468..4b149d61242 100644 --- a/src/core/integrate.cpp +++ b/src/core/integrate.cpp @@ -72,8 +72,6 @@ int integ_switch = INTEG_METHOD_NVT; -int n_verlet_updates = 0; - double time_step = -1.0; double sim_time = 0.0; @@ -206,7 +204,8 @@ int integrate(int n_steps, int reuse_forces) { if (check_runtime_errors(comm_cart)) return 0; - n_verlet_updates = 0; + /* incremented if a Verlet update is done, aka particle resorting. */ + int n_verlet_updates = 0; #ifdef VALGRIND_INSTRUMENTATION CALLGRIND_START_INSTRUMENTATION; @@ -243,6 +242,9 @@ int integrate(int n_steps, int reuse_forces) { virtual_sites()->update(); #endif + if (cell_structure.get_resort_particles() >= Cells::RESORT_LOCAL) + n_verlet_updates++; + // Communication step: distribute ghost positions cells_update_ghosts(global_ghost_flags()); diff --git a/src/core/integrate.hpp b/src/core/integrate.hpp index bc8ea9625b7..ae2998343e8 100644 --- a/src/core/integrate.hpp +++ b/src/core/integrate.hpp @@ -36,9 +36,6 @@ /** Switch determining which integrator to use. */ extern int integ_switch; -/** incremented if a Verlet update is done, aka particle resorting. */ -extern int n_verlet_updates; - /** Time step for the integration. */ extern double time_step; From bfe374168c5d14a7effcaf634840a6cdb910bbb5 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Mon, 13 Jul 2020 11:43:30 +0200 Subject: [PATCH 20/27] core: Inline cells_resort_particles --- src/core/cells.cpp | 13 +++---------- src/core/cells.hpp | 5 ----- src/core/communication.cpp | 4 ++-- 3 files changed, 5 insertions(+), 17 deletions(-) diff --git a/src/core/cells.cpp b/src/core/cells.cpp index a0e695619ef..e7d8086fa00 100644 --- a/src/core/cells.cpp +++ b/src/core/cells.cpp @@ -127,12 +127,6 @@ void cells_re_init(int new_cs) { /*************************************************/ -void cells_resort_particles(int global_flag) { - cell_structure.resort_particles(global_flag); -} - -/*************************************************/ - void check_resort_particles() { const double skin2 = Utils::sqr(skin / 2.0); @@ -153,9 +147,9 @@ void cells_update_ghosts(unsigned data_parts) { auto constexpr resort_only_parts = Cells::DATA_PART_PROPERTIES | Cells::DATA_PART_BONDS; - unsigned int localResort = cell_structure.get_resort_particles(); auto const global_resort = - boost::mpi::all_reduce(comm_cart, localResort, std::bit_or()); + boost::mpi::all_reduce(comm_cart, cell_structure.get_resort_particles(), + std::bit_or()); if (global_resort != Cells::RESORT_NONE) { int global = (global_resort & Cells::RESORT_GLOBAL) @@ -163,8 +157,7 @@ void cells_update_ghosts(unsigned data_parts) { : CELL_NEIGHBOR_EXCHANGE; /* Resort cell system */ - cells_resort_particles(global); - + cell_structure.resort_particles(global); cell_structure.ghosts_update(data_parts); /* Add the ghost particles to the index if we don't already diff --git a/src/core/cells.hpp b/src/core/cells.hpp index f123d5e9cdf..21002fc6016 100644 --- a/src/core/cells.hpp +++ b/src/core/cells.hpp @@ -86,11 +86,6 @@ void cells_re_init(int new_cs); */ void cells_set_use_verlet_lists(bool use_verlet_lists); -/** Sort the particles into the cells and initialize the ghost particle - * structures. - */ -void cells_resort_particles(int global_flag); - /** Update ghost information. If needed, * the particles are also resorted. */ diff --git a/src/core/communication.cpp b/src/core/communication.cpp index 971a2e9c097..8d5426cad1a 100644 --- a/src/core/communication.cpp +++ b/src/core/communication.cpp @@ -649,7 +649,7 @@ void mpi_loop() { std::vector mpi_resort_particles(int global_flag) { mpi_call(mpi_resort_particles_slave, global_flag, 0); - cells_resort_particles(global_flag); + cell_structure.resort_particles(global_flag); clear_particle_node(); @@ -662,7 +662,7 @@ std::vector mpi_resort_particles(int global_flag) { } void mpi_resort_particles_slave(int global_flag, int) { - cells_resort_particles(global_flag); + cell_structure.resort_particles(global_flag); boost::mpi::gather( comm_cart, static_cast(cell_structure.local_particles().size()), 0); From d55462291b72abf7f5a2b93a37ed35d5defdf61e Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Mon, 13 Jul 2020 11:57:22 +0200 Subject: [PATCH 21/27] core: cells: Make CellStructure::local_cells() private --- src/core/CellStructure.hpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/core/CellStructure.hpp b/src/core/CellStructure.hpp index 1b1decdd006..d1ad758af11 100644 --- a/src/core/CellStructure.hpp +++ b/src/core/CellStructure.hpp @@ -243,8 +243,10 @@ struct CellStructure { /** Maximal pair range supported by current cell system. */ Utils::Vector3d max_range() const; - /** Return the global local_cells */ +private: Utils::Span local_cells(); + +public: ParticleRange local_particles(); ParticleRange ghost_particles(); From db8f333d01c3a231ec531f40e0fe61c81ff00ea8 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Mon, 13 Jul 2020 12:24:24 +0200 Subject: [PATCH 22/27] core: coldet: const-correctness --- src/core/collision.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/core/collision.cpp b/src/core/collision.cpp index 526bad17641..30162826885 100644 --- a/src/core/collision.cpp +++ b/src/core/collision.cpp @@ -284,8 +284,9 @@ void bind_at_point_of_collision_calc_vs_pos(const Particle *const p1, } // Considers three particles for three_particle_binding and performs -// the binding if the criteria are met // -void coldet_do_three_particle_bond(Particle &p, Particle &p1, Particle &p2) { +// the binding if the criteria are met +void coldet_do_three_particle_bond(Particle &p, Particle const &p1, + Particle const &p2) { // If p1 and p2 are not closer or equal to the cutoff distance, skip // p1: if (get_mi_vector(p.r.p, p1.r.p, box_geo).norm() > collision_params.distance) From 5d92d135ea401b05bd9d4e96cd8cef09070d0762 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Mon, 13 Jul 2020 16:02:59 +0200 Subject: [PATCH 23/27] core: cells: Move find_current_cell() to CellStructure --- src/core/CellStructure.hpp | 22 ++++++++++++++++++++++ src/core/cells.cpp | 8 +------- 2 files changed, 23 insertions(+), 7 deletions(-) diff --git a/src/core/CellStructure.hpp b/src/core/CellStructure.hpp index d1ad758af11..87cdc879da7 100644 --- a/src/core/CellStructure.hpp +++ b/src/core/CellStructure.hpp @@ -250,6 +250,7 @@ struct CellStructure { ParticleRange local_particles(); ParticleRange ghost_particles(); +private: /** Cell system dependent function to find the right cell for a * particle. * \param p Particle. @@ -258,6 +259,7 @@ struct CellStructure { */ Cell *particle_to_cell(const Particle &p); +public: /** * @brief Add a particle. * @@ -581,6 +583,26 @@ struct CellStructure { * to particles_to_cell. */ void check_particle_sorting(); + +public: + /** + * @brief Find cell a particle is stored in. + * + * For local particles, this returns the cell they + * are stored in, otherwise nullptr is returned. + * + * @param p Particle to find cell for + * @return Cell for particle or nullptr. + */ + Cell *find_current_cell(const Particle &p) { + assert(not get_resort_particles()); + + if (p.l.ghost) { + return nullptr; + } + + return particle_to_cell(p); + } }; #endif // ESPRESSO_CELLSTRUCTURE_HPP diff --git a/src/core/cells.cpp b/src/core/cells.cpp index e7d8086fa00..51d02298dff 100644 --- a/src/core/cells.cpp +++ b/src/core/cells.cpp @@ -177,13 +177,7 @@ void cells_update_ghosts(unsigned data_parts) { } Cell *find_current_cell(const Particle &p) { - assert(not cell_structure.get_resort_particles()); - - if (p.l.ghost) { - return nullptr; - } - - return cell_structure.particle_to_cell(p); + return cell_structure.find_current_cell(p); } const DomainDecomposition *get_domain_decomposition() { From 6654fcdcad8b29e5f95d211178a082751b305142 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 25 Aug 2020 20:53:21 +0200 Subject: [PATCH 24/27] utils: Remove NoOp class --- src/core/dpd.cpp | 1 - .../electrostatics_magnetostatics/icc.cpp | 1 - src/utils/include/utils/NoOp.hpp | 33 ------------------- 3 files changed, 35 deletions(-) delete mode 100644 src/utils/include/utils/NoOp.hpp diff --git a/src/core/dpd.cpp b/src/core/dpd.cpp index a5c7b0bb93d..05b844ee1a0 100644 --- a/src/core/dpd.cpp +++ b/src/core/dpd.cpp @@ -34,7 +34,6 @@ #include "thermostat.hpp" #include -#include #include #include diff --git a/src/core/electrostatics_magnetostatics/icc.cpp b/src/core/electrostatics_magnetostatics/icc.cpp index 73128a7bf07..12351607ebd 100644 --- a/src/core/electrostatics_magnetostatics/icc.cpp +++ b/src/core/electrostatics_magnetostatics/icc.cpp @@ -46,7 +46,6 @@ #include "nonbonded_interactions/nonbonded_interaction_data.hpp" #include "short_range_loop.hpp" -#include #include "electrostatics_magnetostatics/coulomb.hpp" #include "electrostatics_magnetostatics/coulomb_inline.hpp" diff --git a/src/utils/include/utils/NoOp.hpp b/src/utils/include/utils/NoOp.hpp deleted file mode 100644 index 69c6f559420..00000000000 --- a/src/utils/include/utils/NoOp.hpp +++ /dev/null @@ -1,33 +0,0 @@ -/* - * Copyright (C) 2010-2019 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 CORE_UTILS_NO_OP_HPP -#define CORE_UTILS_NO_OP_HPP - -namespace Utils { - -/** - * @brief A NoOp functor that does nothing. - */ -class NoOp { -public: - template void operator()(Args...) const {} -}; - -} // namespace Utils -#endif From 240649a37af8dd3e8ceda1e721b35726a75bbc39 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 25 Aug 2020 21:54:04 +0200 Subject: [PATCH 25/27] core: Header cleanup --- src/core/cells.cpp | 5 ++--- src/core/cells.hpp | 16 ++-------------- src/core/collision.cpp | 2 +- src/core/communication.cpp | 18 ++---------------- src/core/dpd.cpp | 2 +- src/core/electrostatics_magnetostatics/icc.cpp | 12 +++--------- src/core/event.cpp | 1 - .../immersed_boundary/ImmersedBoundaries.cpp | 1 - .../nonbonded_interaction_data.cpp | 5 ----- src/core/particle_data.hpp | 3 --- 10 files changed, 11 insertions(+), 54 deletions(-) diff --git a/src/core/cells.cpp b/src/core/cells.cpp index 51d02298dff..ceaa1dee2f2 100644 --- a/src/core/cells.cpp +++ b/src/core/cells.cpp @@ -27,15 +27,14 @@ #include "cells.hpp" #include "Particle.hpp" #include "communication.hpp" -#include "errorhandling.hpp" #include "event.hpp" #include "grid.hpp" #include "integrate.hpp" -#include "nonbonded_interactions/nonbonded_interaction_data.hpp" -#include "AtomDecomposition.hpp" #include "DomainDecomposition.hpp" +#include "ParticleDecomposition.hpp" +#include #include #include diff --git a/src/core/cells.hpp b/src/core/cells.hpp index 21002fc6016..b190baa304e 100644 --- a/src/core/cells.hpp +++ b/src/core/cells.hpp @@ -38,8 +38,10 @@ * special method like P3M. */ +#include "Cell.hpp" #include "CellStructure.hpp" #include "DomainDecomposition.hpp" +#include "Particle.hpp" #include #include @@ -59,21 +61,9 @@ enum { /*@}*/ -/************************************************************/ -/** \name Exported Variables */ -/************************************************************/ -/*@{*/ - /** Type of cell structure in use. */ extern CellStructure cell_structure; -/*@}*/ - -/************************************************************/ -/** \name Exported Functions */ -/************************************************************/ -/*@{*/ - /** Reinitialize the cell structures. * @param new_cs The new topology to use afterwards. */ @@ -101,8 +91,6 @@ std::vector> mpi_get_pairs(double distance); /** Check if a particle resorting is required. */ void check_resort_particles(); -/*@}*/ - /** * @brief Find the cell in which a particle is stored. * diff --git a/src/core/collision.cpp b/src/core/collision.cpp index 30162826885..43625629abc 100644 --- a/src/core/collision.cpp +++ b/src/core/collision.cpp @@ -28,7 +28,7 @@ #include "event.hpp" #include "grid.hpp" #include "nonbonded_interactions/nonbonded_interaction_data.hpp" -#include "virtual_sites/VirtualSitesRelative.hpp" +#include "virtual_sites.hpp" #include #include diff --git a/src/core/communication.cpp b/src/core/communication.cpp index 8d5426cad1a..c2a7a1eed6e 100644 --- a/src/core/communication.cpp +++ b/src/core/communication.cpp @@ -20,7 +20,6 @@ */ #include #include -#include #include #include #ifdef OPEN_MPI @@ -32,33 +31,23 @@ #include "errorhandling.hpp" +#include "CellStructure.hpp" #include "EspressoSystemInterface.hpp" -#include "Particle.hpp" -#include "bonded_interactions/bonded_tab.hpp" +#include "bonded_interactions/bonded_interaction_data.hpp" #include "cells.hpp" -#include "collision.hpp" #include "cuda_interface.hpp" #include "energy.hpp" #include "event.hpp" -#include "forces.hpp" #include "galilei.hpp" #include "global.hpp" #include "grid.hpp" #include "grid_based_algorithms/lb.hpp" #include "grid_based_algorithms/lb_interface.hpp" -#include "grid_based_algorithms/lb_interpolation.hpp" -#include "grid_based_algorithms/lb_particle_coupling.hpp" #include "integrate.hpp" #include "integrators/steepest_descent.hpp" -#include "io/mpiio/mpiio.hpp" -#include "nonbonded_interactions/nonbonded_tab.hpp" #include "npt.hpp" -#include "partCfg_global.hpp" #include "particle_data.hpp" #include "pressure.hpp" -#include "rotation.hpp" -#include "stokesian_dynamics/sd_interface.hpp" -#include "virtual_sites.hpp" #include "electrostatics_magnetostatics/coulomb.hpp" #include "electrostatics_magnetostatics/dipole.hpp" @@ -67,9 +56,6 @@ #include "serialization/IA_parameters.hpp" -#include -#include - #include #include #include diff --git a/src/core/dpd.cpp b/src/core/dpd.cpp index 05b844ee1a0..33e9d9d1569 100644 --- a/src/core/dpd.cpp +++ b/src/core/dpd.cpp @@ -24,13 +24,13 @@ #include "dpd.hpp" #ifdef DPD +#include "cells.hpp" #include "communication.hpp" #include "event.hpp" #include "grid.hpp" #include "integrate.hpp" #include "nonbonded_interactions/nonbonded_interaction_data.hpp" #include "random.hpp" -#include "short_range_loop.hpp" #include "thermostat.hpp" #include diff --git a/src/core/electrostatics_magnetostatics/icc.cpp b/src/core/electrostatics_magnetostatics/icc.cpp index 12351607ebd..a9b6a7e33f5 100644 --- a/src/core/electrostatics_magnetostatics/icc.cpp +++ b/src/core/electrostatics_magnetostatics/icc.cpp @@ -31,25 +31,19 @@ #ifdef ELECTROSTATICS #include -#include #include -#include "electrostatics_magnetostatics/p3m_gpu.hpp" - #include "Particle.hpp" #include "cells.hpp" #include "communication.hpp" -#include "config.hpp" #include "errorhandling.hpp" #include "event.hpp" -#include "forces.hpp" -#include "nonbonded_interactions/nonbonded_interaction_data.hpp" - -#include "short_range_loop.hpp" #include "electrostatics_magnetostatics/coulomb.hpp" #include "electrostatics_magnetostatics/coulomb_inline.hpp" +#include + iccp3m_struct iccp3m_cfg; void init_forces_iccp3m(const ParticleRange &particles, @@ -101,7 +95,7 @@ int iccp3m_iteration(const ParticleRange &particles, << "ICCP3M: nonpositive dielectric constant is not allowed."; } - auto const pref = 1.0 / (coulomb.prefactor * 6.283185307); + auto const pref = 1.0 / (coulomb.prefactor * 2 * Utils::pi()); iccp3m_cfg.citeration = 0; double globalmax = 1e100; diff --git a/src/core/event.cpp b/src/core/event.cpp index 48dd925d567..e2c1d483320 100644 --- a/src/core/event.cpp +++ b/src/core/event.cpp @@ -25,7 +25,6 @@ */ #include "event.hpp" -#include "Particle.hpp" #include "bonded_interactions/thermalized_bond.hpp" #include "cells.hpp" #include "collision.hpp" diff --git a/src/core/immersed_boundary/ImmersedBoundaries.cpp b/src/core/immersed_boundary/ImmersedBoundaries.cpp index 5dc64c7003a..fab9f2063e2 100644 --- a/src/core/immersed_boundary/ImmersedBoundaries.cpp +++ b/src/core/immersed_boundary/ImmersedBoundaries.cpp @@ -22,7 +22,6 @@ #include "Particle.hpp" #include "bonded_interactions/bonded_interaction_data.hpp" #include "communication.hpp" -#include "errorhandling.hpp" #include "grid.hpp" #include diff --git a/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp b/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp index 2a3dad1e86b..8dd3ee3ca9d 100644 --- a/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp +++ b/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp @@ -25,7 +25,6 @@ #include "bonded_interactions/bonded_interaction_data.hpp" #include "collision.hpp" #include "communication.hpp" -#include "errorhandling.hpp" #include "grid.hpp" #include "serialization/IA_parameters.hpp" @@ -49,10 +48,6 @@ std::vector ia_params; double min_global_cut = INACTIVE_CUTOFF; -/***************************************** - * function prototypes - *****************************************/ - /***************************************** * general low-level functions *****************************************/ diff --git a/src/core/particle_data.hpp b/src/core/particle_data.hpp index 71c2fbc5105..f2301b4b5fe 100644 --- a/src/core/particle_data.hpp +++ b/src/core/particle_data.hpp @@ -74,9 +74,6 @@ enum { * Functions ************************************************/ -/* Other Functions */ -/************************************************/ - /** Invalidate \ref particle_node. This has to be done * at the beginning of the integration. */ From 3e1d418dede4f6be292a30201c6db91a1e5b90c0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 26 Aug 2020 11:19:50 +0200 Subject: [PATCH 26/27] core: cells: Skip non-bonded kernels when inactive Skip non-bonded kernels in the short range loop if no non-bonded interactions are active. --- src/core/CellStructure.hpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/core/CellStructure.hpp b/src/core/CellStructure.hpp index 87cdc879da7..32837dde4c7 100644 --- a/src/core/CellStructure.hpp +++ b/src/core/CellStructure.hpp @@ -33,6 +33,7 @@ #include "algorithm/link_cell.hpp" #include "bond_error.hpp" #include "ghosts.hpp" +#include "nonbonded_interactions/nonbonded_interaction_data.hpp" #include #include @@ -516,7 +517,8 @@ struct CellStructure { * @param pair_kernel Kernel to apply */ template void non_bonded_loop(PairKernel pair_kernel) { - link_cell(pair_kernel); + if (maximal_cutoff() > 0.) + link_cell(pair_kernel); } /** Non-bonded pair loop with potential use @@ -527,6 +529,8 @@ struct CellStructure { template void non_bonded_loop(PairKernel &&pair_kernel, const VerletCriterion &verlet_criterion) { + if (maximal_cutoff() <= 0.) + return; /* In this case the verlet list update is attached to * the pair kernel, and the verlet list is rebuilt as * we go. */ From 66c7a94b2be185a64694c38bcc639db25560916d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 26 Aug 2020 19:18:53 +0200 Subject: [PATCH 27/27] core: cells: Narrow inclusion of header file --- src/core/CellStructure.hpp | 6 +----- src/core/energy.cpp | 4 +++- src/core/forces.cpp | 2 ++ src/core/pressure.cpp | 4 +++- src/core/short_range_loop.hpp | 4 +++- 5 files changed, 12 insertions(+), 8 deletions(-) diff --git a/src/core/CellStructure.hpp b/src/core/CellStructure.hpp index 32837dde4c7..87cdc879da7 100644 --- a/src/core/CellStructure.hpp +++ b/src/core/CellStructure.hpp @@ -33,7 +33,6 @@ #include "algorithm/link_cell.hpp" #include "bond_error.hpp" #include "ghosts.hpp" -#include "nonbonded_interactions/nonbonded_interaction_data.hpp" #include #include @@ -517,8 +516,7 @@ struct CellStructure { * @param pair_kernel Kernel to apply */ template void non_bonded_loop(PairKernel pair_kernel) { - if (maximal_cutoff() > 0.) - link_cell(pair_kernel); + link_cell(pair_kernel); } /** Non-bonded pair loop with potential use @@ -529,8 +527,6 @@ struct CellStructure { template void non_bonded_loop(PairKernel &&pair_kernel, const VerletCriterion &verlet_criterion) { - if (maximal_cutoff() <= 0.) - return; /* In this case the verlet list update is attached to * the pair kernel, and the verlet list is rebuilt as * we go. */ diff --git a/src/core/energy.cpp b/src/core/energy.cpp index 8df39a5b9a3..aa487a7c9ad 100644 --- a/src/core/energy.cpp +++ b/src/core/energy.cpp @@ -31,6 +31,7 @@ #include "energy_inline.hpp" #include "event.hpp" #include "forces.hpp" +#include "nonbonded_interactions/nonbonded_interaction_data.hpp" #include "reduce_observable_stat.hpp" #include "short_range_loop.hpp" @@ -83,7 +84,8 @@ void energy_calc(const double time) { [](Particle const &p1, Particle const &p2, Distance const &d) { add_non_bonded_pair_energy(p1, p2, d.vec21, sqrt(d.dist2), d.dist2, obs_energy); - }); + }, + maximal_cutoff()); calc_long_range_energies(cell_structure.local_particles()); diff --git a/src/core/forces.cpp b/src/core/forces.cpp index b616bdb81b6..4b5f4b34eeb 100644 --- a/src/core/forces.cpp +++ b/src/core/forces.cpp @@ -40,6 +40,7 @@ #include "grid_based_algorithms/lb_interface.hpp" #include "grid_based_algorithms/lb_particle_coupling.hpp" #include "immersed_boundaries.hpp" +#include "nonbonded_interactions/nonbonded_interaction_data.hpp" #include "short_range_loop.hpp" #include @@ -125,6 +126,7 @@ void force_calc(CellStructure &cell_structure) { detect_collision(p1, p2, d.dist2); #endif }, + maximal_cutoff(), VerletCriterion{skin, interaction_range(), coulomb_cutoff, dipole_cutoff, collision_detection_cutoff()}); diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp index d3b1b3027ef..47710bc3d31 100644 --- a/src/core/pressure.cpp +++ b/src/core/pressure.cpp @@ -26,6 +26,7 @@ #include "cells.hpp" #include "communication.hpp" #include "event.hpp" +#include "nonbonded_interactions/nonbonded_interaction_data.hpp" #include "pressure_inline.hpp" #include "reduce_observable_stat.hpp" #include "virtual_sites.hpp" @@ -102,7 +103,8 @@ void pressure_calc() { [](Particle &p1, Particle &p2, Distance const &d) { add_non_bonded_pair_virials(p1, p2, d.vec21, sqrt(d.dist2), obs_pressure); - }); + }, + maximal_cutoff()); calc_long_range_virials(cell_structure.local_particles()); diff --git a/src/core/short_range_loop.hpp b/src/core/short_range_loop.hpp index 937bcc57331..2b903b27cc4 100644 --- a/src/core/short_range_loop.hpp +++ b/src/core/short_range_loop.hpp @@ -38,12 +38,14 @@ struct True { template void short_range_loop(BondKernel bond_kernel, PairKernel pair_kernel, + double distance_cutoff = 1., const VerletCriterion &verlet_criterion = {}) { ESPRESSO_PROFILER_CXX_MARK_FUNCTION; assert(cell_structure.get_resort_particles() == Cells::RESORT_NONE); cell_structure.bond_loop(bond_kernel); - cell_structure.non_bonded_loop(pair_kernel, verlet_criterion); + if (distance_cutoff > 0.) + cell_structure.non_bonded_loop(pair_kernel, verlet_criterion); } #endif