GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/experimental/assembly/multistagemultidomainfvassembler.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 120 165 72.7%
Functions: 29 73 39.7%
Branches: 191 616 31.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
280 //! the number of dof locations of domain i
281 template<std::size_t i>
282 ✗ std::size_t numDofs(Dune::index_constant<i> domainId) const
283 192 { return std::get<domainId>(gridGeometryTuple_)->numDofs(); }
284
285 //! the problem of domain i
286 template<std::size_t i>
287 ✗ const auto& problem(Dune::index_constant<i> domainId) const
288
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_); }
289
290 //! the finite volume grid geometry of domain i
291 template<std::size_t i>
292 ✗ const auto& gridGeometry(Dune::index_constant<i> domainId) const
293
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_); }
294
295 //! the grid view of domain i
296 template<std::size_t i>
297 ✗ const auto& gridView(Dune::index_constant<i> domainId) const
298
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(); }
299
300 //! the grid variables of domain i
301 template<std::size_t i>
302 ✗ GridVariables<i>& gridVariables(Dune::index_constant<i> domainId)
303
15/42
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 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 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_); }
304
305 //! the grid variables of domain i
306 template<std::size_t i>
307 const GridVariables<i>& gridVariables(Dune::index_constant<i> domainId) const
308 { return *std::get<domainId>(gridVariablesTuple_); }
309
310 //! the coupling manager
311 const CouplingManager& couplingManager() const
312 { return *couplingManager_; }
313
314 //! the full Jacobian matrix
315 JacobianMatrix& jacobian()
316
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 54 times.
 216 { return *jacobian_; }
317
318 //! the full residual vector
319 ResidualType& residual()
320
2/4
✓ Branch 3 taken 20 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 20 times.
✗ Branch 7 not taken.
 148 { return *residual_; }
321
322 //! the solution before time integration
323 ✗ const SolutionVector& prevSol() const
324 ✗ { return *prevSol_; }
325
326 /*!
327 * \brief Create a local residual object (used by the local assembler)
328 */
329 template<std::size_t i>
330 ✗ MultiStageFVLocalOperator<LocalResidual<i>> localResidual(Dune::index_constant<i> domainId) const
331
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} }; }
332
333 40 void clearStages()
334 {
335
2/2
✓ Branch 0 taken 20 times.
✓ Branch 1 taken 20 times.
 40 spatialOperatorEvaluations_.clear();
336
2/2
✓ Branch 0 taken 20 times.
✓ Branch 1 taken 20 times.
 40 temporalOperatorEvaluations_.clear();
337
2/2
✓ Branch 0 taken 20 times.
✓ Branch 1 taken 20 times.
 40 stageParams_.reset();
338 40 }
339
340 template<class StageParams>
341 20 void prepareStage(SolutionVector& x, StageParams params)
342 {
343
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 20 times.
 40 stageParams_ = std::move(params);
344
1/2
✓ Branch 0 taken 20 times.
✗ Branch 1 not taken.
 20 const auto curStage = stageParams_->size() - 1;
345
346 // in the first stage, also assemble the old residual
347
1/2
✓ Branch 0 taken 20 times.
✗ Branch 1 not taken.
 20 if (curStage == 1)
348 {
349 // update time in variables?
350 using namespace Dune::Hybrid;
351 20 forEach(std::make_index_sequence<JacobianMatrix::N()>(), [&](const auto domainId)
352 {
353 160 setProblemTime_(*std::get<domainId>(problemTuple_), stageParams_->timeAtStage(curStage));
354 });
355
356 20 resetResidual_(); // residual resized and zero
357 40 spatialOperatorEvaluations_.push_back(*residual_);
358 40 temporalOperatorEvaluations_.push_back(*residual_);
359
360 // assemble stage 0 residuals
361 using namespace Dune::Hybrid;
362 20 forEach(std::make_index_sequence<JacobianMatrix::N()>(), [&](const auto domainId)
363 {
364
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];
365
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];
366 15380 assemble_(domainId, [&](const auto& element)
367 {
368 2560 MultiDomainAssemblerSubDomainView view{*this, domainId};
369
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_);
370
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);
371
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);
372
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);
373 });
374 });
375 }
376
377 // update time in variables?
378 using namespace Dune::Hybrid;
379 20 forEach(std::make_index_sequence<JacobianMatrix::N()>(), [&](const auto domainId)
380 {
381 160 setProblemTime_(*std::get<domainId>(problemTuple_), stageParams_->timeAtStage(curStage));
382 });
383
384 20 resetResidual_(); // residual resized and zero
385 40 spatialOperatorEvaluations_.push_back(*residual_);
386 40 temporalOperatorEvaluations_.push_back(*residual_);
387 20 }
388
389 //! TODO get rid of this (called by Newton but shouldn't be necessary)
390 ✗ bool isStationaryProblem() const
391 ✗ { return false; }
392
393 bool isImplicit() const
394 406016 { return timeSteppingMethod_->implicit(); }
395
396 protected:
397 //! the coupling manager coupling the sub domains
398 std::shared_ptr<CouplingManager> couplingManager_;
399
400 private:
401 /*!
402 * \brief Sets the jacobian build mode
403 */
404 ✗ void setJacobianBuildMode_(JacobianMatrix& jac) const
405 {
406 using namespace Dune::Hybrid;
407
0/2
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 2 forEach(std::make_index_sequence<JacobianMatrix::N()>(), [&](const auto i)
408 {
409
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)
410 {
411 using BlockType = std::decay_t<decltype(jacBlock)>;
412
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)
413 8 jacBlock.setBuildMode(BlockType::BuildMode::random);
414 ✗ else if (jacBlock.buildMode() != BlockType::BuildMode::random)
415 ✗ DUNE_THROW(Dune::NotImplemented, "Only BCRS matrices with random build mode are supported at the moment");
416 });
417 });
418 ✗ }
419
420 /*!
421 * \brief Sets the jacobian sparsity pattern.
422 */
423 void setJacobianPattern_(JacobianMatrix& jac) const
424 {
425 using namespace Dune::Hybrid;
426 2 forEach(std::make_index_sequence<JacobianMatrix::N()>(), [&](const auto domainI)
427 {
428
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)
429 {
430
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);
431
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]);
432 });
433 });
434 }
435
436 /*!
437 * \brief Resizes the residual
438 */
439 void setResidualSize_(ResidualType& res) const
440 {
441 using namespace Dune::Hybrid;
442 194 forEach(integralRange(Dune::Hybrid::size(res)), [&](const auto domainId)
443 192 { res[domainId].resize(this->numDofs(domainId)); });
444 }
445
446 // reset the residual vector to 0.0
447 94 void resetResidual_()
448 {
449
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 94 times.
 94 if(!residual_)
450 {
451 ✗ residual_ = std::make_shared<ResidualType>();
452 ✗ setResidualSize_(*residual_);
453 }
454
455 94 setResidualSize_(constrainedDofs_);
456
457 188 (*residual_) = 0.0;
458 94 constrainedDofs_ = 0.0;
459 94 }
460
461 // reset the jacobian vector to 0.0
462 54 void resetJacobian_()
463 {
464
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 54 times.
 54 if(!jacobian_)
465 {
466 ✗ jacobian_ = std::make_shared<JacobianMatrix>();
467 ✗ setJacobianBuildMode_(*jacobian_);
468 ✗ setJacobianPattern_(*jacobian_);
469 }
470
471 108 (*jacobian_) = 0.0;
472 54 }
473
474 //! Computes the colors
475 ✗ void maybeComputeColors_()
476 {
477 if constexpr (CouplingManagerSupportsMultithreadedAssembly<CouplingManager>::value)
478 if (enableMultithreading_)
479 couplingManager_->computeColorsForAssembly();
480 ✗ }
481
482 template<std::size_t i, class SubRes>
483 void assembleResidual_(Dune::index_constant<i> domainId, SubRes& subRes,
484 const SolutionVector& curSol)
485 {
486 DUNE_THROW(Dune::NotImplemented, "assembleResidual_");
487 }
488
489 /*!
490 * \brief A method assembling something per element
491 * \note Handles exceptions for parallel runs
492 * \throws NumericalProblem on all processes if an exception is thrown during assembly
493 */
494 template<std::size_t i, class AssembleElementFunc>
495 296 void assemble_(Dune::index_constant<i> domainId, AssembleElementFunc&& assembleElement) const
496 {
497 // a state that will be checked on all processes
498 296 bool succeeded = false;
499
500 // try assembling using the local assembly function
501 try
502 {
503 if constexpr (CouplingManagerSupportsMultithreadedAssembly<CouplingManager>::value)
504 {
505 if (enableMultithreading_)
506 {
507 couplingManager_->assembleMultithreaded(
508 domainId, std::forward<AssembleElementFunc>(assembleElement)
509 );
510 return;
511 }
512 }
513
514 // fallback for coupling managers that don't support multithreaded assembly (yet)
515 // or if multithreaded assembly is disabled
516 // let the local assembler add the element contributions
517
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)))
518
1/2
✓ Branch 1 taken 9472 times.
✗ Branch 2 not taken.
 18944 assembleElement(element);
519
520 // if we get here, everything worked well on this process
521 296 succeeded = true;
522 }
523 // throw exception if a problem occurred
524 ✗ catch (NumericalProblem &e)
525 {
526 ✗ std::cout << "rank " << gridView(domainId).comm().rank()
527 ✗ << " caught an exception while assembling:" << e.what()
528 ✗ << "\n";
529 ✗ succeeded = false;
530 }
531
532 // make sure everything worked well on all processes
533
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)
534 ✗ succeeded = gridView(domainId).comm().min(succeeded);
535
536 // if not succeeded rethrow the error on all processes
537
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 148 times.
 296 if (!succeeded)
538 ✗ DUNE_THROW(NumericalProblem, "A process did not succeed in linearizing the system");
539 296 }
540
541 // get diagonal block pattern
542 template<std::size_t i, std::size_t j, typename std::enable_if_t<(i==j), int> = 0>
543 8 Dune::MatrixIndexSet getJacobianPattern_(Dune::index_constant<i> domainI,
544 Dune::index_constant<j> domainJ) const
545 {
546 8 const auto& gg = gridGeometry(domainI);
547
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 16 if (timeSteppingMethod_->implicit())
548 {
549 16 auto pattern = getJacobianPattern<true>(gg);
550 8 couplingManager_->extendJacobianPattern(domainI, pattern);
551 16 return pattern;
552 }
553 else
554 {
555 ✗ auto pattern = getJacobianPattern<false>(gg);
556 ✗ couplingManager_->extendJacobianPattern(domainI, pattern);
557 ✗ return pattern;
558 }
559 }
560
561 // get coupling block pattern
562 template<std::size_t i, std::size_t j, typename std::enable_if_t<(i!=j), int> = 0>
563 4 Dune::MatrixIndexSet getJacobianPattern_(Dune::index_constant<i> domainI,
564 Dune::index_constant<j> domainJ) const
565 {
566
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 8 if (timeSteppingMethod_->implicit())
567 8 return getCouplingJacobianPattern<true>(*couplingManager_,
568 domainI, gridGeometry(domainI),
569 domainJ, gridGeometry(domainJ)
570 16 );
571 else
572 ✗ return getCouplingJacobianPattern<false>(*couplingManager_,
573 domainI, gridGeometry(domainI),
574 domainJ, gridGeometry(domainJ)
575 ✗ );
576 }
577
578 // TODO make this nicer with a is_detected trait in a common location
579 template<class P>
580 ✗ void setProblemTime_(const P& p, const Scalar t)
581 160 { setProblemTimeImpl_(p, t, 0); }
582
583 template<class P>
584 ✗ auto setProblemTimeImpl_(const P& p, const Scalar t, int) -> decltype(p.setTime(0))
585 160 { p.setTime(t); }
586
587 template<class P>
588 void setProblemTimeImpl_(const P& p, const Scalar t, long)
589 {}
590
591 std::shared_ptr<const Experimental::MultiStageMethod<Scalar>> timeSteppingMethod_;
592 std::vector<ResidualType> spatialOperatorEvaluations_;
593 std::vector<ResidualType> temporalOperatorEvaluations_;
594 ResidualType constrainedDofs_;
595 std::shared_ptr<const StageParams> stageParams_;
596
597 //! pointer to the problem to be solved
598 ProblemTuple problemTuple_;
599
600 //! the finite volume geometry of the grid
601 GridGeometryTuple gridGeometryTuple_;
602
603 //! the variables container for the grid
604 GridVariablesTuple gridVariablesTuple_;
605
606 //! Pointer to the previous solution
607 const SolutionVector* prevSol_;
608
609 //! shared pointers to the jacobian matrix and residual
610 std::shared_ptr<JacobianMatrix> jacobian_;
611 std::shared_ptr<ResidualType> residual_;
612
613 //! if multithreaded assembly is enabled
614 bool enableMultithreading_ = false;
615 };
616
617 } // end namespace Dumux
618
619 #endif
620