ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
core/galilei/ComFixed.hpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2017-2022 The ESPResSo project
3 *
4 * This file is part of ESPResSo.
5 *
6 * ESPResSo is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * ESPResSo is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 */
19
20#pragma once
21
22#include "ParticleRange.hpp"
23#include "communication.hpp"
24
25#include <utils/Vector.hpp>
26#include <utils/keys.hpp>
27
28#include <boost/mpi/collectives/all_reduce.hpp>
29#include <boost/mpi/communicator.hpp>
30
31#include <functional>
32#include <unordered_map>
33#include <vector>
34
35class ComFixed {
36private:
37 std::unordered_map<int, int> m_type_index;
38
39 std::vector<Utils::Vector3d>
40 local_type_forces(ParticleRange const &particles) const {
41 std::vector<Utils::Vector3d> ret(m_type_index.size(), Utils::Vector3d{});
42
43 for (auto const &p : particles) {
44 /* Check if type is of interest */
45 auto it = m_type_index.find(p.type());
46 if (it != m_type_index.end()) {
47 ret[it->second] += p.force();
48 }
49 }
50
51 return ret;
52 }
53
54 std::vector<double> local_type_masses(ParticleRange const &particles) const {
55 std::vector<double> ret(m_type_index.size(), 0.);
56
57 for (auto const &p : particles) {
58 /* Check if type is of interest */
59 auto it = m_type_index.find(p.type());
60 if (it != m_type_index.end()) {
61 ret[it->second] += p.mass();
62 }
63 }
64
65 return ret;
66 }
67
68public:
69 ComFixed() = default;
70
71 void set_fixed_types(std::vector<int> const &p_types) {
72 m_type_index.clear();
73
74 int i = 0;
75 for (auto const &p_type : p_types) {
76 m_type_index[p_type] = i++;
77 }
78 }
79
80 std::vector<int> get_fixed_types() const { return Utils::keys(m_type_index); }
81
82 void apply(ParticleRange const &particles) const {
83 /* Bail out early if there is nothing to do. */
84 if (m_type_index.empty())
85 return;
86
87 auto const local_forces = local_type_forces(particles);
88 auto const local_masses = local_type_masses(particles);
89
90 /* Total forces and masses of the types. */
91 std::vector<Utils::Vector3d> forces(m_type_index.size(), Utils::Vector3d{});
92 std::vector<double> masses(m_type_index.size(), 0.0);
93
94 /* Add contributions from all nodes and redistribute them to all. */
95 boost::mpi::all_reduce(::comm_cart, local_forces.data(),
96 static_cast<int>(local_forces.size()), forces.data(),
97 std::plus<Utils::Vector3d>{});
98 boost::mpi::all_reduce(::comm_cart, local_masses.data(),
99 static_cast<int>(local_masses.size()), masses.data(),
100 std::plus<double>{});
101
102 for (auto &p : particles) {
103 /* Check if type is of interest */
104 auto it = m_type_index.find(p.type());
105 if (it != m_type_index.end()) {
106 auto const mass_frac = p.mass() / masses[it->second];
107 auto const &type_force = forces[it->second];
108 for (unsigned int i = 0u; i < 3u; i++) {
109 p.force()[i] -= mass_frac * type_force[i];
110 }
111 }
112 }
113 }
114};
Vector implementation and trait types for boost qvm interoperability.
ComFixed()=default
std::vector< int > get_fixed_types() const
void apply(ParticleRange const &particles) const
void set_fixed_types(std::vector< int > const &p_types)
A range of particles.
boost::mpi::communicator comm_cart
The communicator.
auto keys(Map const &m) -> std::vector< typename Map::key_type >
Return the keys of a map type.
Definition keys.hpp:34