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/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/CMakeLists.txt b/src/core/CMakeLists.txt index 14775881299..a2a807c89e9 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 630e11f2946..a96c1c43f66 100644 --- a/src/core/CellStructure.cpp +++ b/src/core/CellStructure.cpp @@ -26,6 +26,52 @@ #include +#include + +void CellStructure::check_particle_index() { + 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); } @@ -166,6 +212,13 @@ void CellStructure::resort_particles(int global_flag) { for (auto d : diff) { boost::apply_visitor(UpdateParticleIndexVisitor{this}, d); } + + m_rebuild_verlet_list = true; + +#ifdef ADDITIONAL_CHECKS + check_particle_index(); + check_particle_sorting(); +#endif } void CellStructure::set_atom_decomposition(boost::mpi::communicator const &comm, @@ -180,4 +233,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 5d12e08e9a3..87cdc879da7 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 @@ -80,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 @@ -103,6 +131,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. * @@ -212,11 +243,14 @@ 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(); +private: /** Cell system dependent function to find the right cell for a * particle. * \param p Particle. @@ -225,6 +259,7 @@ struct CellStructure { */ Cell *particle_to_cell(const Particle &p); +public: /** * @brief Add a particle. * @@ -293,7 +328,6 @@ struct CellStructure { return assert(m_decomposition), *m_decomposition; } -private: ParticleDecomposition &decomposition() { return assert(m_decomposition), *m_decomposition; } @@ -327,6 +361,7 @@ struct CellStructure { * Cells::DataPart */ void ghosts_update(unsigned data_parts); + /** * @brief Add forces from ghost particles to real particles. */ @@ -355,7 +390,6 @@ struct CellStructure { return partners; } -public: /** * @brief Execute kernel for every bond on particle. * @tparam Handler Callable, which can be invoked with @@ -387,9 +421,10 @@ struct CellStructure { } } -private: - /** 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) { @@ -443,11 +478,130 @@ struct CellStructure { double range, BoxGeometry const &box, LocalBox const &local_geo); +public: + template void bond_loop(BondKernel const &bond_kernel) { + for (auto &p : local_particles()) { + execute_bond_handler(p, bond_kernel); + } + } + +private: /** - * @brief Return true if minimum image convention is - * needed for distance calculation. */ - bool minimum_image_distance() const { - return m_decomposition->minimum_image_distance(); + * @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: + /** 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 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. */ + if (use_verlet_list && m_rebuild_verlet_list) { + m_verlet_list.clear(); + + 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) { + 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 { + /* No verlet lists, just run the kernel with pairs from the cells. */ + link_cell(pair_kernel); + } + } + +private: + /** + * @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(); + +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); } }; 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..6577da19080 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; } @@ -100,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/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/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; }; diff --git a/src/core/algorithm/for_each_pair.hpp b/src/core/algorithm/for_each_pair.hpp deleted file mode 100644 index 20d35f64b12..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 single and pair kernel for each particle (pair) from 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 - * evaluated and @p verlet_criterion is evaluated with the calculated distance. - * Iff true, the pair_kernel is called. - * - * For details see verlet_ia and 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. - * - * 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. - */ -template -void for_each_pair(CellIterator first, CellIterator last, - ParticleKernel &&particle_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), - std::forward(distance_function), - std::forward(verlet_criterion), rebuild); - } else { - link_cell(first, last, std::forward(particle_kernel), - 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 ef4c6bc6e13..28a30f19370 100644 --- a/src/core/algorithm/link_cell.hpp +++ b/src/core/algorithm/link_cell.hpp @@ -26,29 +26,23 @@ namespace Algorithm { * and over all pairs within the cells and with * their neighbors. */ -template +template void link_cell(CellIterator first, CellIterator last, - ParticleKernel &&particle_kernel, PairKernel &&pair_kernel, - DistanceFunction &&distance_function) { + PairKernel &&pair_kernel) { 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); - 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 55c08639cee..00000000000 --- a/src/core/algorithm/verlet_ia.hpp +++ /dev/null @@ -1,111 +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, - ParticleKernel &&particle_kernel, - 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; - - 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); - 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 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, - ParticleKernel &&particle_kernel, 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); - } - } -} -} // 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, - ParticleKernel &&particle_kernel, 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), - std::forward(distance_function)); - } -} -} // namespace Algorithm - -#endif diff --git a/src/core/cells.cpp b/src/core/cells.cpp index a5efb55e2df..ceaa1dee2f2 100644 --- a/src/core/cells.cpp +++ b/src/core/cells.cpp @@ -27,30 +27,22 @@ #include "cells.hpp" #include "Particle.hpp" #include "communication.hpp" -#include "debug.hpp" -#include "errorhandling.hpp" #include "event.hpp" #include "grid.hpp" #include "integrate.hpp" -#include "nonbonded_interactions/nonbonded_interaction_data.hpp" -#include "short_range_loop.hpp" -#include "AtomDecomposition.hpp" #include "DomainDecomposition.hpp" +#include "ParticleDecomposition.hpp" -#include +#include #include #include #include -#include - /** Type of cell structure in use */ CellStructure cell_structure; -bool rebuild_verletlist = true; - /** * @brief Get pairs closer than distance from the cells. * @@ -72,7 +64,7 @@ std::vector> get_pairs(double distance) { ret.emplace_back(p1.p.identity, p2.p.identity); }; - short_range_loop(Utils::NoOp{}, pair_kernel); + cell_structure.non_bonded_loop(pair_kernel); /* Sort pairs */ for (auto &pair : ret) { @@ -134,22 +126,6 @@ void cells_re_init(int new_cs) { /*************************************************/ -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(); - check_particle_sorting(); -#endif - - rebuild_verletlist = true; -} - -/*************************************************/ - void check_resort_particles() { const double skin2 = Utils::sqr(skin / 2.0); @@ -170,9 +146,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) @@ -180,8 +156,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 @@ -201,13 +176,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() { @@ -217,4 +186,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/cells.hpp b/src/core/cells.hpp index d5897a580a2..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,26 +61,9 @@ enum { /*@}*/ -/************************************************************/ -/** \name Exported Variables */ -/************************************************************/ -/*@{*/ - /** 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; - -/*@}*/ - -/************************************************************/ -/** \name Exported Functions */ -/************************************************************/ -/*@{*/ - /** Reinitialize the cell structures. * @param new_cs The new topology to use afterwards. */ @@ -91,11 +76,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. */ @@ -111,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 add90eac340..7e40e947a04 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 @@ -104,12 +104,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 @@ -243,8 +237,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; } @@ -294,8 +287,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) 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/communication.cpp b/src/core/communication.cpp index 971a2e9c097..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 @@ -649,7 +635,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 +648,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); 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/debug.cpp b/src/core/debug.cpp deleted file mode 100644 index fbbb54b6f71..00000000000 --- a/src/core/debug.cpp +++ /dev/null @@ -1,76 +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 "cells.hpp" - -#include - -void check_particle_consistency() { - auto local_particles = cell_structure.local_particles(); - auto const cell_part_cnt = local_particles.size(); - - auto const max_id = cell_structure.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 (cell_structure.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) { - local_part_cnt++; - if (cell_structure.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() { - for (auto cell : cell_structure.local_cells()) { - for (auto const &p : cell->particles()) { - if (cell_structure.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 2d48d12852c..00000000000 --- a/src/core/debug.hpp +++ /dev/null @@ -1,42 +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 - -/** this performs a lot of tests which will very likely detect corruptions of - * the cell structure. - */ -void check_particle_consistency(); - -/** - * @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(); - -#endif diff --git a/src/core/dpd.cpp b/src/core/dpd.cpp index 87846abacfd..33e9d9d1569 100644 --- a/src/core/dpd.cpp +++ b/src/core/dpd.cpp @@ -24,16 +24,16 @@ #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 -#include #include #include @@ -142,8 +142,7 @@ static auto dpd_viscous_stress_local() { on_observable_calc(); Utils::Vector stress{}; - short_range_loop( - Utils::NoOp{}, + cell_structure.non_bonded_loop( [&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..a9b6a7e33f5 100644 --- a/src/core/electrostatics_magnetostatics/icc.cpp +++ b/src/core/electrostatics_magnetostatics/icc.cpp @@ -31,26 +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 #include "electrostatics_magnetostatics/coulomb.hpp" #include "electrostatics_magnetostatics/coulomb_inline.hpp" +#include + iccp3m_struct iccp3m_cfg; void init_forces_iccp3m(const ParticleRange &particles, @@ -102,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; @@ -197,9 +190,9 @@ 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 */ + 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); }); diff --git a/src/core/energy.cpp b/src/core/energy.cpp index af3e3485e15..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" @@ -71,11 +72,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); - }); + }, + maximal_cutoff()); 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/event.cpp b/src/core/event.cpp index 15395fd6418..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" @@ -113,11 +112,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(); @@ -256,6 +250,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. */ @@ -373,4 +369,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/forces.cpp b/src/core/forces.cpp index eaba1e80128..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 @@ -117,7 +118,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 @@ -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/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/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/immersed_boundary/ImmersedBoundaries.cpp b/src/core/immersed_boundary/ImmersedBoundaries.cpp index 6059b54cab0..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 @@ -45,16 +44,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 +72,9 @@ void ImmersedBoundaries::init_volume_conservation(CellStructure &cs) { } } } - } - VolumeInitDone = true; + VolumeInitDone = true; + } } /** Set parameters of volume conservation */ @@ -131,16 +133,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]; @@ -170,13 +168,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; @@ -188,58 +184,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/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; diff --git a/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp b/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp index b774c3aca31..8dd3ee3ca9d 100644 --- a/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp +++ b/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp @@ -23,8 +23,8 @@ */ #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" #include "serialization/IA_parameters.hpp" @@ -48,10 +48,6 @@ std::vector ia_params; double min_global_cut = INACTIVE_CUTOFF; -/***************************************** - * function prototypes - *****************************************/ - /***************************************** * general low-level functions *****************************************/ @@ -198,6 +194,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/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/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/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. */ diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp index db90822ff55..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" @@ -83,11 +84,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); + }, + maximal_cutoff()); 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/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*/ diff --git a/src/core/short_range_loop.hpp b/src/core/short_range_loop.hpp index f4d5c569b74..2b903b27cc4 100644 --- a/src/core/short_range_loop.hpp +++ b/src/core/short_range_loop.hpp @@ -19,99 +19,33 @@ #ifndef CORE_SHORT_RANGE_HPP #define CORE_SHORT_RANGE_HPP -#include "algorithm/for_each_pair.hpp" #include "cells.hpp" -#include "grid.hpp" -#include "integrate.hpp" -#include #include -#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; -}; +#include 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 Decide which distance function to use depending on the - * cell system, and call the pair code. - */ -template -void decide_distance(CellIterator first, CellIterator last, - ParticleKernel &&particle_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); - } else { - Algorithm::for_each_pair( - first, last, std::forward(particle_kernel), - std::forward(pair_kernel), EuclidianDistance{}, - std::forward(verlet_criterion), - cell_structure.use_verlet_list, rebuild_verletlist); - } -} - /** - * @brief Functor that returns true for any argument. + * @brief Functor that returns true for + * any arguments. */ struct True { template bool operator()(T...) const { return true; } }; } // namespace detail -template -void short_range_loop(ParticleKernel &&particle_kernel, - PairKernel &&pair_kernel, +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); - 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); - - rebuild_verletlist = false; - } else { - for (auto &p : cell_structure.local_particles()) { - particle_kernel(p); - } - } + cell_structure.bond_loop(bond_kernel); + if (distance_cutoff > 0.) + cell_structure.non_bonded_loop(pair_kernel, verlet_criterion); } - #endif diff --git a/src/core/unit_tests/CMakeLists.txt b/src/core/unit_tests/CMakeLists.txt index 72cc7e12d1d..e5a8bb4e837 100644 --- a/src/core/unit_tests/CMakeLists.txt +++ b/src/core/unit_tests/CMakeLists.txt @@ -31,7 +31,6 @@ unit_test(NAME ParticleIterator_test SRC ParticleIterator_test.cpp DEPENDS EspressoUtils) unit_test(NAME p3m_test SRC p3m_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 field_coupling_couplings SRC field_coupling_couplings_test.cpp diff --git a/src/core/unit_tests/link_cell_test.cpp b/src/core/unit_tests/link_cell_test.cpp index 55456bcf32b..8f711f14d98 100644 --- a/src/core/unit_tests/link_cell_test.cpp +++ b/src/core/unit_tests/link_cell_test.cpp @@ -54,27 +54,12 @@ 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 - * 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); - }); - - /* 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; })); + Algorithm::link_cell(cells.begin(), cells.end(), + [&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); + }); 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 c1360eec77f..00000000000 --- a/src/core/unit_tests/verlet_ia_test.cpp +++ /dev/null @@ -1,147 +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); - 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); - }, - [](Particle const &p1, Particle const &p2) { - return Distance{p1.p.identity <= p2.p.identity}; - }, - 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); - }, - [](Particle const &p1, Particle const &p2) { - return Distance{p1.p.identity <= p2.p.identity}; - }, - 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); - }, - [](Particle const &p1, Particle const &p2) { - return Distance{p1.p.identity <= p2.p.identity}; - }, - 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); -} 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