GCC Code Coverage Report

Directory: ./
File: src/interaction.cpp
Date: 2025-11-05 16:49:12
Exec Total Coverage
Lines: 188 188 100.0%
Functions: 22 22 100.0%
Branches: 198 258 76.7%
Line Branch Exec Source
1 #include <algorithm>
2 #include <iostream>
3 #include <iterator>
4 #ifndef NO_OMP
5 #include <omp.h>
6 #endif
7
8 #include "cells.hpp"
9 #include "integrator.hpp"
10 #include "interaction.hpp"
11 #include "particle.hpp"
12 #include "potentials.hpp"
13 #include "system.hpp"
14
15 107 Interaction::Interaction(System *system, const YAML::Node &config)
16 107 : system(system), particles(system->particles),
17
1/2
✓ Branch 1 taken 107 times.
✗ Branch 2 not taken.
 107 externals(PotentialFactory::createExternals(system, config["external"])),
18
2/4
✓ Branch 1 taken 107 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 107 times.
✗ Branch 5 not taken.
 107 internal(PotentialFactory::createInternal(system, config["internal"],
19
8/14
✓ Branch 1 taken 107 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 50 times.
✓ Branch 4 taken 57 times.
✓ Branch 6 taken 50 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 50 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 50 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 50 times.
✓ Branch 15 taken 57 times.
✗ Branch 20 not taken.
✗ Branch 23 not taken.
 214 config["strategy"] && config["strategy"]["newton3"]
20
13/19
✓ Branch 0 taken 50 times.
✓ Branch 1 taken 57 times.
✓ Branch 3 taken 50 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 50 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 50 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 35 times.
✓ Branch 12 taken 15 times.
✓ Branch 13 taken 50 times.
✓ Branch 14 taken 57 times.
✓ Branch 16 taken 50 times.
✓ Branch 17 taken 57 times.
✓ Branch 19 taken 50 times.
✓ Branch 20 taken 57 times.
✗ Branch 22 not taken.
✗ Branch 25 not taken.
✗ Branch 28 not taken.
 264 ? config["strategy"]["newton3"].as<bool>()
21 : true)),
22 321 allDoE([this]() {
23 107 bool doesE = true;
24
2/2
✓ Branch 0 taken 52 times.
✓ Branch 1 taken 107 times.
 159 for (External *e : externals) {
25 52 doesE &= e->doesE();
26 }
27 107 doesE &= internal->doesE();
28 107 return doesE;
29
1/2
✓ Branch 1 taken 107 times.
✗ Branch 2 not taken.
 107 }()),
30 321 allDoF([this]() {
31 107 bool doesF = true;
32
2/2
✓ Branch 0 taken 52 times.
✓ Branch 1 taken 107 times.
 159 for (External *e : externals) {
33 52 doesF &= e->doesF();
34 }
35 107 doesF &= internal->doesF();
36 107 return doesF;
37
1/2
✓ Branch 2 taken 107 times.
✗ Branch 3 not taken.
 214 }()) {
38
1/2
✓ Branch 1 taken 107 times.
✗ Branch 2 not taken.
 107 std::clog << "\nInteraction method:";
39
2/4
✓ Branch 1 taken 107 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 107 times.
✗ Branch 5 not taken.
 107 std::clog << "\n\tNewton III: " << this->internal->newton3;
40
2/2
✓ Branch 1 taken 19 times.
✓ Branch 2 taken 88 times.
 107 if ((system->L.array() < internal->range).any()) {
41 19 std::clog << "\nWARNING: The range of the internal interactions is larger than some length of "
42
1/2
✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
 19 "the system size.";
43 }
44 107 }
45
46 210 Interaction::~Interaction() {
47
2/2
✓ Branch 0 taken 52 times.
✓ Branch 1 taken 105 times.
 314 for (External *e : externals) {
48
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 104 delete e;
49 }
50
1/2
✓ Branch 0 taken 105 times.
✗ Branch 1 not taken.
 210 delete internal;
51 }
52
53 16 AllInteraction::AllInteraction(System *system, const YAML::Node &config)
54 16 : Interaction(system, config) {
55
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 std::clog << "\n\tstrategy: all";
56 16 }
57
58 2222 double AllInteraction::iteratePairs(Particle &p, PairFunction &func) {
59 2222 int n = particles.size();
60 2222 double sum = 0.;
61
2/2
✓ Branch 0 taken 355886 times.
✓ Branch 1 taken 2222 times.
 358108 for (int j = 0; j < n; ++j) {
62
2/2
✓ Branch 0 taken 353664 times.
✓ Branch 1 taken 2222 times.
 355886 Particle &pj = particles[j];
63
2/2
✓ Branch 0 taken 353664 times.
✓ Branch 1 taken 2222 times.
 355886 if (&p != &pj) {
64 353664 sum += func(p, pj);
65 }
66 }
67 2222 return sum;
68 }
69
70 1488 double AllInteraction::iterateTriples(Particle &p, TripleFunction &func) {
71 1488 int n = particles.size();
72 1488 double sum = 0.;
73
2/2
✓ Branch 0 taken 237294 times.
✓ Branch 1 taken 1488 times.
 238782 for (int j = 0; j < n; ++j) {
74
2/2
✓ Branch 0 taken 1488 times.
✓ Branch 1 taken 235806 times.
 237294 Particle &pj = particles[j];
75
2/2
✓ Branch 0 taken 1488 times.
✓ Branch 1 taken 235806 times.
 237294 if (&p == &pj) {
76 1488 continue;
77 }
78
2/2
✓ Branch 0 taken 26911706 times.
✓ Branch 1 taken 235806 times.
 27147512 for (int k = j + 1; k < n; ++k) {
79
2/2
✓ Branch 0 taken 26675900 times.
✓ Branch 1 taken 235806 times.
 26911706 Particle &pk = particles[k];
80
2/2
✓ Branch 0 taken 26675900 times.
✓ Branch 1 taken 235806 times.
 26911706 if (&p != &pk) {
81 26675900 sum += func(p, pj, pk);
82 }
83 }
84 }
85 1488 return sum;
86 }
87
88 17 double AllInteraction::iteratePairs(PairFunction &func, bool newton3, bool parallel) {
89 17 int n = particles.size();
90 17 double sum = 0.;
91 17 #pragma omp parallel for if (!newton3 && parallel) reduction(+ : sum)
92 for (int i = 0; i < n; ++i) {
93 Particle &pi = particles[i];
94 int jStart = newton3 ? i + 1 : 0;
95 for (int j = jStart; j < n; ++j) {
96 Particle &pj = particles[j];
97 if (&pi != &pj) {
98 sum += func(pi, pj);
99 }
100 }
101 }
102 17 return sum;
103 }
104
105 11 double AllInteraction::iterateTriples(TripleFunction &func, bool newton3, bool parallel) {
106 11 int n = particles.size();
107 11 double sum = 0.;
108 11 #pragma omp parallel for if (!newton3 && parallel) reduction(+ : sum)
109 for (int i = 0; i < n; ++i) {
110 Particle &pi = particles[i];
111 for (int j = 0; j < n; ++j) {
112 Particle &pj = particles[j];
113 if (&pi == &pj) {
114 continue;
115 }
116 for (int k = j + 1; k < n; ++k) {
117 Particle &pk = particles[k];
118 if (&pi != &pk) {
119 sum += func(pi, pj, pk);
120 }
121 }
122 }
123 }
124 11 return sum;
125 }
126
127 83 CellNeighborInteraction::CellNeighborInteraction(System *system, const YAML::Node &config)
128 249 : Interaction(system, config), cellCount([this]() -> ArrayI {
129 83 ArrayI count = floor(this->system->L.array() / internal->range).max(1).cast<int>();
130 #ifndef NO_OMP
131
2/2
✓ Branch 0 taken 39 times.
✓ Branch 1 taken 36 times.
 75 if (internal->newton3 &&
132
17/18
✓ Branch 0 taken 60 times.
✓ Branch 1 taken 15 times.
✓ Branch 2 taken 23 times.
✓ Branch 3 taken 37 times.
✓ Branch 4 taken 23 times.
✓ Branch 5 taken 15 times.
✓ Branch 6 taken 1 times.
✓ Branch 7 taken 22 times.
✓ Branch 8 taken 22 times.
✓ Branch 9 taken 15 times.
✓ Branch 10 taken 1 times.
✓ Branch 11 taken 21 times.
✓ Branch 12 taken 75 times.
✓ Branch 13 taken 8 times.
✓ Branch 14 taken 39 times.
✓ Branch 15 taken 36 times.
✓ Branch 16 taken 39 times.
✗ Branch 17 not taken.
 234 count.unaryExpr([](int x) -> bool { return x > 1 && x % 2; }).any() &&
133 39 omp_get_max_threads() > 1) {
134 39 std::clog << "\n\tWARNING: forcing even cell count in all directions to prevent race "
135 39 "conditions with newton3.";
136
3/4
✓ Branch 0 taken 117 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 77 times.
✓ Branch 3 taken 40 times.
 156 count = count.unaryExpr([](int x) -> int { return x > 1 && x % 2 ? --x : x; });
137 }
138 #endif
139 83 return count;
140 }()),
141 83 cellL(system->L.array() / cellCount.cast<double>()),
142
1/2
✓ Branch 1 taken 83 times.
✗ Branch 2 not taken.
 83 cells(cellL, cellCount,
143
3/9
✓ Branch 1 taken 83 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 83 times.
✗ Branch 12 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 83 times.
✗ Branch 20 not taken.
✗ Branch 23 not taken.
 166 config["interaction"] && config["interaction"]["cellColoring"]
144
5/13
✗ Branch 0 not taken.
✓ Branch 1 taken 83 times.
✓ Branch 12 taken 83 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 83 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 83 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 83 times.
✗ Branch 23 not taken.
✗ Branch 26 not taken.
✗ Branch 29 not taken.
 166 ? config["interaction"]["cellColoring"].as<std::string>()
145 : "",
146
2/4
✓ Branch 2 taken 83 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 83 times.
✗ Branch 6 not taken.
 166 internal->threebody()) {
147
1/2
✓ Branch 1 taken 83 times.
✗ Branch 2 not taken.
 83 std::clog << "\n\tstrategy: cells";
148 83 }
149
150 4762 void CellNeighborInteraction::updateMeta() {
151
2/2
✓ Branch 0 taken 982414 times.
✓ Branch 1 taken 4762 times.
 987176 for (auto &c : cells) {
152 982414 c.particleIndices.clear();
153 }
154 4762 size_t n = particles.size();
155
2/2
✓ Branch 0 taken 16671 times.
✓ Branch 1 taken 4762 times.
 21433 for (size_t i = 0; i < n; ++i) {
156 16671 int cIndex = cells.getIndex(particles[i].r);
157 16671 cells[cIndex].particleIndices.push_back(i);
158 16671 cells.particle2cell.insert_or_assign(i, cIndex);
159 }
160 4762 }
161
162 2730 void CellNeighborInteraction::updateParticleMeta(Particle &p,
163 ParticleUpdateType particleUpdateType) {
164
2/2
✓ Branch 0 taken 66 times.
✓ Branch 1 taken 2664 times.
 2730 int i = std::distance(&particles.front(), &p);
165
2/2
✓ Branch 0 taken 66 times.
✓ Branch 1 taken 2664 times.
 2730 if (particleUpdateType != ParticleUpdateType::INSERT) { // p was deleted or moved
166
1/2
✓ Branch 1 taken 66 times.
✗ Branch 2 not taken.
 66 Cell &cOld = cells[cells.particle2cell[i]];
167
1/2
✓ Branch 0 taken 66 times.
✗ Branch 1 not taken.
 66 if (auto found = std::find(cOld.particleIndices.begin(), cOld.particleIndices.end(), i);
168
1/2
✓ Branch 0 taken 66 times.
✗ Branch 1 not taken.
 66 found != cOld.particleIndices.end()) {
169 66 std::iter_swap(found, cOld.particleIndices.end() - 1);
170 66 cOld.particleIndices.pop_back();
171 }
172 }
173
2/2
✓ Branch 0 taken 65 times.
✓ Branch 1 taken 1 times.
 66 if (particleUpdateType != ParticleUpdateType::DELETE) { // p was inserted or moved
174 2729 int cIndex = cells.getIndex(p.r);
175 2729 cells[cIndex].particleIndices.push_back(i);
176 2729 cells.particle2cell.insert_or_assign(i, cIndex);
177 }
178 2730 }
179
180 2723 double CellNeighborInteraction::iteratePairs(Particle &p, PairFunction &func) {
181 2723 double sum = 0.;
182 2723 Cell &c = cells[p.r];
183
2/2
✓ Branch 0 taken 39184 times.
✓ Branch 1 taken 2723 times.
 41907 for (auto &[_, icj] : c.pairsAll) {
184 39184 Cell &cj = cells[icj];
185 39184 int nj = cj.particleIndices.size();
186
2/2
✓ Branch 0 taken 293974 times.
✓ Branch 1 taken 39184 times.
 333158 for (int j = 0; j < nj; ++j) {
187
2/2
✓ Branch 0 taken 291251 times.
✓ Branch 1 taken 2723 times.
 293974 Particle &pj = particles[cj.particleIndices[j]];
188
2/2
✓ Branch 0 taken 291251 times.
✓ Branch 1 taken 2723 times.
 293974 if (&p != &pj) {
189 291251 sum += func(p, pj);
190 }
191 }
192 }
193 2723 return sum;
194 }
195
196 1489 double CellNeighborInteraction::iterateTriples(Particle &p, TripleFunction &func) {
197 1489 double sum = 0.;
198 1489 Cell &c = cells[p.r];
199
2/2
✓ Branch 0 taken 49509 times.
✓ Branch 1 taken 1489 times.
 50998 for (auto &[_, icj, ick] : c.triplesAll) {
200 49509 Cell &cj = cells[icj];
201 49509 Cell &ck = cells[ick];
202 49509 int nj = cj.particleIndices.size();
203 49509 int nk = ck.particleIndices.size();
204
2/2
✓ Branch 0 taken 1260077 times.
✓ Branch 1 taken 49509 times.
 1309586 for (int j = 0; j < nj; ++j) {
205
2/2
✓ Branch 0 taken 6977 times.
✓ Branch 1 taken 1253100 times.
 1260077 Particle &pj = particles[cj.particleIndices[j]];
206
2/2
✓ Branch 0 taken 6977 times.
✓ Branch 1 taken 1253100 times.
 1260077 if (&p == &pj) {
207 6977 continue;
208 }
209
2/2
✓ Branch 0 taken 235860 times.
✓ Branch 1 taken 1017240 times.
 1253100 int kStart = icj == ick ? j + 1 : 0;
210
2/2
✓ Branch 0 taken 26830388 times.
✓ Branch 1 taken 1253100 times.
 28083488 for (int k = kStart; k < nk; ++k) {
211
2/2
✓ Branch 0 taken 26676176 times.
✓ Branch 1 taken 154212 times.
 26830388 Particle &pk = particles[ck.particleIndices[k]];
212
2/2
✓ Branch 0 taken 26676176 times.
✓ Branch 1 taken 154212 times.
 26830388 if (&p != &pk) {
213 26676176 sum += func(p, pj, pk);
214 }
215 }
216 }
217 }
218 1489 return sum;
219 }
220
221 1024206 double CellNeighborInteraction::doCell(Cell &c, PairFunction &func, bool newton3) {
222 1024206 double sum = 0.;
223
2/2
✓ Branch 0 taken 1024096 times.
✓ Branch 1 taken 110 times.
 1024206 auto &pairs = newton3 ? c.pairsNewton3 : c.pairsAll;
224
2/2
✓ Branch 0 taken 14231079 times.
✓ Branch 1 taken 1024206 times.
 15255285 for (auto &[ici, icj] : pairs) {
225 14231079 Cell &ci = cells[ici];
226 14231079 Cell &cj = cells[icj];
227 14231079 int ni = ci.particleIndices.size();
228 14231079 int nj = cj.particleIndices.size();
229
2/2
✓ Branch 0 taken 227849 times.
✓ Branch 1 taken 14231079 times.
 14458928 for (int i = 0; i < ni; ++i) {
230
2/2
✓ Branch 0 taken 212716 times.
✓ Branch 1 taken 15133 times.
 227849 Particle &pi = particles[ci.particleIndices[i]];
231
4/4
✓ Branch 0 taken 212716 times.
✓ Branch 1 taken 15133 times.
✓ Branch 2 taken 17499 times.
✓ Branch 3 taken 195217 times.
 227849 int jStart = newton3 && ici == icj ? i + 1 : 0;
232
2/2
✓ Branch 0 taken 445090 times.
✓ Branch 1 taken 227849 times.
 672939 for (int j = jStart; j < nj; ++j) {
233
2/2
✓ Branch 0 taken 443959 times.
✓ Branch 1 taken 1131 times.
 445090 Particle &pj = particles[cj.particleIndices[j]];
234
2/2
✓ Branch 0 taken 443959 times.
✓ Branch 1 taken 1131 times.
 445090 if (&pi != &pj) {
235 443959 sum += func(pi, pj);
236 }
237 }
238 }
239 }
240 1024206 return sum;
241 }
242
243 5298 double CellNeighborInteraction::doCell(Cell &c, TripleFunction &func, bool newton3) {
244 5298 double sum = 0.;
245
2/2
✓ Branch 0 taken 5280 times.
✓ Branch 1 taken 18 times.
 5298 auto &triples = newton3 ? c.triplesNewton3 : c.triplesAll;
246
2/2
✓ Branch 0 taken 190098 times.
✓ Branch 1 taken 5298 times.
 195396 for (auto &[ici, icj, ick] : triples) {
247 190098 Cell &ci = cells[ici];
248 190098 Cell &cj = cells[icj];
249 190098 Cell &ck = cells[ick];
250 190098 int ni = ci.particleIndices.size();
251 190098 int nj = cj.particleIndices.size();
252 190098 int nk = ck.particleIndices.size();
253
2/2
✓ Branch 0 taken 120385 times.
✓ Branch 1 taken 190098 times.
 310483 for (int i = 0; i < ni; ++i) {
254 120385 Particle &pi = particles[ci.particleIndices[i]];
255
2/2
✓ Branch 0 taken 2390766 times.
✓ Branch 1 taken 120385 times.
 2511151 for (int j = 0; j < nj; ++j) {
256
2/2
✓ Branch 0 taken 18760 times.
✓ Branch 1 taken 2372006 times.
 2390766 Particle &pj = particles[cj.particleIndices[j]];
257
2/2
✓ Branch 0 taken 18760 times.
✓ Branch 1 taken 2372006 times.
 2390766 if (&pi == &pj) {
258 18760 continue;
259 }
260
2/2
✓ Branch 0 taken 475536 times.
✓ Branch 1 taken 1896470 times.
 2372006 int kStart = icj == ick ? j + 1 : 0;
261
2/2
✓ Branch 0 taken 80267115 times.
✓ Branch 1 taken 2372006 times.
 82639121 for (int k = kStart; k < nk; ++k) {
262
2/2
✓ Branch 0 taken 80029347 times.
✓ Branch 1 taken 237768 times.
 80267115 Particle &pk = particles[ck.particleIndices[k]];
263
2/2
✓ Branch 0 taken 80029347 times.
✓ Branch 1 taken 237768 times.
 80267115 if (&pi != &pk) {
264 80029347 sum += func(pi, pj, pk);
265 }
266 }
267 }
268 }
269 }
270 5298 return sum;
271 }
272
273 template <typename T>
274 12186 double CellNeighborInteraction::iterate(T &func, bool newton3, bool parallel) {
275 12186 double sum = 0.;
276
2/2
✓ Branch 0 taken 6080 times.
✓ Branch 1 taken 13 times.
 12186 if (newton3) {
277
2/2
✓ Branch 0 taken 97280 times.
✓ Branch 1 taken 6080 times.
 206720 for (auto &sameColorIndices : cells.sameColorIndices) {
278 194560 #pragma omp parallel for schedule(dynamic) if (parallel) reduction(+ : sum)
279 for (int ic : sameColorIndices) {
280 Cell &c = cells[ic];
281 sum += doCell(c, func, newton3);
282 }
283 }
284 } else {
285 26 #pragma omp parallel for schedule(dynamic) if (parallel) reduction(+ : sum)
286 for (Cell &c : cells) {
287 sum += doCell(c, func, newton3);
288 }
289 }
290 12186 return sum;
291 }
292
293 5415 double CellNeighborInteraction::iteratePairs(PairFunction &func, bool newton3, bool parallel) {
294 5415 return iterate(func, newton3, parallel);
295 }
296
297 678 double CellNeighborInteraction::iterateTriples(TripleFunction &func, bool newton3, bool parallel) {
298 678 return iterate(func, newton3, parallel);
299 }
300
301 99 Interaction *InteractionFactory::create(System *system, const YAML::Node &configInteraction) {
302 99 Interaction *interaction = nullptr;
303
11/16
✓ Branch 1 taken 99 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 50 times.
✓ Branch 4 taken 49 times.
✓ Branch 6 taken 50 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 50 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 49 times.
✓ Branch 12 taken 1 times.
✓ Branch 14 taken 50 times.
✓ Branch 15 taken 49 times.
✓ Branch 18 taken 16 times.
✓ Branch 19 taken 83 times.
✗ Branch 22 not taken.
✗ Branch 25 not taken.
 198 if (configInteraction["strategy"] && configInteraction["strategy"]["neighbors"] &&
304
13/19
✓ Branch 1 taken 49 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 49 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 49 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 33 times.
✓ Branch 10 taken 16 times.
✓ Branch 11 taken 49 times.
✓ Branch 12 taken 50 times.
✓ Branch 13 taken 49 times.
✓ Branch 14 taken 50 times.
✓ Branch 16 taken 49 times.
✓ Branch 17 taken 50 times.
✓ Branch 19 taken 50 times.
✓ Branch 20 taken 49 times.
✗ Branch 22 not taken.
✗ Branch 25 not taken.
✗ Branch 28 not taken.
 181 configInteraction["strategy"]["neighbors"].as<std::string>() == "all") {
305
1/2
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
 16 interaction = new AllInteraction(system, configInteraction);
306 } else {
307
1/2
✓ Branch 2 taken 83 times.
✗ Branch 3 not taken.
 83 interaction = new CellNeighborInteraction(system, configInteraction);
308 }
309 99 return interaction;
310 }
311