GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/assembly/fvassembler.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 165 173 95.4%
Functions: 1999 2479 80.6%
Branches: 168 359 46.8%
Line Branch Exec Source
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 //
4 // SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5 // SPDX-License-Identifier: GPL-3.0-or-later
6 //
7 /*!
8 * \file
9 * \ingroup Assembly
10 * \brief A linear system assembler (residual and Jacobian) for finite volume schemes
11 */
12 #ifndef DUMUX_FV_ASSEMBLER_HH
13 #define DUMUX_FV_ASSEMBLER_HH
14
15 #include <vector>
16 #include <deque>
17 #include <type_traits>
18 #include <memory>
19
20 #include <dune/istl/matrixindexset.hh>
21
22 #include <dumux/common/properties.hh>
23 #include <dumux/common/timeloop.hh>
24 #include <dumux/common/gridcapabilities.hh>
25 #include <dumux/common/typetraits/vector.hh>
26 #include <dumux/common/typetraits/periodic.hh>
27
28 #include <dumux/discretization/method.hh>
29 #include <dumux/linear/parallelhelpers.hh>
30 #include <dumux/linear/dunevectors.hh>
31
32 #include <dumux/assembly/coloring.hh>
33 #include <dumux/assembly/jacobianpattern.hh>
34 #include <dumux/assembly/diffmethod.hh>
35
36 #include <dumux/parallel/multithreading.hh>
37 #include <dumux/parallel/parallel_for.hh>
38
39 #include "cvfelocalassembler.hh"
40 #include "cclocalassembler.hh"
41 #include "fclocalassembler.hh"
42
43 namespace Dumux::Detail {
44
45 template<class DiscretizationMethod>
46 struct LocalAssemblerChooser;
47
48 template<class DM>
49 struct LocalAssemblerChooser<DiscretizationMethods::CVFE<DM>>
50 {
51 template<class TypeTag, class Impl, DiffMethod diffMethod, bool isImplicit>
52 using type = CVFELocalAssembler<TypeTag, Impl, diffMethod, isImplicit>;
53 };
54
55 template<>
56 struct LocalAssemblerChooser<DiscretizationMethods::CCMpfa>
57 {
58 template<class TypeTag, class Impl, DiffMethod diffMethod, bool isImplicit>
59 using type = CCLocalAssembler<TypeTag, Impl, diffMethod, isImplicit>;
60 };
61
62 template<>
63 struct LocalAssemblerChooser<DiscretizationMethods::CCTpfa>
64 {
65 template<class TypeTag, class Impl, DiffMethod diffMethod, bool isImplicit>
66 using type = CCLocalAssembler<TypeTag, Impl, diffMethod, isImplicit>;
67 };
68
69 template<>
70 struct LocalAssemblerChooser<DiscretizationMethods::FCStaggered>
71 {
72 template<class TypeTag, class Impl, DiffMethod diffMethod, bool isImplicit>
73 using type = FaceCenteredLocalAssembler<TypeTag, Impl, diffMethod, isImplicit>;
74 };
75
76 template<class TypeTag, class Impl, DiffMethod diffMethod, bool isImplicit>
77 using LocalAssemblerChooser_t = typename LocalAssemblerChooser<
78 typename GetPropType<TypeTag, Properties::GridGeometry>::DiscretizationMethod
79 >::template type<TypeTag, Impl, diffMethod, isImplicit>;
80
81 } // end namespace Dumux::Detail
82
83 namespace Dumux {
84
85 /*!
86 * \ingroup Assembly
87 * \brief A linear system assembler (residual and Jacobian) for finite volume schemes (box, tpfa, mpfa, ...)
88 * \tparam TypeTag The TypeTag
89 * \tparam diffMethod The differentiation method to residual compute derivatives
90 * \tparam isImplicit Specifies whether the time discretization is implicit or not not (i.e. explicit)
91 * \tparam LocalResidual The local residual (local operator) of the chosen model
92 */
93 template<class TypeTag, DiffMethod diffMethod, bool isImplicit = true, class LocalResidual = GetPropType<TypeTag, Properties::LocalResidual>>
94 class FVAssembler
95 {
96 using GridGeo = GetPropType<TypeTag, Properties::GridGeometry>;
97 using GridView = typename GridGeo::GridView;
98 using Element = typename GridView::template Codim<0>::Entity;
99 using ElementSeed = typename GridView::Grid::template Codim<0>::EntitySeed;
100 using TimeLoop = TimeLoopBase<GetPropType<TypeTag, Properties::Scalar>>;
101
102 static constexpr bool isBox = GridGeo::discMethod == DiscretizationMethods::box;
103
104 using ThisType = FVAssembler<TypeTag, diffMethod, isImplicit, LocalResidual>;
105 using LocalAssembler = typename Detail::LocalAssemblerChooser_t<TypeTag, ThisType, diffMethod, isImplicit>;
106
107 public:
108 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
109 using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>;
110 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
111 using ResidualType = typename Detail::NativeDuneVectorType<SolutionVector>::type;
112
113 using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
114
115 using GridGeometry = GridGeo;
116 using Problem = GetPropType<TypeTag, Properties::Problem>;
117
118 /*!
119 * \brief The constructor for stationary problems
120 * \note the grid variables might be temporarily changed during assembly (if caching is enabled)
121 * it is however guaranteed that the state after assembly will be the same as before
122 */
123 131 FVAssembler(std::shared_ptr<const Problem> problem,
124 std::shared_ptr<const GridGeometry> gridGeometry,
125 std::shared_ptr<GridVariables> gridVariables)
126 131 : problem_(problem)
127 131 , gridGeometry_(gridGeometry)
128 131 , gridVariables_(gridVariables)
129
1/2
✓ Branch 1 taken 129 times.
✗ Branch 2 not taken.
 131 , timeLoop_()
130
1/2
✓ Branch 1 taken 129 times.
✗ Branch 2 not taken.
 131 , isStationaryProblem_(true)
131 {
132 static_assert(isImplicit, "Explicit assembler for stationary problem doesn't make sense!");
133 131 enableMultithreading_ = SupportsColoring<typename GridGeometry::DiscretizationMethod>::value
134
1/3
✓ Branch 1 taken 129 times.
✗ Branch 2 not taken.
✗ Branch 0 not taken.
 131 && Grid::Capabilities::supportsMultithreading(gridGeometry_->gridView())
135 && !Multithreading::isSerial()
136
4/7
✓ Branch 0 taken 28 times.
✓ Branch 1 taken 100 times.
✓ Branch 3 taken 28 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 28 times.
✗ Branch 2 not taken.
 130 && getParam<bool>("Assembly.Multithreading", true);
137
138
1/2
✓ Branch 1 taken 129 times.
✗ Branch 2 not taken.
 131 maybeComputeColors_();
139 131 }
140
141 /*!
142 * \brief The constructor for instationary problems
143 * \note the grid variables might be temporarily changed during assembly (if caching is enabled)
144 * it is however guaranteed that the state after assembly will be the same as before
145 */
146 253 FVAssembler(std::shared_ptr<const Problem> problem,
147 std::shared_ptr<const GridGeometry> gridGeometry,
148 std::shared_ptr<GridVariables> gridVariables,
149 std::shared_ptr<const TimeLoop> timeLoop,
150 const SolutionVector& prevSol)
151 253 : problem_(problem)
152 253 , gridGeometry_(gridGeometry)
153 253 , gridVariables_(gridVariables)
154 253 , timeLoop_(timeLoop)
155 253 , prevSol_(&prevSol)
156
1/2
✓ Branch 1 taken 251 times.
✗ Branch 2 not taken.
 253 , isStationaryProblem_(!timeLoop)
157 {
158 253 enableMultithreading_ = SupportsColoring<typename GridGeometry::DiscretizationMethod>::value
159
1/3
✓ Branch 1 taken 251 times.
✗ Branch 2 not taken.
✗ Branch 0 not taken.
 253 && Grid::Capabilities::supportsMultithreading(gridGeometry_->gridView())
160 && !Multithreading::isSerial()
161
4/7
✓ Branch 0 taken 33 times.
✓ Branch 1 taken 218 times.
✓ Branch 3 taken 33 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 33 times.
✗ Branch 2 not taken.
 253 && getParam<bool>("Assembly.Multithreading", true);
162
163
1/2
✓ Branch 1 taken 251 times.
✗ Branch 2 not taken.
 253 maybeComputeColors_();
164 253 }
165
166 /*!
167 * \brief Assembles the global Jacobian of the residual
168 * and the residual for the current solution.
169 */
170 template<class PartialReassembler = DefaultPartialReassembler>
171 57218 void assembleJacobianAndResidual(const SolutionVector& curSol, const PartialReassembler* partialReassembler = nullptr)
172 {
173 57218 checkAssemblerState_();
174 57218 resetJacobian_(partialReassembler);
175 57218 resetResidual_();
176
177 5376890 assemble_([&](const Element& element)
178 {
179
1/2
✓ Branch 1 taken 926074 times.
✗ Branch 2 not taken.
 37592848 LocalAssembler localAssembler(*this, element, curSol);
180
4/7
✓ Branch 1 taken 37204668 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1756942 times.
✓ Branch 4 taken 80500 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 307680 times.
✗ Branch 7 not taken.
 37592848 localAssembler.assembleJacobianAndResidual(*jacobian_, *residual_, *gridVariables_, partialReassembler);
181 37592848 });
182
183 57218 enforcePeriodicConstraints_(*jacobian_, *residual_, curSol, *gridGeometry_);
184 57218 }
185
186 /*!
187 * \brief Assembles only the global Jacobian of the residual.
188 */
189 23 void assembleJacobian(const SolutionVector& curSol)
190 {
191 23 checkAssemblerState_();
192 23 resetJacobian_();
193
194 23 assemble_([&](const Element& element)
195 {
196 181900 LocalAssembler localAssembler(*this, element, curSol);
197
1/2
✓ Branch 1 taken 181900 times.
✗ Branch 2 not taken.
 362664 localAssembler.assembleJacobian(*jacobian_, *gridVariables_);
198 181900 });
199 23 }
200
201 //! compute the residuals using the internal residual
202 1741 void assembleResidual(const SolutionVector& curSol)
203 {
204 1741 resetResidual_();
205 1741 assembleResidual(*residual_, curSol);
206 1741 }
207
208 //! assemble a residual r
209 1741 void assembleResidual(ResidualType& r, const SolutionVector& curSol) const
210 {
211 1741 checkAssemblerState_();
212
213 1741 assemble_([&](const Element& element)
214 {
215
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 2307240 LocalAssembler localAssembler(*this, element, curSol);
216
1/4
✓ Branch 1 taken 2307240 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 2307240 localAssembler.assembleResidual(r);
217 2307240 });
218 }
219
220 /*!
221 * \brief Tells the assembler which jacobian and residual to use.
222 * This also resizes the containers to the required sizes and sets the
223 * sparsity pattern of the jacobian matrix.
224 */
225 27 void setLinearSystem(std::shared_ptr<JacobianMatrix> A,
226 std::shared_ptr<ResidualType> r)
227 {
228 27 jacobian_ = A;
229 27 residual_ = r;
230
231 // check and/or set the BCRS matrix's build mode
232
2/2
✓ Branch 0 taken 24 times.
✓ Branch 1 taken 1 times.
 27 if (jacobian_->buildMode() == JacobianMatrix::BuildMode::unknown)
233 25 jacobian_->setBuildMode(JacobianMatrix::random);
234
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 2 else if (jacobian_->buildMode() != JacobianMatrix::BuildMode::random)
235 ✗ DUNE_THROW(Dune::NotImplemented, "Only BCRS matrices with random build mode are supported at the moment");
236
237 27 setResidualSize_();
238 27 setJacobianPattern_();
239 27 }
240
241 /*!
242 * \brief The version without arguments uses the default constructor to create
243 * the jacobian and residual objects in this assembler if you don't need them outside this class
244 */
245 342 void setLinearSystem()
246 {
247
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 337 times.
 345 jacobian_ = std::make_shared<JacobianMatrix>();
248 342 jacobian_->setBuildMode(JacobianMatrix::random);
249
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 337 times.
 345 residual_ = std::make_shared<ResidualType>();
250
251 342 setResidualSize_();
252 342 setJacobianPattern_();
253 342 }
254
255 /*!
256 * \brief Resizes jacobian and residual and recomputes colors
257 */
258 13 void updateAfterGridAdaption()
259 {
260 13 setResidualSize_();
261 13 setJacobianPattern_();
262 13 maybeComputeColors_();
263 13 }
264
265 //! Returns the number of degrees of freedom
266 792 std::size_t numDofs() const
267
1/3
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 396 { return gridGeometry_->numDofs(); }
268
269 //! The problem
270 508352582 const Problem& problem() const
271
19/20
✓ Branch 1 taken 630385 times.
✓ Branch 2 taken 93236 times.
✓ Branch 4 taken 26590 times.
✓ Branch 5 taken 576 times.
✓ Branch 7 taken 62 times.
✓ Branch 8 taken 38 times.
✓ Branch 11 taken 39 times.
✓ Branch 12 taken 4 times.
✓ Branch 14 taken 3 times.
✓ Branch 15 taken 1 times.
✓ Branch 17 taken 1 times.
✓ Branch 18 taken 1 times.
✓ Branch 3 taken 169 times.
✓ Branch 6 taken 200 times.
✓ Branch 9 taken 68 times.
✓ Branch 10 taken 34 times.
✓ Branch 13 taken 34 times.
✓ Branch 0 taken 167792 times.
✓ Branch 16 taken 1 times.
✗ Branch 19 not taken.
 516079009 { return *problem_; }
272
273 //! The global finite volume geometry
274 132166015 const GridGeometry& gridGeometry() const
275
15/18
✓ Branch 1 taken 791 times.
✓ Branch 2 taken 11584 times.
✓ Branch 13 taken 3200 times.
✗ Branch 14 not taken.
✓ Branch 3 taken 981674 times.
✓ Branch 4 taken 7218357 times.
✓ Branch 9 taken 2260005 times.
✓ Branch 10 taken 1761308 times.
✓ Branch 15 taken 2501 times.
✗ Branch 16 not taken.
✓ Branch 6 taken 478734 times.
✓ Branch 7 taken 12187155 times.
✓ Branch 5 taken 2484278 times.
✓ Branch 11 taken 8175508 times.
✓ Branch 12 taken 1619 times.
✓ Branch 8 taken 333872 times.
✓ Branch 20 taken 76800 times.
✗ Branch 21 not taken.
 90307971 { return *gridGeometry_; }
276
277 //! The gridview
278
6/14
✓ Branch 1 taken 386 times.
✓ Branch 2 taken 4267 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 48993 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 779 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✓ Branch 14 taken 61 times.
✓ Branch 5 taken 20 times.
✗ Branch 10 not taken.
 58545 const GridView& gridView() const
279
13/29
✓ Branch 1 taken 386 times.
✓ Branch 2 taken 4267 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 48993 times.
✓ Branch 6 taken 4413977 times.
✓ Branch 7 taken 8604743 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 26163472 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 307680 times.
✓ Branch 14 taken 779 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✓ Branch 19 taken 335024 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✓ Branch 23 taken 400 times.
✓ Branch 24 taken 61 times.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✓ Branch 29 taken 77672 times.
✗ Branch 10 not taken.
✓ Branch 15 taken 27648 times.
✗ Branch 5 not taken.
✗ Branch 20 not taken.
✗ Branch 25 not taken.
 40035599 { return gridGeometry().gridView(); }
280
281 //! The global grid variables
282 32537 GridVariables& gridVariables()
283 61319 { return *gridVariables_; }
284
285 //! The global grid variables
286 40207027 const GridVariables& gridVariables() const
287
2/4
✓ Branch 1 taken 35887749 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 79300 times.
✗ Branch 5 not taken.
 40207027 { return *gridVariables_; }
288
289 //! The jacobian matrix
290 55835 JacobianMatrix& jacobian()
291
2/3
✓ Branch 2 taken 11 times.
✗ Branch 3 not taken.
✓ Branch 1 taken 18 times.
 60949 { return *jacobian_; }
292
293 //! The residual vector (rhs)
294 70598 ResidualType& residual()
295
7/18
✓ Branch 5 taken 10 times.
✓ Branch 6 taken 31 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 104 times.
✓ Branch 11 taken 10 times.
✗ Branch 12 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✓ Branch 3 taken 13020 times.
✓ Branch 4 taken 18 times.
✗ Branch 7 not taken.
✗ Branch 10 not taken.
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 70598 { return *residual_; }
296
297 //! The solution of the previous time step
298 27070452 const SolutionVector& prevSol() const
299
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 3470 times.
 70206560 { return *prevSol_; }
300
301 /*!
302 * \brief Set time loop for instationary problems
303 * \note calling this turns this into a stationary assembler
304 */
305 void setTimeLoop(std::shared_ptr<const TimeLoop> timeLoop)
306 { timeLoop_ = timeLoop; isStationaryProblem_ = !static_cast<bool>(timeLoop); }
307
308 /*!
309 * \brief Sets the solution from which to start the time integration. Has to be
310 * called prior to assembly for time-dependent problems.
311 */
312 6 void setPreviousSolution(const SolutionVector& u)
313
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 { prevSol_ = &u; }
314
315 /*!
316 * \brief Whether we are assembling a stationary or instationary problem
317 */
318 202945648 bool isStationaryProblem() const
319
9/12
✓ Branch 0 taken 149306017 times.
✓ Branch 1 taken 466717 times.
✓ Branch 2 taken 33405665 times.
✓ Branch 3 taken 3920895 times.
✓ Branch 4 taken 9490848 times.
✓ Branch 5 taken 4239764 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 884972 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 923090 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 307680 times.
 202945648 { return isStationaryProblem_; }
320
321 /*!
322 * \brief Create a local residual object (used by the local assembler)
323 */
324 40081993 LocalResidual localResidual() const
325
2/4
✓ Branch 1 taken 35887753 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 79301 times.
✗ Branch 5 not taken.
 40081993 { return LocalResidual(problem_.get(), timeLoop_.get()); }
326
327 /*!
328 * \brief Update the grid variables
329 */
330 13366 void updateGridVariables(const SolutionVector &cursol)
331 {
332
1/2
✓ Branch 1 taken 950 times.
✗ Branch 2 not taken.
 13811 gridVariables().update(cursol);
333 950 }
334
335 /*!
336 * \brief Reset the gridVariables
337 */
338
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 279 void resetTimeStep(const SolutionVector &cursol)
339 {
340
2/5
✓ Branch 0 taken 248 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 32 times.
✗ Branch 4 not taken.
 280 gridVariables().resetTimeStep(cursol);
341 }
342
343 private:
344 /*!
345 * \brief Resizes the jacobian and sets the jacobian' sparsity pattern.
346 */
347
1/3
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 402 void setJacobianPattern_()
348 {
349 // resize the jacobian and the residual
350 402 const auto numDofs = this->numDofs();
351 402 jacobian_->setSize(numDofs, numDofs);
352
353 // create occupation pattern of the jacobian
354 402 const auto occupationPattern = getJacobianPattern<isImplicit>(gridGeometry());
355
356 // export pattern to jacobian
357
1/2
✓ Branch 1 taken 396 times.
✗ Branch 2 not taken.
 402 occupationPattern.exportIdx(*jacobian_);
358 402 }
359
360 //! Resizes the residual
361 400 void setResidualSize_()
362
1/3
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 406 { residual_->resize(numDofs()); }
363
364 //! Computes the colors
365 399 void maybeComputeColors_()
366 {
367
2/2
✓ Branch 0 taken 386 times.
✓ Branch 1 taken 7 times.
 399 if (enableMultithreading_)
368 782 elementSets_ = computeColoring(gridGeometry()).sets;
369 399 }
370
371 // reset the residual vector to 0.0
372 58911 void resetResidual_()
373 {
374
2/2
✓ Branch 0 taken 18 times.
✓ Branch 1 taken 58119 times.
 58911 if(!residual_)
375 {
376
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18 times.
 18 residual_ = std::make_shared<ResidualType>();
377 18 setResidualSize_();
378 }
379
380 58911 (*residual_) = 0.0;
381 58911 }
382
383 // reset the Jacobian matrix to 0.0
384 template <class PartialReassembler = DefaultPartialReassembler>
385 57241 void resetJacobian_(const PartialReassembler *partialReassembler = nullptr)
386 {
387
2/2
✓ Branch 0 taken 18 times.
✓ Branch 1 taken 56401 times.
 57241 if(!jacobian_)
388 {
389
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18 times.
 18 jacobian_ = std::make_shared<JacobianMatrix>();
390 18 jacobian_->setBuildMode(JacobianMatrix::random);
391 18 setJacobianPattern_();
392 }
393
394
2/2
✓ Branch 0 taken 5837 times.
✓ Branch 1 taken 50582 times.
 57241 if (partialReassembler)
395 5255 partialReassembler->resetJacobian(*this);
396 else
397 51986 *jacobian_ = 0.0;
398 57241 }
399
400 // check if the assembler is in a correct state for assembly
401 58934 void checkAssemblerState_() const
402 {
403
3/4
✓ Branch 0 taken 57804 times.
✓ Branch 1 taken 356 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 57804 times.
 58934 if (!isStationaryProblem_ && !prevSol_)
404 ✗ DUNE_THROW(Dune::InvalidStateException, "Assembling instationary problem but previous solution was not set!");
405 58934 }
406
407 /*!
408 * \brief A method assembling something per element
409 * \note Handles exceptions for parallel runs
410 * \throws NumericalProblem on all processes if an exception is thrown during assembly
411 */
412 template<typename AssembleElementFunc>
413 116310 void assemble_(AssembleElementFunc&& assembleElement) const
414 {
415 // a state that will be checked on all processes
416 116310 bool succeeded = false;
417
418 // try assembling using the local assembly function
419 try
420 {
421
2/2
✓ Branch 0 taken 57775 times.
✓ Branch 1 taken 385 times.
 116310 if (enableMultithreading_)
422 {
423
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 57775 times.
 115541 assert(elementSets_.size() > 0);
424
425 // make this element loop run in parallel
426 // for this we have to color the elements so that we don't get
427 // race conditions when writing into the global matrix
428 // each color can be assembled using multiple threads
429
2/2
✓ Branch 0 taken 261208 times.
✓ Branch 1 taken 57775 times.
 637876 for (const auto& elements : elementSets_)
430 {
431
2/2
✓ Branch 1 taken 1034 times.
✓ Branch 2 taken 260174 times.
 80022317 Dumux::parallelFor(elements.size(), [&](const std::size_t i)
432 {
433
3/8
✗ Branch 0 not taken.
✓ Branch 1 taken 39182172 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 670352 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 78072 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 39970996 const auto element = gridView().grid().entity(elements[i]);
434
3/6
✓ Branch 1 taken 8604743 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 27648 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 400 times.
✗ Branch 8 not taken.
 39970996 assembleElement(element);
435 8632791 });
436 }
437 }
438 else
439
7/13
✓ Branch 1 taken 384 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 110992 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 384 times.
✓ Branch 8 taken 10001 times.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 3 not taken.
✓ Branch 6 taken 100992 times.
✓ Branch 9 taken 10000 times.
 424739 for (const auto& element : elements(gridView()))
440
1/2
✓ Branch 1 taken 110992 times.
✗ Branch 2 not taken.
 211984 assembleElement(element);
441
442 // if we get here, everything worked well on this process
443 107745 succeeded = true;
444 }
445 // throw exception if a problem occurred
446 ✗ catch (NumericalProblem &e)
447 {
448 ✗ std::cout << "rank " << gridView().comm().rank()
449 ✗ << " caught an exception while assembling:" << e.what()
450 ✗ << "\n";
451 ✗ succeeded = false;
452 }
453
454 // make sure everything worked well on all processes
455
4/4
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 55606 times.
✓ Branch 2 taken 7124 times.
✓ Branch 3 taken 45261 times.
 116310 if (gridView().comm().size() > 1)
456 12118 succeeded = gridView().comm().min(succeeded);
457
458 // if not succeeded rethrow the error on all processes
459
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 58160 times.
 116310 if (!succeeded)
460 ✗ DUNE_THROW(NumericalProblem, "A process did not succeed in linearizing the system");
461 116310 }
462
463 template<class GG>
464 17566 void enforcePeriodicConstraints_(JacobianMatrix& jac, ResidualType& res, const SolutionVector& curSol, const GG& gridGeometry)
465 {
466 if constexpr (Detail::hasPeriodicDofMap<GG>())
467 {
468
2/2
✓ Branch 0 taken 202 times.
✓ Branch 1 taken 17566 times.
 17768 for (const auto& m : gridGeometry.periodicDofMap())
469 {
470
2/2
✓ Branch 0 taken 101 times.
✓ Branch 1 taken 101 times.
 202 if (m.first < m.second)
471 {
472 // add the second row to the first
473 101 res[m.first] += res[m.second];
474 101 const auto end = jac[m.second].end();
475
2/2
✓ Branch 0 taken 703 times.
✓ Branch 1 taken 101 times.
 804 for (auto it = jac[m.second].begin(); it != end; ++it)
476 703 jac[m.first][it.index()] += (*it);
477
478 // enforce constraint in second row
479 101 res[m.second] = curSol[m.second] - curSol[m.first];
480
481 // set derivatives accordingly in jacobian, i.e. id for m.second and -id for m.first
482 auto setMatrixBlock = [] (auto& matrixBlock, double diagValue)
483 {
484
0/4
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 202 for (int eIdx = 0; eIdx < matrixBlock.N(); ++eIdx)
485 202 matrixBlock[eIdx][eIdx] = diagValue;
486 };
487
488
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 804 times.
✓ Branch 2 taken 703 times.
✓ Branch 3 taken 101 times.
 804 for (auto it = jac[m.second].begin(); it != end; ++it)
489 {
490 703 auto& matrixBlock = *it;
491
2/2
✓ Branch 0 taken 101 times.
✓ Branch 1 taken 602 times.
 703 matrixBlock = 0.0;
492
493 assert(matrixBlock.N() == matrixBlock.M());
494
2/2
✓ Branch 0 taken 101 times.
✓ Branch 1 taken 602 times.
 703 if(it.index() == m.second)
495 703 setMatrixBlock(matrixBlock, 1.0);
496
497
2/2
✓ Branch 0 taken 101 times.
✓ Branch 1 taken 602 times.
 703 if(it.index() == m.first)
498 703 setMatrixBlock(matrixBlock, -1.0);
499
500 }
501 }
502 }
503 }
504 17566 }
505
506 //! pointer to the problem to be solved
507 std::shared_ptr<const Problem> problem_;
508
509 //! the finite volume geometry of the grid
510 std::shared_ptr<const GridGeometry> gridGeometry_;
511
512 //! the variables container for the grid
513 std::shared_ptr<GridVariables> gridVariables_;
514
515 //! the time loop for instationary problem assembly
516 std::shared_ptr<const TimeLoop> timeLoop_;
517
518 //! an observing pointer to the previous solution for instationary problems
519 const SolutionVector* prevSol_ = nullptr;
520
521 //! if this assembler is assembling an instationary problem
522 bool isStationaryProblem_;
523
524 //! shared pointers to the jacobian matrix and residual
525 std::shared_ptr<JacobianMatrix> jacobian_;
526 std::shared_ptr<ResidualType> residual_;
527
528 //! element sets for parallel assembly
529 bool enableMultithreading_ = false;
530 std::deque<std::vector<ElementSeed>> elementSets_;
531 };
532
533 } // namespace Dumux
534
535 #endif
536