GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/multidomain/multibinarycouplingmanager.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 52 52 100.0%
Functions: 81 81 100.0%
Branches: 85 124 68.5%
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 * \copydoc Dumux::MultiBinaryCouplingManager
11 */
12
13 #ifndef DUMUX_MULTIDOMAIN_MULTIBINARY_COUPLINGMANAGER_HH
14 #define DUMUX_MULTIDOMAIN_MULTIBINARY_COUPLINGMANAGER_HH
15
16 #include <utility>
17 #include <memory>
18 #include <dune/common/hybridutilities.hh>
19 #include <dumux/common/properties.hh>
20 #include <dumux/multidomain/traits.hh>
21
22 namespace Dumux {
23
24 namespace Detail {
25
26 template <std::size_t, typename Tuple>
27 struct HasIndex;
28
29 template <std::size_t i, typename... Indices>
30 struct HasIndex<i, std::tuple<Indices...>>
31 : std::disjunction<std::is_same<Dune::index_constant<i>, Indices>...>
32 {};
33
34 } // end namespace Detail
35
36 /*!
37 * \ingroup MultiDomain
38 * \brief Coupling manager that combines an arbitrary number of binary coupling manager (coupling two domains each)
39 * \tparam MDTraits the multidomain traits
40 * \tparam CouplingMap a coupling policy class
41 * \tparam CouplingMgrs the binary sub-coupling manager types
42 *
43 * The coupling policy has to provide the interfaces
44 * - CouplingMap::coupledDomains(i): returns a tuple of Dune::index_constants with the coupled domains
45 * - CouplingMap::globalToLocal(i, j): maps the indices i, j to the local index pair of the responsible sub coupling manager
46 * - CouplingMap::managerMap(): returns a two-dimensional array mapping two indices to the coupling manager index
47 */
48 template<class MDTraits, class CouplingMap, class ...CouplingMgrs>
49 class MultiBinaryCouplingManager
50 {
51 // the sub domain type tags
52 template<std::size_t id>
53 using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
54
55 template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
56 template<std::size_t id> using GridView = typename GridGeometry<id>::GridView;
57 template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
58 template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView;
59 template<std::size_t id> using SubControlVolumeFace = typename FVElementGeometry<id>::SubControlVolumeFace;
60 template<std::size_t id> using SubControlVolume = typename FVElementGeometry<id>::SubControlVolume;
61 using CouplingStencil = std::vector<std::size_t>;
62
63 template<std::size_t id>
64 using SubCouplingManagerT = typename std::tuple_element_t<id, std::tuple<CouplingMgrs...>>;
65
66 using CMIndices = std::make_index_sequence<sizeof...(CouplingMgrs)>;
67 using CouplingManagers = typename Detail::MultiDomainTupleSharedPtr<SubCouplingManagerT, CMIndices>::type;
68
69 template<std::size_t id>
70 using SubSolutionVector = std::decay_t<decltype(std::declval<typename MDTraits::SolutionVector>()[Dune::index_constant<id>()])>;
71 using SolutionVectors = typename MDTraits::template TupleOfSharedPtr<SubSolutionVector>;
72
73 static constexpr auto couplingManagerMap_ = CouplingMap::managerMap();
74
75 //! Returns the local domain indices used within the binary sub coupling managers.
76 template<std::size_t i, std::size_t j>
77 static constexpr auto globalToLocal_(Dune::index_constant<i> domainI, Dune::index_constant<j> domainJ)
78 {
79 static_assert(i <= MDTraits::numSubDomains && j <= MDTraits::numSubDomains);
80 return CouplingMap::globalToLocal(domainI, domainJ);
81 }
82
83 //! If two domain are coupled
84 template<class Map, std::size_t i, std::size_t j>
85 static constexpr bool isCoupled_(Dune::index_constant<i> domainI, Dune::index_constant<j>)
86 { return Detail::HasIndex<j, std::decay_t<decltype(Map::coupledDomains(domainI))>>::value; }
87
88 //! Returns the coupling manager index for a given domain combination
89 template<std::size_t i, std::size_t j>
90 static constexpr auto subCouplingManagerIndex_(Dune::index_constant<i> domainI, Dune::index_constant<j> domainJ)
91 {
92 static_assert(
93 isCoupled_<CouplingMap>(domainI, domainJ),
94 "Sub-coupling manager only exists for coupled domains."
95 );
96 return couplingManagerMap_[i][j];
97 }
98
99 public:
100 template<std::size_t i, std::size_t j>
101 using SubCouplingManager = SubCouplingManagerT<couplingManagerMap_[i][j]>;
102
103 21 MultiBinaryCouplingManager()
104 21 {
105 using namespace Dune::Hybrid;
106
1/2
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
 147 forEach(couplingManagers_, [&](auto&& couplingManager)
107 {
108 63 couplingManager = std::make_shared<typename std::decay_t<decltype(couplingManager)>::element_type>();
109 });
110
111
2/3
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 17 times.
✗ Branch 3 not taken.
 96 forEach(solutionVectors_, [&](auto&& solutionVector)
112 {
113 63 solutionVector = std::make_shared<typename std::decay_t<decltype(solutionVector)>::element_type>();
114 });
115 21 }
116
117 //! return the binary sub-coupling manager
118 template<std::size_t i, std::size_t j>
119 49620461 auto& subCouplingManager(Dune::index_constant<i> domainI,
120 Dune::index_constant<j> domainJ)
121 {
122 57133293 constexpr auto idx = subCouplingManagerIndex_(domainI, domainJ);
123 63 return *std::get<idx>(couplingManagers_);
124 }
125
126 //! return the binary sub-coupling manager
127 template<std::size_t i, std::size_t j>
128 188824801 const auto& subCouplingManager(Dune::index_constant<i> domainI,
129 Dune::index_constant<j> domainJ) const
130 {
131 192558175 constexpr auto idx = subCouplingManagerIndex_(domainI, domainJ);
132
7/12
✓ Branch 1 taken 47442 times.
✓ Branch 2 taken 29422 times.
✓ Branch 10 taken 43571 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 57022 times.
✗ Branch 13 not taken.
✓ Branch 3 taken 16764 times.
✓ Branch 14 taken 11703 times.
✗ Branch 15 not taken.
✓ Branch 8 taken 5130 times.
✗ Branch 9 not taken.
✗ Branch 4 not taken.
 173369104 return *std::get<idx>(couplingManagers_);
133 }
134
135 //! apply a function to the domainI-domainJ sub coupling manager using its local indices
136 template<std::size_t i, std::size_t j, class Apply>
137 71073082 decltype(auto) subApply(Dune::index_constant<i> domainI,
138 Dune::index_constant<j> domainJ,
139 Apply&& apply)
140 {
141 constexpr auto localIndices = globalToLocal_(domainI, domainJ);
142 100259283 return apply(subCouplingManager(domainI, domainJ), localIndices.first, localIndices.second);
143 }
144
145 //! apply a function to the domainI-domainJ sub coupling manager using its local indices
146 template<std::size_t i, std::size_t j, class Apply>
147
12/16
✓ Branch 0 taken 11147 times.
✓ Branch 1 taken 3555 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 5187 times.
✓ Branch 9 taken 8412 times.
✓ Branch 10 taken 537836 times.
✓ Branch 11 taken 458 times.
✓ Branch 12 taken 155 times.
✓ Branch 14 taken 3916 times.
✓ Branch 15 taken 3886 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 2408 times.
✓ Branch 13 taken 69 times.
✓ Branch 16 taken 18 times.
✗ Branch 5 not taken.
✗ Branch 8 not taken.
 7419811 decltype(auto) subApply(Dune::index_constant<i> domainI,
148 Dune::index_constant<j> domainJ,
149 const Apply& apply) const
150 {
151 constexpr auto localIndices = globalToLocal_(domainI, domainJ);
152
33/34
✓ Branch 0 taken 11147 times.
✓ Branch 1 taken 208323 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 129652 times.
✓ Branch 6 taken 107328 times.
✓ Branch 7 taken 111115 times.
✓ Branch 9 taken 140216 times.
✓ Branch 10 taken 298682 times.
✓ Branch 14 taken 33347 times.
✓ Branch 15 taken 36520 times.
✓ Branch 16 taken 61836 times.
✓ Branch 17 taken 113227 times.
✓ Branch 18 taken 2396 times.
✓ Branch 19 taken 37351 times.
✓ Branch 21 taken 8480 times.
✓ Branch 22 taken 3840 times.
✓ Branch 24 taken 594409 times.
✓ Branch 25 taken 109016 times.
✓ Branch 27 taken 31945 times.
✓ Branch 28 taken 70149 times.
✓ Branch 12 taken 69916 times.
✓ Branch 13 taken 210132 times.
✓ Branch 20 taken 92973 times.
✓ Branch 23 taken 236 times.
✓ Branch 26 taken 27200 times.
✓ Branch 29 taken 82844 times.
✓ Branch 2 taken 214894 times.
✓ Branch 5 taken 78052 times.
✓ Branch 8 taken 7112 times.
✓ Branch 30 taken 726 times.
✓ Branch 31 taken 240 times.
✓ Branch 33 taken 4164 times.
✓ Branch 34 taken 4164 times.
✓ Branch 11 taken 20720 times.
 7763073 return apply(subCouplingManager(domainI, domainJ), localIndices.first, localIndices.second);
153 }
154
155 //! apply a function to a sub coupling manager containing this domain
156 template<std::size_t i, class Apply>
157 decltype(auto) subApply(Dune::index_constant<i> domainI, Apply&& apply)
158 {
159 constexpr auto dm = CouplingMap::coupledDomains(domainI);
160 static_assert(std::tuple_size_v<std::decay_t<decltype(dm)>> != 0, "subApply on uncoupled domain!");
161 constexpr auto domainJ = std::get<0>(dm);
162 constexpr auto localIndices = globalToLocal_(domainI, domainJ);
163 return apply(subCouplingManager(domainI, domainJ), localIndices.first);
164 }
165
166 //! apply a function to a sub coupling manager containing this domain
167 template<std::size_t i, class Apply>
168
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 178360 times.
 178396 decltype(auto) subApply(Dune::index_constant<i> domainI, const Apply& apply) const
169 {
170 constexpr auto dm = CouplingMap::coupledDomains(domainI);
171 static_assert(std::tuple_size_v<std::decay_t<decltype(dm)>> != 0, "subApply on uncoupled domain!");
172 constexpr auto domainJ = std::get<0>(dm);
173 constexpr auto localIndices = globalToLocal_(domainI, domainJ);
174
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 178360 times.
 178432 return apply(subCouplingManager(domainI, domainJ), localIndices.first);
175 }
176
177 //! Update the solution vector before assembly
178 732 void updateSolution(const typename MDTraits::SolutionVector& curSol)
179 {
180 using namespace Dune::Hybrid;
181 732 forEach(integralRange(Dune::Hybrid::size(solutionVectors_)), [&](const auto id)
182 {
183 732 *std::get<id>(solutionVectors_) = curSol[id];
184 });
185 }
186
187 /*!
188 * \brief extend the jacobian pattern of the diagonal block of domain i
189 * by those entries that are not already in the uncoupled pattern
190 * \note per default we do not add such additional dependencies
191 * \note Such additional dependencies can arise from the coupling, e.g. if a coupling source
192 * term depends on a non-local average of a quantity of the same domain
193 * \warning if you overload this also implement evalAdditionalDomainDerivatives
194 */
195 template<std::size_t id, class JacobianPattern>
196 void extendJacobianPattern(Dune::index_constant<id> domainI, JacobianPattern& pattern) const
197 {}
198
199 /*!
200 * \brief Return the coupling element stencil for a given bulk domain element
201 */
202 template<std::size_t i, class Entity, std::size_t j>
203 2240092 const auto& couplingStencil(Dune::index_constant<i> domainI,
204 const Entity& entity,
205 Dune::index_constant<j> domainJ) const
206 {
207 // if the domains are coupled according to the map, forward to sub-coupling manager
208 if constexpr (isCoupled_<CouplingMap>(domainI, domainJ))
209
7/11
✓ Branch 5 taken 15 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 30 times.
✓ Branch 8 taken 135615 times.
✓ Branch 10 taken 15 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 135630 times.
✓ Branch 13 taken 15 times.
✓ Branch 16 taken 137888 times.
✗ Branch 17 not taken.
✗ Branch 9 not taken.
 2240152 return subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj) -> const auto& {
210
3/6
✓ Branch 5 taken 135615 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 135615 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 137888 times.
✗ Branch 12 not taken.
 2240092 return cm.couplingStencil(ii, entity, jj);
211 });
212 else
213 return emptyStencil_;
214 }
215
216 // ! evaluate coupling residual for the derivative low dim DOF with respect to bulk DOF
217 // ! we only need to evaluate the part of the residual that will be influence by the bulk DOF
218 template<std::size_t i, class LocalAssemblerI, std::size_t j>
219 decltype(auto) evalCouplingResidual(Dune::index_constant<i> domainI,
220 const SubControlVolumeFace<i>& scvfI,
221 const LocalAssemblerI& localAssemblerI,
222 Dune::index_constant<j> domainJ,
223 std::size_t dofIdxGlobalJ) const
224 {
225 if constexpr (i == j || isCoupled_<CouplingMap>(domainI, domainJ))
226 return subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj) -> decltype(auto) {
227 return cm.evalCouplingResidual(ii, scvfI, localAssemblerI, jj, dofIdxGlobalJ);
228 });
229 else
230 {
231 DUNE_THROW(Dune::InvalidStateException,
232 "Calling evalCouplingResidual for uncoupled domains " << i << " and " << j
233 );
234 return localAssemblerI.evalLocalResidual();
235 }
236 }
237
238 //! evaluate coupling residual for the derivative low dim DOF with respect to bulk DOF
239 //! we only need to evaluate the part of the residual that will be influence by the bulk DOF
240 template<std::size_t i, class LocalAssemblerI, std::size_t j>
241 3707788 decltype(auto) evalCouplingResidual(Dune::index_constant<i> domainI,
242 const LocalAssemblerI& localAssemblerI,
243 Dune::index_constant<j> domainJ,
244 std::size_t dofIdxGlobalJ) const
245 {
246 if constexpr (i == j || isCoupled_<CouplingMap>(domainI, domainJ))
247
1/2
✓ Branch 7 taken 1839816 times.
✗ Branch 8 not taken.
 3707788 return subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj) -> decltype(auto) {
248
1/2
✓ Branch 7 taken 1839816 times.
✗ Branch 8 not taken.
 3996508 return cm.evalCouplingResidual(ii, localAssemblerI, jj, dofIdxGlobalJ);
249 288720 });
250 else
251 {
252 DUNE_THROW(Dune::InvalidStateException,
253 "Calling evalCouplingResidual for uncoupled domains " << i << " and " << j
254 );
255 return localAssemblerI.evalLocalResidual();
256 }
257 }
258
259 //! evaluate coupling residual for the derivative low dim DOF with respect to bulk DOF
260 //! we only need to evaluate the part of the residual that will be influence by the bulk DOF
261 template<std::size_t i, class LocalAssemblerI, std::size_t j>
262 5644860 decltype(auto) evalCouplingResidual(Dune::index_constant<i> domainI,
263 const LocalAssemblerI& localAssemblerI,
264 const SubControlVolume<i>& scvI,
265 Dune::index_constant<j> domainJ,
266 std::size_t dofIdxGlobalJ) const
267 {
268 if constexpr (i == j || isCoupled_<CouplingMap>(domainI, domainJ))
269
1/2
✓ Branch 3 taken 2306380 times.
✗ Branch 4 not taken.
 7951240 return subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj) -> decltype(auto) {
270
1/2
✓ Branch 3 taken 2306380 times.
✗ Branch 4 not taken.
 5644860 return cm.evalCouplingResidual(ii, localAssemblerI, scvI, jj, dofIdxGlobalJ);
271 2306380 });
272 else
273 {
274 DUNE_THROW(Dune::InvalidStateException,
275 "Calling evalCouplingResidual for uncoupled domains " << i << " and " << j
276 );
277 return localAssemblerI.evalLocalResidual();
278 }
279 }
280
281 /*!
282 * \brief Update the coupling context for the bulk face residual w.r.t to the lowDim dofs
283 */
284 template<std::size_t i, class LocalAssemblerI, std::size_t j, class PrimaryVariables>
285 41449280 void updateCouplingContext(Dune::index_constant<i> domainI,
286 const LocalAssemblerI& localAssemblerI,
287 Dune::index_constant<j> domainJ,
288 const std::size_t dofIdxGlobalJ,
289 const PrimaryVariables& priVars,
290 int pvIdxJ)
291 {
292 // only one other manager needs to take action if i != j
293 if constexpr (i != j)
294 {
295 if constexpr (isCoupled_<CouplingMap>(domainI, domainJ))
296
2/4
✓ Branch 10 taken 1767420 times.
✗ Branch 11 not taken.
✓ Branch 4 taken 1545600 times.
✗ Branch 5 not taken.
 10367088 subApply(domainI, domainJ, [&](auto& cm, auto&& ii, auto&& jj){
297 10367088 cm.updateCouplingContext(ii, localAssemblerI, jj, dofIdxGlobalJ, priVars, pvIdxJ);
298 });
299 }
300 else
301 {
302 // for i == j, we need to call all relevant managers where domainI is involved and make sure that the
303 // context is updated with respect to its own domain (I-I coupling context)
304 using namespace Dune::Hybrid;
305 60864576 forEach(CouplingMap::coupledDomains(domainI), [&](const auto domainJ){
306 static_assert(domainI != domainJ);
307 37560224 subApply(domainI, domainJ, [&](auto& cm, auto&& ii, auto&& jj){
308 37560224 cm.updateCouplingContext(ii, localAssemblerI, ii, dofIdxGlobalJ, priVars, pvIdxJ);
309 });
310 });
311 }
312 }
313
314 //! Bind the coupling context for a low dim element TODO remove Assembler
315 template<std::size_t i, class Assembler = int>
316 1306497 void bindCouplingContext(Dune::index_constant<i> domainI, const Element<i>& element, const Assembler& assembler = 0)
317 {
318 // for the coupling blocks
319 using namespace Dune::Hybrid;
320 1693086 forEach(CouplingMap::coupledDomains(domainI), [&](const auto domainJ){
321 1306497 subApply(domainI, domainJ, [&](auto& cm, auto&& ii, auto&& jj) -> void {
322 1306497 cm.bindCouplingContext(ii, element, assembler);
323 });
324 });
325 }
326
327 /*!
328 * \brief return the numeric epsilon used for deflecting primary variables of coupled domain i.
329 * \note specialization for free-flow schemes
330 */
331 template<std::size_t i>
332 126 decltype(auto) numericEpsilon(Dune::index_constant<i> domainI,
333 const std::string& paramGroup) const
334 {
335
6/12
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 21 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 21 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 21 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 21 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 21 times.
✗ Branch 17 not taken.
 126 return subApply(domainI, [&](const auto& cm, auto&& ii) -> decltype(auto) {
336
6/12
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 21 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 21 times.
✗ Branch 10 not taken.
✓ Branch 13 taken 21 times.
✗ Branch 14 not taken.
✓ Branch 17 taken 21 times.
✗ Branch 18 not taken.
✓ Branch 21 taken 21 times.
✗ Branch 22 not taken.
 142 return cm.numericEpsilon(ii, paramGroup);
337 16 });
338 }
339
340 /*!
341 * \brief evaluate additional derivatives of the element residual of a domain with respect
342 * to dofs in the same domain that are not in the regular stencil (see CouplingManager::extendJacobianPattern)
343 * \note Such additional dependencies can arise from the coupling, e.g. if a coupling source
344 * term depends on a non-local average of a quantity of the same domain
345 */
346 template<std::size_t i, class LocalAssemblerI, class JacobianMatrixDiagBlock, class GridVariables>
347 void evalAdditionalDomainDerivatives(Dune::index_constant<i> domainI,
348 const LocalAssemblerI& localAssemblerI,
349 const typename LocalAssemblerI::LocalResidual::ElementResidualVector& origResiduals,
350 JacobianMatrixDiagBlock& A,
351 GridVariables& gridVariables)
352 {}
353
354 template<std::size_t i, class LocalAssemblerI, class UpdatableElementVolVars, class UpdatableFluxVarCache>
355 9352648 void updateCoupledVariables(Dune::index_constant<i> domainI,
356 const LocalAssemblerI& localAssemblerI,
357 UpdatableElementVolVars& elemVolVars,
358 UpdatableFluxVarCache& elemFluxVarsCache)
359 {
360 // for the coupling blocks
361 using namespace Dune::Hybrid;
362 9352648 forEach(CouplingMap::coupledDomains(domainI), [&](const auto domainJ){
363 7512832 subApply(domainI, domainJ, [&](auto& cm, auto&& ii, auto&& jj){
364 9352648 cm.updateCoupledVariables(ii, localAssemblerI, elemVolVars, elemFluxVarsCache);
365 });
366 });
367 }
368
369 protected:
370 SolutionVectors& curSol()
371 63 { return solutionVectors_; }
372
373 const SolutionVectors& curSol() const
374 { return solutionVectors_; }
375
376 private:
377 CouplingManagers couplingManagers_;
378 SolutionVectors solutionVectors_;
379
380 CouplingStencil emptyStencil_;
381 };
382
383 } // end namespace Dumux
384
385 #endif
386