GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/assembly/fvassembler.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 131 149 87.9%
Functions: 1960 3206 61.1%
Branches: 260 545 47.7%
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 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 */
92 template<class TypeTag, DiffMethod diffMethod, bool isImplicit = true>
93 class FVAssembler
94 {
95 using GridGeo = GetPropType<TypeTag, Properties::GridGeometry>;
96 using GridView = typename GridGeo::GridView;
97 using LocalResidual = GetPropType<TypeTag, Properties::LocalResidual>;
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>;
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 120 FVAssembler(std::shared_ptr<const Problem> problem,
124 std::shared_ptr<const GridGeometry> gridGeometry,
125 std::shared_ptr<GridVariables> gridVariables)
126 : problem_(problem)
127 , gridGeometry_(gridGeometry)
128 , gridVariables_(gridVariables)
129 , timeLoop_()
130
4/20
✓ Branch 4 taken 118 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 118 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 118 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 118 times.
✗ 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.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
 120 , isStationaryProblem_(true)
131 {
132 static_assert(isImplicit, "Explicit assembler for stationary problem doesn't make sense!");
133 120 enableMultithreading_ = SupportsColoring<typename GridGeometry::DiscretizationMethod>::value
134
5/8
✗ Branch 0 not taken.
✓ Branch 1 taken 118 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 24 times.
✓ Branch 4 taken 94 times.
✓ Branch 5 taken 24 times.
✓ Branch 7 taken 94 times.
✗ Branch 8 not taken.
 360 && Grid::Capabilities::supportsMultithreading(gridGeometry_->gridView())
135 && !Multithreading::isSerial()
136
5/7
✓ Branch 0 taken 24 times.
✓ Branch 1 taken 94 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 24 times.
✓ Branch 4 taken 94 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 24 times.
 120 && getParam<bool>("Assembly.Multithreading", true);
137
138
1/2
✓ Branch 1 taken 118 times.
✗ Branch 2 not taken.
 120 maybeComputeColors_();
139 120 }
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 234 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 : problem_(problem)
152 , gridGeometry_(gridGeometry)
153 , gridVariables_(gridVariables)
154 , timeLoop_(timeLoop)
155 , prevSol_(&prevSol)
156
3/19
✓ Branch 5 taken 232 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 232 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 232 times.
✗ 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.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
 234 , isStationaryProblem_(!timeLoop)
157 {
158 234 enableMultithreading_ = SupportsColoring<typename GridGeometry::DiscretizationMethod>::value
159
5/8
✗ Branch 0 not taken.
✓ Branch 1 taken 232 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 34 times.
✓ Branch 4 taken 198 times.
✓ Branch 5 taken 34 times.
✓ Branch 7 taken 198 times.
✗ Branch 8 not taken.
 702 && Grid::Capabilities::supportsMultithreading(gridGeometry_->gridView())
160 && !Multithreading::isSerial()
161
5/7
✓ Branch 0 taken 28 times.
✓ Branch 1 taken 201 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 28 times.
✓ Branch 4 taken 195 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 28 times.
 231 && getParam<bool>("Assembly.Multithreading", true);
162
163
1/2
✓ Branch 1 taken 232 times.
✗ Branch 2 not taken.
 234 maybeComputeColors_();
164 234 }
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 56699 void assembleJacobianAndResidual(const SolutionVector& curSol, const PartialReassembler* partialReassembler = nullptr)
172 {
173 56699 checkAssemblerState_();
174 56699 resetJacobian_(partialReassembler);
175 56699 resetResidual_();
176
177 10676843 assemble_([&](const Element& element)
178 {
179
1/2
✓ Branch 1 taken 926074 times.
✗ Branch 2 not taken.
 77407946 LocalAssembler localAssembler(*this, element, curSol);
180
13/22
✓ Branch 1 taken 36927082 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36927082 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 36927082 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 36927082 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 1756942 times.
✓ Branch 13 taken 79300 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 307680 times.
✓ Branch 16 taken 79300 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 307680 times.
✓ Branch 19 taken 79300 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 307680 times.
✓ Branch 22 taken 79300 times.
✗ Branch 23 not taken.
✓ Branch 24 taken 307680 times.
✗ Branch 25 not taken.
 149256248 localAssembler.assembleJacobianAndResidual(*jacobian_, *residual_, *gridVariables_, partialReassembler);
181 });
182
183 226796 enforcePeriodicConstraints_(*jacobian_, *residual_, curSol, *gridGeometry_);
184 56699 }
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 15223 assemble_([&](const Element& element)
195 {
196
1/11
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 363872 LocalAssembler localAssembler(*this, element, curSol);
197
3/15
✓ Branch 1 taken 181900 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 181900 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 181900 times.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ 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.
 726464 localAssembler.assembleJacobian(*jacobian_, *gridVariables_);
198 });
199 23 }
200
201 //! compute the residuals using the internal residual
202 1642 void assembleResidual(const SolutionVector& curSol)
203 {
204 1642 resetResidual_();
205 1642 assembleResidual(*residual_, curSol);
206 1642 }
207
208 //! assemble a residual r
209 void assembleResidual(ResidualType& r, const SolutionVector& curSol) const
210 {
211 1642 checkAssemblerState_();
212
213 1642 assemble_([&](const Element& element)
214 {
215
1/6
✓ Branch 1 taken 1003600 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 6538880 LocalAssembler localAssembler(*this, element, curSol);
216
1/4
✓ Branch 1 taken 2267640 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 2267640 localAssembler.assembleResidual(r);
217 });
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 14 void setLinearSystem(std::shared_ptr<JacobianMatrix> A,
226 std::shared_ptr<ResidualType> r)
227 {
228 14 jacobian_ = A;
229 14 residual_ = r;
230
231 // check and/or set the BCRS matrix's build mode
232
4/4
✓ Branch 0 taken 11 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 11 times.
✓ Branch 3 taken 1 times.
 28 if (jacobian_->buildMode() == JacobianMatrix::BuildMode::unknown)
233 24 jacobian_->setBuildMode(JacobianMatrix::random);
234
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
 4 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 21 setResidualSize_();
238 14 setJacobianPattern_();
239 14 }
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 325 void setLinearSystem()
246 {
247
3/4
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 320 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 323 times.
 328 jacobian_ = std::make_shared<JacobianMatrix>();
248 650 jacobian_->setBuildMode(JacobianMatrix::random);
249
3/4
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 320 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 323 times.
 328 residual_ = std::make_shared<ResidualType>();
250
251 462 setResidualSize_();
252 325 setJacobianPattern_();
253 325 }
254
255 /*!
256 * \brief Resizes jacobian and residual and recomputes colors
257 */
258 13 void updateAfterGridAdaption()
259 {
260 21 setResidualSize_();
261 13 setJacobianPattern_();
262 13 maybeComputeColors_();
263 13 }
264
265 //! Returns the number of degrees of freedom
266 std::size_t numDofs() const
267
2/6
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
 732 { return gridGeometry_->numDofs(); }
268
269 //! The problem
270 const Problem& problem() const
271
40/48
✓ Branch 0 taken 175697 times.
✓ Branch 1 taken 494353 times.
✓ Branch 2 taken 175697 times.
✓ Branch 3 taken 506696 times.
✓ Branch 4 taken 12994 times.
✓ Branch 5 taken 86353 times.
✓ Branch 6 taken 13433 times.
✓ Branch 7 taken 87409 times.
✓ Branch 8 taken 176 times.
✓ Branch 9 taken 160904 times.
✓ Branch 10 taken 240 times.
✓ Branch 11 taken 147934 times.
✓ Branch 12 taken 116 times.
✓ Branch 13 taken 37 times.
✓ Branch 14 taken 42 times.
✓ Branch 15 taken 82 times.
✓ Branch 16 taken 37 times.
✓ Branch 17 taken 42 times.
✓ Branch 18 taken 82 times.
✓ Branch 19 taken 22 times.
✓ Branch 20 taken 42 times.
✓ Branch 21 taken 73 times.
✓ Branch 22 taken 22 times.
✗ Branch 23 not taken.
✓ Branch 24 taken 73 times.
✓ Branch 25 taken 19 times.
✗ Branch 26 not taken.
✓ Branch 27 taken 3 times.
✓ Branch 28 taken 19 times.
✓ Branch 29 taken 1 times.
✓ Branch 30 taken 3 times.
✓ Branch 31 taken 4 times.
✓ Branch 32 taken 1 times.
✓ Branch 33 taken 65904 times.
✓ Branch 34 taken 4 times.
✓ Branch 35 taken 1 times.
✓ Branch 36 taken 65904 times.
✓ Branch 37 taken 3600 times.
✓ Branch 38 taken 1 times.
✗ Branch 39 not taken.
✓ Branch 40 taken 3600 times.
✓ Branch 41 taken 1 times.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✓ Branch 44 taken 1 times.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
 1027488262 { return *problem_; }
272
273 //! The global finite volume geometry
274 const GridGeometry& gridGeometry() const
275
23/31
✗ Branch 0 not taken.
✓ Branch 1 taken 7362144 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1115592 times.
✓ Branch 4 taken 7356096 times.
✓ Branch 5 taken 3773 times.
✓ Branch 6 taken 1106074 times.
✓ Branch 7 taken 15353996 times.
✓ Branch 8 taken 303 times.
✓ Branch 9 taken 1111878 times.
✓ Branch 10 taken 15353996 times.
✓ Branch 11 taken 1250005 times.
✓ Branch 12 taken 1111878 times.
✓ Branch 13 taken 1762399 times.
✓ Branch 14 taken 1250002 times.
✓ Branch 15 taken 8175508 times.
✓ Branch 16 taken 1762396 times.
✓ Branch 17 taken 18 times.
✓ Branch 18 taken 8175508 times.
✓ Branch 19 taken 2518 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 1 times.
✓ Branch 22 taken 2500 times.
✓ Branch 23 taken 1 times.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✓ Branch 27 taken 76800 times.
✗ Branch 28 not taken.
✓ Branch 30 taken 76800 times.
✗ Branch 31 not taken.
 183922340 { return *gridGeometry_; }
276
277 //! The gridview
278 const GridView& gridView() const
279
12/29
✓ Branch 1 taken 1095 times.
✓ Branch 2 taken 4266 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 48743 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 4413977 times.
✓ Branch 7 taken 8227075 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 25044212 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 307680 times.
✓ Branch 14 taken 680 times.
✓ Branch 15 taken 27648 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✓ Branch 19 taken 295424 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✓ Branch 24 taken 53 times.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✓ Branch 29 taken 76872 times.
 38497965 { return gridGeometry().gridView(); }
280
281 //! The global grid variables
282 GridVariables& gridVariables()
283 121712 { return *gridVariables_; }
284
285 //! The global grid variables
286 const GridVariables& gridVariables() const
287
12/24
✓ Branch 1 taken 36115495 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36115495 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 36115495 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 36115495 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 36115495 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 36115495 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 79300 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 79300 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 79300 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 79300 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 79300 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 79300 times.
✗ Branch 35 not taken.
 238581612 { return *gridVariables_; }
288
289 //! The jacobian matrix
290 JacobianMatrix& jacobian()
291
4/7
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 18 times.
✓ Branch 5 taken 11 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 11 times.
✗ Branch 9 not taken.
 231672 { return *jacobian_; }
292
293 //! The residual vector (rhs)
294 ResidualType& residual()
295
14/35
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 18 times.
✓ Branch 5 taken 12888 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 18 times.
✓ Branch 8 taken 12888 times.
✓ Branch 9 taken 10 times.
✓ Branch 10 taken 18 times.
✓ Branch 11 taken 31 times.
✓ Branch 12 taken 10 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 31 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 104 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✓ Branch 20 taken 104 times.
✓ Branch 21 taken 10 times.
✗ Branch 22 not taken.
✓ Branch 24 taken 10 times.
✗ Branch 25 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
 139818 { return *residual_; }
296
297 //! The solution of the previous time step
298 ✗ const SolutionVector& prevSol() const
299 ✗ { 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 void setPreviousSolution(const SolutionVector& u)
313 { prevSol_ = &u; }
314
315 /*!
316 * \brief Whether we are assembling a stationary or instationary problem
317 */
318 ✗ bool isStationaryProblem() const
319 ✗ { return isStationaryProblem_; }
320
321 /*!
322 * \brief Create a local residual object (used by the local assembler)
323 */
324 LocalResidual localResidual() const
325
4/8
✓ Branch 1 taken 36115499 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36115499 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 79301 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 79301 times.
✗ Branch 11 not taken.
 79527214 { return LocalResidual(problem_.get(), timeLoop_.get()); }
326
327 /*!
328 * \brief Update the grid variables
329 */
330 ✗ void updateGridVariables(const SolutionVector &cursol)
331 {
332
2/4
✓ Branch 1 taken 950 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 950 times.
✗ Branch 5 not taken.
 155666 gridVariables().update(cursol);
333 ✗ }
334
335 /*!
336 * \brief Reset the gridVariables
337 */
338 ✗ void resetTimeStep(const SolutionVector &cursol)
339 {
340
2/8
✓ Branch 0 taken 270 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 270 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 604 gridVariables().resetTimeStep(cursol);
341 ✗ }
342
343 private:
344 /*!
345 * \brief Resizes the jacobian and sets the jacobian' sparsity pattern.
346 */
347 372 void setJacobianPattern_()
348 {
349 // resize the jacobian and the residual
350
1/3
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 372 const auto numDofs = this->numDofs();
351 744 jacobian_->setSize(numDofs, numDofs);
352
353 // create occupation pattern of the jacobian
354 1116 const auto occupationPattern = getJacobianPattern<isImplicit>(gridGeometry());
355
356 // export pattern to jacobian
357
2/4
✓ Branch 1 taken 366 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 366 times.
✗ Branch 5 not taken.
 744 occupationPattern.exportIdx(*jacobian_);
358 372 }
359
360 //! Resizes the residual
361 212 void setResidualSize_()
362
2/6
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
 896 { residual_->resize(numDofs()); }
363
364 //! Computes the colors
365 369 void maybeComputeColors_()
366 {
367
2/2
✓ Branch 0 taken 354 times.
✓ Branch 1 taken 9 times.
 369 if (enableMultithreading_)
368 720 elementSets_ = computeColoring(gridGeometry()).sets;
369 369 }
370
371 // reset the residual vector to 0.0
372 58341 void resetResidual_()
373 {
374
2/2
✓ Branch 0 taken 18 times.
✓ Branch 1 taken 57549 times.
 58341 if(!residual_)
375 {
376
2/4
✗ Branch 1 not taken.
✓ Branch 2 taken 18 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 18 times.
 18 residual_ = std::make_shared<ResidualType>();
377 18 setResidualSize_();
378 }
379
380 116682 (*residual_) = 0.0;
381 58341 }
382
383 // reset the Jacobian matrix to 0.0
384 template <class PartialReassembler = DefaultPartialReassembler>
385 56722 void resetJacobian_(const PartialReassembler *partialReassembler = nullptr)
386 {
387
2/2
✓ Branch 0 taken 18 times.
✓ Branch 1 taken 55930 times.
 56722 if(!jacobian_)
388 {
389
2/4
✗ Branch 1 not taken.
✓ Branch 2 taken 18 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 18 times.
 18 jacobian_ = std::make_shared<JacobianMatrix>();
390 36 jacobian_->setBuildMode(JacobianMatrix::random);
391 18 setJacobianPattern_();
392 }
393
394
2/2
✓ Branch 0 taken 5810 times.
✓ Branch 1 taken 50138 times.
 56722 if (partialReassembler)
395 5241 partialReassembler->resetJacobian(*this);
396 else
397 102963 *jacobian_ = 0.0;
398 56722 }
399
400 // check if the assembler is in a correct state for assembly
401 58364 void checkAssemblerState_() const
402 {
403
3/4
✓ Branch 0 taken 57251 times.
✓ Branch 1 taken 339 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 57251 times.
 58364 if (!isStationaryProblem_ && !prevSol_)
404 ✗ DUNE_THROW(Dune::InvalidStateException, "Assembling instationary problem but previous solution was not set!");
405 58364 }
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 115171 void assemble_(AssembleElementFunc&& assembleElement) const
414 {
415 // a state that will be checked on all processes
416 115171 bool succeeded = false;
417
418 // try assembling using the local assembly function
419 try
420 {
421
2/2
✓ Branch 0 taken 56497 times.
✓ Branch 1 taken 1093 times.
 115171 if (enableMultithreading_)
422 {
423
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 56497 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 56497 times.
 225970 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
4/4
✓ Branch 0 taken 249835 times.
✓ Branch 1 taken 56497 times.
✓ Branch 2 taken 249835 times.
✓ Branch 3 taken 56497 times.
 1225148 for (const auto& elements : elementSets_)
430 {
431
2/2
✓ Branch 2 taken 248895 times.
✓ Branch 3 taken 940 times.
 77423704 Dumux::parallelFor(elements.size(), [&](const std::size_t i)
432 {
433
7/18
✗ Branch 0 not taken.
✓ Branch 1 taken 37685244 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 37685244 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 603104 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 630752 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 27648 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 76872 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 76872 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
 85202059 const auto element = gridView().grid().entity(elements[i]);
434
2/4
✓ Branch 1 taken 8227075 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 27648 times.
✗ Branch 5 not taken.
 38433268 assembleElement(element);
435 });
436 }
437 }
438 else
439
8/19
✓ Branch 1 taken 1093 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1093 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1330334 times.
✓ Branch 7 taken 1093 times.
✓ Branch 8 taken 1330334 times.
✓ Branch 9 taken 1093 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 1330334 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✓ Branch 14 taken 1330334 times.
✗ Branch 15 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
 2665040 for (const auto& element : elements(gridView()))
440
1/2
✓ Branch 1 taken 1330334 times.
✗ Branch 2 not taken.
 2660668 assembleElement(element);
441
442 // if we get here, everything worked well on this process
443 106607 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
8/8
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 53762 times.
✓ Branch 2 taken 1488 times.
✓ Branch 3 taken 56102 times.
✓ Branch 4 taken 6016 times.
✓ Branch 5 taken 47291 times.
✓ Branch 6 taken 4530 times.
✓ Branch 7 taken 44951 times.
 230352 if (gridView().comm().size() > 1)
456 36090 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 57590 times.
 115171 if (!succeeded)
460 ✗ DUNE_THROW(NumericalProblem, "A process did not succeed in linearizing the system");
461 115171 }
462
463 template<class GG>
464 17061 void enforcePeriodicConstraints_(JacobianMatrix& jac, ResidualType& res, const SolutionVector& curSol, const GG& gridGeometry)
465 {
466 if constexpr (Detail::hasPeriodicDofMap<GG>())
467 {
468
4/4
✓ Branch 0 taken 202 times.
✓ Branch 1 taken 17061 times.
✓ Branch 2 taken 101 times.
✓ Branch 3 taken 101 times.
 68446 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 303 res[m.first] += res[m.second];
474 202 const auto end = jac[m.second].end();
475
4/4
✓ Branch 0 taken 703 times.
✓ Branch 1 taken 101 times.
✓ Branch 2 taken 703 times.
✓ Branch 3 taken 101 times.
 1507 for (auto it = jac[m.second].begin(); it != end; ++it)
476 3515 jac[m.first][it.index()] += (*it);
477
478 // enforce constraint in second row
479 404 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 ✗ for (int eIdx = 0; eIdx < matrixBlock.N(); ++eIdx)
485 404 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.
 905 for (auto it = jac[m.second].begin(); it != end; ++it)
489 {
490
2/2
✓ Branch 0 taken 101 times.
✓ Branch 1 taken 602 times.
 703 auto& matrixBlock = *it;
491
2/2
✓ Branch 0 taken 101 times.
✓ Branch 1 taken 602 times.
 703 matrixBlock = 0.0;
492
493
4/4
✓ Branch 0 taken 101 times.
✓ Branch 1 taken 602 times.
✓ Branch 2 taken 101 times.
✓ Branch 3 taken 602 times.
 1406 assert(matrixBlock.N() == matrixBlock.M());
494
4/4
✓ Branch 0 taken 101 times.
✓ Branch 1 taken 602 times.
✓ Branch 2 taken 101 times.
✓ Branch 3 taken 602 times.
 1406 if(it.index() == m.second)
495 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 setMatrixBlock(matrixBlock, -1.0);
499
500 }
501 }
502 }
503 }
504 17061 }
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