GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/multidomain/fvassembler.hh
Date: 2025-07-12 19:18:49
Exec Total Coverage
Lines: 200 205 97.6%
Functions: 2768 3287 84.2%
Branches: 303 675 44.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 MultiDomain
10 * \ingroup Assembly
11 * \brief A linear system assembler (residual and Jacobian) for finite volume schemes
12 * with multiple domains
13 */
14 #ifndef DUMUX_MULTIDOMAIN_FV_ASSEMBLER_HH
15 #define DUMUX_MULTIDOMAIN_FV_ASSEMBLER_HH
16
17 #include <type_traits>
18 #include <tuple>
19
20 #include <dune/common/hybridutilities.hh>
21 #include <dune/istl/matrixindexset.hh>
22
23 #include <dumux/common/exceptions.hh>
24 #include <dumux/common/properties.hh>
25 #include <dumux/common/timeloop.hh>
26 #include <dumux/common/typetraits/utility.hh>
27 #include <dumux/common/typetraits/periodic.hh>
28 #include <dumux/common/gridcapabilities.hh>
29 #include <dumux/discretization/method.hh>
30 #include <dumux/assembly/diffmethod.hh>
31 #include <dumux/assembly/jacobianpattern.hh>
32 #include <dumux/linear/parallelhelpers.hh>
33 #include <dumux/parallel/multithreading.hh>
34
35 #include "couplingjacobianpattern.hh"
36 #include "subdomaincclocalassembler.hh"
37 #include "subdomaincvfelocalassembler.hh"
38 #include "subdomainstaggeredlocalassembler.hh"
39 #include "subdomainfclocalassembler.hh"
40 #include "assemblerview.hh"
41
42 #include <dumux/discretization/method.hh>
43
44 namespace Dumux {
45
46 namespace Grid::Capabilities {
47
48 namespace Detail {
49 // helper for multi-domain models
50 template<class T, std::size_t... I>
51
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 bool allGridsSupportsMultithreadingImpl(const T& gridGeometries, std::index_sequence<I...>)
52 {
53
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 26 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 26 return (... && supportsMultithreading(std::get<I>(gridGeometries)->gridView()));
54 }
55 } // end namespace Detail
56
57 // helper for multi-domain models (all grids have to support multithreading)
58 template<class... GG>
59 26 bool allGridsSupportsMultithreading(const std::tuple<GG...>& gridGeometries)
60 {
61 1 return Detail::allGridsSupportsMultithreadingImpl<std::tuple<GG...>>(gridGeometries, std::make_index_sequence<sizeof...(GG)>());
62 }
63
64 } // end namespace Grid::Capabilities
65
66 namespace Detail {
67
68 //! helper struct detecting if sub-problem has a constraints() function
69 template<class P>
70 using SubProblemConstraintsDetector = decltype(std::declval<P>().constraints());
71
72 template<class P>
73 constexpr inline bool hasSubProblemGlobalConstraints()
74 { return Dune::Std::is_detected<SubProblemConstraintsDetector, P>::value; }
75
76 } // end namespace Detail
77
78 /*!
79 * \ingroup MultiDomain
80 * \ingroup Assembly
81 * \brief trait that is specialized for coupling manager supporting multithreaded assembly
82 */
83 template<class CM>
84 struct CouplingManagerSupportsMultithreadedAssembly : public std::false_type
85 {};
86
87 /*!
88 * \ingroup MultiDomain
89 * \ingroup Assembly
90 * \brief A linear system assembler (residual and Jacobian) for finite volume schemes (box, tpfa, mpfa, ...)
91 * with multiple domains
92 * \tparam MDTraits the multidimensional traits
93 * \tparam diffMethod the differentiation method to residual compute derivatives
94 * \tparam useImplicitAssembly if to use an implicit or explicit time discretization
95 */
96 template<class MDTraits, class CMType, DiffMethod diffMethod, bool useImplicitAssembly = true>
97 class MultiDomainFVAssembler
98 {
99 template<std::size_t id>
100 using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
101
102 public:
103 using Traits = MDTraits;
104
105 using Scalar = typename MDTraits::Scalar;
106
107 template<std::size_t id>
108 using LocalResidual = typename MDTraits::template SubDomain<id>::LocalResidual;
109
110 template<std::size_t id>
111 using GridVariables = typename MDTraits::template SubDomain<id>::GridVariables;
112
113 template<std::size_t id>
114 using GridGeometry = typename MDTraits::template SubDomain<id>::GridGeometry;
115
116 template<std::size_t id>
117 using Problem = typename MDTraits::template SubDomain<id>::Problem;
118
119 using JacobianMatrix = typename MDTraits::JacobianMatrix;
120 using SolutionVector = typename MDTraits::SolutionVector;
121 using ResidualType = typename MDTraits::ResidualVector;
122
123 using CouplingManager = CMType;
124
125 /*!
126 * \brief Returns true if the assembler considers implicit assembly.
127 */
128 static constexpr bool isImplicit()
129 { return useImplicitAssembly; }
130
131 private:
132
133 using ProblemTuple = typename MDTraits::template TupleOfSharedPtrConst<Problem>;
134 using GridGeometryTuple = typename MDTraits::template TupleOfSharedPtrConst<GridGeometry>;
135 using GridVariablesTuple = typename MDTraits::template TupleOfSharedPtr<GridVariables>;
136
137 using TimeLoop = TimeLoopBase<Scalar>;
138 using ThisType = MultiDomainFVAssembler<MDTraits, CouplingManager, diffMethod, isImplicit()>;
139
140 template<std::size_t id>
141 using SubDomainAssemblerView = MultiDomainAssemblerSubDomainView<ThisType, id>;
142
143 template<class DiscretizationMethod, std::size_t id>
144 struct SubDomainAssemblerType;
145
146 template<std::size_t id>
147 struct SubDomainAssemblerType<DiscretizationMethods::CCTpfa, id>
148 {
149 using type = SubDomainCCLocalAssembler<id, SubDomainTypeTag<id>, SubDomainAssemblerView<id>, diffMethod, isImplicit()>;
150 };
151
152 template<std::size_t id>
153 struct SubDomainAssemblerType<DiscretizationMethods::CCMpfa, id>
154 {
155 using type = SubDomainCCLocalAssembler<id, SubDomainTypeTag<id>, SubDomainAssemblerView<id>, diffMethod, isImplicit()>;
156 };
157
158 template<std::size_t id, class DM>
159 struct SubDomainAssemblerType<DiscretizationMethods::CVFE<DM>, id>
160 {
161 using type = SubDomainCVFELocalAssembler<id, SubDomainTypeTag<id>, SubDomainAssemblerView<id>, diffMethod, isImplicit()>;
162 };
163
164 template<std::size_t id>
165 struct SubDomainAssemblerType<DiscretizationMethods::Staggered, id>
166 {
167 using type = SubDomainStaggeredLocalAssembler<id, SubDomainTypeTag<id>, SubDomainAssemblerView<id>, diffMethod, isImplicit()>;
168 };
169
170 template<std::size_t id>
171 struct SubDomainAssemblerType<DiscretizationMethods::FCStaggered, id>
172 {
173 using type = SubDomainFaceCenteredLocalAssembler<id, SubDomainTypeTag<id>, SubDomainAssemblerView<id>, diffMethod, isImplicit()>;
174 };
175
176 template<std::size_t id>
177 using SubDomainAssembler = typename SubDomainAssemblerType<typename GridGeometry<id>::DiscretizationMethod, id>::type;
178
179 public:
180
181
182 /*!
183 * \brief The constructor for stationary problems
184 * \note the grid variables might be temporarily changed during assembly (if caching is enabled)
185 * it is however guaranteed that the state after assembly will be the same as before
186 */
187 129 MultiDomainFVAssembler(ProblemTuple problem,
188 GridGeometryTuple gridGeometry,
189 GridVariablesTuple gridVariables,
190 std::shared_ptr<CouplingManager> couplingManager)
191 129 : couplingManager_(couplingManager)
192
1/2
✓ Branch 1 taken 129 times.
✗ Branch 2 not taken.
 129 , problemTuple_(std::move(problem))
193 129 , gridGeometryTuple_(std::move(gridGeometry))
194 129 , gridVariablesTuple_(std::move(gridVariables))
195
1/2
✓ Branch 1 taken 129 times.
✗ Branch 2 not taken.
 129 , timeLoop_()
196 129 , isStationaryProblem_(true)
197
1/2
✓ Branch 1 taken 129 times.
✗ Branch 2 not taken.
 129 , warningIssued_(false)
198 {
199 static_assert(isImplicit(), "Explicit assembler for stationary problem doesn't make sense!");
200
1/2
✓ Branch 1 taken 129 times.
✗ Branch 2 not taken.
 129 std::cout << "Instantiated assembler for a stationary problem." << std::endl;
201
202 129 enableMultithreading_ = CouplingManagerSupportsMultithreadedAssembly<CouplingManager>::value
203
1/2
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
 21 && Grid::Capabilities::allGridsSupportsMultithreading(gridGeometryTuple_)
204 && !Multithreading::isSerial()
205
1/2
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
 21 && getParam<bool>("Assembly.Multithreading", true);
206
207 129 maybeComputeColors_();
208 21 }
209
210 /*!
211 * \brief The constructor for instationary problems
212 * \note the grid variables might be temporarily changed during assembly (if caching is enabled)
213 * it is however guaranteed that the state after assembly will be the same as before
214 */
215 92 MultiDomainFVAssembler(ProblemTuple problem,
216 GridGeometryTuple gridGeometry,
217 GridVariablesTuple gridVariables,
218 std::shared_ptr<CouplingManager> couplingManager,
219 std::shared_ptr<const TimeLoop> timeLoop,
220 const SolutionVector& prevSol)
221 92 : couplingManager_(couplingManager)
222 92 , problemTuple_(std::move(problem))
223 92 , gridGeometryTuple_(std::move(gridGeometry))
224 92 , gridVariablesTuple_(std::move(gridVariables))
225 92 , timeLoop_(timeLoop)
226 92 , prevSol_(&prevSol)
227 92 , isStationaryProblem_(false)
228
1/2
✓ Branch 1 taken 92 times.
✗ Branch 2 not taken.
 92 , warningIssued_(false)
229 {
230
1/2
✓ Branch 1 taken 92 times.
✗ Branch 2 not taken.
 92 std::cout << "Instantiated assembler for an instationary problem." << std::endl;
231
232 92 enableMultithreading_ = CouplingManagerSupportsMultithreadedAssembly<CouplingManager>::value
233
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 5 && Grid::Capabilities::allGridsSupportsMultithreading(gridGeometryTuple_)
234 && !Multithreading::isSerial()
235
4/7
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 2 not taken.
 5 && getParam<bool>("Assembly.Multithreading", true);
236
237 92 maybeComputeColors_();
238 5 }
239
240 /*!
241 * \brief Assembles the global Jacobian of the residual
242 * and the residual for the current solution.
243 */
244 10582 void assembleJacobianAndResidual(const SolutionVector& curSol)
245 {
246 10582 checkAssemblerState_();
247 10582 resetJacobian_();
248 10582 resetResidual_();
249
250 using namespace Dune::Hybrid;
251 80786 forEach(std::make_index_sequence<JacobianMatrix::N()>(), [&](const auto domainId)
252 {
253 19061 auto& jacRow = (*jacobian_)[domainId];
254 19061 auto& subRes = (*residual_)[domainId];
255 19061 const auto& subSol = curSol[domainId];
256 38122 this->assembleJacobianAndResidual_(domainId, jacRow, subRes, curSol);
257
258 38122 const auto gridGeometry = std::get<domainId>(gridGeometryTuple_);
259
8/10
✓ Branch 1 taken 287 times.
✓ Branch 2 taken 6609 times.
✓ Branch 3 taken 945 times.
✓ Branch 4 taken 1694 times.
✓ Branch 6 taken 1500 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 500 times.
✗ Branch 9 not taken.
✓ Branch 0 taken 7285 times.
✓ Branch 5 taken 241 times.
 19061 enforcePeriodicConstraints_(domainId, jacRow, subRes, *gridGeometry, subSol);
260
4/8
✓ Branch 0 taken 7572 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 7572 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2417 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1500 times.
✗ Branch 7 not taken.
 19061 enforceGlobalDirichletConstraints_(domainId, jacRow, subRes, *gridGeometry, subSol);
261 17138 });
262 10582 }
263
264 //! compute the residuals using the internal residual
265 1510 void assembleResidual(const SolutionVector& curSol)
266 {
267
1/8
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 1 taken 99 times.
✗ Branch 2 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 1510 resetResidual_();
268
1/8
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 1 taken 99 times.
✗ Branch 2 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 1510 assembleResidual(*residual_, curSol);
269 99 }
270
271 //! assemble a residual r
272 1510 void assembleResidual(ResidualType& r, const SolutionVector& curSol)
273 {
274 1510 r = 0.0;
275
276 1510 checkAssemblerState_();
277
278 // update the grid variables for the case of active caching
279 1510 updateGridVariables(curSol);
280
281 using namespace Dune::Hybrid;
282 6038 forEach(integralRange(Dune::Hybrid::size(r)), [&](const auto domainId)
283 {
284 1510 auto& subRes = r[domainId];
285 4528 this->assembleResidual_(domainId, subRes, curSol);
286 });
287 1509 }
288
289 /*!
290 * \brief Tells the assembler which jacobian and residual to use.
291 * This also resizes the containers to the required sizes and sets the
292 * sparsity pattern of the jacobian matrix.
293 */
294 3 void setLinearSystem(std::shared_ptr<JacobianMatrix> A,
295 std::shared_ptr<ResidualType> r)
296 {
297 3 jacobian_ = A;
298 3 residual_ = r;
299
300 3 setJacobianBuildMode(*jacobian_);
301 3 setJacobianPattern_(*jacobian_);
302 3 setResidualSize_(*residual_);
303 3 }
304
305 /*!
306 * \brief The version without arguments uses the default constructor to create
307 * the jacobian and residual objects in this assembler if you don't need them outside this class
308 */
309 204 void setLinearSystem()
310 {
311
2/2
✓ Branch 0 taken 7 times.
✓ Branch 1 taken 197 times.
 211 jacobian_ = std::make_shared<JacobianMatrix>();
312
2/2
✓ Branch 0 taken 7 times.
✓ Branch 1 taken 197 times.
 211 residual_ = std::make_shared<ResidualType>();
313
314 204 setJacobianBuildMode(*jacobian_);
315 244 setJacobianPattern_(*jacobian_);
316 244 setResidualSize_(*residual_);
317 204 }
318
319 /*!
320 * \brief Sets the jacobian build mode
321 */
322 228 void setJacobianBuildMode(JacobianMatrix& jac) const
323 {
324 using namespace Dune::Hybrid;
325 456 forEach(std::make_index_sequence<JacobianMatrix::N()>(), [&](const auto i)
326 {
327
1/2
✓ Branch 0 taken 1112 times.
✗ Branch 1 not taken.
 2640 forEach(jac[i], [&](auto& jacBlock)
328 {
329 using BlockType = std::decay_t<decltype(jacBlock)>;
330
6/12
✓ Branch 0 taken 422 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 411 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 166 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 74 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 25 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 14 times.
✗ Branch 11 not taken.
 1112 if (jacBlock.buildMode() == BlockType::BuildMode::unknown)
331 1112 jacBlock.setBuildMode(BlockType::BuildMode::random);
332 ✗ else if (jacBlock.buildMode() != BlockType::BuildMode::random)
333
0/120
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
✗ Branch 47 not taken.
✗ Branch 48 not taken.
✗ Branch 50 not taken.
✗ Branch 51 not taken.
✗ Branch 53 not taken.
✗ Branch 54 not taken.
✗ Branch 56 not taken.
✗ Branch 57 not taken.
✗ Branch 59 not taken.
✗ Branch 60 not taken.
✗ Branch 62 not taken.
✗ Branch 63 not taken.
✗ Branch 69 not taken.
✗ Branch 70 not taken.
✗ Branch 72 not taken.
✗ Branch 73 not taken.
✗ Branch 75 not taken.
✗ Branch 76 not taken.
✗ Branch 78 not taken.
✗ Branch 79 not taken.
✗ Branch 81 not taken.
✗ Branch 82 not taken.
✗ Branch 84 not taken.
✗ Branch 85 not taken.
✗ Branch 87 not taken.
✗ Branch 88 not taken.
✗ Branch 90 not taken.
✗ Branch 91 not taken.
✗ Branch 93 not taken.
✗ Branch 94 not taken.
✗ Branch 96 not taken.
✗ Branch 97 not taken.
✗ Branch 103 not taken.
✗ Branch 104 not taken.
✗ Branch 106 not taken.
✗ Branch 107 not taken.
✗ Branch 109 not taken.
✗ Branch 110 not taken.
✗ Branch 112 not taken.
✗ Branch 113 not taken.
✗ Branch 115 not taken.
✗ Branch 116 not taken.
✗ Branch 118 not taken.
✗ Branch 119 not taken.
✗ Branch 121 not taken.
✗ Branch 122 not taken.
✗ Branch 124 not taken.
✗ Branch 125 not taken.
✗ Branch 127 not taken.
✗ Branch 128 not taken.
✗ Branch 130 not taken.
✗ Branch 131 not taken.
✗ Branch 137 not taken.
✗ Branch 138 not taken.
✗ Branch 140 not taken.
✗ Branch 141 not taken.
✗ Branch 143 not taken.
✗ Branch 144 not taken.
✗ Branch 146 not taken.
✗ Branch 147 not taken.
✗ Branch 149 not taken.
✗ Branch 150 not taken.
✗ Branch 152 not taken.
✗ Branch 153 not taken.
✗ Branch 155 not taken.
✗ Branch 156 not taken.
✗ Branch 158 not taken.
✗ Branch 159 not taken.
✗ Branch 161 not taken.
✗ Branch 162 not taken.
✗ Branch 164 not taken.
✗ Branch 165 not taken.
✗ Branch 171 not taken.
✗ Branch 172 not taken.
✗ Branch 174 not taken.
✗ Branch 175 not taken.
✗ Branch 177 not taken.
✗ Branch 178 not taken.
✗ Branch 180 not taken.
✗ Branch 181 not taken.
✗ Branch 183 not taken.
✗ Branch 184 not taken.
✗ Branch 186 not taken.
✗ Branch 187 not taken.
✗ Branch 189 not taken.
✗ Branch 190 not taken.
✗ Branch 192 not taken.
✗ Branch 193 not taken.
✗ Branch 195 not taken.
✗ Branch 196 not taken.
✗ Branch 198 not taken.
✗ Branch 199 not taken.
 1112 DUNE_THROW(Dune::NotImplemented, "Only BCRS matrices with random build mode are supported at the moment");
334 });
335 });
336 }
337
338 /*!
339 * \brief Resizes jacobian and residual and recomputes colors
340 */
341 void updateAfterGridAdaption()
342 {
343 setJacobianPattern_(*jacobian_);
344 setResidualSize_(*residual_);
345 maybeComputeColors_();
346 }
347
348 /*!
349 * \brief Updates the grid variables with the given solution
350 */
351 10222 void updateGridVariables(const SolutionVector& curSol)
352 {
353 using namespace Dune::Hybrid;
354
3/6
✓ Branch 4 taken 1058 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1057 times.
✗ Branch 8 not taken.
✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
 9285 forEach(integralRange(Dune::Hybrid::size(gridVariablesTuple_)), [&](const auto domainId)
355 13257 { this->gridVariables(domainId).update(curSol[domainId]); });
356 }
357
358 /*!
359 * \brief Resets the grid variables to the last time step
360 */
361 19 void resetTimeStep(const SolutionVector& curSol)
362 {
363 using namespace Dune::Hybrid;
364
3/7
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✗ Branch 0 not taken.
 19 forEach(integralRange(Dune::Hybrid::size(gridVariablesTuple_)), [&](const auto domainId)
365 23 { this->gridVariables(domainId).resetTimeStep(curSol[domainId]); });
366 }
367
368 //! the number of dof locations of domain i
369 template<std::size_t i>
370 496 std::size_t numDofs(Dune::index_constant<i> domainId) const
371 465 { return std::get<domainId>(gridGeometryTuple_)->numDofs(); }
372
373 //! the problem of domain i
374 template<std::size_t i>
375 283270316 const auto& problem(Dune::index_constant<i> domainId) const
376
49/53
✓ Branch 0 taken 531280 times.
✓ Branch 1 taken 52060359 times.
✓ Branch 2 taken 21006 times.
✓ Branch 3 taken 10397990 times.
✓ Branch 4 taken 6851916 times.
✓ Branch 5 taken 6822006 times.
✓ Branch 6 taken 1638936 times.
✓ Branch 7 taken 60450441 times.
✓ Branch 9 taken 8495379 times.
✓ Branch 10 taken 6147 times.
✓ Branch 12 taken 5433556 times.
✓ Branch 13 taken 1502569 times.
✓ Branch 15 taken 20538044 times.
✓ Branch 16 taken 49 times.
✓ Branch 18 taken 10030787 times.
✓ Branch 19 taken 21 times.
✓ Branch 21 taken 1308538 times.
✓ Branch 22 taken 19 times.
✓ Branch 24 taken 712255 times.
✓ Branch 25 taken 858527 times.
✓ Branch 29 taken 2441134 times.
✓ Branch 30 taken 542897 times.
✓ Branch 32 taken 39156 times.
✓ Branch 33 taken 790179 times.
✓ Branch 8 taken 1100258 times.
✓ Branch 20 taken 4275218 times.
✓ Branch 23 taken 310731 times.
✓ Branch 11 taken 2769790 times.
✓ Branch 14 taken 466870 times.
✓ Branch 17 taken 48 times.
✓ Branch 26 taken 5234904 times.
✓ Branch 28 taken 646182 times.
✓ Branch 31 taken 574665 times.
✓ Branch 27 taken 1381599 times.
✓ Branch 34 taken 3 times.
✓ Branch 36 taken 3 times.
✓ Branch 37 taken 374004 times.
✓ Branch 39 taken 3 times.
✓ Branch 40 taken 374004 times.
✓ Branch 41 taken 39043 times.
✓ Branch 43 taken 3 times.
✓ Branch 44 taken 1 times.
✓ Branch 35 taken 15 times.
✓ Branch 38 taken 39056 times.
✓ Branch 42 taken 1 times.
✓ Branch 45 taken 2 times.
✗ Branch 46 not taken.
✓ Branch 48 taken 1 times.
✗ Branch 49 not taken.
✓ Branch 50 taken 1 times.
✗ Branch 51 not taken.
✓ Branch 53 taken 1 times.
✗ Branch 54 not taken.
 696492547 { return *std::get<domainId>(problemTuple_); }
377
378 //! the finite volume grid geometry of domain i
379 template<std::size_t i>
380 39062556 const auto& gridGeometry(Dune::index_constant<i> domainId) const
381
8/12
✓ Branch 1 taken 7468405 times.
✓ Branch 2 taken 193914 times.
✓ Branch 4 taken 24428846 times.
✓ Branch 5 taken 851040 times.
✗ Branch 3 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 773889 times.
✓ Branch 8 taken 39040 times.
✓ Branch 9 taken 150000 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 2892000 times.
✗ Branch 13 not taken.
 37540948 { return *std::get<domainId>(gridGeometryTuple_); }
382
383 //! the grid view of domain i
384 template<std::size_t i>
385
9/16
✓ Branch 1 taken 1370 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 807 times.
✓ Branch 5 taken 5679 times.
✓ Branch 6 taken 7 times.
✓ Branch 7 taken 8089 times.
✓ Branch 9 taken 910 times.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 0 not taken.
✓ Branch 3 taken 1371 times.
✓ Branch 8 taken 7 times.
✓ Branch 11 taken 910 times.
 21496 const auto& gridView(Dune::index_constant<i> domainId) const
386
26/46
✓ Branch 1 taken 131 times.
✓ Branch 2 taken 1239 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 131 times.
✓ Branch 6 taken 1240 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 1371 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 146 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 1283 times.
✓ Branch 16 taken 2 times.
✓ Branch 17 taken 5201 times.
✓ Branch 19 taken 1168 times.
✓ Branch 20 taken 4380 times.
✓ Branch 22 taken 1283 times.
✓ Branch 23 taken 171 times.
✓ Branch 24 taken 1500 times.
✓ Branch 25 taken 5259 times.
✓ Branch 27 taken 169 times.
✓ Branch 28 taken 6628 times.
✗ Branch 29 not taken.
✓ Branch 30 taken 162 times.
✗ Branch 32 not taken.
✓ Branch 33 taken 748 times.
✓ Branch 35 taken 162 times.
✓ Branch 36 taken 748 times.
✗ Branch 37 not taken.
✓ Branch 38 taken 162 times.
✗ Branch 40 not taken.
✓ Branch 41 taken 748 times.
✗ Branch 10 not taken.
✓ Branch 18 taken 138 times.
✗ Branch 26 not taken.
✗ Branch 31 not taken.
✗ Branch 21 not taken.
✓ Branch 13 taken 823 times.
✗ Branch 5 not taken.
✗ Branch 34 not taken.
✗ Branch 39 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✓ Branch 44 taken 748 times.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
 42993 { return gridGeometry(domainId).gridView(); }
387
388 //! the grid variables of domain i
389 template<std::size_t i>
390 41827293 GridVariables<i>& gridVariables(Dune::index_constant<i> domainId)
391
15/23
✓ Branch 1 taken 894377 times.
✓ Branch 2 taken 5308770 times.
✓ Branch 4 taken 2554442 times.
✓ Branch 5 taken 21462027 times.
✓ Branch 7 taken 848114 times.
✓ Branch 8 taken 201045 times.
✓ Branch 10 taken 668890 times.
✓ Branch 11 taken 352352 times.
✓ Branch 13 taken 1171704 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 100012 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 100264 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 105781 times.
✗ Branch 23 not taken.
✓ Branch 3 taken 347500 times.
✓ Branch 6 taken 974000 times.
✗ Branch 9 not taken.
✗ Branch 12 not taken.
✗ Branch 15 not taken.
✓ Branch 25 taken 1928000 times.
✗ Branch 26 not taken.
 40520192 { return *std::get<domainId>(gridVariablesTuple_); }
392
393 //! the grid variables of domain i
394 template<std::size_t i>
395 500 const GridVariables<i>& gridVariables(Dune::index_constant<i> domainId) const
396
2/4
✓ Branch 1 taken 200 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 200 times.
✗ Branch 5 not taken.
 500 { return *std::get<domainId>(gridVariablesTuple_); }
397
398 //! the coupling manager
399 const CouplingManager& couplingManager() const
400 { return *couplingManager_; }
401
402 //! the full Jacobian matrix
403 7572 JacobianMatrix& jacobian()
404
2/6
✗ Branch 0 not taken.
✓ Branch 1 taken 6642 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 13 times.
✗ Branch 2 not taken.
✗ Branch 5 not taken.
 15123 { return *jacobian_; }
405
406 //! the full residual vector
407 11506 ResidualType& residual()
408
5/19
✓ Branch 5 taken 15 times.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 11 taken 6 times.
✓ Branch 12 taken 1 times.
✗ 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 2404 times.
✗ Branch 4 not taken.
✗ Branch 7 not taken.
✗ Branch 10 not taken.
✓ Branch 1 taken 119 times.
✗ Branch 2 not taken.
✗ Branch 13 not taken.
 11506 { return *residual_; }
409
410 //! the solution of the previous time step
411 5138804 const SolutionVector& prevSol() const
412
6/10
✓ Branch 4 taken 264000 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 100000 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 100000 times.
✗ Branch 11 not taken.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 100000 times.
✓ Branch 9 taken 617000 times.
 5845108 { return *prevSol_; }
413
414 /*!
415 * \brief Set time loop for instationary problems
416 * \note calling this turns this into a stationary assembler
417 */
418 void setTimeManager(std::shared_ptr<const TimeLoop> timeLoop)
419 { timeLoop_ = timeLoop; isStationaryProblem_ = !(static_cast<bool>(timeLoop)); }
420
421 /*!
422 * \brief Sets the solution from which to start the time integration. Has to be
423 * called prior to assembly for time-dependent problems.
424 */
425 void setPreviousSolution(const SolutionVector& u)
426 { prevSol_ = &u; }
427
428 /*!
429 * \brief Whether we are assembling a stationary or instationary problem
430 */
431 381011627 bool isStationaryProblem() const
432
18/22
✓ Branch 0 taken 58579817 times.
✓ Branch 1 taken 3640910 times.
✓ Branch 2 taken 158890249 times.
✓ Branch 3 taken 14550531 times.
✓ Branch 4 taken 30769995 times.
✓ Branch 5 taken 25484776 times.
✓ Branch 6 taken 33405966 times.
✓ Branch 7 taken 25181614 times.
✓ Branch 8 taken 2757828 times.
✓ Branch 9 taken 6618831 times.
✓ Branch 10 taken 192855 times.
✓ Branch 11 taken 12587902 times.
✓ Branch 12 taken 2892000 times.
✓ Branch 13 taken 388353 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 150000 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 2892000 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 100000 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 1928000 times.
 381011627 { return isStationaryProblem_; }
433
434 /*!
435 * \brief Create a local residual object (used by the local assembler)
436 */
437 template<std::size_t i>
438 39183620 LocalResidual<i> localResidual(Dune::index_constant<i> domainId) const
439
7/15
✓ Branch 1 taken 7662908 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 24465960 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1627331 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 150618 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 163489 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 2892000 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 300 times.
✗ Branch 20 not taken.
✗ Branch 6 not taken.
 38750336 { return LocalResidual<i>(std::get<domainId>(problemTuple_).get(), timeLoop_.get()); }
440
441 protected:
442 //! the coupling manager coupling the sub domains
443 std::shared_ptr<CouplingManager> couplingManager_;
444
445 private:
446 /*!
447 * \brief Sets the jacobian sparsity pattern.
448 */
449 234 void setJacobianPattern_(JacobianMatrix& jac) const
450 {
451 using namespace Dune::Hybrid;
452 1754 forEach(std::make_index_sequence<JacobianMatrix::N()>(), [&](const auto domainI)
453 {
454 4816 forEach(integralRange(Dune::Hybrid::size(jac[domainI])), [&](const auto domainJ)
455 {
456 1981 const auto pattern = this->getJacobianPattern_(domainI, domainJ);
457
9/18
✓ Branch 1 taken 225 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 225 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 225 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 225 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 43 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 43 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 42 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 42 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 40 times.
✗ Branch 26 not taken.
 1110 pattern.exportIdx(jac[domainI][domainJ]);
458 1112 });
459 });
460 194 }
461
462 /*!
463 * \brief Resizes the residual
464 */
465 232 void setResidualSize_(ResidualType& res) const
466 {
467 using namespace Dune::Hybrid;
468 1876 forEach(integralRange(Dune::Hybrid::size(res)), [&](const auto domainId)
469 465 { res[domainId].resize(this->numDofs(domainId)); });
470 159 }
471
472 // reset the residual vector to 0.0
473 12092 void resetResidual_()
474 {
475
2/2
✓ Branch 0 taken 21 times.
✓ Branch 1 taken 10561 times.
 12092 if(!residual_)
476 {
477
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 21 times.
 21 residual_ = std::make_shared<ResidualType>();
478 21 setResidualSize_(*residual_);
479 }
480
481 12092 (*residual_) = 0.0;
482 12092 }
483
484 // reset the jacobian vector to 0.0
485 10582 void resetJacobian_()
486 {
487
2/2
✓ Branch 0 taken 21 times.
✓ Branch 1 taken 9051 times.
 10582 if(!jacobian_)
488 {
489
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 21 times.
 21 jacobian_ = std::make_shared<JacobianMatrix>();
490 21 setJacobianBuildMode(*jacobian_);
491 21 setJacobianPattern_(*jacobian_);
492 }
493
494 10582 (*jacobian_) = 0.0;
495 10582 }
496
497 //! Computes the colors
498 221 void maybeComputeColors_()
499 {
500 if constexpr (CouplingManagerSupportsMultithreadedAssembly<CouplingManager>::value)
501
1/2
✓ Branch 0 taken 26 times.
✗ Branch 1 not taken.
 26 if (enableMultithreading_)
502
1/2
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
 26 couplingManager_->computeColorsForAssembly();
503 }
504
505 // check if the assembler is in a correct state for assembly
506 12092 void checkAssemblerState_() const
507 {
508
3/4
✓ Branch 0 taken 10226 times.
✓ Branch 1 taken 356 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 10226 times.
 12092 if (!isStationaryProblem_ && !prevSol_)
509 ✗ DUNE_THROW(Dune::InvalidStateException, "Assembling instationary problem but previous solution was not set!");
510
511
3/4
✓ Branch 0 taken 356 times.
✓ Branch 1 taken 10226 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 356 times.
 12092 if (isStationaryProblem_ && prevSol_)
512 ✗ DUNE_THROW(Dune::InvalidStateException, "Assembling stationary problem but a previous solution was set."
513 << " Did you forget to set the timeLoop to make this problem instationary?");
514 12092 }
515
516 template<std::size_t i, class JacRow, class SubRes>
517 19061 void assembleJacobianAndResidual_(Dune::index_constant<i> domainId, JacRow& jacRow, SubRes& subRes,
518 const SolutionVector& curSol)
519 {
520 97841685 assemble_(domainId, [&](const auto& element)
521 {
522
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 24455656 MultiDomainAssemblerSubDomainView view{*this, domainId};
523
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 24455656 SubDomainAssembler<i> subDomainAssembler(view, element, curSol, *couplingManager_);
524
10/13
✓ Branch 1 taken 5354528 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 261096 times.
✓ Branch 4 taken 13592648 times.
✓ Branch 6 taken 518162 times.
✓ Branch 7 taken 492181 times.
✓ Branch 8 taken 754509 times.
✓ Branch 9 taken 57708 times.
✓ Branch 5 taken 1304082 times.
✓ Branch 10 taken 1928000 times.
✓ Branch 13 taken 964000 times.
✗ Branch 14 not taken.
✗ Branch 11 not taken.
 27146800 subDomainAssembler.assembleJacobianAndResidual(jacRow, subRes, gridVariablesTuple_);
525 24455656 });
526 }
527
528 template<std::size_t i, class SubRes>
529 1509 void assembleResidual_(Dune::index_constant<i> domainId, SubRes& subRes,
530 const SolutionVector& curSol)
531 {
532 58251476 assemble_(domainId, [&](const auto& element)
533 {
534 14562492 MultiDomainAssemblerSubDomainView view{*this, domainId};
535 14562492 SubDomainAssembler<i> subDomainAssembler(view, element, curSol, *couplingManager_);
536
3/6
✓ Branch 1 taken 2773386 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 11789105 times.
✓ Branch 5 taken 1 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 14562492 subDomainAssembler.assembleResidual(subRes);
537 14562491 });
538 }
539
540 /*!
541 * \brief A method assembling something per element
542 * \note Handles exceptions for parallel runs
543 * \throws NumericalProblem on all processes if an exception is thrown during assembly
544 */
545 template<std::size_t i, class AssembleElementFunc>
546 44160 void assemble_(Dune::index_constant<i> domainId, AssembleElementFunc&& assembleElement) const
547 {
548 // a state that will be checked on all processes
549
1/2
✓ Branch 1 taken 16590 times.
✗ Branch 2 not taken.
 44160 bool succeeded = false;
550
551 // try assembling using the local assembly function
552 try
553 {
554 if constexpr (CouplingManagerSupportsMultithreadedAssembly<CouplingManager>::value)
555 {
556
1/2
✓ Branch 0 taken 584 times.
✗ Branch 1 not taken.
 1168 if (enableMultithreading_)
557 {
558
1/2
✓ Branch 1 taken 584 times.
✗ Branch 2 not taken.
 1168 couplingManager_->assembleMultithreaded(
559 domainId, std::forward<AssembleElementFunc>(assembleElement)
560 );
561 1168 return;
562 }
563 }
564
565 // fallback for coupling managers that don't support multithreaded assembly (yet)
566 // or if multithreaded assembly is disabled
567 // let the local assembler add the element contributions
568
12/13
✓ Branch 1 taken 16590 times.
✓ Branch 2 taken 2962815 times.
✓ Branch 4 taken 29341154 times.
✓ Branch 5 taken 4 times.
✓ Branch 7 taken 27206912 times.
✓ Branch 8 taken 50713 times.
✓ Branch 10 taken 3789113 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 3789085 times.
✓ Branch 13 taken 2346 times.
✓ Branch 3 taken 5278 times.
✓ Branch 6 taken 2069606 times.
✓ Branch 9 taken 326348 times.
 144365246 for (const auto& element : elements(gridView(domainId)))
569
2/2
✓ Branch 1 taken 36076810 times.
✓ Branch 2 taken 1 times.
 72153622 assembleElement(element);
570
571 // if we get here, everything worked well on this process
572 42990 succeeded = true;
573 }
574 // throw exception if a problem occurred
575
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 2 catch (NumericalProblem &e)
576 {
577
2/5
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 2 not taken.
 2 std::cout << "rank " << gridView(domainId).comm().rank()
578 2 << " caught an exception while assembling:" << e.what()
579
4/8
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 1 times.
✗ Branch 12 not taken.
 2 << "\n";
580 2 succeeded = false;
581 }
582
583 // make sure everything worked well on all processes
584
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 19150 times.
✓ Branch 2 taken 2346 times.
✓ Branch 3 taken 13982 times.
 42992 if (gridView(domainId).comm().size() > 1)
585 ✗ succeeded = gridView(domainId).comm().min(succeeded);
586
587 // if not succeeded rethrow the error on all processes
588
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 21495 times.
 42992 if (!succeeded)
589
11/22
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 1 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 1 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 1 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 1 times.
✗ Branch 29 not taken.
✓ Branch 32 taken 1 times.
✗ Branch 33 not taken.
 8 DUNE_THROW(NumericalProblem, "A process did not succeed in linearizing the system");
590 42990 }
591
592 // get diagonal block pattern
593 template<std::size_t i, std::size_t j, typename std::enable_if_t<(i==j), int> = 0>
594 548 Dune::MatrixIndexSet getJacobianPattern_(Dune::index_constant<i> domainI,
595 Dune::index_constant<j> domainJ) const
596 {
597 548 const auto& gg = gridGeometry(domainI);
598 548 auto pattern = getJacobianPattern<isImplicit()>(gg);
599
4/8
✓ Branch 1 taken 252 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 200 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 42 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
 548 couplingManager_->extendJacobianPattern(domainI, pattern);
600 115 return pattern;
601 ✗ }
602
603 // get coupling block pattern
604 template<std::size_t i, std::size_t j, typename std::enable_if_t<(i!=j), int> = 0>
605
2/4
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
 616 Dune::MatrixIndexSet getJacobianPattern_(Dune::index_constant<i> domainI,
606 Dune::index_constant<j> domainJ) const
607 {
608
2/4
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
 536 return getCouplingJacobianPattern<isImplicit()>(*couplingManager_,
609 domainI, gridGeometry(domainI),
610
2/4
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
 616 domainJ, gridGeometry(domainJ));
611 }
612
613 // build periodic constraints into the system matrix
614 template<std::size_t i, class JacRow, class Res, class GG, class Sol>
615 2258 void enforcePeriodicConstraints_(Dune::index_constant<i> domainI, JacRow& jacRow, Res& res, const GG& gridGeometry, const Sol& curSol)
616 {
617 if constexpr (Detail::hasPeriodicDofMap<GG>())
618 {
619
2/2
✓ Branch 0 taken 644 times.
✓ Branch 1 taken 1718 times.
 2902 for (const auto& m : gridGeometry.periodicDofMap())
620 {
621
2/2
✓ Branch 0 taken 322 times.
✓ Branch 1 taken 322 times.
 644 if (m.first < m.second)
622 {
623 322 auto& jac = jacRow[domainI];
624
625 // add the second row to the first
626 322 res[m.first] += res[m.second];
627
628 322 const auto end = jac[m.second].end();
629
2/2
✓ Branch 0 taken 2838 times.
✓ Branch 1 taken 322 times.
 3160 for (auto it = jac[m.second].begin(); it != end; ++it)
630 2838 jac[m.first][it.index()] += (*it);
631
632
633 // enforce the solution of the first periodic DOF to the second one
634 322 res[m.second] = curSol[m.second] - curSol[m.first];
635
636 // set derivatives accordingly in jacobian, i.e. id for m.second and -id for m.first
637 auto setMatrixBlock = [] (auto& matrixBlock, double diagValue)
638 {
639
0/4
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 644 for (int eIdx = 0; eIdx < matrixBlock.N(); ++eIdx)
640 644 matrixBlock[eIdx][eIdx] = diagValue;
641 };
642
643
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 3160 times.
✓ Branch 2 taken 2838 times.
✓ Branch 3 taken 322 times.
 3160 for (auto it = jac[m.second].begin(); it != end; ++it)
644 {
645 2838 auto& matrixBlock = *it;
646
2/2
✓ Branch 0 taken 322 times.
✓ Branch 1 taken 2516 times.
 2838 matrixBlock = 0.0;
647
648 assert(matrixBlock.N() == matrixBlock.M());
649
2/2
✓ Branch 0 taken 322 times.
✓ Branch 1 taken 2516 times.
 2838 if(it.index() == m.second)
650 2838 setMatrixBlock(matrixBlock, 1.0);
651
652
2/2
✓ Branch 0 taken 322 times.
✓ Branch 1 taken 2516 times.
 2838 if(it.index() == m.first)
653 2838 setMatrixBlock(matrixBlock, -1.0);
654
655 }
656
657 using namespace Dune::Hybrid;
658 322 forEach(makeIncompleteIntegerSequence<JacRow::size(), domainI>(), [&](const auto couplingDomainId)
659 {
660 322 auto& jacCoupling = jacRow[couplingDomainId];
661
662
3/8
✗ Branch 0 not taken.
✓ Branch 1 taken 644 times.
✓ Branch 2 taken 322 times.
✓ Branch 3 taken 322 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 644 for (auto it = jacCoupling[m.second].begin(); it != jacCoupling[m.second].end(); ++it)
663 322 jacCoupling[m.first][it.index()] += (*it);
664
665
2/4
✓ Branch 0 taken 322 times.
✓ Branch 1 taken 322 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 644 for (auto it = jacCoupling[m.second].begin(); it != jacCoupling[m.second].end(); ++it)
666 322 (*it) = 0.0;
667 });
668 }
669 }
670 }
671 2258 }
672
673 // enforce global constraints into the system matrix
674 template<std::size_t i, class JacRow, class Res, class GG, class Sol>
675 void enforceGlobalDirichletConstraints_(Dune::index_constant<i> domainI, JacRow& jacRow, Res& res, const GG& gridGeometry, const Sol& curSol)
676 {
677 if constexpr (Detail::hasSubProblemGlobalConstraints<Problem<domainI>>())
678 {
679 auto& jac = jacRow[domainI];
680 auto applyDirichletConstraint = [&] (const auto& dofIdx,
681 const auto& values,
682 const auto eqIdx,
683 const auto pvIdx)
684 {
685 (res)[dofIdx][eqIdx] = curSol[dofIdx][pvIdx] - values[pvIdx];
686
687 auto& row = jac[dofIdx];
688 for (auto col = row.begin(); col != row.end(); ++col)
689 row[col.index()][eqIdx] = 0.0;
690
691 jac[dofIdx][dofIdx][eqIdx][pvIdx] = 1.0;
692
693 using namespace Dune::Hybrid;
694 forEach(makeIncompleteIntegerSequence<JacRow::size(), domainI>(), [&](const auto couplingDomainId)
695 {
696 auto& jacCoupling = jacRow[couplingDomainId];
697
698 auto& rowCoupling = jacCoupling[dofIdx];
699 for (auto colCoupling = rowCoupling.begin(); colCoupling != rowCoupling.end(); ++colCoupling)
700 rowCoupling[colCoupling.index()][eqIdx] = 0.0;
701 });
702
703 };
704
705 for (const auto& constraintData : this->problem(domainI).constraints())
706 {
707 const auto& constraintInfo = constraintData.constraintInfo();
708 const auto& values = constraintData.values();
709 const auto dofIdx = constraintData.dofIndex();
710 // set the constraint in residual and jacobian
711 for (int eqIdx = 0; eqIdx < constraintInfo.size(); ++eqIdx)
712 {
713 if (constraintInfo.isConstraintEquation(eqIdx))
714 {
715 const auto pvIdx = constraintInfo.eqToPriVarIndex(eqIdx);
716 assert(0 <= pvIdx && pvIdx < constraintInfo.size());
717 applyDirichletConstraint(dofIdx, values, eqIdx, pvIdx);
718 }
719 }
720 }
721 }
722 }
723
724 //! pointer to the problem to be solved
725 ProblemTuple problemTuple_;
726
727 //! the finite volume geometry of the grid
728 GridGeometryTuple gridGeometryTuple_;
729
730 //! the variables container for the grid
731 GridVariablesTuple gridVariablesTuple_;
732
733 //! the time loop for instationary problem assembly
734 std::shared_ptr<const TimeLoop> timeLoop_;
735
736 //! an observing pointer to the previous solution for instationary problems
737 const SolutionVector* prevSol_ = nullptr;
738
739 //! if this assembler is assembling an instationary problem
740 bool isStationaryProblem_;
741
742 //! shared pointers to the jacobian matrix and residual
743 std::shared_ptr<JacobianMatrix> jacobian_;
744 std::shared_ptr<ResidualType> residual_;
745
746 //! Issue a warning if the calculation is used in parallel with overlap. This could be a static local variable if it wasn't for g++7 yielding a linker error.
747 mutable bool warningIssued_;
748
749 //! if multithreaded assembly is enabled
750 bool enableMultithreading_ = false;
751 };
752
753 } // end namespace Dumux
754
755 #endif
756