GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/experimental/assembly/multistagemultidomainfvassembler.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 120 166 72.3%
Functions: 29 73 39.7%
Branches: 190 594 32.0%
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 Experimental
10 * \ingroup MultiDomain
11 * \ingroup Assembly
12 * \brief A linear system assembler (residual and Jacobian) for finite volume schemes
13 * with multiple domains
14 */
15 #ifndef DUMUX_EXPERIMENTAL_MULTISTAGE_MULTIDOMAIN_FV_ASSEMBLER_HH
16 #define DUMUX_EXPERIMENTAL_MULTISTAGE_MULTIDOMAIN_FV_ASSEMBLER_HH
17
18 #include <vector>
19 #include <type_traits>
20 #include <tuple>
21 #include <memory>
22
23 #include <dune/common/hybridutilities.hh>
24 #include <dune/istl/matrixindexset.hh>
25
26 #include <dumux/common/exceptions.hh>
27 #include <dumux/common/properties.hh>
28 #include <dumux/common/typetraits/utility.hh>
29 #include <dumux/common/gridcapabilities.hh>
30
31 #include <dumux/discretization/method.hh>
32 #include <dumux/assembly/diffmethod.hh>
33 #include <dumux/assembly/jacobianpattern.hh>
34 #include <dumux/linear/parallelhelpers.hh>
35 #include <dumux/parallel/multithreading.hh>
36
37 #include <dumux/multidomain/couplingjacobianpattern.hh>
38 #include <dumux/multidomain/assemblerview.hh>
39
40 #include <dumux/experimental/assembly/subdomaincclocalassembler.hh>
41 #include <dumux/experimental/assembly/subdomaincvfelocalassembler.hh>
42 #include <dumux/experimental/assembly/multistagefvlocaloperator.hh>
43
44 #include <dumux/experimental/timestepping/multistagemethods.hh>
45 #include <dumux/experimental/timestepping/multistagetimestepper.hh>
46
47 namespace Dumux::Grid::Capabilities {
48
49 namespace Detail {
50 // helper for multi-domain models
51 template<class T, std::size_t... I>
52 bool allGridsSupportsMultithreadingImpl(const T& gridGeometries, std::index_sequence<I...>)
53 {
54 return (... && supportsMultithreading(std::get<I>(gridGeometries)->gridView()));
55 }
56 } // end namespace Detail
57
58 // helper for multi-domain models (all grids have to support multithreading)
59 template<class... GG>
60 bool allGridsSupportsMultithreading(const std::tuple<GG...>& gridGeometries)
61 {
62 return Detail::allGridsSupportsMultithreadingImpl<std::tuple<GG...>>(gridGeometries, std::make_index_sequence<sizeof...(GG)>());
63 }
64
65 } // end namespace Dumux::Grid::Capabilities
66
67 namespace Dumux {
68
69 /*!
70 * \ingroup Experimental
71 * \ingroup MultiDomain
72 * \ingroup Assembly
73 * \brief Type trait that is specialized for coupling manager supporting multithreaded assembly
74 * \note A coupling manager implementation that wants to enable multithreaded assembly has to specialize this trait
75 */
76 template<class CM>
77 struct CouplingManagerSupportsMultithreadedAssembly : public std::false_type
78 {};
79
80 } // end namespace Dumux
81
82 namespace Dumux::Experimental {
83
84 /*!
85 * \ingroup Experimental
86 * \ingroup MultiDomain
87 * \ingroup Assembly
88 * \brief A linear system assembler (residual and Jacobian) for finite volume schemes (box, tpfa, mpfa, ...)
89 * with multiple domains
90 * \tparam MDTraits the multidimensional traits
91 * \tparam diffMethod the differentiation method to residual compute derivatives
92 */
93 template<class MDTraits, class CMType, DiffMethod diffMethod>
94 class MultiStageMultiDomainFVAssembler
95 {
96 template<std::size_t id>
97 using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
98
99 public:
100 using Traits = MDTraits;
101
102 using Scalar = typename MDTraits::Scalar;
103 using StageParams = Experimental::MultiStageParams<Scalar>;
104
105 //! TODO get rid of this GetPropType
106 template<std::size_t id>
107 using LocalResidual = GetPropType<SubDomainTypeTag<id>, Properties::LocalResidual>;
108
109 template<std::size_t id>
110 using GridVariables = typename MDTraits::template SubDomain<id>::GridVariables;
111
112 template<std::size_t id>
113 using GridGeometry = typename MDTraits::template SubDomain<id>::GridGeometry;
114
115 template<std::size_t id>
116 using Problem = typename MDTraits::template SubDomain<id>::Problem;
117
118 using JacobianMatrix = typename MDTraits::JacobianMatrix;
119 using SolutionVector = typename MDTraits::SolutionVector;
120 using ResidualType = typename MDTraits::ResidualVector;
121
122 using CouplingManager = CMType;
123
124 private:
125
126 using ProblemTuple = typename MDTraits::template TupleOfSharedPtrConst<Problem>;
127 using GridGeometryTuple = typename MDTraits::template TupleOfSharedPtrConst<GridGeometry>;
128 using GridVariablesTuple = typename MDTraits::template TupleOfSharedPtr<GridVariables>;
129
130 using ThisType = MultiStageMultiDomainFVAssembler<MDTraits, CouplingManager, diffMethod>;
131
132 template<std::size_t id>
133 using SubDomainAssemblerView = MultiDomainAssemblerSubDomainView<ThisType, id>;
134
135 template<class DiscretizationMethod, std::size_t id>
136 struct SubDomainAssemblerType;
137
138 template<std::size_t id>
139 struct SubDomainAssemblerType<DiscretizationMethods::CCTpfa, id>
140 {
141 using type = Experimental::SubDomainCCLocalAssembler<id, SubDomainTypeTag<id>, SubDomainAssemblerView<id>, diffMethod>;
142 };
143
144 template<std::size_t id>
145 struct SubDomainAssemblerType<DiscretizationMethods::CCMpfa, id>
146 {
147 using type = Experimental::SubDomainCCLocalAssembler<id, SubDomainTypeTag<id>, SubDomainAssemblerView<id>, diffMethod>;
148 };
149
150 template<std::size_t id, class DM>
151 struct SubDomainAssemblerType<DiscretizationMethods::CVFE<DM>, id>
152 {
153 using type = Experimental::SubDomainCVFELocalAssembler<id, SubDomainTypeTag<id>, SubDomainAssemblerView<id>, diffMethod>;
154 };
155
156 template<std::size_t id>
157 using SubDomainAssembler = typename SubDomainAssemblerType<typename GridGeometry<id>::DiscretizationMethod, id>::type;
158
159 public:
160 /*!
161 * \brief The constructor for instationary problems
162 * \note the grid variables might be temporarily changed during assembly (if caching is enabled)
163 * it is however guaranteed that the state after assembly will be the same as before
164 */
165 2 MultiStageMultiDomainFVAssembler(ProblemTuple problem,
166 GridGeometryTuple gridGeometry,
167 GridVariablesTuple gridVariables,
168 std::shared_ptr<CouplingManager> couplingManager,
169 std::shared_ptr<const Experimental::MultiStageMethod<Scalar>> msMethod,
170 const SolutionVector& prevSol)
171 : couplingManager_(couplingManager)
172 , timeSteppingMethod_(msMethod)
173 2 , problemTuple_(std::move(problem))
174 2 , gridGeometryTuple_(std::move(gridGeometry))
175 2 , gridVariablesTuple_(std::move(gridVariables))
176
0/10
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
 2 , prevSol_(&prevSol)
177 {
178
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 4 std::cout << "Instantiated assembler for an instationary problem." << std::endl;
179
180 2 enableMultithreading_ = CouplingManagerSupportsMultithreadedAssembly<CouplingManager>::value
181 && Grid::Capabilities::allGridsSupportsMultithreading(gridGeometryTuple_)
182 && !Multithreading::isSerial()
183 && getParam<bool>("Assembly.Multithreading", true);
184
185 2 maybeComputeColors_();
186 2 }
187
188 /*!
189 * \brief Assembles the global Jacobian of the residual
190 * and the residual for the current solution.
191 */
192 54 void assembleJacobianAndResidual(const SolutionVector& curSol)
193 {
194 54 resetJacobian_();
195
196 54 resetResidual_();
197 108 spatialOperatorEvaluations_.back() = 0.0;
198 108 temporalOperatorEvaluations_.back() = 0.0;
199
200
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 54 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 54 times.
 162 if (stageParams_->size() != spatialOperatorEvaluations_.size())
201 ✗ DUNE_THROW(Dune::InvalidStateException, "Wrong number of residuals");
202
203 using namespace Dune::Hybrid;
204
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 1134 forEach(std::make_index_sequence<JacobianMatrix::N()>(), [&](const auto domainId)
205 {
206 // assemble the spatial and temporal residual of the current time step and the Jacobian
207 // w.r.t to the current solution (the current solution on the current stage)
208 216 auto& jacRow = (*jacobian_)[domainId];
209 216 auto& spatial = spatialOperatorEvaluations_.back()[domainId];
210 216 auto& temporal = temporalOperatorEvaluations_.back()[domainId];
211
212 55620 assemble_(domainId, [&](const auto& element)
213 {
214 6912 MultiDomainAssemblerSubDomainView view{*this, domainId};
215 20844 SubDomainAssembler<domainId()> subDomainAssembler(view, element, curSol, *couplingManager_);
216
6/12
✓ Branch 1 taken 3456 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3456 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3456 times.
✗ Branch 8 not taken.
✓ Branch 11 taken 3456 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 3456 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 3456 times.
✗ Branch 18 not taken.
 20736 subDomainAssembler.assembleJacobianAndResidual(
217
4/8
✓ Branch 1 taken 3456 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3456 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3456 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3456 times.
✗ Branch 11 not taken.
 13824 jacRow, (*residual_)[domainId],
218 gridVariablesTuple_,
219 *stageParams_, temporal, spatial,
220 constrainedDofs_[domainId]
221 );
222 });
223
224 // assemble the full residual for the time integration stage
225
2/4
✓ Branch 3 taken 54 times.
✗ Branch 4 not taken.
✓ Branch 8 taken 54 times.
✗ Branch 9 not taken.
 432 auto constantResidualComponent = (*residual_)[domainId];
226 constantResidualComponent = 0.0;
227
8/8
✓ Branch 0 taken 54 times.
✓ Branch 1 taken 54 times.
✓ Branch 2 taken 54 times.
✓ Branch 3 taken 54 times.
✓ Branch 4 taken 54 times.
✓ Branch 5 taken 54 times.
✓ Branch 6 taken 54 times.
✓ Branch 7 taken 54 times.
 324 for (std::size_t k = 0; k < stageParams_->size()-1; ++k)
228 {
229
6/12
✓ Branch 0 taken 54 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 54 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 54 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 54 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 54 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 54 times.
✗ Branch 11 not taken.
 324 if (!stageParams_->skipTemporal(k))
230 540 constantResidualComponent.axpy(stageParams_->temporalWeight(k), temporalOperatorEvaluations_[k][domainId]);
231
6/12
✗ Branch 0 not taken.
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 54 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 54 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 54 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 54 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 54 times.
 324 if (!stageParams_->skipSpatial(k))
232 ✗ constantResidualComponent.axpy(stageParams_->spatialWeight(k), spatialOperatorEvaluations_[k][domainId]);
233 }
234
235 // masked summation of constant residual component onto this stage's resiudal component
236
4/4
✓ Branch 0 taken 6750 times.
✓ Branch 1 taken 54 times.
✓ Branch 2 taken 3456 times.
✓ Branch 3 taken 54 times.
 10314 for (std::size_t i = 0; i < constantResidualComponent.size(); ++i)
237
8/8
✓ Branch 0 taken 20250 times.
✓ Branch 1 taken 6750 times.
✓ Branch 2 taken 20250 times.
✓ Branch 3 taken 6750 times.
✓ Branch 4 taken 6912 times.
✓ Branch 5 taken 3456 times.
✓ Branch 6 taken 6912 times.
✓ Branch 7 taken 3456 times.
 74736 for (std::size_t ii = 0; ii < constantResidualComponent[i].size(); ++ii)
238
12/16
✓ Branch 0 taken 4374 times.
✓ Branch 1 taken 15876 times.
✓ Branch 2 taken 4374 times.
✓ Branch 3 taken 15876 times.
✓ Branch 4 taken 4374 times.
✓ Branch 5 taken 15876 times.
✓ Branch 6 taken 4374 times.
✓ Branch 7 taken 15876 times.
✓ Branch 8 taken 6912 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 6912 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 6912 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 6912 times.
✗ Branch 15 not taken.
 108648 (*residual_)[domainId][i][ii] += constrainedDofs_[domainId][i][ii] > 0.5 ? 0.0 : constantResidualComponent[i][ii];
239 });
240 54 }
241
242 //! compute the residuals using the internal residual
243 ✗ void assembleResidual(const SolutionVector& curSol)
244 ✗ { DUNE_THROW(Dune::NotImplemented, "residual"); }
245
246 /*!
247 * \brief The version without arguments uses the default constructor to create
248 * the jacobian and residual objects in this assembler if you don't need them outside this class
249 */
250 2 void setLinearSystem()
251 {
252
2/4
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
 2 jacobian_ = std::make_shared<JacobianMatrix>();
253
2/4
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
 2 residual_ = std::make_shared<ResidualType>();
254
255 2 setJacobianBuildMode_(*jacobian_);
256 6 setJacobianPattern_(*jacobian_);
257 6 setResidualSize_(*residual_);
258 2 }
259
260 /*!
261 * \brief Updates the grid variables with the given solution
262 */
263 ✗ void updateGridVariables(const SolutionVector& curSol)
264 {
265 using namespace Dune::Hybrid;
266 162 forEach(integralRange(Dune::Hybrid::size(gridVariablesTuple_)), [&](const auto domainId)
267 432 { this->gridVariables(domainId).update(curSol[domainId]); });
268 ✗ }
269
270 /*!
271 * \brief Resets the grid variables to the last time step
272 */
273 ✗ void resetTimeStep(const SolutionVector& curSol)
274 {
275 using namespace Dune::Hybrid;
276 ✗ forEach(integralRange(Dune::Hybrid::size(gridVariablesTuple_)), [&](const auto domainId)
277 ✗ { this->gridVariables(domainId).resetTimeStep(curSol[domainId]); });
278
279 ✗ this->clearStages();
280 ✗ }
281
282 //! the number of dof locations of domain i
283 template<std::size_t i>
284 ✗ std::size_t numDofs(Dune::index_constant<i> domainId) const
285 192 { return std::get<domainId>(gridGeometryTuple_)->numDofs(); }
286
287 //! the problem of domain i
288 template<std::size_t i>
289 ✗ const auto& problem(Dune::index_constant<i> domainId) const
290
27/66
✗ 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 31 taken 2 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 2 times.
✗ Branch 35 not taken.
✓ Branch 37 taken 2 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 2 times.
✗ Branch 41 not taken.
✓ Branch 43 taken 2 times.
✗ Branch 44 not taken.
✓ Branch 46 taken 2 times.
✗ Branch 47 not taken.
✓ Branch 49 taken 2 times.
✗ Branch 50 not taken.
✓ Branch 52 taken 2 times.
✗ Branch 53 not taken.
✓ Branch 55 taken 2 times.
✗ Branch 56 not taken.
✓ Branch 60 taken 2 times.
✓ Branch 61 taken 6910 times.
✓ Branch 62 taken 2 times.
✓ Branch 63 taken 6910 times.
✓ Branch 64 taken 2 times.
✓ Branch 65 taken 6910 times.
✗ Branch 67 not taken.
✗ Branch 68 not taken.
✗ Branch 70 not taken.
✗ Branch 71 not taken.
✗ Branch 73 not taken.
✗ Branch 74 not taken.
✗ Branch 76 not taken.
✗ Branch 77 not taken.
✗ Branch 79 not taken.
✗ Branch 80 not taken.
✗ Branch 82 not taken.
✗ Branch 83 not taken.
✓ Branch 88 taken 2 times.
✗ Branch 89 not taken.
✓ Branch 91 taken 2 times.
✗ Branch 92 not taken.
✓ Branch 94 taken 2 times.
✗ Branch 95 not taken.
✓ Branch 97 taken 2 times.
✗ Branch 98 not taken.
✓ Branch 100 taken 2 times.
✗ Branch 101 not taken.
✓ Branch 103 taken 2 times.
✗ Branch 104 not taken.
✓ Branch 111 taken 2 times.
✓ Branch 112 taken 269566 times.
✓ Branch 113 taken 2 times.
✓ Branch 114 taken 269566 times.
✓ Branch 115 taken 2 times.
✓ Branch 116 taken 269566 times.
 1455390 { return *std::get<domainId>(problemTuple_); }
291
292 //! the finite volume grid geometry of domain i
293 template<std::size_t i>
294 ✗ const auto& gridGeometry(Dune::index_constant<i> domainId) const
295
14/36
✓ Branch 9 taken 54 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 54 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 17 taken 54 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✓ Branch 20 taken 54 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✓ Branch 25 taken 4736 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 4736 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 4736 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 4736 times.
✗ Branch 35 not taken.
✓ Branch 37 taken 4736 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 4736 times.
✗ Branch 41 not taken.
✓ Branch 43 taken 20 times.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
✓ Branch 46 taken 20 times.
✗ Branch 48 not taken.
✗ Branch 49 not taken.
✓ Branch 51 taken 20 times.
✗ Branch 52 not taken.
✗ Branch 53 not taken.
✓ Branch 54 taken 20 times.
✗ Branch 56 not taken.
✗ Branch 57 not taken.
 28720 { return *std::get<domainId>(gridGeometryTuple_); }
296
297 //! the grid view of domain i
298 template<std::size_t i>
299 ✗ const auto& gridView(Dune::index_constant<i> domainId) const
300
8/24
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 54 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 9 taken 54 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 54 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 17 taken 20 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✓ Branch 20 taken 20 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✓ Branch 25 taken 20 times.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✓ Branch 28 taken 20 times.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
 296 { return gridGeometry(domainId).gridView(); }
301
302 //! the grid variables of domain i
303 template<std::size_t i>
304 ✗ GridVariables<i>& gridVariables(Dune::index_constant<i> domainId)
305
15/30
✓ Branch 19 taken 4736 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 4736 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 4736 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 4736 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 4736 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 4736 times.
✗ Branch 35 not taken.
✓ Branch 37 taken 4736 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 4736 times.
✗ Branch 41 not taken.
✓ Branch 43 taken 4736 times.
✗ Branch 44 not taken.
✓ Branch 46 taken 4736 times.
✗ Branch 47 not taken.
✓ Branch 49 taken 4736 times.
✗ Branch 50 not taken.
✓ Branch 52 taken 4736 times.
✗ Branch 53 not taken.
✓ Branch 55 taken 4736 times.
✗ Branch 56 not taken.
✓ Branch 58 taken 4736 times.
✗ Branch 59 not taken.
✓ Branch 61 taken 4736 times.
✗ Branch 62 not taken.
 71364 { return *std::get<domainId>(gridVariablesTuple_); }
306
307 //! the grid variables of domain i
308 template<std::size_t i>
309 const GridVariables<i>& gridVariables(Dune::index_constant<i> domainId) const
310 { return *std::get<domainId>(gridVariablesTuple_); }
311
312 //! the coupling manager
313 const CouplingManager& couplingManager() const
314 { return *couplingManager_; }
315
316 //! the full Jacobian matrix
317 JacobianMatrix& jacobian()
318
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 54 times.
 216 { return *jacobian_; }
319
320 //! the full residual vector
321 ResidualType& residual()
322
2/4
✓ Branch 3 taken 20 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 20 times.
✗ Branch 7 not taken.
 148 { return *residual_; }
323
324 //! the solution before time integration
325 ✗ const SolutionVector& prevSol() const
326 ✗ { return *prevSol_; }
327
328 /*!
329 * \brief Create a local residual object (used by the local assembler)
330 */
331 template<std::size_t i>
332 ✗ MultiStageFVLocalOperator<LocalResidual<i>> localResidual(Dune::index_constant<i> domainId) const
333
12/24
✓ Branch 1 taken 4736 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4736 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4736 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 4736 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 4736 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 4736 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 4736 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 4736 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 20 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 20 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 20 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 20 times.
✗ Branch 35 not taken.
 37968 { return { LocalResidual<i>{std::get<domainId>(problemTuple_).get(), nullptr} }; }
334
335 40 void clearStages()
336 {
337
2/2
✓ Branch 0 taken 20 times.
✓ Branch 1 taken 20 times.
 40 spatialOperatorEvaluations_.clear();
338
2/2
✓ Branch 0 taken 20 times.
✓ Branch 1 taken 20 times.
 40 temporalOperatorEvaluations_.clear();
339
2/2
✓ Branch 0 taken 20 times.
✓ Branch 1 taken 20 times.
 40 stageParams_.reset();
340 40 }
341
342 template<class StageParams>
343 20 void prepareStage(SolutionVector& x, StageParams params)
344 {
345
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 20 times.
 20 stageParams_ = std::move(params);
346
1/2
✓ Branch 0 taken 20 times.
✗ Branch 1 not taken.
 20 const auto curStage = stageParams_->size() - 1;
347
348 // in the first stage, also assemble the old residual
349
1/2
✓ Branch 0 taken 20 times.
✗ Branch 1 not taken.
 20 if (curStage == 1)
350 {
351 // update time in variables?
352 using namespace Dune::Hybrid;
353 20 forEach(std::make_index_sequence<JacobianMatrix::N()>(), [&](const auto domainId)
354 {
355 160 setProblemTime_(*std::get<domainId>(problemTuple_), stageParams_->timeAtStage(curStage));
356 });
357
358 20 resetResidual_(); // residual resized and zero
359 40 spatialOperatorEvaluations_.push_back(*residual_);
360 40 temporalOperatorEvaluations_.push_back(*residual_);
361
362 // assemble stage 0 residuals
363 using namespace Dune::Hybrid;
364 20 forEach(std::make_index_sequence<JacobianMatrix::N()>(), [&](const auto domainId)
365 {
366
0/8
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
 40 auto& spatial = spatialOperatorEvaluations_.back()[domainId];
367
0/8
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
 40 auto& temporal = temporalOperatorEvaluations_.back()[domainId];
368 15380 assemble_(domainId, [&](const auto& element)
369 {
370 2560 MultiDomainAssemblerSubDomainView view{*this, domainId};
371
0/4
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 7680 SubDomainAssembler<domainId()> subDomainAssembler(view, element, x, *couplingManager_);
372
2/4
✓ Branch 1 taken 1280 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1280 times.
✗ Branch 5 not taken.
 2560 subDomainAssembler.localResidual().spatialWeight(1.0);
373
2/4
✓ Branch 1 taken 1280 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1280 times.
✗ Branch 5 not taken.
 2560 subDomainAssembler.localResidual().temporalWeight(1.0);
374
2/4
✓ Branch 1 taken 1280 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1280 times.
✗ Branch 5 not taken.
 2560 subDomainAssembler.assembleCurrentResidual(spatial, temporal);
375 });
376 });
377 }
378
379 // update time in variables?
380 using namespace Dune::Hybrid;
381 20 forEach(std::make_index_sequence<JacobianMatrix::N()>(), [&](const auto domainId)
382 {
383 160 setProblemTime_(*std::get<domainId>(problemTuple_), stageParams_->timeAtStage(curStage));
384 });
385
386 20 resetResidual_(); // residual resized and zero
387 40 spatialOperatorEvaluations_.push_back(*residual_);
388 40 temporalOperatorEvaluations_.push_back(*residual_);
389 20 }
390
391 //! TODO get rid of this (called by Newton but shouldn't be necessary)
392 ✗ bool isStationaryProblem() const
393 ✗ { return false; }
394
395 bool isImplicit() const
396 406016 { return timeSteppingMethod_->implicit(); }
397
398 protected:
399 //! the coupling manager coupling the sub domains
400 std::shared_ptr<CouplingManager> couplingManager_;
401
402 private:
403 /*!
404 * \brief Sets the jacobian build mode
405 */
406 ✗ void setJacobianBuildMode_(JacobianMatrix& jac) const
407 {
408 using namespace Dune::Hybrid;
409
0/2
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 2 forEach(std::make_index_sequence<JacobianMatrix::N()>(), [&](const auto i)
410 {
411
0/8
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
 16 forEach(jac[i], [&](auto& jacBlock)
412 {
413 using BlockType = std::decay_t<decltype(jacBlock)>;
414
4/8
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
 8 if (jacBlock.buildMode() == BlockType::BuildMode::unknown)
415 8 jacBlock.setBuildMode(BlockType::BuildMode::random);
416 ✗ else if (jacBlock.buildMode() != BlockType::BuildMode::random)
417 ✗ DUNE_THROW(Dune::NotImplemented, "Only BCRS matrices with random build mode are supported at the moment");
418 });
419 });
420 ✗ }
421
422 /*!
423 * \brief Sets the jacobian sparsity pattern.
424 */
425 void setJacobianPattern_(JacobianMatrix& jac) const
426 {
427 using namespace Dune::Hybrid;
428 2 forEach(std::make_index_sequence<JacobianMatrix::N()>(), [&](const auto domainI)
429 {
430
0/8
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
 48 forEach(integralRange(Dune::Hybrid::size(jac[domainI])), [&](const auto domainJ)
431 {
432
0/8
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
 16 const auto pattern = this->getJacobianPattern_(domainI, domainJ);
433
12/48
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 2 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 2 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 2 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 2 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 2 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 2 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 2 times.
✗ Branch 35 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✗ Branch 52 not taken.
✗ Branch 53 not taken.
✗ Branch 55 not taken.
✗ Branch 56 not taken.
✗ Branch 58 not taken.
✗ Branch 59 not taken.
✗ Branch 61 not taken.
✗ Branch 62 not taken.
✗ Branch 64 not taken.
✗ Branch 65 not taken.
✗ Branch 67 not taken.
✗ Branch 68 not taken.
✗ Branch 70 not taken.
✗ Branch 71 not taken.
 24 pattern.exportIdx(jac[domainI][domainJ]);
434 });
435 });
436 }
437
438 /*!
439 * \brief Resizes the residual
440 */
441 void setResidualSize_(ResidualType& res) const
442 {
443 using namespace Dune::Hybrid;
444 194 forEach(integralRange(Dune::Hybrid::size(res)), [&](const auto domainId)
445 192 { res[domainId].resize(this->numDofs(domainId)); });
446 }
447
448 // reset the residual vector to 0.0
449 94 void resetResidual_()
450 {
451
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 94 times.
 94 if(!residual_)
452 {
453 ✗ residual_ = std::make_shared<ResidualType>();
454 ✗ setResidualSize_(*residual_);
455 }
456
457 94 setResidualSize_(constrainedDofs_);
458
459 188 (*residual_) = 0.0;
460 94 constrainedDofs_ = 0.0;
461 94 }
462
463 // reset the jacobian vector to 0.0
464 54 void resetJacobian_()
465 {
466
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 54 times.
 54 if(!jacobian_)
467 {
468 ✗ jacobian_ = std::make_shared<JacobianMatrix>();
469 ✗ setJacobianBuildMode_(*jacobian_);
470 ✗ setJacobianPattern_(*jacobian_);
471 }
472
473 108 (*jacobian_) = 0.0;
474 54 }
475
476 //! Computes the colors
477 ✗ void maybeComputeColors_()
478 {
479 if constexpr (CouplingManagerSupportsMultithreadedAssembly<CouplingManager>::value)
480 if (enableMultithreading_)
481 couplingManager_->computeColorsForAssembly();
482 ✗ }
483
484 template<std::size_t i, class SubRes>
485 void assembleResidual_(Dune::index_constant<i> domainId, SubRes& subRes,
486 const SolutionVector& curSol)
487 {
488 DUNE_THROW(Dune::NotImplemented, "assembleResidual_");
489 }
490
491 /*!
492 * \brief A method assembling something per element
493 * \note Handles exceptions for parallel runs
494 * \throws NumericalProblem on all processes if an exception is thrown during assembly
495 */
496 template<std::size_t i, class AssembleElementFunc>
497 296 void assemble_(Dune::index_constant<i> domainId, AssembleElementFunc&& assembleElement) const
498 {
499 // a state that will be checked on all processes
500 296 bool succeeded = false;
501
502 // try assembling using the local assembly function
503 try
504 {
505 if constexpr (CouplingManagerSupportsMultithreadedAssembly<CouplingManager>::value)
506 {
507 if (enableMultithreading_)
508 {
509 couplingManager_->assembleMultithreaded(
510 domainId, std::forward<AssembleElementFunc>(assembleElement)
511 );
512 return;
513 }
514 }
515
516 // fallback for coupling managers that don't support multithreaded assembly (yet)
517 // or if multithreaded assembly is disabled
518 // let the local assembler add the element contributions
519
4/8
✓ Branch 1 taken 148 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 148 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 9620 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 9472 times.
✗ Branch 10 not taken.
 19536 for (const auto& element : elements(gridView(domainId)))
520
1/2
✓ Branch 1 taken 9472 times.
✗ Branch 2 not taken.
 18944 assembleElement(element);
521
522 // if we get here, everything worked well on this process
523 296 succeeded = true;
524 }
525 // throw exception if a problem occurred
526 ✗ catch (NumericalProblem &e)
527 {
528 ✗ std::cout << "rank " << gridView(domainId).comm().rank()
529 ✗ << " caught an exception while assembling:" << e.what()
530 ✗ << "\n";
531 ✗ succeeded = false;
532 }
533
534 // make sure everything worked well on all processes
535
4/8
✗ Branch 0 not taken.
✓ Branch 1 taken 148 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 148 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 148 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 148 times.
 592 if (gridView(domainId).comm().size() > 1)
536 ✗ succeeded = gridView(domainId).comm().min(succeeded);
537
538 // if not succeeded rethrow the error on all processes
539
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 148 times.
 296 if (!succeeded)
540 ✗ DUNE_THROW(NumericalProblem, "A process did not succeed in linearizing the system");
541 296 }
542
543 // get diagonal block pattern
544 template<std::size_t i, std::size_t j, typename std::enable_if_t<(i==j), int> = 0>
545 8 Dune::MatrixIndexSet getJacobianPattern_(Dune::index_constant<i> domainI,
546 Dune::index_constant<j> domainJ) const
547 {
548 8 const auto& gg = gridGeometry(domainI);
549
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 16 if (timeSteppingMethod_->implicit())
550 {
551 16 auto pattern = getJacobianPattern<true>(gg);
552 8 couplingManager_->extendJacobianPattern(domainI, pattern);
553 16 return pattern;
554 }
555 else
556 {
557 ✗ auto pattern = getJacobianPattern<false>(gg);
558 ✗ couplingManager_->extendJacobianPattern(domainI, pattern);
559 ✗ return pattern;
560 }
561 }
562
563 // get coupling block pattern
564 template<std::size_t i, std::size_t j, typename std::enable_if_t<(i!=j), int> = 0>
565 4 Dune::MatrixIndexSet getJacobianPattern_(Dune::index_constant<i> domainI,
566 Dune::index_constant<j> domainJ) const
567 {
568
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 8 if (timeSteppingMethod_->implicit())
569 8 return getCouplingJacobianPattern<true>(*couplingManager_,
570 domainI, gridGeometry(domainI),
571 domainJ, gridGeometry(domainJ)
572 16 );
573 else
574 ✗ return getCouplingJacobianPattern<false>(*couplingManager_,
575 domainI, gridGeometry(domainI),
576 domainJ, gridGeometry(domainJ)
577 ✗ );
578 }
579
580 // TODO make this nicer with a is_detected trait in a common location
581 template<class P>
582 ✗ void setProblemTime_(const P& p, const Scalar t)
583 160 { setProblemTimeImpl_(p, t, 0); }
584
585 template<class P>
586 ✗ auto setProblemTimeImpl_(const P& p, const Scalar t, int) -> decltype(p.setTime(0))
587 160 { p.setTime(t); }
588
589 template<class P>
590 void setProblemTimeImpl_(const P& p, const Scalar t, long)
591 {}
592
593 std::shared_ptr<const Experimental::MultiStageMethod<Scalar>> timeSteppingMethod_;
594 std::vector<ResidualType> spatialOperatorEvaluations_;
595 std::vector<ResidualType> temporalOperatorEvaluations_;
596 ResidualType constrainedDofs_;
597 std::shared_ptr<const StageParams> stageParams_;
598
599 //! pointer to the problem to be solved
600 ProblemTuple problemTuple_;
601
602 //! the finite volume geometry of the grid
603 GridGeometryTuple gridGeometryTuple_;
604
605 //! the variables container for the grid
606 GridVariablesTuple gridVariablesTuple_;
607
608 //! Pointer to the previous solution
609 const SolutionVector* prevSol_;
610
611 //! shared pointers to the jacobian matrix and residual
612 std::shared_ptr<JacobianMatrix> jacobian_;
613 std::shared_ptr<ResidualType> residual_;
614
615 //! if multithreaded assembly is enabled
616 bool enableMultithreading_ = false;
617 };
618
619 } // end namespace Dumux
620
621 #endif
622