GCC Code Coverage Report

Directory: ./
File: src/cells.cpp
Date: 2024-04-18 12:22:13
Exec Total Coverage
Lines: 124 128 96.9%
Functions: 5 5 100.0%
Branches: 123 168 73.2%
Line Branch Exec Source
1 #include <algorithm>
2 #include <iostream>
3 #include <stdexcept>
4 #include <unordered_set>
5
6 #include "cells.hpp"
7 #include "particle.hpp"
8
9 88 static int getNumColors(const std::string &mode) {
10 88 int colors;
11
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 88 times.
 88 if (mode == "c08") {
12 colors = 8;
13 ✗ } else if (mode == "c04") {
14 colors = 4;
15 } else {
16 ✗ throw std::invalid_argument("Cell color scheme " + mode + " not implemented");
17 }
18 88 return colors;
19 }
20
21 88 Cells::Cells(const Vec &cellL, const ArrayI &count, const std::string &colorScheme, bool triples)
22
2/4
✓ Branch 4 taken 88 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 88 times.
✗ Branch 7 not taken.
 88 : cellLinv(1. / cellL.array()), multiindex(count), colorScheme(colorScheme) {
23
1/2
✓ Branch 0 taken 88 times.
✗ Branch 1 not taken.
 88 if (this->colorScheme.empty()) {
24 // Set default color scheme if none given
25 88 if constexpr (DIM == 3) {
26
1/2
✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
 88 this->colorScheme = "c08";
27 }
28 if constexpr (DIM == 2) {
29 this->colorScheme = "c04";
30 }
31 }
32
1/2
✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
 88 this->cells.resize(multiindex.linsize);
33
2/2
✓ Branch 0 taken 7550 times.
✓ Branch 1 taken 88 times.
 7638 for (int i = 0; i < multiindex.linsize; ++i) {
34
1/2
✓ Branch 1 taken 7550 times.
✗ Branch 2 not taken.
 7550 this->cells[i].particleIndices = {};
35
3/6
✓ Branch 1 taken 7550 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7550 times.
✗ Branch 5 not taken.
✓ Branch 10 taken 7550 times.
✗ Branch 11 not taken.
 15100 this->cells[i].pairsAll = getPairs(i, "all");
36
1/2
✓ Branch 1 taken 7550 times.
✗ Branch 2 not taken.
 7550 this->cells[i].pairsNewton3 = getPairs(i, this->colorScheme);
37
2/2
✓ Branch 0 taken 54 times.
✓ Branch 1 taken 7496 times.
 7550 if (triples) {
38
3/6
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 54 times.
✗ Branch 5 not taken.
✓ Branch 10 taken 54 times.
✗ Branch 11 not taken.
 108 this->cells[i].triplesAll = getTriples(i, "all");
39
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 54 this->cells[i].triplesNewton3 = getTriples(i, this->colorScheme);
40 }
41 }
42
1/2
✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
 88 int colors = getNumColors(this->colorScheme);
43
1/2
✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
 88 this->sameColorIndices.resize(colors);
44
2/2
✓ Branch 0 taken 704 times.
✓ Branch 1 taken 88 times.
 792 for (int c = 0; c < colors; c++) {
45
1/2
✓ Branch 1 taken 704 times.
✗ Branch 2 not taken.
 1408 this->sameColorIndices.push_back(getSameColorIndices(c, this->colorScheme));
46 }
47
2/4
✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 88 times.
✗ Branch 5 not taken.
 88 std::clog << "\n\tcell color scheme: " << this->colorScheme;
48
1/2
✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
 88 std::clog << "\n\tcell count:";
49
2/2
✓ Branch 0 taken 264 times.
✓ Branch 1 taken 88 times.
 352 for (int d = 0; d < DIM; ++d) {
50
2/4
✓ Branch 1 taken 264 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 264 times.
✗ Branch 6 not taken.
 264 std::clog << " " << this->multiindex.shape[d];
51 }
52
1/2
✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
 88 std::clog << "\n\tcell size:";
53
2/2
✓ Branch 0 taken 264 times.
✓ Branch 1 taken 88 times.
 352 for (int d = 0; d < DIM; ++d) {
54
2/4
✓ Branch 1 taken 264 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 264 times.
✗ Branch 6 not taken.
 264 std::clog << " " << cellL[d];
55 }
56 88 }
57
58 15103 std::vector<std::pair<int, int>> Cells::getPairs(int cellIndex, const std::string &mode) {
59 15103 ArrayI cellMulti = multiindex.lin2multi(cellIndex);
60
2/2
✓ Branch 0 taken 7550 times.
✓ Branch 1 taken 7553 times.
 15103 std::vector<std::pair<ArrayI, ArrayI>> shiftPairs;
61
2/2
✓ Branch 0 taken 7550 times.
✓ Branch 1 taken 7553 times.
 15103 if (mode == "all") {
62 7550 ArrayI shift1 = ArrayI::Zero();
63
2/2
✓ Branch 2 taken 203850 times.
✓ Branch 3 taken 7550 times.
 211400 for (ArrayI shift2 : Multiindex<DIM>(ArrayI(3))) {
64 203850 shift2 -= ArrayI(1);
65
1/2
✓ Branch 1 taken 203850 times.
✗ Branch 2 not taken.
 203850 shiftPairs.push_back({shift1, shift2});
66 }
67
1/3
✗ Branch 0 not taken.
✓ Branch 1 taken 7553 times.
✗ Branch 3 not taken.
 7553 } else if (mode == "c08" || mode == "c04") {
68 7553 Multiindex<DIM> range(ArrayI(2));
69
2/2
✓ Branch 0 taken 60424 times.
✓ Branch 1 taken 7553 times.
 67977 for (ArrayI shift1 : range) {
70
2/2
✓ Branch 0 taken 483392 times.
✓ Branch 1 taken 60424 times.
 543816 for (ArrayI shift2 : range) {
71
2/2
✓ Branch 1 taken 279461 times.
✓ Branch 2 taken 203931 times.
 483392 if ((shift1 * shift2 == 1).any()) {
72 279461 continue;
73 }
74
1/2
✓ Branch 1 taken 203931 times.
✗ Branch 2 not taken.
 203931 shiftPairs.push_back({shift1, shift2});
75 }
76 }
77 } else if (mode == "c04") {
78 throw std::invalid_argument("Cell mode " + mode + " not implemented");
79 }
80 // Build up the cell pair list. This is quite involved because we aim for a general solution that
81 // also works in narrow systems, where periodic boundary conditions make our lives hard...
82 15103 std::set<std::pair<int, int>> pairsSet; // Use set for construction to ensure uniqueness
83 // Now go through all *potential* cell pairs
84 15103 std::pair<int, int> pair;
85
2/2
✓ Branch 1 taken 407781 times.
✓ Branch 2 taken 15103 times.
 422884 for (auto &[shift1, shift2] : shiftPairs) {
86 407781 ArrayI multi1 = cellMulti + shift1;
87 407781 ArrayI multi2 = cellMulti + shift2;
88
4/4
✓ Branch 3 taken 406515 times.
✓ Branch 4 taken 1266 times.
✓ Branch 5 taken 401268 times.
✓ Branch 6 taken 5247 times.
 814296 if ((multiindex.shape == 1 && (multi1 == 1 || multi2 == 1)).any() ||
89
2/2
✓ Branch 3 taken 401268 times.
✓ Branch 4 taken 5247 times.
 406515 (multiindex.shape == 2 && (multi1 == 2 || multi2 == 2)).any()) {
90 // If there are less than 3 cells in a direction, ignore the shifts that invoke the periodic
91 // boundary conditions. This way, double counting is prevented.
92 6513 continue;
93 }
94 // Transform to linear index
95
2/2
✓ Branch 2 taken 201389 times.
✓ Branch 3 taken 199879 times.
 401268 pair = {multiindex.multi2lin(multi1, true), multiindex.multi2lin(multi2, true)};
96
3/4
✓ Branch 0 taken 201389 times.
✓ Branch 1 taken 199879 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 201389 times.
 401268 if (mode == "c08" || mode == "c04") {
97 // For newton3 modes: sort to avoid pair duplicates that only differ in order
98
2/2
✓ Branch 0 taken 96163 times.
✓ Branch 1 taken 103716 times.
 199879 auto [c1, c2] = pair;
99
2/2
✓ Branch 0 taken 96163 times.
✓ Branch 1 taken 103716 times.
 199879 if (c1 > c2) {
100 96163 pair = {c2, c1};
101 }
102 }
103
1/2
✓ Branch 1 taken 401268 times.
✗ Branch 2 not taken.
 401268 pairsSet.insert(pair);
104 }
105 // Go back to vector for random access performance
106 15103 std::vector<std::pair<int, int>> pairs;
107
2/2
✓ Branch 0 taken 303556 times.
✓ Branch 1 taken 15103 times.
 318659 for (auto &cp : pairsSet) {
108
1/2
✓ Branch 1 taken 303556 times.
✗ Branch 2 not taken.
 303556 pairs.push_back(cp);
109 }
110 30206 return pairs;
111 15103 }
112
113 108 std::vector<std::tuple<int, int, int>> Cells::getTriples(int cellIndex, const std::string &mode) {
114 108 ArrayI cellMulti = multiindex.lin2multi(cellIndex);
115
2/2
✓ Branch 0 taken 54 times.
✓ Branch 1 taken 54 times.
 108 std::vector<std::tuple<ArrayI, ArrayI, ArrayI>> shiftTriples;
116
2/2
✓ Branch 0 taken 54 times.
✓ Branch 1 taken 54 times.
 108 if (mode == "all") {
117 54 ArrayI shift1 = ArrayI::Zero();
118 54 ArrayI start(-1);
119 54 Multiindex<DIM> range(ArrayI(3));
120
2/2
✓ Branch 0 taken 1458 times.
✓ Branch 1 taken 54 times.
 1512 for (ArrayI shift2 : range) {
121 1458 shift2 += start;
122
2/2
✓ Branch 0 taken 39366 times.
✓ Branch 1 taken 1458 times.
 40824 for (ArrayI shift3 : range) {
123 39366 shift3 += start;
124
2/2
✓ Branch 1 taken 20844 times.
✓ Branch 2 taken 18522 times.
 39366 if (((shift2 - shift3).abs() > 1).any()) {
125 // cell 2 and 3 are not neighbors, ignore those shifts
126 20844 continue;
127 }
128
1/2
✓ Branch 1 taken 18522 times.
✗ Branch 2 not taken.
 18522 shiftTriples.push_back({shift1, shift2, shift3});
129 }
130 }
131
1/3
✗ Branch 0 not taken.
✓ Branch 1 taken 54 times.
✗ Branch 3 not taken.
 54 } else if (mode == "c08" || mode == "c04") {
132 54 Multiindex<DIM> range(ArrayI(2));
133
2/2
✓ Branch 0 taken 432 times.
✓ Branch 1 taken 54 times.
 486 for (ArrayI shift1 : range) {
134
2/2
✓ Branch 0 taken 3456 times.
✓ Branch 1 taken 432 times.
 3888 for (ArrayI shift2 : range) {
135
2/2
✓ Branch 0 taken 27648 times.
✓ Branch 1 taken 3456 times.
 31104 for (ArrayI shift3 : range) {
136
2/2
✓ Branch 1 taken 9126 times.
✓ Branch 2 taken 18522 times.
 27648 if ((shift1 * shift2 * shift3 == 1).any()) {
137 9126 continue;
138 }
139
1/2
✓ Branch 1 taken 18522 times.
✗ Branch 2 not taken.
 18522 shiftTriples.push_back({shift1, shift2, shift3});
140 }
141 }
142 }
143 } else {
144 ✗ throw std::invalid_argument("Cell mode " + mode + " not implemented");
145 }
146 108 std::set<std::tuple<int, int, int>> triplesSet; // Use set for construction to ensure uniqueness
147 108 std::tuple<int, int, int> triple;
148
2/2
✓ Branch 1 taken 37044 times.
✓ Branch 2 taken 108 times.
 37152 for (auto &[shift1, shift2, shift3] : shiftTriples) {
149 37044 ArrayI multi1 = cellMulti + shift1;
150 37044 ArrayI multi2 = cellMulti + shift2;
151 37044 ArrayI multi3 = cellMulti + shift3;
152
4/4
✓ Branch 4 taken 28350 times.
✓ Branch 5 taken 8694 times.
✓ Branch 6 taken 10125 times.
✓ Branch 7 taken 18225 times.
 65394 if ((multiindex.shape == 1 && (multi1 == 1 || multi2 == 1 || multi3 == 1)).any() ||
153
2/2
✓ Branch 4 taken 10125 times.
✓ Branch 5 taken 18225 times.
 28350 (multiindex.shape == 2 && (multi1 == 2 || multi2 == 2 || multi3 == 2)).any()) {
154 // If there are less than 3 cells in a direction, ignore the shifts that invoke the periodic
155 // boundary conditions. This way, double counting is prevented.
156 26919 continue;
157 }
158 // Transform to linear index
159 10125 triple = {multiindex.multi2lin(multi1, true), multiindex.multi2lin(multi2, true),
160
2/2
✓ Branch 1 taken 4135 times.
✓ Branch 2 taken 5990 times.
 10125 multiindex.multi2lin(multi3, true)};
161 // sort to avoid duplicates that only differ in order of second and third cell
162
2/2
✓ Branch 0 taken 4135 times.
✓ Branch 1 taken 5990 times.
 10125 auto [c1, c2, c3] = triple;
163
2/2
✓ Branch 0 taken 4135 times.
✓ Branch 1 taken 5990 times.
 10125 if (c2 > c3) {
164 4135 triple = {c1, c3, c2};
165 }
166
1/2
✓ Branch 1 taken 10125 times.
✗ Branch 2 not taken.
 10125 triplesSet.insert(triple);
167 }
168 // Go back to vector for random access performance
169 108 std::vector<std::tuple<int, int, int>> triples;
170
2/2
✓ Branch 0 taken 2908 times.
✓ Branch 1 taken 108 times.
 3016 for (auto &ct : triplesSet) {
171
1/2
✓ Branch 1 taken 2908 times.
✗ Branch 2 not taken.
 2908 triples.push_back(ct);
172 }
173 216 return triples;
174 108 }
175
176 704 std::vector<int> Cells::getSameColorIndices(int color, const std::string &mode) {
177
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 704 times.
 704 std::vector<int> sameColorCells;
178
1/3
✗ Branch 0 not taken.
✓ Branch 1 taken 704 times.
✗ Branch 3 not taken.
 704 if (mode == "c08" || mode == "c04") {
179 704 std::vector<ArrayI> startingMulti;
180 704 Multiindex<DIM> rangeStartingMulti(ArrayI(2));
181
2/2
✓ Branch 0 taken 5632 times.
✓ Branch 1 taken 704 times.
 6336 for (ArrayI multi : rangeStartingMulti) {
182
1/2
✓ Branch 1 taken 5632 times.
✗ Branch 2 not taken.
 5632 startingMulti.push_back(multi);
183 }
184 704 ArrayI start = startingMulti[color];
185
2/2
✓ Branch 2 taken 130 times.
✓ Branch 3 taken 574 times.
 704 if ((multiindex.shape == 1 && start == 1).any()) {
186 // We would go around the PBC, so skip this color (i.e. leave the corresponding index list
187 // empty)
188 130 return {};
189 }
190 574 Multiindex<DIM> rangeSameColorCells(multiindex.shape - start);
191
2/2
✓ Branch 0 taken 44900 times.
✓ Branch 1 taken 574 times.
 45474 for (ArrayI multi : rangeSameColorCells) {
192
2/2
✓ Branch 0 taken 37350 times.
✓ Branch 1 taken 7550 times.
 44900 if (multi.unaryExpr([](int x) -> bool { return x % 2; }).any()) {
193 // Ignore odd differences
194 37350 continue;
195 }
196 7550 int linind = multiindex.multi2lin(start + multi, true);
197
1/2
✓ Branch 1 taken 7550 times.
✗ Branch 2 not taken.
 7550 sameColorCells.push_back(linind);
198 }
199 704 } else {
200 ✗ throw std::invalid_argument("Cell mode " + mode + " not implemented");
201 }
202 574 return sameColorCells;
203 704 }
204