GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/linear/parallelhelpers.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 221 256 86.3%
Functions: 225 4202 5.4%
Branches: 312 572 54.5%
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-FileCopyrightInfo: 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 Linear
10 * \brief Provides a helper class for nonoverlapping decomposition
11 */
12 #ifndef DUMUX_LINEAR_PARALLELHELPERS_HH
13 #define DUMUX_LINEAR_PARALLELHELPERS_HH
14
15 #include <dune/common/exceptions.hh>
16 #include <dune/geometry/dimension.hh>
17 #include <dune/grid/common/datahandleif.hh>
18 #include <dune/grid/common/partitionset.hh>
19 #include <dune/istl/owneroverlapcopy.hh>
20 #include <dune/istl/paamg/pinfo.hh>
21 #include <dune/istl/bvector.hh>
22 #include <dune/istl/multitypeblockvector.hh>
23 #include <dumux/parallel/vectorcommdatahandle.hh>
24 #include <dumux/common/gridcapabilities.hh>
25
26 namespace Dumux::Detail {
27
28 template<class LinearSolverTraits, bool canCommunicate = false>
29 class ParallelISTLHelperImpl {};
30
31 template<class LinearSolverTraits>
32 class ParallelISTLHelperImpl<LinearSolverTraits, true>
33 {
34 using GridView = typename LinearSolverTraits::GridView;
35 using DofMapper = typename LinearSolverTraits::DofMapper;
36 static constexpr int dofCodim = LinearSolverTraits::dofCodim;
37
38 // TODO: this is some large number (replace by limits?)
39 static constexpr std::size_t ghostMarker_ = 1<<24;
40
41 class BaseGatherScatter
42 {
43 public:
44 250 BaseGatherScatter(const DofMapper& mapper)
45 250 : mapper_(mapper) {}
46
47 template<class EntityType>
48 ✗ int index(const EntityType& e) const
49 62660 { return mapper_.index(e); }
50
51 ✗ bool contains(int dim, int codim) const
52 916 { return dofCodim == codim; }
53
54 //! returns true if size per entity of given dim and codim is a constant
55 ✗ bool fixedSize(int dim, int codim) const
56 ✗ { return true; }
57
58 template<class EntityType>
59 ✗ std::size_t size(EntityType& e) const
60 ✗ { return 1; }
61
62 template<class EntityType>
63 ✗ bool isNeitherInteriorNorBorderEntity(EntityType& e) const
64
19/72
✓ Branch 0 taken 610 times.
✓ Branch 1 taken 888 times.
✓ Branch 2 taken 786 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 838 times.
✓ Branch 5 taken 1596 times.
✓ Branch 6 taken 1016 times.
✓ Branch 7 taken 160 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 176 times.
✓ Branch 13 taken 772 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 176 times.
✓ Branch 16 taken 554 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 554 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 782 times.
✗ Branch 21 not taken.
✓ Branch 22 taken 678 times.
✓ Branch 23 taken 104 times.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✓ Branch 28 taken 104 times.
✓ Branch 29 taken 112 times.
✗ Branch 30 not taken.
✓ Branch 31 taken 104 times.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✗ Branch 48 not taken.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✗ Branch 51 not taken.
✗ Branch 52 not taken.
✗ Branch 53 not taken.
✗ Branch 54 not taken.
✗ Branch 55 not taken.
✗ Branch 56 not taken.
✗ Branch 57 not taken.
✗ Branch 58 not taken.
✗ Branch 59 not taken.
✗ Branch 60 not taken.
✗ Branch 61 not taken.
✗ Branch 62 not taken.
✗ Branch 63 not taken.
✗ Branch 64 not taken.
✗ Branch 65 not taken.
✗ Branch 66 not taken.
✓ Branch 67 taken 400 times.
✗ Branch 68 not taken.
✗ Branch 69 not taken.
✗ Branch 70 not taken.
✗ Branch 71 not taken.
 9602 { return e.partitionType() != Dune::InteriorEntity && e.partitionType() != Dune::BorderEntity; }
65
66 private:
67 const DofMapper& mapper_;
68 };
69
70 /*!
71 * \brief Writes ghostMarker_ to each data item (of the container) that is gathered or scattered
72 * and is neither interior nor border.
73 *
74 * Can be used to mark ghost cells.
75 */
76 class GhostGatherScatter
77 : public BaseGatherScatter
78 , public Dune::CommDataHandleIF<GhostGatherScatter, std::size_t>
79 {
80 public:
81 using DataType = std::size_t;
82 using BaseGatherScatter::contains;
83 using BaseGatherScatter::fixedSize;
84 using BaseGatherScatter::size;
85
86 50 GhostGatherScatter(std::vector<std::size_t>& ranks, const DofMapper& mapper)
87 : BaseGatherScatter(mapper)
88 100 , ranks_(ranks)
89 {}
90
91 template<class MessageBuffer, class EntityType>
92 11056 void gather(MessageBuffer& buff, const EntityType& e) const
93 {
94
2/2
✓ Branch 1 taken 3942 times.
✓ Branch 2 taken 440 times.
 11056 auto& data = ranks_[this->index(e)];
95
3/3
✓ Branch 0 taken 3942 times.
✓ Branch 1 taken 440 times.
✓ Branch 2 taken 302 times.
 11056 if (this->isNeitherInteriorNorBorderEntity(e))
96 ✗ data = ghostMarker_;
97 11056 buff.write(data);
98 11056 }
99
100 template<class MessageBuffer, class EntityType>
101 7224 void scatter(MessageBuffer& buff, const EntityType& e, std::size_t n)
102 {
103 std::size_t x;
104
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 7224 buff.read(x);
105
2/5
✓ Branch 1 taken 1128 times.
✓ Branch 2 taken 2484 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 7224 auto& data = ranks_[this->index(e)];
106
2/5
✓ Branch 0 taken 1128 times.
✓ Branch 1 taken 2484 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 7224 if (this->isNeitherInteriorNorBorderEntity(e))
107 4968 data = ghostMarker_;
108 7224 }
109 private:
110 std::vector<std::size_t>& ranks_;
111 };
112
113 /*!
114 * \brief GatherScatter handle that sets ghostMarker_ for data items neither associated to
115 * the interior or border and take the minimum when scattering.
116 *
117 * Used to compute an owner rank for each unknown.
118 */
119 class InteriorBorderGatherScatter
120 : public BaseGatherScatter
121 , public Dune::CommDataHandleIF<InteriorBorderGatherScatter, std::size_t>
122 {
123 public:
124 using DataType = std::size_t;
125 using BaseGatherScatter::contains;
126 using BaseGatherScatter::fixedSize;
127 using BaseGatherScatter::size;
128
129 50 InteriorBorderGatherScatter(std::vector<std::size_t>& ranks, const DofMapper& mapper)
130 : BaseGatherScatter(mapper)
131 100 , ranks_(ranks)
132 {}
133
134 /*!
135 * \brief Pack data (to be send to other processes) from this process to message buffer
136 */
137 template<class MessageBuffer, class EntityType>
138 2852 void gather(MessageBuffer& buff, const EntityType& e) const
139 {
140
1/2
✓ Branch 1 taken 1392 times.
✗ Branch 2 not taken.
 2852 auto& data = ranks_[this->index(e)];
141
2/3
✓ Branch 0 taken 1392 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 34 times.
 2852 if (this->isNeitherInteriorNorBorderEntity(e))
142 ✗ data = ghostMarker_;
143 2852 buff.write(data);
144 2852 }
145
146 /*!
147 * \brief Unpack data (received from other process) from message buffer to this process
148 */
149 template<class MessageBuffer, class EntityType>
150 2784 void scatter(MessageBuffer& buff, const EntityType& e, std::size_t n)
151 {
152 std::size_t x;
153
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 2784 buff.read(x);
154
1/4
✓ Branch 1 taken 1392 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 2784 auto& data = ranks_[this->index(e)];
155
156 // we leave ghost unchanged
157 // for other dofs, the process with the lowest rank
158 // is assigned to be the (unique) owner
159 using std::min;
160
3/7
✓ Branch 0 taken 1392 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 696 times.
✓ Branch 3 taken 696 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 6960 data = this->isNeitherInteriorNorBorderEntity(e) ? x : min(data, x);
161 2784 }
162 private:
163 std::vector<std::size_t>& ranks_;
164 };
165
166 /*!
167 * \brief GatherScatter handle for finding out about neighbouring processor ranks.
168 *
169 */
170 struct NeighbourGatherScatter
171 : public BaseGatherScatter
172 , public Dune::CommDataHandleIF<NeighbourGatherScatter, int>
173 {
174 using DataType = int;
175 using BaseGatherScatter::contains;
176 using BaseGatherScatter::fixedSize;
177 using BaseGatherScatter::size;
178
179 50 NeighbourGatherScatter(const DofMapper& mapper, int rank, std::set<int>& neighbours)
180 : BaseGatherScatter(mapper)
181 , rank_(rank)
182 100 , neighbours_(neighbours)
183 {}
184
185 template<class MessageBuffer, class EntityType>
186 ✗ void gather(MessageBuffer& buff, const EntityType& e) const
187 {
188
1/8
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 800 times.
 9888 buff.write(rank_);
189 ✗ }
190
191 template<class MessageBuffer, class EntityType>
192 ✗ void scatter(MessageBuffer& buff, const EntityType& e, std::size_t n)
193 {
194 int x;
195
0/13
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
 202 buff.read(x);
196
0/13
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
 202 neighbours_.insert(x);
197 ✗ }
198 private:
199 int rank_;
200 std::set<int>& neighbours_;
201 };
202
203
204 /*!
205 * \brief GatherScatter handle for finding out about neighbouring processor ranks.
206 *
207 */
208 struct SharedGatherScatter
209 : public BaseGatherScatter
210 , public Dune::CommDataHandleIF<SharedGatherScatter, int>
211 {
212 using DataType = int;
213 using BaseGatherScatter::contains;
214 using BaseGatherScatter::fixedSize;
215 using BaseGatherScatter::size;
216
217 50 SharedGatherScatter(std::vector<int>& shared, const DofMapper& mapper)
218 : BaseGatherScatter(mapper)
219 100 , shared_(shared)
220 {}
221
222 template<class MessageBuffer, class EntityType>
223 ✗ void gather(MessageBuffer& buff, EntityType& e) const
224 {
225 9686 int data = true;
226
1/8
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 800 times.
 9888 buff.write(data);
227 ✗ }
228
229 template<class MessageBuffer, class EntityType>
230 ✗ void scatter(MessageBuffer& buff, const EntityType &e, std::size_t n)
231 {
232 int x;
233 ✗ buff.read(x);
234 ✗ auto& data = shared_[this->index(e)];
235 ✗ data = data || x;
236 ✗ }
237 private:
238 std::vector<int>& shared_;
239 };
240
241 /*!
242 * \brief GatherScatter handle for finding out about neighbouring processor ranks.
243 *
244 */
245 template<typename GlobalIndex>
246 struct GlobalIndexGatherScatter
247 : public BaseGatherScatter
248 , public Dune::CommDataHandleIF<GlobalIndexGatherScatter<GlobalIndex>, GlobalIndex>
249 {
250 using DataType = GlobalIndex;
251 using BaseGatherScatter::contains;
252 using BaseGatherScatter::fixedSize;
253 using BaseGatherScatter::size;
254
255 50 GlobalIndexGatherScatter(std::vector<GlobalIndex>& globalIndices, const DofMapper& mapper)
256 : BaseGatherScatter(mapper)
257 100 , globalIndices_(globalIndices)
258 {}
259
260 template<class MessageBuffer, class EntityType>
261 19372 void gather(MessageBuffer& buff, const EntityType& e) const
262 {
263 19372 buff.write(globalIndices_[this->index(e)]);
264 19372 }
265
266 template<class MessageBuffer, class EntityType>
267 19372 void scatter(MessageBuffer& buff, const EntityType& e, std::size_t n)
268 {
269 19372 DataType x;
270 19372 buff.read(x);
271 using std::min;
272 19372 DataType& data = globalIndices_[this->index(e)];
273 19372 data = min(data, x);
274 19372 }
275 private:
276 std::vector<GlobalIndex>& globalIndices_;
277 };
278
279 public:
280
281 50 ParallelISTLHelperImpl(const GridView& gridView, const DofMapper& mapper)
282
5/12
✓ Branch 1 taken 48 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 48 times.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 48 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 146 : gridView_(gridView), mapper_(mapper)
283 {
284 if constexpr (Detail::canCommunicate<typename GridView::Traits::Grid, dofCodim>)
285
1/2
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
 50 initGhostsAndOwners_();
286 else
287 DUNE_THROW(Dune::InvalidStateException,
288 "Cannot initialize parallel helper for a grid that cannot communicate codim-" << dofCodim << "-entities."
289 );
290 50 }
291
292 bool isGhost(std::size_t i) const
293 554974 { return isGhost_[i] == ghostMarker_; }
294
295 bool isOwned(std::size_t i) const
296 { return isOwned_[i] == 1; }
297
298 /*!
299 * \brief Creates a parallel index set for a communicator
300 * \param comm The OwnerOverlapCopyCommunication communicators
301 *
302 * The parallel index set contains for each dof index,
303 * the triplet (globalIndex, (localIndex, attribute))
304 * where attribute is owner, overlap or copy. Each dof
305 * is uniquely owned by exactly one process. This allows
306 * to apply local operators additively to get a global operator
307 * without communication.
308 */
309 template<class Comm>
310 50 void createParallelIndexSet(Comm& comm) const
311 {
312 if constexpr (Detail::canCommunicate<typename GridView::Traits::Grid, dofCodim>)
313 {
314
5/6
✗ Branch 0 not taken.
✓ Branch 1 taken 30 times.
✓ Branch 2 taken 20 times.
✓ Branch 3 taken 30 times.
✓ Branch 4 taken 20 times.
✓ Branch 5 taken 30 times.
 54 if (gridView_.comm().size() <= 1)
315 {
316 ✗ comm.remoteIndices().template rebuild<false>();
317 ✗ return;
318 }
319
320 // First find out which dofs we share with other processors
321
5/10
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 22 times.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 4 times.
✗ Branch 11 not taken.
 174 std::vector<int> isShared(mapper_.size(), false);
322 100 SharedGatherScatter sgs(isShared, mapper_);
323
2/4
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
 50 gridView_.communicate(sgs, Dune::All_All_Interface, Dune::ForwardCommunication);
324
325 // Count shared dofs that we own
326 using GlobalIndex = typename Comm::ParallelIndexSet::GlobalIndex;
327
1/2
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
 50 GlobalIndex count = 0;
328
4/4
✓ Branch 0 taken 49665 times.
✓ Branch 1 taken 50 times.
✓ Branch 2 taken 49665 times.
✓ Branch 3 taken 50 times.
 49765 for (std::size_t i = 0; i < isShared.size(); ++i)
329
8/8
✓ Branch 0 taken 8229 times.
✓ Branch 1 taken 41436 times.
✓ Branch 2 taken 8229 times.
✓ Branch 3 taken 41436 times.
✓ Branch 4 taken 3422 times.
✓ Branch 5 taken 4807 times.
✓ Branch 6 taken 3422 times.
✓ Branch 7 taken 4807 times.
 99330 if (isShared[i] && isOwned_[i] == 1)
330 ++count;
331
332
9/16
✗ Branch 0 not taken.
✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 48 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 28 times.
✓ Branch 6 taken 20 times.
✓ Branch 7 taken 2 times.
✓ Branch 8 taken 48 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 30 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
 154 std::vector<GlobalIndex> counts(gridView_.comm().size());
333
7/9
✗ Branch 0 not taken.
✓ Branch 1 taken 30 times.
✓ Branch 2 taken 20 times.
✓ Branch 3 taken 28 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 20 times.
✓ Branch 6 taken 28 times.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 54 gridView_.comm().allgather(&count, 1, counts.data());
334
335 // Compute start index start_p = \sum_{i=0}^{i<p} counts_i
336
4/6
✗ Branch 0 not taken.
✓ Branch 1 taken 30 times.
✓ Branch 2 taken 20 times.
✓ Branch 3 taken 28 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 52 const int rank = gridView_.comm().rank();
337
1/2
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
 50 auto start = std::accumulate(counts.begin(), counts.begin() + rank, GlobalIndex(0));
338
339 // starting from start (offset), number global indices owned by this process consecutively
340
4/11
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 22 times.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 28 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 224 std::vector<GlobalIndex> globalIndices(mapper_.size(), std::numeric_limits<GlobalIndex>::max());
341
4/4
✓ Branch 0 taken 49665 times.
✓ Branch 1 taken 50 times.
✓ Branch 2 taken 49665 times.
✓ Branch 3 taken 50 times.
 49715 for (std::size_t i = 0; i < globalIndices.size(); ++i)
342 {
343
8/8
✓ Branch 0 taken 44858 times.
✓ Branch 1 taken 4807 times.
✓ Branch 2 taken 44858 times.
✓ Branch 3 taken 4807 times.
✓ Branch 4 taken 3422 times.
✓ Branch 5 taken 41436 times.
✓ Branch 6 taken 3422 times.
✓ Branch 7 taken 41436 times.
 99330 if (isOwned_[i] == 1 && isShared[i])
344 {
345 3422 globalIndices[i] = start;
346 3422 ++start;
347 }
348 }
349
350 // publish global indices for the shared DOFS to other processors.
351 100 GlobalIndexGatherScatter<GlobalIndex> gigs(globalIndices, mapper_);
352
2/4
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
 50 gridView_.communicate(gigs, Dune::All_All_Interface, Dune::ForwardCommunication);
353
354 // put the information into the parallel index set
355
1/2
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
 50 resizeAndFillIndexSet_(comm, globalIndices);
356
357 // Compute neighbours using communication
358
4/6
✗ Branch 0 not taken.
✓ Branch 1 taken 30 times.
✓ Branch 2 taken 20 times.
✓ Branch 3 taken 28 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 100 std::set<int> neighbours;
359
4/6
✗ Branch 0 not taken.
✓ Branch 1 taken 30 times.
✓ Branch 2 taken 20 times.
✓ Branch 3 taken 28 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 52 NeighbourGatherScatter ngs(mapper_, gridView_.comm().rank(), neighbours);
360
2/4
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
 50 gridView_.communicate(ngs, Dune::All_All_Interface, Dune::ForwardCommunication);
361
2/4
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 50 times.
✗ Branch 5 not taken.
 100 comm.remoteIndices().setNeighbours(neighbours);
362
363
1/2
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
 50 comm.remoteIndices().template rebuild<false>();
364 }
365 else
366 DUNE_THROW(Dune::InvalidStateException,
367 "Cannot build parallel index set for a grid that cannot communicate codim-" << dofCodim << "-entities."
368 );
369 }
370
371 //! Return the dofMapper
372 ✗ const DofMapper& dofMapper() const
373 ✗ { return mapper_; }
374
375 //! Return the gridView
376 const GridView& gridView() const
377
3/5
✗ Branch 0 not taken.
✓ Branch 1 taken 218 times.
✓ Branch 2 taken 28 times.
✓ Branch 3 taken 22 times.
✗ Branch 4 not taken.
 844 { return gridView_; }
378
379 private:
380 /*!
381 * \brief Initialize ghost and owner flags
382 *
383 * All dofs on ghost entities are marked.
384 * All dofs are assigned a unique owner process. This will help us to create a unique partition
385 * and unique representations of a vector according to Definition 2.5 in Blatt and Bastian (2009)
386 * https://doi.org/10.1504/IJCSE.2008.021112
387 */
388 50 void initGhostsAndOwners_()
389 {
390
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 28 times.
 52 const auto rank = gridView_.comm().rank();
391 74 isOwned_.assign(mapper_.size(), rank);
392
393 // find out about ghosts
394 100 GhostGatherScatter ggs(isOwned_, mapper_);
395
396
5/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 48 times.
✓ Branch 2 taken 30 times.
✓ Branch 3 taken 20 times.
✓ Branch 4 taken 30 times.
✗ Branch 5 not taken.
 54 if (gridView_.comm().size() > 1)
397 50 gridView_.communicate(ggs, Dune::InteriorBorder_All_Interface, Dune::ForwardCommunication);
398
399 // isGhost_ contains rank when the dof is not on a ghost entity (interior or border)
400 // and ghostMarker_ when the dof is on a ghost
401 50 isGhost_ = isOwned_;
402
403 // partition interior/border uniquely
404 // after the communication each dof is assigned to one unique owner process rank
405 100 InteriorBorderGatherScatter dh(isOwned_, mapper_);
406
407
5/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 48 times.
✓ Branch 2 taken 30 times.
✓ Branch 3 taken 20 times.
✓ Branch 4 taken 30 times.
✗ Branch 5 not taken.
 54 if (gridView_.comm().size() > 1)
408 50 gridView_.communicate(dh, Dune::InteriorBorder_InteriorBorder_Interface, Dune::ForwardCommunication);
409
410 // convert vector into mask vector
411
4/4
✓ Branch 0 taken 49665 times.
✓ Branch 1 taken 50 times.
✓ Branch 2 taken 49665 times.
✓ Branch 3 taken 50 times.
 99480 for (auto& v : isOwned_)
412
2/2
✓ Branch 0 taken 4807 times.
✓ Branch 1 taken 44858 times.
 54472 v = (v == rank) ? 1 : 0;
413 50 }
414
415 template<class Comm, class GlobalIndices>
416 50 void resizeAndFillIndexSet_(Comm& comm, const GlobalIndices& globalIndices) const
417 {
418 50 comm.indexSet().beginResize();
419
420 // add triplets characterizing each dof in the parallel index set
421 // (globalIndex, (localIndex, attribute)) where attribute is owner, overlap or copy
422
4/4
✓ Branch 1 taken 49665 times.
✓ Branch 2 taken 50 times.
✓ Branch 3 taken 49665 times.
✓ Branch 4 taken 50 times.
 49715 for (std::size_t localIndex = 0; localIndex < globalIndices.size(); ++localIndex)
423 {
424 49665 const auto globalIndex = globalIndices[localIndex];
425
2/2
✓ Branch 0 taken 8229 times.
✓ Branch 1 taken 41436 times.
 148995 if (globalIndex != std::numeric_limits<typename GlobalIndices::value_type>::max())
426 {
427
2/2
✓ Branch 0 taken 4807 times.
✓ Branch 1 taken 3422 times.
 8229 const bool isOwned = isOwned_[localIndex] > 0;
428
2/2
✓ Branch 0 taken 4807 times.
✓ Branch 1 taken 3422 times.
 8229 const auto attr = getAttribute_(comm, isOwned, isGhost(localIndex));
429 using LocalIndex = typename Comm::ParallelIndexSet::LocalIndex;
430 16458 comm.indexSet().add(globalIndex, LocalIndex{localIndex, attr});
431 }
432 }
433
434 50 comm.indexSet().endResize();
435 50 }
436
437 template<class Comm>
438 Dune::OwnerOverlapCopyAttributeSet::AttributeSet
439 getAttribute_(const Comm& comm, const bool isOwned, const bool isGhost) const
440 {
441
2/2
✓ Branch 0 taken 4807 times.
✓ Branch 1 taken 3422 times.
 8229 if (isOwned) // this process owns the dof (uniquely)
442 return Dune::OwnerOverlapCopyAttributeSet::owner;
443
4/4
✓ Branch 0 taken 4094 times.
✓ Branch 1 taken 713 times.
✓ Branch 2 taken 1576 times.
✓ Branch 3 taken 2518 times.
 4807 else if (isGhost && (comm.category() == Dune::SolverCategory::nonoverlapping) )
444 return Dune::OwnerOverlapCopyAttributeSet::overlap;
445 else
446 return Dune::OwnerOverlapCopyAttributeSet::copy;
447 }
448
449 const GridView gridView_; //!< the grid view
450 const DofMapper& mapper_; //!< the dof mapper
451 //! vector to identify unique decomposition (1: this rank owns the dof; 0: owned by remote process)
452 std::vector<std::size_t> isOwned_;
453 //!< vector to identify ghost dofs (ghostMarker_: this dof is on a ghost entity; contains and unspecified value otherwise)
454 std::vector<std::size_t> isGhost_;
455
456 };
457
458 } // end namespace Dumux::Detail
459
460 namespace Dumux {
461
462 /*!
463 * \ingroup Linear
464 * \brief A parallel helper class providing a parallel
465 * decomposition of all degrees of freedom
466 */
467 template<class LinearSolverTraits>
468 using ParallelISTLHelper =
469 Detail::ParallelISTLHelperImpl<
470 LinearSolverTraits, LinearSolverTraits::canCommunicate
471 >;
472
473
474 template<class GridView, class DofMapper, int dofCodim>
475 ✗ class ParallelVectorHelper
476 {
477 public:
478 448 ParallelVectorHelper(const GridView& gridView, const DofMapper& mapper)
479
2/4
✓ Branch 1 taken 162 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 162 times.
✗ Branch 5 not taken.
 324 : gridView_(gridView), mapper_(mapper)
480 {}
481
482 //! \brief Make a vector consistent for non-overlapping domain decomposition methods
483 template<class Block, class Alloc>
484 1200 void makeNonOverlappingConsistent(Dune::BlockVector<Block, Alloc>& v) const
485 {
486 if constexpr (Detail::canCommunicate<typename GridView::Traits::Grid, dofCodim>)
487 {
488
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 308 times.
 1200 VectorCommDataHandleSum<DofMapper, Dune::BlockVector<Block, Alloc>, dofCodim, Block> gs(mapper_, v);
489
4/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1200 times.
✓ Branch 2 taken 308 times.
✓ Branch 3 taken 892 times.
✓ Branch 4 taken 308 times.
✗ Branch 5 not taken.
 1200 if (gridView_.comm().size() > 1)
490 1200 gridView_.communicate(gs, Dune::InteriorBorder_InteriorBorder_Interface,
491 Dune::ForwardCommunication);
492 }
493 else
494 DUNE_THROW(Dune::InvalidStateException, "Cannot call makeNonOverlappingConsistent for a grid that cannot communicate codim-" << dofCodim << "-entities.");
495 1200 }
496
497 //! \brief Make a vector consistent for non-overlapping domain decomposition methods
498 template<class... Blocks>
499 void makeNonOverlappingConsistent(Dune::MultiTypeBlockVector<Blocks...>& v) const
500 {
501 DUNE_THROW(Dune::NotImplemented, "makeNonOverlappingConsistent for Dune::MultiTypeBlockVector");
502 }
503
504 private:
505 const GridView gridView_; //!< the grid view
506 const DofMapper& mapper_; //!< the dof mapper
507 };
508
509 /*!
510 * \ingroup Linear
511 * \brief Helper class for adding up matrix entries for border entities
512 *
513 * Border means all degrees of freedom located
514 * on lower-dimensional entities (faces, edges, vertices)
515 * that form the the process boundary
516 */
517 template<class Matrix, class GridView,
518 class RowDofMapper, int rowDofCodim>
519 class ParallelMatrixHelper
520 {
521 using IdType = typename GridView::Traits::Grid::Traits::GlobalIdSet::IdType;
522
523 /*!
524 * \brief A DataHandle class to exchange matrix sparsity patterns.
525 *
526 * We look at a 2D example with a nonoverlapping grid,
527 * two processes and no ghosts with Q1 discretization.
528 * Process 0 has the left part of the domain
529 * with three cells and eight vertices (1-8),
530 * Process 1 the right part with three cells
531 * and eight vertices (2,4,7-12).
532 * <pre>
533 * 1 _ 2 2 _ 9 _ 10
534 * | | | | |
535 * 3 _ 4 _ 7 4 _ 7 _ 11
536 * | | | | |
537 * 5 _ 6 _ 8 8 _ 12
538 * </pre>
539 * If we look at vertex 7 and the corresponding entries in the matrix for P0,
540 * there will be entries for (7,4) and (7,8), but not for (7,2).
541 * The MatrixPatternExchange class will find these entries and fills a vector "sparsity",
542 * that contains all missing connections.
543 */
544 template<class ColIsGhostFunc>
545 struct MatrixPatternExchange
546 : public Dune::CommDataHandleIF<MatrixPatternExchange<ColIsGhostFunc>, IdType>
547 {
548 //! Export type of data for message buffer
549 using DataType = IdType;
550
551 /*!
552 * \brief Construct a new sparsity pattern exchange handle
553 *
554 * \param rowEntityMapper an index mapper for entities associated with row entities
555 * \param globalToLocal map (column dof) from global ID to processor-local index
556 * \param localToGlobal map (column dof) from processor-local index to global ID
557 * \param A the matrix for which we want to extend the pattern
558 * \param[out] sparsityPattern the extended sparsity pattern
559 * \param isGhostColumDof a function that returns whether a given column index is associated with a ghost entity
560 */
561 940 MatrixPatternExchange(const RowDofMapper& rowEntityMapper,
562 const std::map<IdType,int>& globalToLocal,
563 const std::map<int,IdType>& localToGlobal,
564 Matrix& A,
565 std::vector<std::set<int>>& sparsityPattern,
566 const ColIsGhostFunc& isGhostColumDof)
567 : rowEntityMapper_(rowEntityMapper), idToIndex_(globalToLocal), indexToID_(localToGlobal)
568 940 , sparsityPattern_(sparsityPattern), A_(A), isGhostColumDof_(isGhostColumDof)
569 {
570 940 sparsityPattern_.resize(A.N());
571 }
572
573 /*!
574 * \brief Returns true if data for given valid codim should be communicated
575 */
576 ✗ bool contains (int dim, int codim) const
577 48 { return (codim == rowDofCodim); }
578
579 //! returns true if size per entity of given dim and codim is a constant
580 ✗ bool fixedSize(int dim, int codim) const
581 ✗ { return false; }
582
583 /*!
584 * \brief How many objects of type DataType have to be sent for a given entity
585 */
586 template<class EntityType>
587 186600 std::size_t size(EntityType& e) const
588 {
589 186600 const auto rowIdx = rowEntityMapper_.index(e);
590
591 // all column entity indices of this row that are in the index set
592 186600 std::size_t n = 0;
593
6/6
✓ Branch 0 taken 908884 times.
✓ Branch 1 taken 93300 times.
✓ Branch 2 taken 908884 times.
✓ Branch 3 taken 93300 times.
✓ Branch 4 taken 908884 times.
✓ Branch 5 taken 93300 times.
 2190968 for (auto colIt = A_[rowIdx].begin(); colIt != A_[rowIdx].end(); ++colIt)
594 3635536 if (indexToID_.count(colIt.index()))
595 565392 n++;
596
597 186600 return n;
598 }
599
600 /*!
601 * \brief Pack data from user to message buffer
602 */
603 template<class MessageBuffer, class EntityType>
604 143760 void gather(MessageBuffer& buff, const EntityType& e) const
605 {
606 143760 const auto rowIdx = rowEntityMapper_.index(e);
607
608 // send all column entity ids of this row that are in the index set
609
8/8
✓ Branch 0 taken 8896 times.
✓ Branch 1 taken 795468 times.
✓ Branch 2 taken 8896 times.
✓ Branch 3 taken 795468 times.
✓ Branch 4 taken 8896 times.
✓ Branch 5 taken 795468 times.
✓ Branch 6 taken 723588 times.
✓ Branch 7 taken 70336 times.
 1752488 for (auto colIt = A_[rowIdx].begin(); colIt != A_[rowIdx].end(); ++colIt)
610
4/4
✓ Branch 0 taken 269466 times.
✓ Branch 1 taken 463018 times.
✓ Branch 2 taken 269466 times.
✓ Branch 3 taken 463018 times.
 5859872 if (auto it = indexToID_.find(colIt.index()); it != indexToID_.end())
611 1077864 buff.write(it->second);
612 143760 }
613
614 /*!
615 * \brief Unpack data from message buffer to user
616 */
617 template<class MessageBuffer, class EntityType>
618 143760 void scatter(MessageBuffer& buff, const EntityType& e, std::size_t n)
619 {
620 143760 const auto rowIdx = rowEntityMapper_.index(e);
621
2/2
✓ Branch 1 taken 269466 times.
✓ Branch 2 taken 71880 times.
 682692 for (std::size_t k = 0; k < n; k++)
622 {
623 // receive n column entity IDs
624 530112 IdType id;
625 538932 buff.read(id);
626
627 // only add entries that are contained in the index set
628
5/5
✓ Branch 0 taken 4410 times.
✓ Branch 1 taken 264848 times.
✓ Branch 2 taken 4618 times.
✓ Branch 3 taken 264848 times.
✓ Branch 4 taken 208 times.
 556572 if (const auto it = idToIndex_.find(id); it != idToIndex_.end())
629 {
630
1/2
✓ Branch 0 taken 269258 times.
✗ Branch 1 not taken.
 538516 const auto colIdx = it->second; // get local column index
631 // add this entry (if it doesn't exist yet)
632 // (and only if the column dof is not associated with a ghost entity)
633 1077032 if (!isGhostColumDof_(colIdx))
634 // std::set takes care that each index only is inserted once
635 1077032 sparsityPattern_[rowIdx].insert(colIdx);
636 }
637 }
638 143760 }
639
640 private:
641 const RowDofMapper& rowEntityMapper_;
642 const std::map<IdType, int>& idToIndex_;
643 const std::map<int, IdType>& indexToID_;
644 std::vector<std::set<int>>& sparsityPattern_;
645 Matrix& A_;
646 const ColIsGhostFunc& isGhostColumDof_;
647
648 }; // class MatrixPatternExchange
649
650 //! Local matrix blocks associated with the global id set
651 1345874 struct MatrixEntry
652 {
653 IdType first;
654 typename Matrix::block_type second;
655 };
656
657 //! A DataHandle class to exchange matrix entries
658 struct MatrixEntryExchange
659 : public Dune::CommDataHandleIF<MatrixEntryExchange, MatrixEntry>
660 {
661 //! Export type of data for message buffer
662 using DataType = MatrixEntry;
663
664 940 MatrixEntryExchange(const RowDofMapper& rowEntityMapper,
665 const std::map<IdType, int>& globalToLocal,
666 const std::map<int, IdType>& localToGlobal,
667 Matrix& A)
668 940 : rowEntityMapper_(rowEntityMapper), idToIndex_(globalToLocal), indexToID_(localToGlobal), A_(A)
669 {}
670
671 /*!
672 * \brief Returns true if data for given valid codim should be communicated
673 */
674 ✗ bool contains(int dim, int codim) const
675 48 { return (codim == rowDofCodim); }
676
677 //! returns true if size per entity of given dim and codim is a constant
678 ✗ bool fixedSize(int dim, int codim) const
679 ✗ { return false; }
680
681 /*!
682 * \brief How many objects of type DataType have to be sent for a given entity
683 */
684 template<class EntityType>
685 186600 std::size_t size(EntityType& e) const
686 {
687 186600 const auto rowIdx = rowEntityMapper_.index(e);
688 186600 std::size_t n = 0;
689
6/6
✓ Branch 0 taken 909386 times.
✓ Branch 1 taken 93300 times.
✓ Branch 2 taken 909386 times.
✓ Branch 3 taken 93300 times.
✓ Branch 4 taken 909386 times.
✓ Branch 5 taken 93300 times.
 2191972 for (auto colIt = A_[rowIdx].begin(); colIt != A_[rowIdx].end(); ++colIt)
690 3637544 if (indexToID_.count(colIt.index()))
691 566396 n++;
692
693 186600 return n;
694 }
695
696 /*!
697 * \brief Pack data from user to message buffer
698 */
699 template<class MessageBuffer, class EntityType>
700 143760 void gather(MessageBuffer& buff, const EntityType& e) const
701 {
702 143760 const auto rowIdx = rowEntityMapper_.index(e);
703 // send all matrix entries for which the column index is in the index set
704
8/8
✓ Branch 0 taken 8896 times.
✓ Branch 1 taken 795970 times.
✓ Branch 2 taken 8896 times.
✓ Branch 3 taken 795970 times.
✓ Branch 4 taken 8896 times.
✓ Branch 5 taken 795970 times.
✓ Branch 6 taken 724090 times.
✓ Branch 7 taken 70336 times.
 1753492 for (auto colIt = A_[rowIdx].begin(); colIt != A_[rowIdx].end(); ++colIt)
705
4/4
✓ Branch 0 taken 269968 times.
✓ Branch 1 taken 463018 times.
✓ Branch 2 taken 269968 times.
✓ Branch 3 taken 463018 times.
 5863888 if (auto it = indexToID_.find(colIt.index()); it != indexToID_.end())
706 2159744 buff.write(MatrixEntry{ it->second, *colIt });
707 143760 }
708
709 /*!
710 * \brief Unpack data from message buffer to user
711 */
712 template<class MessageBuffer, class EntityType>
713 143760 void scatter(MessageBuffer& buff, const EntityType& e, std::size_t n)
714 {
715 143760 const auto rowIdx = rowEntityMapper_.index(e);
716
2/2
✓ Branch 1 taken 269968 times.
✓ Branch 2 taken 71880 times.
 683696 for (std::size_t k = 0; k < n; k++)
717 {
718 539936 MatrixEntry m;
719 539936 buff.read(m);
720 539936 const auto& [colDofID, matrixBlock] = m; // unpack
721
722 // only add entries in the index set and for which A has allocated memory
723
5/5
✓ Branch 0 taken 4410 times.
✓ Branch 1 taken 265350 times.
✓ Branch 2 taken 4618 times.
✓ Branch 3 taken 265350 times.
✓ Branch 4 taken 208 times.
 557576 if (auto it = idToIndex_.find(colDofID); it != idToIndex_.end())
724
2/4
✗ Branch 4 not taken.
✓ Branch 5 taken 269760 times.
✓ Branch 6 taken 269760 times.
✗ Branch 7 not taken.
 2158080 if (A_[rowIdx].find(it->second) != A_[rowIdx].end())
725 1649416 A_[rowIdx][it->second] += matrixBlock;
726 }
727 143760 }
728
729 private:
730 const RowDofMapper& rowEntityMapper_;
731 const std::map<IdType, int>& idToIndex_;
732 const std::map<int, IdType>& indexToID_;
733 Matrix& A_;
734
735 }; // class MatrixEntryExchange
736
737 public:
738
739 940 ParallelMatrixHelper(const GridView& gridView, const RowDofMapper& mapper)
740
0/2
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 2820 : gridView_(gridView), mapper_(mapper)
741 {
742 940 idToIndex_.clear();
743 940 indexToID_.clear();
744
745
6/10
✓ Branch 1 taken 850 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 90 times.
✓ Branch 4 taken 850 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 90 times.
✓ Branch 7 taken 850 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 850 times.
✗ Branch 11 not taken.
 1970 std::vector<bool> handledDof(gridView_.size(rowDofCodim), false);
746
13/16
✓ Branch 1 taken 940 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 36128 times.
✓ Branch 4 taken 848 times.
✓ Branch 5 taken 18720 times.
✓ Branch 6 taken 90 times.
✓ Branch 7 taken 758 times.
✓ Branch 8 taken 18720 times.
✓ Branch 9 taken 1290386 times.
✓ Branch 10 taken 758 times.
✓ Branch 11 taken 1290386 times.
✓ Branch 12 taken 758 times.
✓ Branch 14 taken 1290386 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 1290386 times.
✗ Branch 18 not taken.
 1328212 for (const auto& element : elements(gridView_))
747 {
748
5/5
✓ Branch 0 taken 69264 times.
✓ Branch 1 taken 6555570 times.
✓ Branch 2 taken 69264 times.
✓ Branch 3 taken 5246464 times.
✓ Branch 4 taken 1309106 times.
 6642150 for (int i = 0; i < element.subEntities(rowDofCodim); ++i)
749 {
750
1/6
✓ Branch 1 taken 5229148 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 10452680 const auto entity = element.template subEntity<rowDofCodim>(i);
751
5/5
✓ Branch 0 taken 5024050 times.
✓ Branch 1 taken 274362 times.
✓ Branch 2 taken 6390 times.
✓ Branch 3 taken 68634 times.
✓ Branch 4 taken 69120 times.
 5367676 if (entity.partitionType() == Dune::BorderEntity)
752 {
753
1/2
✓ Branch 1 taken 5760 times.
✗ Branch 2 not taken.
 284302 const auto localRowIdx = mapper_.index(entity);
754
3/6
✓ Branch 0 taken 284302 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 284302 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 284302 times.
✗ Branch 5 not taken.
 852906 if (!handledDof[localRowIdx])
755 {
756
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 284302 times.
✓ Branch 3 taken 284302 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 284302 times.
✗ Branch 7 not taken.
 284302 IdType dofIdxGlobal = gridView_.grid().globalIdSet().id(entity);
757
1/2
✓ Branch 1 taken 284302 times.
✗ Branch 2 not taken.
 284302 idToIndex_.emplace(dofIdxGlobal, localRowIdx);
758
1/2
✓ Branch 1 taken 284302 times.
✗ Branch 2 not taken.
 284302 indexToID_.emplace(localRowIdx, dofIdxGlobal);
759 }
760 }
761 }
762 }
763 940 }
764
765 /*!
766 * \brief communicates values for the sparsity pattern of the new matrix.
767 * \param A Matrix to operate on.
768 * \param isGhost Function returning if a column dof index is on a ghost entity
769 */
770 template<class IsGhostFunc>
771 940 void extendMatrix(Matrix& A, const IsGhostFunc& isGhost)
772 {
773 if constexpr (Detail::canCommunicate<typename GridView::Grid, rowDofCodim>)
774 {
775
5/6
✗ Branch 0 not taken.
✓ Branch 1 taken 182 times.
✓ Branch 2 taken 758 times.
✓ Branch 3 taken 182 times.
✓ Branch 4 taken 758 times.
✓ Branch 5 taken 182 times.
 940 if (gridView_.comm().size() <= 1)
776 ✗ return;
777
778 // make a copy of the matrix as originally assembled
779 1880 Matrix matrixAsAssembled(A);
780
781 // first get entries in sparsity pattern from other processes
782 940 std::size_t numNonZeroEntries = 0;
783
1/2
✓ Branch 1 taken 940 times.
✗ Branch 2 not taken.
 940 std::vector<std::set<int>> sparsityPattern; // column indices for every row
784
1/2
✓ Branch 1 taken 940 times.
✗ Branch 2 not taken.
 940 MatrixPatternExchange<IsGhostFunc> dataHandle(mapper_, idToIndex_, indexToID_, A, sparsityPattern, isGhost);
785
1/2
✓ Branch 1 taken 940 times.
✗ Branch 2 not taken.
 940 gridView_.communicate(
786 dataHandle, Dune::InteriorBorder_InteriorBorder_Interface, Dune::ForwardCommunication
787 );
788
789 // add own entries to the sparsity pattern and count number of non-zeros
790
4/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1485606 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1485606 times.
✓ Branch 4 taken 1484666 times.
✓ Branch 5 taken 940 times.
 2970272 for (auto rowIt = A.begin(); rowIt != A.end(); ++rowIt)
791 {
792 2969332 const auto colEndIt = A[rowIt.index()].end();
793
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 13905886 times.
✓ Branch 2 taken 12421220 times.
✓ Branch 3 taken 1484666 times.
 15390552 for (auto colIt = rowIt->begin(); colIt != colEndIt; ++colIt)
794 37263660 if (!sparsityPattern[rowIt.index()].count(colIt.index()))
795
1/2
✓ Branch 1 taken 12152464 times.
✗ Branch 2 not taken.
 12152464 sparsityPattern[rowIt.index()].insert(colIt.index());
796
797 4453998 numNonZeroEntries += sparsityPattern[rowIt.index()].size();
798 }
799
800 // insert any additional entries into the matrix
801
1/2
✓ Branch 1 taken 940 times.
✗ Branch 2 not taken.
 940 A.setSize(matrixAsAssembled.N(), matrixAsAssembled.M(), numNonZeroEntries);
802
1/2
✓ Branch 1 taken 940 times.
✗ Branch 2 not taken.
 940 A.setBuildMode(Matrix::row_wise);
803
1/2
✓ Branch 1 taken 940 times.
✗ Branch 2 not taken.
 1880 auto citer = A.createbegin();
804
6/8
✓ Branch 0 taken 1484666 times.
✓ Branch 1 taken 940 times.
✓ Branch 2 taken 1484666 times.
✓ Branch 3 taken 940 times.
✓ Branch 5 taken 1484666 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 1484666 times.
✗ Branch 9 not taken.
 2972152 for (auto i = sparsityPattern.begin(), end = sparsityPattern.end(); i!=end; ++i, ++citer)
805
2/2
✓ Branch 0 taken 12421722 times.
✓ Branch 1 taken 1484666 times.
 16875720 for (auto si = i->begin(), send = i->end(); si!=send; ++si)
806
2/4
✓ Branch 1 taken 12421722 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 12421722 times.
✗ Branch 5 not taken.
 24843444 citer.insert(*si);
807
808 // reset matrix to contain original values
809
1/4
✓ Branch 1 taken 940 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 940 A = 0;
810 940 const auto rowEndIt = matrixAsAssembled.end();
811
4/4
✓ Branch 0 taken 1484666 times.
✓ Branch 1 taken 940 times.
✓ Branch 2 taken 1484666 times.
✓ Branch 3 taken 940 times.
 1485606 for (auto rowIt = matrixAsAssembled.begin(); rowIt != rowEndIt; ++rowIt)
812
5/8
✗ Branch 0 not taken.
✓ Branch 1 taken 13905886 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 13905886 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 13905886 times.
✓ Branch 6 taken 12421220 times.
✓ Branch 7 taken 1484666 times.
 40232992 for (auto colIt = matrixAsAssembled[rowIt.index()].begin(); colIt != matrixAsAssembled[rowIt.index()].end(); ++colIt)
813
4/8
✓ Branch 1 taken 12421220 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 12421220 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 12421220 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 12421220 times.
✗ Branch 11 not taken.
 49684880 A[rowIt.index()][colIt.index()] = *colIt;
814
815 // The matrix A has now a possible extended sparsity pattern but any additional entries are
816 // initialized to zero and have to be filled by the sumEntries function
817 }
818 else
819 DUNE_THROW(Dune::InvalidStateException, "Cannot call extendMatrix for a grid that cannot communicate codim-" << rowDofCodim << "-entities.");
820 }
821
822 /*!
823 * \brief Sums up the entries corresponding to border entities (usually vertices or faces)
824 * \param A Matrix to operate on
825 *
826 * The idea is as follows: (Blatt and Bastian (2009) https://doi.org/10.1504/IJCSE.2008.021112)
827 * The local matrix operator stores for each row that corresponds to a dof (uniquely) owned by this process
828 * the full row of the global operator with the same entries as the global operator.
829
830 * This, together with some masking procedure (end of Section 2.4.1) allows to compute a matrix-vector product
831 * given a consistent vector representation such that the result is in a valid representation. Each valid
832 * representation can be transformed to a additive unique representation by setting all entries for
833 * dofs that are not (uniquely) owned by the process to zero.
834 */
835 940 void sumEntries(Matrix& A)
836 {
837 if constexpr (Detail::canCommunicate<typename GridView::Grid, rowDofCodim>)
838 {
839
5/6
✗ Branch 0 not taken.
✓ Branch 1 taken 182 times.
✓ Branch 2 taken 758 times.
✓ Branch 3 taken 182 times.
✓ Branch 4 taken 758 times.
✓ Branch 5 taken 182 times.
 940 if (gridView_.comm().size() <= 1)
840 ✗ return;
841
842 1880 MatrixEntryExchange dataHandle(mapper_, idToIndex_, indexToID_, A);
843 940 gridView_.communicate(
844 dataHandle, Dune::InteriorBorder_InteriorBorder_Interface, Dune::ForwardCommunication
845 );
846 }
847 else
848 DUNE_THROW(Dune::InvalidStateException,
849 "Cannot call sumEntries for a grid that cannot communicate codim-" << rowDofCodim << "-entities."
850 );
851 }
852
853 private:
854 const GridView gridView_;
855 const RowDofMapper& mapper_;
856 std::map<IdType, int> idToIndex_;
857 std::map<int, IdType> indexToID_;
858
859 };
860
861 /*!
862 * \brief Prepare a matrix for parallel solvers
863 */
864 template<class LinearSolverTraits, class ParallelTraits,
865 class Matrix, class ParallelHelper>
866 1880 void prepareMatrixParallel(Matrix& A, ParallelHelper& pHelper)
867 {
868 if constexpr (ParallelTraits::isNonOverlapping)
869 {
870 // extend the matrix pattern such that it is usable for a parallel solver
871 // and make right-hand side consistent
872 using GridView = typename LinearSolverTraits::GridView;
873 using DofMapper = typename LinearSolverTraits::DofMapper;
874 static constexpr int dofCodim = LinearSolverTraits::dofCodim;
875 3760 ParallelMatrixHelper<Matrix, GridView, DofMapper, dofCodim> matrixHelper(pHelper.gridView(), pHelper.dofMapper());
876
7/24
✓ Branch 0 taken 262476 times.
✓ Branch 1 taken 940 times.
✓ Branch 2 taken 262476 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2334 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2334 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 4448 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 4448 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
 540396 matrixHelper.extendMatrix(A, [&pHelper](auto idx){ return pHelper.isGhost(idx); });
877
1/2
✓ Branch 1 taken 940 times.
✗ Branch 2 not taken.
 1880 matrixHelper.sumEntries(A);
878 }
879 1880 }
880
881 /*!
882 * \brief Prepare a vector for parallel solvers
883 */
884 template<class LinearSolverTraits, class ParallelTraits,
885 class Vector, class ParallelHelper>
886 ✗ void prepareVectorParallel(Vector& b, ParallelHelper& pHelper)
887 {
888 if constexpr (ParallelTraits::isNonOverlapping)
889 {
890 // extend the matrix pattern such that it is usable for a parallel solver
891 // and make right-hand side consistent
892 using GridView = typename LinearSolverTraits::GridView;
893 using DofMapper = typename LinearSolverTraits::DofMapper;
894 static constexpr int dofCodim = LinearSolverTraits::dofCodim;
895
0/5
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 858 ParallelVectorHelper<GridView, DofMapper, dofCodim> vectorHelper(pHelper.gridView(), pHelper.dofMapper());
896
0/5
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 286 vectorHelper.makeNonOverlappingConsistent(b);
897 }
898 ✗ }
899
900 /*!
901 * \brief Prepare linear algebra variables for parallel solvers
902 */
903 template<class LinearSolverTraits, class ParallelTraits,
904 class Matrix, class Vector, class ParallelHelper>
905 572 void prepareLinearAlgebraParallel(Matrix& A, Vector& b, ParallelHelper& pHelper)
906 {
907 966 prepareMatrixParallel<LinearSolverTraits, ParallelTraits>(A, pHelper);
908 966 prepareVectorParallel<LinearSolverTraits, ParallelTraits>(b, pHelper);
909 572 }
910
911 /*!
912 * \brief Prepare linear algebra variables for parallel solvers
913 */
914 template<class LinearSolverTraits, class ParallelTraits,
915 class Matrix, class Vector, class ParallelHelper>
916 void prepareLinearAlgebraParallel(Matrix& A, Vector& b,
917 std::shared_ptr<typename ParallelTraits::Comm>& comm,
918 std::shared_ptr<typename ParallelTraits::LinearOperator>& fop,
919 std::shared_ptr<typename ParallelTraits::ScalarProduct>& sp,
920 ParallelHelper& pHelper)
921 {
922 prepareLinearAlgebraParallel<LinearSolverTraits, ParallelTraits>(A, b, pHelper);
923 const auto category = ParallelTraits::isNonOverlapping ?
924 Dune::SolverCategory::nonoverlapping : Dune::SolverCategory::overlapping;
925
926 comm = std::make_shared<typename ParallelTraits::Comm>(pHelper.gridView().comm(), category);
927 pHelper.createParallelIndexSet(*comm);
928 fop = std::make_shared<typename ParallelTraits::LinearOperator>(A, *comm);
929 sp = std::make_shared<typename ParallelTraits::ScalarProduct>(*comm);
930 }
931
932 } // end namespace Dumux
933
934 #endif
935