GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/linear/parallelhelpers.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 249 258 96.5%
Functions: 261 1685 15.5%
Branches: 226 338 66.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 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 48 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 260 BaseGatherScatter(const DofMapper& mapper)
45 260 : mapper_(mapper) {}
46
47 template<class EntityType>
48 43222 int index(const EntityType& e) const
49 86444 { return mapper_.index(e); }
50
51 916 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 13972 bool isNeitherInteriorNorBorderEntity(EntityType& e) const
64
14/22
✓ Branch 0 taken 1144 times.
✓ Branch 1 taken 300 times.
✓ Branch 2 taken 844 times.
✓ Branch 3 taken 300 times.
✓ Branch 4 taken 264 times.
✓ Branch 5 taken 444 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 264 times.
✓ Branch 8 taken 1345 times.
✓ Branch 9 taken 225 times.
✓ Branch 10 taken 1195 times.
✓ Branch 11 taken 603 times.
✓ Branch 12 taken 280 times.
✓ Branch 13 taken 840 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 280 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
 8590 { 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 52 GhostGatherScatter(std::vector<std::size_t>& ranks, const DofMapper& mapper)
87 : BaseGatherScatter(mapper)
88 52 , ranks_(ranks)
89 {}
90
91 template<class MessageBuffer, class EntityType>
92 9616 void gather(MessageBuffer& buff, const EntityType& e) const
93 {
94
2/2
✓ Branch 1 taken 330 times.
✓ Branch 2 taken 4052 times.
 9616 auto& data = ranks_[this->index(e)];
95
3/3
✓ Branch 1 taken 4052 times.
✓ Branch 2 taken 334 times.
✓ Branch 0 taken 330 times.
 10904 if (this->isNeitherInteriorNorBorderEntity(e))
96 ✗ data = ghostMarker_;
97 9616 buff.write(data);
98 9616 }
99
100 template<class MessageBuffer, class EntityType>
101 9616 void scatter(MessageBuffer& buff, const EntityType& e, std::size_t n)
102 {
103 std::size_t x;
104 9616 buff.read(x);
105
2/2
✓ Branch 1 taken 3204 times.
✓ Branch 2 taken 1178 times.
 9616 auto& data = ranks_[this->index(e)];
106
3/3
✓ Branch 1 taken 1478 times.
✓ Branch 2 taken 34 times.
✓ Branch 0 taken 3204 times.
 9616 if (this->isNeitherInteriorNorBorderEntity(e))
107 7062 data = ghostMarker_;
108 9616 }
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 52 InteriorBorderGatherScatter(std::vector<std::size_t>& ranks, const DofMapper& mapper)
130 : BaseGatherScatter(mapper)
131 52 , 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 2554 void gather(MessageBuffer& buff, const EntityType& e) const
139 {
140
2/2
✓ Branch 1 taken 264 times.
✓ Branch 2 taken 1128 times.
 2554 auto& data = ranks_[this->index(e)];
141
3/3
✓ Branch 1 taken 1128 times.
✓ Branch 2 taken 34 times.
✓ Branch 0 taken 264 times.
 2554 if (this->isNeitherInteriorNorBorderEntity(e))
142 ✗ data = ghostMarker_;
143 2554 buff.write(data);
144 2554 }
145
146 /*!
147 * \brief Unpack data (received from other process) from message buffer to this process
148 */
149 template<class MessageBuffer, class EntityType>
150 2554 void scatter(MessageBuffer& buff, const EntityType& e, std::size_t n)
151 {
152 std::size_t x;
153 2554 buff.read(x);
154
2/2
✓ Branch 1 taken 264 times.
✓ Branch 2 taken 1128 times.
 2554 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
4/4
✓ Branch 0 taken 264 times.
✓ Branch 1 taken 1162 times.
✓ Branch 2 taken 696 times.
✓ Branch 3 taken 696 times.
 4112 data = this->isNeitherInteriorNorBorderEntity(e) ? x : min(data, x);
161 2554 }
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 52 NeighbourGatherScatter(const DofMapper& mapper, int rank, std::set<int>& neighbours)
180 : BaseGatherScatter(mapper)
181 52 , rank_(rank)
182
1/2
✓ Branch 1 taken 52 times.
✗ Branch 2 not taken.
 52 , neighbours_(neighbours)
183 {}
184
185 template<class MessageBuffer, class EntityType>
186 9750 void gather(MessageBuffer& buff, const EntityType& e) const
187 {
188
2/4
✓ Branch 1 taken 800 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 800 times.
 10550 buff.write(rank_);
189 }
190
191 template<class MessageBuffer, class EntityType>
192
1/2
✓ Branch 1 taken 2564 times.
✗ Branch 2 not taken.
 9750 void scatter(MessageBuffer& buff, const EntityType& e, std::size_t n)
193 {
194 int x;
195 9750 buff.read(x);
196
3/4
✓ Branch 1 taken 800 times.
✓ Branch 2 taken 1764 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 800 times.
 9750 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 52 SharedGatherScatter(std::vector<int>& shared, const DofMapper& mapper)
218 : BaseGatherScatter(mapper)
219 52 , shared_(shared)
220 {}
221
222 template<class MessageBuffer, class EntityType>
223 9750 void gather(MessageBuffer& buff, EntityType& e) const
224 {
225
1/2
✓ Branch 1 taken 800 times.
✗ Branch 2 not taken.
 9750 int data = true;
226
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 800 times.
 9750 buff.write(data);
227 }
228
229 template<class MessageBuffer, class EntityType>
230 16734 void scatter(MessageBuffer& buff, const EntityType &e, std::size_t n)
231 {
232 int x;
233 16734 buff.read(x);
234
2/2
✓ Branch 1 taken 8293 times.
✓ Branch 2 taken 1457 times.
 16734 auto& data = shared_[this->index(e)];
235
3/4
✓ Branch 0 taken 8293 times.
✓ Branch 1 taken 1457 times.
✓ Branch 2 taken 8293 times.
✗ Branch 3 not taken.
 16734 data = data || x;
236 16734 }
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 52 GlobalIndexGatherScatter(std::vector<GlobalIndex>& globalIndices, const DofMapper& mapper)
256 : BaseGatherScatter(mapper)
257 52 , globalIndices_(globalIndices)
258 {}
259
260 template<class MessageBuffer, class EntityType>
261 16734 void gather(MessageBuffer& buff, const EntityType& e) const
262 {
263 16734 buff.write(globalIndices_[this->index(e)]);
264 16734 }
265
266 template<class MessageBuffer, class EntityType>
267 16734 void scatter(MessageBuffer& buff, const EntityType& e, std::size_t n)
268 {
269 16734 DataType x;
270 16734 buff.read(x);
271 using std::min;
272 16734 DataType& data = globalIndices_[this->index(e)];
273 16734 data = min(data, x);
274 16734 }
275 private:
276 std::vector<GlobalIndex>& globalIndices_;
277 };
278
279 public:
280
281 52 ParallelISTLHelperImpl(const GridView& gridView, const DofMapper& mapper)
282
2/5
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 1 taken 50 times.
 52 : gridView_(gridView), mapper_(mapper)
283 {
284 if constexpr (Detail::canCommunicate<typename GridView::Traits::Grid, dofCodim>)
285
1/2
✓ Branch 1 taken 52 times.
✗ Branch 2 not taken.
 52 initGhostsAndOwners_();
286 else
287 DUNE_THROW(Dune::InvalidStateException,
288 "Cannot initialize parallel helper for a grid that cannot communicate codim-" << dofCodim << "-entities."
289 );
290 52 }
291
292 277551 bool isGhost(std::size_t i) const
293 277551 { 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 52 void createParallelIndexSet(Comm& comm) const
311 {
312 if constexpr (Detail::canCommunicate<typename GridView::Traits::Grid, dofCodim>)
313 {
314
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
✓ Branch 2 taken 20 times.
✓ Branch 3 taken 30 times.
 52 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
2/3
✓ Branch 2 taken 48 times.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
 52 std::vector<int> isShared(mapper_.size(), false);
322
1/2
✓ Branch 1 taken 52 times.
✗ Branch 2 not taken.
 52 SharedGatherScatter sgs(isShared, mapper_);
323
2/4
✓ Branch 1 taken 52 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 30 times.
✗ Branch 5 not taken.
 74 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 52 times.
✗ Branch 2 not taken.
 52 GlobalIndex count = 0;
328
2/2
✓ Branch 0 taken 50081 times.
✓ Branch 1 taken 52 times.
 50133 for (std::size_t i = 0; i < isShared.size(); ++i)
329
4/4
✓ Branch 0 taken 8293 times.
✓ Branch 1 taken 41788 times.
✓ Branch 2 taken 3454 times.
✓ Branch 3 taken 4839 times.
 50081 if (isShared[i] && isOwned_[i] == 1)
330 50081 ++count;
331
332
3/5
✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
✓ Branch 3 taken 30 times.
✗ Branch 4 not taken.
✓ Branch 2 taken 20 times.
 52 std::vector<GlobalIndex> counts(gridView_.comm().size());
333
3/5
✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
✓ Branch 3 taken 30 times.
✗ Branch 4 not taken.
✓ Branch 2 taken 20 times.
 52 gridView_.comm().allgather(&count, 1, counts.data());
334
335 // Compute start index start_p = \sum_{i=0}^{i<p} counts_i
336
3/5
✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
✓ Branch 3 taken 30 times.
✗ Branch 4 not taken.
✓ Branch 2 taken 20 times.
 52 const int rank = gridView_.comm().rank();
337
2/4
✓ Branch 1 taken 52 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 52 times.
✗ Branch 5 not taken.
 52 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
2/7
✓ Branch 1 taken 52 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 3 not taken.
 52 std::vector<GlobalIndex> globalIndices(mapper_.size(), std::numeric_limits<GlobalIndex>::max());
341
2/2
✓ Branch 0 taken 50081 times.
✓ Branch 1 taken 52 times.
 50133 for (std::size_t i = 0; i < globalIndices.size(); ++i)
342 {
343
4/4
✓ Branch 0 taken 45242 times.
✓ Branch 1 taken 4839 times.
✓ Branch 2 taken 3454 times.
✓ Branch 3 taken 41788 times.
 50081 if (isOwned_[i] == 1 && isShared[i])
344 {
345 3454 globalIndices[i] = start;
346 50081 ++start;
347 }
348 }
349
350 // publish global indices for the shared DOFS to other processors.
351
1/2
✓ Branch 1 taken 52 times.
✗ Branch 2 not taken.
 52 GlobalIndexGatherScatter<GlobalIndex> gigs(globalIndices, mapper_);
352
2/4
✓ Branch 1 taken 52 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 30 times.
✗ Branch 5 not taken.
 74 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 52 times.
✗ Branch 2 not taken.
 52 resizeAndFillIndexSet_(comm, globalIndices);
356
357 // Compute neighbours using communication
358 52 std::set<int> neighbours;
359
3/5
✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
✓ Branch 3 taken 30 times.
✗ Branch 4 not taken.
✓ Branch 2 taken 20 times.
 52 NeighbourGatherScatter ngs(mapper_, gridView_.comm().rank(), neighbours);
360
3/5
✓ Branch 1 taken 52 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 22 times.
✗ Branch 6 not taken.
✓ Branch 4 taken 30 times.
 74 gridView_.communicate(ngs, Dune::All_All_Interface, Dune::ForwardCommunication);
361
1/2
✓ Branch 1 taken 52 times.
✗ Branch 2 not taken.
 52 comm.remoteIndices().setNeighbours(neighbours);
362
363
1/2
✓ Branch 1 taken 52 times.
✗ Branch 2 not taken.
 52 comm.remoteIndices().template rebuild<false>();
364 156 }
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 2140 const DofMapper& dofMapper() const
373 2140 { return mapper_; }
374
375 //! Return the gridView
376 948 const GridView& gridView() const
377
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 106 times.
 2242 { 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 52 void initGhostsAndOwners_()
389 {
390
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 30 times.
 52 const auto rank = gridView_.comm().rank();
391 52 isOwned_.assign(mapper_.size(), rank);
392
393 // find out about ghosts
394
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 30 times.
 52 GhostGatherScatter ggs(isOwned_, mapper_);
395
396
3/4
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 50 times.
✓ Branch 2 taken 30 times.
✗ Branch 3 not taken.
 52 if (gridView_.comm().size() > 1)
397 52 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 52 isGhost_ = isOwned_;
402
403 // partition interior/border uniquely
404 // after the communication each dof is assigned to one unique owner process rank
405
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 30 times.
 52 InteriorBorderGatherScatter dh(isOwned_, mapper_);
406
407
3/4
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 50 times.
✓ Branch 2 taken 30 times.
✗ Branch 3 not taken.
 52 if (gridView_.comm().size() > 1)
408 52 gridView_.communicate(dh, Dune::InteriorBorder_InteriorBorder_Interface, Dune::ForwardCommunication);
409
410 // convert vector into mask vector
411
2/2
✓ Branch 0 taken 50081 times.
✓ Branch 1 taken 52 times.
 50133 for (auto& v : isOwned_)
412
2/2
✓ Branch 0 taken 4839 times.
✓ Branch 1 taken 45242 times.
 54920 v = (v == rank) ? 1 : 0;
413 52 }
414
415 template<class Comm, class GlobalIndices>
416 52 void resizeAndFillIndexSet_(Comm& comm, const GlobalIndices& globalIndices) const
417 {
418 52 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
2/2
✓ Branch 1 taken 50081 times.
✓ Branch 2 taken 52 times.
 50133 for (std::size_t localIndex = 0; localIndex < globalIndices.size(); ++localIndex)
423 {
424 50081 const auto globalIndex = globalIndices[localIndex];
425
2/2
✓ Branch 0 taken 8293 times.
✓ Branch 1 taken 41788 times.
 100162 if (globalIndex != std::numeric_limits<typename GlobalIndices::value_type>::max())
426 {
427
2/2
✓ Branch 0 taken 4839 times.
✓ Branch 1 taken 3454 times.
 8293 const bool isOwned = isOwned_[localIndex] > 0;
428 16586 const auto attr = getAttribute_(comm, isOwned, isGhost(localIndex));
429 using LocalIndex = typename Comm::ParallelIndexSet::LocalIndex;
430 8293 comm.indexSet().add(globalIndex, LocalIndex{localIndex, attr});
431 }
432 }
433
434 52 comm.indexSet().endResize();
435 52 }
436
437 template<class Comm>
438 Dune::OwnerOverlapCopyAttributeSet::AttributeSet
439 8293 getAttribute_(const Comm& comm, const bool isOwned, const bool isGhost) const
440 {
441
2/2
✓ Branch 0 taken 4839 times.
✓ Branch 1 taken 3454 times.
 8293 if (isOwned) // this process owns the dof (uniquely)
442 return Dune::OwnerOverlapCopyAttributeSet::owner;
443
4/4
✓ Branch 0 taken 4126 times.
✓ Branch 1 taken 713 times.
✓ Branch 2 taken 1608 times.
✓ Branch 3 taken 2518 times.
 4839 else if (isGhost && (comm.category() == Dune::SolverCategory::nonoverlapping) )
444 return Dune::OwnerOverlapCopyAttributeSet::overlap;
445 else
446 2321 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 1200 ParallelVectorHelper(const GridView& gridView, const DofMapper& mapper)
479
1/2
✓ Branch 1 taken 162 times.
✗ Branch 2 not taken.
 262 : 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
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1200 times.
✓ Branch 2 taken 308 times.
✗ Branch 3 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 940 : rowEntityMapper_(rowEntityMapper), idToIndex_(globalToLocal), indexToID_(localToGlobal)
568 940 , sparsityPattern_(sparsityPattern), A_(A), isGhostColumDof_(isGhostColumDof)
569 {
570
1/4
✓ Branch 1 taken 940 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 940 sparsityPattern_.resize(A.N());
571 940 }
572
573 /*!
574 * \brief Returns true if data for given valid codim should be communicated
575 */
576 48 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 162106 std::size_t size(EntityType& e) const
588 {
589 162106 const auto rowIdx = rowEntityMapper_.index(e);
590
591 // all column entity indices of this row that are in the index set
592 162106 std::size_t n = 0;
593
2/2
✓ Branch 0 taken 908884 times.
✓ Branch 1 taken 93300 times.
 1781348 for (auto colIt = A_[rowIdx].begin(); colIt != A_[rowIdx].end(); ++colIt)
594 3238484 if (indexToID_.count(colIt.index()))
595 543120 n++;
596
597 162106 return n;
598 }
599
600 /*!
601 * \brief Pack data from user to message buffer
602 */
603 template<class MessageBuffer, class EntityType>
604 140686 void gather(MessageBuffer& buff, const EntityType& e) const
605 {
606 140686 const auto rowIdx = rowEntityMapper_.index(e);
607
608 // send all column entity ids of this row that are in the index set
609
4/4
✓ Branch 0 taken 8896 times.
✓ Branch 1 taken 795468 times.
✓ Branch 2 taken 723588 times.
✓ Branch 3 taken 70336 times.
 1583528 for (auto colIt = A_[rowIdx].begin(); colIt != A_[rowIdx].end(); ++colIt)
610
2/2
✓ Branch 0 taken 269466 times.
✓ Branch 1 taken 463018 times.
 2885684 if (auto it = indexToID_.find(colIt.index()); it != indexToID_.end())
611 1950606 buff.write(it->second);
612 140686 }
613
614 /*!
615 * \brief Unpack data from message buffer to user
616 */
617 template<class MessageBuffer, class EntityType>
618 140686 void scatter(MessageBuffer& buff, const EntityType& e, std::size_t n)
619 {
620 140686 const auto rowIdx = rowEntityMapper_.index(e);
621
2/2
✓ Branch 1 taken 269466 times.
✓ Branch 2 taken 71880 times.
 670576 for (std::size_t k = 0; k < n; k++)
622 {
623 // receive n column entity IDs
624
2/2
✓ Branch 0 taken 4448 times.
✓ Branch 1 taken 184 times.
 525480 IdType id;
625 529890 buff.read(id);
626
627 // only add entries that are contained in the index set
628
3/3
✓ Branch 1 taken 260584 times.
✓ Branch 2 taken 24 times.
✓ Branch 0 taken 8858 times.
 534300 if (const auto it = idToIndex_.find(id); it != idToIndex_.end())
629 {
630 529658 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 529658 if (!isGhostColumDof_(colIdx))
634 // std::set takes care that each index only is inserted once
635 529658 sparsityPattern_[rowIdx].insert(colIdx);
636 }
637 }
638 140686 }
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
2/4
✓ Branch 0 taken 4448 times.
✓ Branch 1 taken 184 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 814980 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 48 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 162106 std::size_t size(EntityType& e) const
686 {
687 162106 const auto rowIdx = rowEntityMapper_.index(e);
688 162106 std::size_t n = 0;
689
2/2
✓ Branch 0 taken 909386 times.
✓ Branch 1 taken 93300 times.
 1782352 for (auto colIt = A_[rowIdx].begin(); colIt != A_[rowIdx].end(); ++colIt)
690 3240492 if (indexToID_.count(colIt.index()))
691 544124 n++;
692
693 162106 return n;
694 }
695
696 /*!
697 * \brief Pack data from user to message buffer
698 */
699 template<class MessageBuffer, class EntityType>
700 140686 void gather(MessageBuffer& buff, const EntityType& e) const
701 {
702 140686 const auto rowIdx = rowEntityMapper_.index(e);
703 // send all matrix entries for which the column index is in the index set
704
4/4
✓ Branch 0 taken 8896 times.
✓ Branch 1 taken 795970 times.
✓ Branch 2 taken 724090 times.
✓ Branch 3 taken 70336 times.
 1584532 for (auto colIt = A_[rowIdx].begin(); colIt != A_[rowIdx].end(); ++colIt)
705
2/2
✓ Branch 0 taken 269968 times.
✓ Branch 1 taken 463018 times.
 2887692 if (auto it = indexToID_.find(colIt.index()); it != indexToID_.end())
706 530894 buff.write(MatrixEntry{ it->second, *colIt });
707 140686 }
708
709 /*!
710 * \brief Unpack data from message buffer to user
711 */
712 template<class MessageBuffer, class EntityType>
713 140686 void scatter(MessageBuffer& buff, const EntityType& e, std::size_t n)
714 {
715 140686 const auto rowIdx = rowEntityMapper_.index(e);
716
2/2
✓ Branch 1 taken 269968 times.
✓ Branch 2 taken 71880 times.
 671580 for (std::size_t k = 0; k < n; k++)
717 {
718
2/2
✓ Branch 0 taken 4448 times.
✓ Branch 1 taken 184 times.
 530894 MatrixEntry m;
719 530894 buff.read(m);
720 530894 const auto& [colDofID, matrixBlock] = m; // unpack
721
722 // only add entries in the index set and for which A has allocated memory
723
3/3
✓ Branch 1 taken 261086 times.
✓ Branch 2 taken 24 times.
✓ Branch 0 taken 8858 times.
 535304 if (auto it = idToIndex_.find(colDofID); it != idToIndex_.end())
724
2/4
✗ Branch 1 not taken.
✓ Branch 2 taken 269760 times.
✓ Branch 3 taken 269760 times.
✗ Branch 4 not taken.
 530662 if (A_[rowIdx].find(it->second) != A_[rowIdx].end())
725 1039374 A_[rowIdx][it->second] += matrixBlock;
726 }
727 140686 }
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 2 not taken.
✗ Branch 3 not taken.
 940 : gridView_(gridView), mapper_(mapper)
741 {
742 940 idToIndex_.clear();
743
1/2
✓ Branch 2 taken 940 times.
✗ Branch 3 not taken.
 940 indexToID_.clear();
744
745
3/6
✓ Branch 1 taken 940 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 940 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 940 times.
✗ Branch 8 not taken.
 940 std::vector<bool> handledDof(gridView_.size(rowDofCodim), false);
746
7/12
✓ Branch 1 taken 940 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36886 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 848 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1290386 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 1290386 times.
✓ Branch 13 taken 758 times.
✓ Branch 6 taken 18720 times.
✗ Branch 3 not taken.
 2656908 for (const auto& element : elements(gridView_))
747 {
748
4/5
✓ Branch 1 taken 6555570 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 5229148 times.
✓ Branch 4 taken 1309106 times.
✓ Branch 0 taken 69264 times.
 13069488 for (int i = 0; i < element.subEntities(rowDofCodim); ++i)
749 {
750
2/3
✓ Branch 1 taken 274362 times.
✗ Branch 2 not taken.
✓ Branch 0 taken 5024050 times.
 5298412 const auto entity = element.template subEntity<rowDofCodim>(i);
751
2/2
✓ Branch 0 taken 12150 times.
✓ Branch 1 taken 131994 times.
 5298412 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
1/2
✓ Branch 0 taken 284302 times.
✗ Branch 1 not taken.
 284302 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 272152 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
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 182 times.
✓ Branch 2 taken 758 times.
✓ Branch 3 taken 182 times.
 940 if (gridView_.comm().size() <= 1)
776 ✗ return;
777
778 // make a copy of the matrix as originally assembled
779 940 Matrix matrixAsAssembled(A);
780
781 // first get entries in sparsity pattern from other processes
782 940 std::size_t numNonZeroEntries = 0;
783 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.
 1698 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
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1485606 times.
✓ Branch 2 taken 1484666 times.
✓ Branch 3 taken 940 times.
 1485606 for (auto rowIt = A.begin(); rowIt != A.end(); ++rowIt)
791 {
792 1484666 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.
 13905886 for (auto colIt = rowIt->begin(); colIt != colEndIt; ++colIt)
794 24842440 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 1484666 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 940 auto citer = A.createbegin();
804
3/4
✓ Branch 1 taken 1484666 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1484666 times.
✓ Branch 4 taken 940 times.
 1485606 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.
 13906388 for (auto si = i->begin(), send = i->end(); si!=send; ++si)
806
1/2
✓ Branch 1 taken 12421722 times.
✗ Branch 2 not taken.
 24843444 citer.insert(*si);
807
808 // reset matrix to contain original values
809
1/2
✓ Branch 1 taken 940 times.
✗ Branch 2 not taken.
 940 A = 0;
810 940 const auto rowEndIt = matrixAsAssembled.end();
811
2/2
✓ Branch 0 taken 1484666 times.
✓ Branch 1 taken 940 times.
 1485606 for (auto rowIt = matrixAsAssembled.begin(); rowIt != rowEndIt; ++rowIt)
812
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 13905886 times.
✓ Branch 2 taken 12421220 times.
✓ Branch 3 taken 1484666 times.
 13905886 for (auto colIt = matrixAsAssembled[rowIt.index()].begin(); colIt != matrixAsAssembled[rowIt.index()].end(); ++colIt)
813
1/2
✓ Branch 1 taken 12421220 times.
✗ Branch 2 not taken.
 12421220 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 940 }
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
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 182 times.
✓ Branch 2 taken 758 times.
✓ Branch 3 taken 182 times.
 940 if (gridView_.comm().size() <= 1)
840 ✗ return;
841
842 940 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 940 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 940 ParallelMatrixHelper<Matrix, GridView, DofMapper, dofCodim> matrixHelper(pHelper.gridView(), pHelper.dofMapper());
876
4/11
✓ Branch 0 taken 266924 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2334 times.
✓ Branch 3 taken 182 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 758 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 8 not taken.
 270198 matrixHelper.extendMatrix(A, [&pHelper](auto idx){ return pHelper.isGhost(idx); });
877
1/2
✓ Branch 1 taken 940 times.
✗ Branch 2 not taken.
 940 matrixHelper.sumEntries(A);
878 940 }
879 940 }
880
881 /*!
882 * \brief Prepare a vector for parallel solvers
883 */
884 template<class LinearSolverTraits, class ParallelTraits,
885 class Vector, class ParallelHelper>
886 1038 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 1038 ParallelVectorHelper<GridView, DofMapper, dofCodim> vectorHelper(pHelper.gridView(), pHelper.dofMapper());
896
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 1038 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 286 void prepareLinearAlgebraParallel(Matrix& A, Vector& b, ParallelHelper& pHelper)
906 {
907 286 prepareMatrixParallel<LinearSolverTraits, ParallelTraits>(A, pHelper);
908 286 prepareVectorParallel<LinearSolverTraits, ParallelTraits>(b, pHelper);
909 286 }
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