GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/assembly/coloring.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 81 82 98.8%
Functions: 268 280 95.7%
Branches: 190 257 73.9%
Line Branch Exec Source
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 //
4 // SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5 // SPDX-License-Identifier: GPL-3.0-or-later
6 //
7 /*!
8 * \file
9 * \ingroup Assembly
10 * \brief Coloring schemes for shared-memory-parallel assembly
11 */
12 #ifndef DUMUX_ASSEMBLY_COLORING_HH
13 #define DUMUX_ASSEMBLY_COLORING_HH
14
15 #include <vector>
16 #include <deque>
17 #include <iostream>
18 #include <tuple>
19 #include <algorithm>
20
21 #include <dune/common/timer.hh>
22 #include <dune/common/exceptions.hh>
23
24 #include <dumux/io/format.hh>
25 #include <dumux/discretization/method.hh>
26
27 #ifndef DOXYGEN // hide from doxygen
28 namespace Dumux::Detail {
29
30 //! Compute a map from dof indices to element indices (helper data for coloring algorithm)
31 template <class GridGeometry>
32 std::vector<std::vector<std::size_t>>
33 500 computeConnectedElements(const GridGeometry& gg)
34 {
35 500 std::vector<std::vector<std::size_t>> connectedElements;
36
37 if constexpr (GridGeometry::discMethod == DiscretizationMethods::cctpfa)
38 {
39
3/5
✓ Branch 1 taken 188 times.
✓ Branch 2 taken 34 times.
✓ Branch 4 taken 184 times.
✗ Branch 5 not taken.
✗ Branch 3 not taken.
 269 connectedElements.resize(gg.gridView().size(0));
40
1/2
✓ Branch 1 taken 184 times.
✗ Branch 2 not taken.
 269 const auto& eMapper = gg.elementMapper();
41
8/10
✓ Branch 1 taken 13341 times.
✓ Branch 2 taken 5512 times.
✓ Branch 4 taken 35932 times.
✓ Branch 5 taken 12 times.
✓ Branch 7 taken 25 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 62076 times.
✓ Branch 10 taken 25 times.
✓ Branch 3 taken 1281404 times.
✗ Branch 6 not taken.
 4538864 for (const auto& element : elements(gg.gridView()))
42 {
43
1/2
✓ Branch 1 taken 80514 times.
✗ Branch 2 not taken.
 2310056 const auto eIdx = eMapper.index(element);
44
17/21
✓ Branch 1 taken 6886829 times.
✓ Branch 2 taken 121600 times.
✓ Branch 7 taken 267427 times.
✓ Branch 8 taken 65729 times.
✓ Branch 9 taken 311065 times.
✓ Branch 10 taken 131168 times.
✓ Branch 11 taken 1468 times.
✓ Branch 12 taken 103988 times.
✓ Branch 17 taken 62076 times.
✓ Branch 18 taken 18438 times.
✓ Branch 0 taken 1232424 times.
✓ Branch 4 taken 6754539 times.
✓ Branch 5 taken 69511 times.
✓ Branch 6 taken 77486 times.
✗ Branch 13 not taken.
✗ Branch 19 not taken.
✓ Branch 20 taken 81504 times.
✓ Branch 21 taken 4046 times.
✓ Branch 3 taken 134060 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
 18832081 for (const auto& intersection : intersections(gg.gridView(), element))
45
2/4
✓ Branch 0 taken 493634 times.
✓ Branch 1 taken 9591 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 25083583 if (intersection.neighbor())
46
3/5
✓ Branch 1 taken 338864 times.
✓ Branch 2 taken 6803380 times.
✓ Branch 4 taken 327424 times.
✗ Branch 5 not taken.
✗ Branch 3 not taken.
 12673464 connectedElements[eMapper.index(intersection.outside())].push_back(eIdx);
47 }
48 }
49
50 else if constexpr (GridGeometry::discMethod == DiscretizationMethods::box
51 || GridGeometry::discMethod == DiscretizationMethods::pq1bubble)
52 {
53 static constexpr int dim = GridGeometry::GridView::dimension;
54
3/5
✓ Branch 1 taken 146 times.
✓ Branch 2 taken 27 times.
✓ Branch 4 taken 145 times.
✗ Branch 5 not taken.
✗ Branch 3 not taken.
 179 connectedElements.resize(gg.gridView().size(dim));
55
1/2
✓ Branch 1 taken 145 times.
✗ Branch 2 not taken.
 179 const auto& vMapper = gg.vertexMapper();
56
10/12
✓ Branch 1 taken 145 times.
✓ Branch 2 taken 42447 times.
✓ Branch 4 taken 159739 times.
✓ Branch 5 taken 1 times.
✓ Branch 7 taken 45 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 35631 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 35631 times.
✓ Branch 13 taken 16 times.
✓ Branch 3 taken 10028 times.
✓ Branch 6 taken 32603 times.
 665521 for (const auto& element : elements(gg.gridView()))
57 {
58
1/2
✓ Branch 1 taken 68234 times.
✗ Branch 2 not taken.
 261411 const auto eIdx = gg.elementMapper().index(element);
59
5/5
✓ Branch 1 taken 548868 times.
✓ Branch 2 taken 82447 times.
✓ Branch 3 taken 281742 times.
✓ Branch 4 taken 68234 times.
✓ Branch 0 taken 538008 times.
 1571039 for (int i = 0; i < element.subEntities(dim); i++)
60
2/4
✓ Branch 1 taken 941622 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 937590 times.
✗ Branch 5 not taken.
 1098166 connectedElements[vMapper.subIndex(element, i, dim)].push_back(eIdx);
61 }
62 }
63
64 else if constexpr (
65 GridGeometry::discMethod == DiscretizationMethods::fcstaggered
66 || GridGeometry::discMethod == DiscretizationMethods::ccmpfa
67 )
68 {
69 // for MPFA-O schemes the assembly of each element residual touches all vertex neighbors
70 // for face-centered staggered it is all codim-2 neighbors (vertex neighbors in 2D, edge neighbors in 3D)
71 // but we use vertex neighbors also in 3D for simplicity
72 37 std::vector<std::vector<std::size_t>> vToElements;
73 static constexpr int dim = GridGeometry::GridView::dimension;
74
3/5
✓ Branch 1 taken 34 times.
✓ Branch 2 taken 3 times.
✓ Branch 4 taken 34 times.
✗ Branch 5 not taken.
✗ Branch 3 not taken.
 37 vToElements.resize(gg.gridView().size(dim));
75
1/2
✓ Branch 1 taken 34 times.
✗ Branch 2 not taken.
 37 const auto& vMapper = gg.vertexMapper();
76
9/12
✓ Branch 1 taken 34 times.
✓ Branch 2 taken 6978 times.
✓ Branch 4 taken 23788 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 17 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 12782 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 12782 times.
✓ Branch 13 taken 3 times.
✓ Branch 6 taken 14484 times.
✓ Branch 3 taken 3 times.
 96386 for (const auto& element : elements(gg.gridView()))
77 {
78
1/2
✓ Branch 1 taken 27266 times.
✗ Branch 2 not taken.
 43528 const auto eIdx = gg.elementMapper().index(element);
79
5/5
✓ Branch 1 taken 166548 times.
✓ Branch 2 taken 6978 times.
✓ Branch 3 taken 109064 times.
✓ Branch 4 taken 27266 times.
✓ Branch 0 taken 37136 times.
 274572 for (int i = 0; i < element.subEntities(dim); i++)
80
2/4
✓ Branch 1 taken 167134 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 167134 times.
✗ Branch 5 not taken.
 167134 vToElements[vMapper.subIndex(element, i, dim)].push_back(eIdx);
81 }
82
83
3/5
✓ Branch 1 taken 34 times.
✓ Branch 2 taken 3 times.
✓ Branch 4 taken 34 times.
✗ Branch 5 not taken.
✗ Branch 3 not taken.
 37 connectedElements.resize(gg.gridView().size(0));
84
7/9
✓ Branch 1 taken 7012 times.
✓ Branch 2 taken 3 times.
✓ Branch 4 taken 17 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 12782 times.
✓ Branch 10 taken 3 times.
✓ Branch 3 taken 23785 times.
 52858 for (const auto& element : elements(gg.gridView()))
85 {
86
1/2
✓ Branch 1 taken 27266 times.
✗ Branch 2 not taken.
 43528 const auto eIdx = gg.elementMapper().index(element);
87
5/5
✓ Branch 1 taken 166548 times.
✓ Branch 2 taken 6978 times.
✓ Branch 3 taken 109064 times.
✓ Branch 4 taken 27266 times.
✓ Branch 0 taken 37136 times.
 274572 for (int i = 0; i < element.subEntities(dim); i++)
88 {
89
2/4
✓ Branch 1 taken 167134 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 167134 times.
✗ Branch 5 not taken.
 167134 const auto& e = vToElements[vMapper.subIndex(element, i, dim)];
90
1/2
✓ Branch 1 taken 167134 times.
✗ Branch 2 not taken.
 167134 connectedElements[eIdx].insert(connectedElements[eIdx].end(), e.begin(), e.end());
91 }
92
93 // make unique
94 43528 std::sort(connectedElements[eIdx].begin(), connectedElements[eIdx].end());
95
1/2
✓ Branch 3 taken 27266 times.
✗ Branch 4 not taken.
 43528 connectedElements[eIdx].erase(
96 43528 std::unique(connectedElements[eIdx].begin(), connectedElements[eIdx].end()),
97 43528 connectedElements[eIdx].end()
98 );
99 }
100 37 }
101
102 // nothing has to be precomputed here as only immediate face neighbors are connected
103 else if constexpr (GridGeometry::discMethod == DiscretizationMethods::fcdiamond)
104 return connectedElements;
105
106 else
107 DUNE_THROW(Dune::NotImplemented,
108 "Missing coloring scheme implementation for this discretization method"
109 );
110
111 485 return connectedElements;
112 ✗ }
113
114 /*!
115 * \brief Compute the colors of neighboring nodes in the dependency graph
116 *
117 * Neighboring nodes are those elements that manipulate the
118 * same data structures (e.g. system matrix, volvars, flux cache) in the same places
119 *
120 * \param gridGeometry the grid geometry
121 * \param element the element we want to color
122 * \param colors a vector of current colors for each element (not assigned: -1)
123 * \param connectedElements a map from implementation-defined indices to element indices
124 * \param neighborColors a vector to add the colors of neighbor nodes to
125 */
126 template<class GridGeometry, class ConnectedElements>
127 2646473 void addNeighborColors(const GridGeometry& gg,
128 const typename GridGeometry::LocalView::Element& element,
129 const std::vector<int>& colors,
130 const ConnectedElements& connectedElements,
131 std::vector<int>& neighborColors)
132 {
133 if constexpr (
134 GridGeometry::discMethod == DiscretizationMethods::cctpfa
135 || GridGeometry::discMethod == DiscretizationMethods::ccmpfa
136 || GridGeometry::discMethod == DiscretizationMethods::fcstaggered
137 )
138 {
139 // we modify neighbor elements during the assembly
140 // check who else modifies these neighbor elements
141 2353584 const auto& eMapper = gg.elementMapper();
142
12/15
✓ Branch 2 taken 6979330 times.
✓ Branch 3 taken 121600 times.
✓ Branch 5 taken 333359 times.
✓ Branch 6 taken 124607 times.
✓ Branch 8 taken 209970 times.
✓ Branch 9 taken 2102 times.
✓ Branch 10 taken 176408 times.
✗ Branch 11 not taken.
✓ Branch 15 taken 137396 times.
✓ Branch 16 taken 6090 times.
✓ Branch 7 taken 375295 times.
✓ Branch 1 taken 1241708 times.
✓ Branch 4 taken 98986 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 17106035 for (const auto& intersection : intersections(gg.gridView(), element))
143 {
144
2/4
✓ Branch 0 taken 621144 times.
✓ Branch 1 taken 12965 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 25287005 if (intersection.neighbor())
145 {
146 // direct face neighbors
147
1/2
✓ Branch 1 taken 445508 times.
✗ Branch 2 not taken.
 12836376 const auto nIdx = eMapper.index(intersection.outside());
148
1/2
✓ Branch 1 taken 621144 times.
✗ Branch 2 not taken.
 12836376 neighborColors.push_back(colors[nIdx]);
149
150 // neighbor-neighbors
151
8/8
✓ Branch 1 taken 459212 times.
✓ Branch 2 taken 38423604 times.
✓ Branch 5 taken 1615109 times.
✓ Branch 6 taken 875740 times.
✓ Branch 7 taken 1576869 times.
✓ Branch 8 taken 296672 times.
✓ Branch 3 taken 6709156 times.
✓ Branch 4 taken 1748846 times.
 86480891 for (const auto nnIdx : connectedElements[eMapper.index(intersection.outside())])
152
1/2
✓ Branch 1 taken 3202023 times.
✗ Branch 2 not taken.
 73644515 neighborColors.push_back(colors[nnIdx]);
153 }
154 }
155 }
156
157 else if constexpr (GridGeometry::discMethod == DiscretizationMethods::box
158 || GridGeometry::discMethod == DiscretizationMethods::pq1bubble)
159 {
160 // we modify the vertex dofs of our element during the assembly
161 // check who else modifies these vertex dofs
162 261411 const auto& vMapper = gg.vertexMapper();
163 static constexpr int dim = GridGeometry::GridView::dimension;
164 // direct vertex neighbors
165
3/3
✓ Branch 0 taken 629880 times.
✓ Branch 1 taken 471209 times.
✓ Branch 2 taken 78234 times.
 1571039 for (int i = 0; i < element.subEntities(dim); i++)
166
3/3
✓ Branch 2 taken 5555210 times.
✓ Branch 3 taken 937590 times.
✓ Branch 1 taken 8124 times.
 8325168 for (auto eIdx : connectedElements[vMapper.subIndex(element, i, dim)])
167 7227002 neighborColors.push_back(colors[eIdx]);
168 }
169
170 else if constexpr (GridGeometry::discMethod == DiscretizationMethods::fcdiamond)
171 {
172 // we modify neighbor faces during the assembly
173 // check who else modifies these neighbor elements
174 31478 const auto& eMapper = gg.elementMapper();
175
8/12
✓ Branch 2 taken 84500 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 32778 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 400 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 500 times.
✗ Branch 11 not taken.
✓ Branch 15 taken 360 times.
✓ Branch 16 taken 40 times.
✓ Branch 7 taken 42956 times.
✓ Branch 1 taken 21100 times.
 235928 for (const auto& intersection : intersections(gg.gridView(), element))
176
2/2
✓ Branch 0 taken 32398 times.
✓ Branch 1 taken 680 times.
 200298 if (intersection.neighbor())
177
2/4
✓ Branch 1 taken 32398 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 32398 times.
✗ Branch 5 not taken.
 147256 neighborColors.push_back(colors[eMapper.index(intersection.outside())]);
178 }
179
180 else
181 DUNE_THROW(Dune::NotImplemented,
182 "Missing coloring scheme implementation for this discretization method"
183 );
184 2646473 }
185
186 /*!
187 * \brief Find the smallest color (integer >= 0) _not_ present in the given list of colors
188 * \param colors list of colors which are already taken
189 * \param notAssigned container to store which colors are not yet taken (is resized as required)
190 */
191 2175767 int smallestAvailableColor(const std::vector<int>& colors,
192 std::vector<bool>& colorUsed)
193 {
194 2175767 const int numColors = colors.size();
195 2175767 colorUsed.assign(numColors, false);
196
197 // The worst case for e.g. numColors=3 is colors={0, 1, 2}
198 // in which case we return 3 as smallest available color
199 // That means, we only track candidates in the (half-open) interval [0, numColors)
200 // Mark candidate colors which are present in colors
201
2/2
✓ Branch 0 taken 60283027 times.
✓ Branch 1 taken 2175767 times.
 62458794 for (int i = 0; i < numColors; i++)
202
4/4
✓ Branch 0 taken 25696774 times.
✓ Branch 1 taken 34586253 times.
✓ Branch 2 taken 25682946 times.
✓ Branch 3 taken 13828 times.
 60283027 if (colors[i] >= 0 && colors[i] < numColors)
203 25682946 colorUsed[colors[i]] = true;
204
205 // return smallest color not in colors
206
2/2
✓ Branch 0 taken 8477073 times.
✓ Branch 1 taken 450 times.
 8477523 for (int i = 0; i < numColors; i++)
207
2/2
✓ Branch 0 taken 2175317 times.
✓ Branch 1 taken 6301756 times.
 8477073 if (!colorUsed[i])
208 2175317 return i;
209
210 return numColors;
211 }
212
213 } // end namespace Dumux::Detail
214 #endif // DOXYGEN
215
216 namespace Dumux {
217
218 /*!
219 * \brief Compute iterable lists of element seeds partitioned by color
220 *
221 * Splits up the elements of a grid view into partitions such that
222 * all elements in one partition do not modify global data structures
223 * at the same place during assembly. This is used to allow for
224 * lock-free thread-parallel (shared memory) assembly routines.
225 *
226 * Implements a simply greedy graph coloring algorithm:
227 * For each node (element), assign the smallest available color
228 * not used by any of the neighboring nodes (element with conflicting memory access)
229 * The greedy algorithm doesn't necessarily return the smallest
230 * possible number of colors (that's a hard problem) but is fast
231 *
232 * Returns a struct with access to the colors of each element (member colors)
233 * and vector of element seed sets of the same color (member sets)
234 *
235 * \param gg the grid geometry
236 * \param verbosity the verbosity level
237 */
238 template<class GridGeometry>
239 395 auto computeColoring(const GridGeometry& gg, int verbosity = 1)
240 {
241 395 Dune::Timer timer;
242
243 using ElementSeed = typename GridGeometry::GridView::Grid::template Codim<0>::EntitySeed;
244 393 struct Coloring
245 {
246 using Sets = std::deque<std::vector<ElementSeed>>;
247 using Colors = std::vector<int>;
248
249
1/2
✓ Branch 2 taken 395 times.
✗ Branch 3 not taken.
 395 Coloring(std::size_t size) : sets{}, colors(size, -1) {}
250
251 Sets sets;
252 Colors colors;
253 };
254
255 395 Coloring coloring(gg.gridView().size(0));
256
257 // pre-reserve some memory for helper arrays to avoid reallocation
258
1/2
✓ Branch 1 taken 395 times.
✗ Branch 2 not taken.
 395 std::vector<int> neighborColors; neighborColors.reserve(30);
259
1/2
✓ Branch 1 taken 395 times.
✗ Branch 2 not taken.
 395 std::vector<bool> colorUsed; colorUsed.reserve(30);
260
261 // dof to element map to speed up neighbor search
262
1/2
✓ Branch 1 taken 380 times.
✗ Branch 2 not taken.
 395 const auto connectedElements = Detail::computeConnectedElements(gg);
263
264
13/14
✓ Branch 1 taken 100823 times.
✓ Branch 2 taken 43 times.
✓ Branch 4 taken 489620 times.
✓ Branch 5 taken 40443 times.
✓ Branch 7 taken 293 times.
✓ Branch 8 taken 62427 times.
✓ Branch 9 taken 112554 times.
✓ Branch 10 taken 53 times.
✓ Branch 12 taken 112543 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 112543 times.
✓ Branch 15 taken 53 times.
✓ Branch 6 taken 531746 times.
✓ Branch 3 taken 62363 times.
 2032342 for (const auto& element : elements(gg.gridView()))
265 {
266 // compute neighbor colors based on discretization-dependent stencil
267 704638 neighborColors.clear();
268
1/2
✓ Branch 1 taken 704638 times.
✗ Branch 2 not taken.
 704638 Detail::addNeighborColors(gg, element, coloring.colors, connectedElements, neighborColors);
269
270 // find smallest color (positive integer) not in neighborColors
271
1/2
✓ Branch 1 taken 704638 times.
✗ Branch 2 not taken.
 704638 const auto color = Detail::smallestAvailableColor(neighborColors, colorUsed);
272
273 // assign color to element
274
5/5
✓ Branch 1 taken 699332 times.
✓ Branch 2 taken 1598 times.
✓ Branch 3 taken 174130 times.
✓ Branch 4 taken 840 times.
✓ Branch 0 taken 3708 times.
 704638 coloring.colors[gg.elementMapper().index(element)] = color;
275
276 // add element to the set of elements with the same color
277
2/2
✓ Branch 0 taken 702157 times.
✓ Branch 1 taken 2481 times.
 704638 if (color < coloring.sets.size())
278
4/5
✓ Branch 1 taken 1566 times.
✓ Branch 2 taken 660225 times.
✓ Branch 4 taken 112141 times.
✗ Branch 5 not taken.
✓ Branch 3 taken 40366 times.
 814298 coloring.sets[color].push_back(element.seed());
279 else
280
4/5
✓ Branch 1 taken 2045 times.
✓ Branch 2 taken 34 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2045 times.
✓ Branch 5 taken 34 times.
 4560 coloring.sets.push_back(std::vector<ElementSeed>{ element.seed() });
281 }
282
283
1/2
✓ Branch 0 taken 395 times.
✗ Branch 1 not taken.
 395 if (verbosity > 0)
284
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 400 std::cout << Fmt::format("Colored {} elements with {} colors in {} seconds.\n",
285
3/5
✓ Branch 2 taken 357 times.
✓ Branch 3 taken 38 times.
✓ Branch 5 taken 352 times.
✗ Branch 6 not taken.
✗ Branch 4 not taken.
 790 gg.gridView().size(0), coloring.sets.size(), timer.elapsed());
286
287 395 return coloring;
288 1201 }
289
290 //! Traits specifying if a given discretization tag supports coloring
291 template<class DiscretizationMethod>
292 struct SupportsColoring : public std::false_type {};
293
294 template<> struct SupportsColoring<DiscretizationMethods::CCTpfa> : public std::true_type {};
295 template<> struct SupportsColoring<DiscretizationMethods::CCMpfa> : public std::true_type {};
296 template<> struct SupportsColoring<DiscretizationMethods::Box> : public std::true_type {};
297 template<> struct SupportsColoring<DiscretizationMethods::FCStaggered> : public std::true_type {};
298 template<> struct SupportsColoring<DiscretizationMethods::FCDiamond> : public std::true_type {};
299 template<> struct SupportsColoring<DiscretizationMethods::PQ1Bubble> : public std::true_type {};
300
301 } // end namespace Dumux
302
303 #endif
304