GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/couplingmapper_staggered_cctpfa.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 62 81 76.5%
Functions: 1 14 7.1%
Branches: 85 177 48.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 FreeFlowPorousMediumCoupling
10 * \copydoc Dumux::FFMomentumPMCouplingMapperStaggeredCCTpfa
11 */
12
13 #ifndef DUMUX_MULTIDOMAIN_FREEFLOWMOMENTUM_POROUSMEDIUM_COUPLINGMAPPER_STAGGERED_TPFA_HH
14 #define DUMUX_MULTIDOMAIN_FREEFLOWMOMENTUM_POROUSMEDIUM_COUPLINGMAPPER_STAGGERED_TPFA_HH
15
16 #include <iostream>
17 #include <unordered_map>
18 #include <tuple>
19 #include <vector>
20
21 #include <dune/common/timer.hh>
22 #include <dune/common/exceptions.hh>
23 #include <dune/common/indices.hh>
24
25 #include <dumux/geometry/intersectingentities.hh>
26 #include <dumux/discretization/method.hh>
27
28 namespace Dumux {
29
30 /*!
31 * \ingroup FreeFlowPorousMediumCoupling
32 * \brief the default mapper for conforming equal dimension boundary coupling between two domains (box or cc)
33 * \todo how to extend to arbitrary number of domains?
34 */
35 template<class MDTraits, class CouplingManager>
36 class FFMomentumPMCouplingMapperStaggeredCCTpfa
37 {
38 using Scalar = typename MDTraits::Scalar;
39
40 template<std::size_t i> using GridGeometry = typename MDTraits::template SubDomain<i>::GridGeometry;
41 template<std::size_t i> using SubControlVolume = typename GridGeometry<i>::SubControlVolume;
42 template<std::size_t i> using SubControlVolumeFace = typename GridGeometry<i>::SubControlVolumeFace;
43 template<std::size_t i> using GridView = typename GridGeometry<i>::GridView;
44 template<std::size_t i> using Element = typename GridView<i>::template Codim<0>::Entity;
45
46 template<std::size_t i>
47 static constexpr auto domainIdx()
48 { return typename MDTraits::template SubDomain<i>::Index{}; }
49
50 template<std::size_t i>
51 static constexpr bool isFcStaggered()
52 { return GridGeometry<i>::discMethod == DiscretizationMethods::fcstaggered; }
53
54 template<std::size_t i>
55 static constexpr bool isCCTpfa()
56 { return GridGeometry<i>::discMethod == DiscretizationMethods::cctpfa; }
57
58 struct ScvfInfoPM
59 {
60 std::size_t eIdxOutside;
61 std::size_t flipScvfIdx;
62 };
63
64 struct ScvfInfoFF
65 {
66 std::size_t eIdxOutside;
67 std::size_t flipScvfIdx;
68 std::size_t dofIdxOutside;
69 };
70
71 using FlipScvfMapTypePM = std::unordered_map<std::size_t, ScvfInfoPM>;
72 using FlipScvfMapTypeFF = std::unordered_map<std::size_t, ScvfInfoFF>;
73 using MapType = std::unordered_map<std::size_t, std::vector<std::size_t>>;
74
75 static constexpr std::size_t numSD = MDTraits::numSubDomains;
76
77 public:
78 /*!
79 * \brief Main update routine
80 */
81 14 void update(const CouplingManager& couplingManager)
82 {
83 static_assert(numSD == 2, "More than two subdomains not implemented!");
84 static_assert(isFcStaggered<0>() && isCCTpfa<1>(), "Only coupling between fcstaggered and cctpfa implemented!");
85
86 14 Dune::Timer watch;
87 28 std::cout << "Initializing the coupling map..." << std::endl;
88
89
2/2
✓ Branch 0 taken 28 times.
✓ Branch 1 taken 14 times.
 42 for (std::size_t domIdx = 0; domIdx < numSD; ++domIdx)
90 56 stencils_[domIdx].clear();
91
92 28 std::get<CouplingManager::freeFlowMomentumIndex>(scvfInfo_).clear();
93 28 std::get<CouplingManager::porousMediumIndex>(scvfInfo_).clear();
94
95
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 14 times.
 14 const auto& freeFlowMomentumProblem = couplingManager.problem(CouplingManager::freeFlowMomentumIndex);
96
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 14 times.
 14 const auto& porousMediumProblem = couplingManager.problem(CouplingManager::porousMediumIndex);
97 14 const auto& freeFlowMomentumGG = freeFlowMomentumProblem.gridGeometry();
98 14 const auto& porousMediumGG = porousMediumProblem.gridGeometry();
99
100 28 isCoupledFFDof_.resize(freeFlowMomentumGG.numScvf(), false);
101 28 isCoupledFFElement_.resize(freeFlowMomentumGG.gridView().size(0), false);
102 42 isCoupledScvf_[CouplingManager::freeFlowMomentumIndex].resize(freeFlowMomentumGG.numScvf(), false);
103 28 isCoupledScvf_[CouplingManager::porousMediumIndex].resize(porousMediumGG.numScvf(), false);
104
105
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 14 auto pmFvGeometry = localView(porousMediumGG);
106
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 14 auto ffFvGeometry = localView(freeFlowMomentumGG);
107
108
4/8
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 135214 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 135200 times.
✗ Branch 10 not taken.
 270428 for (const auto& pmElement : elements(porousMediumGG.gridView()))
109 {
110
1/2
✓ Branch 1 taken 135200 times.
✗ Branch 2 not taken.
 135200 pmFvGeometry.bindElement(pmElement);
111
112
4/4
✓ Branch 0 taken 540800 times.
✓ Branch 1 taken 135200 times.
✓ Branch 2 taken 540800 times.
✓ Branch 3 taken 135200 times.
 811200 for (const auto& pmScvf : scvfs(pmFvGeometry))
113 {
114 // skip all non-boundaries
115
2/2
✓ Branch 0 taken 4640 times.
✓ Branch 1 taken 536160 times.
 540800 if (!pmScvf.boundary())
116 539640 continue;
117
118 // get elements intersecting with the scvf center
119 // for robustness add epsilon in unit outer normal direction
120
1/2
✓ Branch 1 taken 4640 times.
✗ Branch 2 not taken.
 9280 const auto eps = (pmScvf.ipGlobal() - pmFvGeometry.geometry(pmScvf).corner(0)).two_norm()*1e-8;
121 9280 auto globalPos = pmScvf.ipGlobal(); globalPos.axpy(eps, pmScvf.unitOuterNormal());
122
3/8
✓ Branch 1 taken 4640 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4640 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1160 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 11600 const auto indices = intersectingEntities(globalPos, freeFlowMomentumGG.boundingBoxTree());
123
124 // skip if no intersection was found
125
4/4
✓ Branch 0 taken 3480 times.
✓ Branch 1 taken 1160 times.
✓ Branch 2 taken 3480 times.
✓ Branch 3 taken 1160 times.
 9280 if (indices.empty())
126
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3480 times.
 3480 continue;
127
128 // sanity check
129
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1160 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1160 times.
 2320 if (indices.size() > 1)
130 ✗ DUNE_THROW(Dune::InvalidStateException, "Are you sure your sub-domain grids are conformingly discretized on the common interface?");
131
132 // add the pair to the multimap
133 2320 const auto pmElemIdx = porousMediumGG.elementMapper().index(pmElement);
134 1160 const auto ffElemIdx = indices[0];
135
1/2
✓ Branch 1 taken 1160 times.
✗ Branch 2 not taken.
 1160 const auto& ffElement = freeFlowMomentumGG.element(ffElemIdx);
136 1160 ffFvGeometry.bindElement(ffElement);
137
138
4/4
✓ Branch 1 taken 15108 times.
✓ Branch 2 taken 1160 times.
✓ Branch 3 taken 15108 times.
✓ Branch 4 taken 1160 times.
 16268 for (const auto& ffScvf : scvfs(ffFvGeometry)) // TODO: maybe better iterate over all scvs
139 {
140 // TODO this only takes the scv directly at the interface, maybe extend
141
4/4
✓ Branch 0 taken 3564 times.
✓ Branch 1 taken 11544 times.
✓ Branch 2 taken 1188 times.
✓ Branch 3 taken 2376 times.
 15108 if (!ffScvf.boundary() || !ffScvf.isFrontal())
142 continue;
143
144 4752 const auto dist = (pmScvf.ipGlobal() - ffScvf.ipGlobal()).two_norm();
145
2/2
✓ Branch 0 taken 1160 times.
✓ Branch 1 taken 28 times.
 1188 if (dist > eps)
146 continue;
147
148
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1160 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1160 times.
 2320 const auto& ffScv = ffFvGeometry.scv(ffScvf.insideScvIdx());
149
3/6
✓ Branch 1 taken 1160 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1160 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1160 times.
✗ Branch 8 not taken.
 2320 stencils_[CouplingManager::porousMediumIndex][pmElemIdx].push_back(ffScv.dofIndex());
150
3/6
✓ Branch 1 taken 1160 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1160 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1160 times.
✗ Branch 8 not taken.
 2320 stencils_[CouplingManager::freeFlowMomentumIndex][ffScv.dofIndex()].push_back(pmElemIdx);
151
152 // mark the scvf and find and mark the flip scvf
153 3480 isCoupledScvf_[CouplingManager::porousMediumIndex][pmScvf.index()] = true;
154 3480 isCoupledScvf_[CouplingManager::freeFlowMomentumIndex][ffScvf.index()] = true;
155
156 // add all lateral free-flow scvfs touching the coupling interface ("standing" or "lying")
157
4/4
✓ Branch 1 taken 4640 times.
✓ Branch 2 taken 1160 times.
✓ Branch 3 taken 4640 times.
✓ Branch 4 taken 1160 times.
 5800 for (const auto& otherFfScvf : scvfs(ffFvGeometry, ffScv))
158 {
159
2/2
✓ Branch 0 taken 2320 times.
✓ Branch 1 taken 2320 times.
 4640 if (otherFfScvf.isLateral())
160 {
161 // the orthogonal scvf has to lie on the coupling interface
162 2320 const auto& lateralOrthogonalScvf = ffFvGeometry.lateralOrthogonalScvf(otherFfScvf);
163
1/2
✓ Branch 1 taken 2320 times.
✗ Branch 2 not taken.
 2320 isCoupledLateralScvf_[lateralOrthogonalScvf.index()] = true;
164
165 // the "standing" scvf is marked as coupled if it does not lie on a boundary or if that boundary itself is a coupling interface
166
2/2
✓ Branch 0 taken 2292 times.
✓ Branch 1 taken 28 times.
 2320 if (!otherFfScvf.boundary())
167
1/2
✓ Branch 1 taken 2292 times.
✗ Branch 2 not taken.
 2292 isCoupledLateralScvf_[otherFfScvf.index()] = true;
168 else
169 {
170 // for robustness add epsilon in unit outer normal direction
171
1/2
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
 112 const auto otherScvfEps = (otherFfScvf.ipGlobal() - ffFvGeometry.geometry(otherFfScvf).corner(0)).two_norm()*1e-8;
172 56 auto otherScvfGlobalPos = otherFfScvf.center(); otherScvfGlobalPos.axpy(otherScvfEps, otherFfScvf.unitOuterNormal());
173
174
5/10
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 28 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 28 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 28 times.
 56 if (!intersectingEntities(otherScvfGlobalPos, porousMediumGG.boundingBoxTree()).empty())
175 ✗ isCoupledLateralScvf_[otherFfScvf.index()] = true;
176 }
177 }
178 }
179
2/4
✓ Branch 1 taken 1160 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1160 times.
✗ Branch 5 not taken.
 2320 isCoupledFFDof_[ffScv.dofIndex()] = true;
180
2/4
✓ Branch 1 taken 1160 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1160 times.
✗ Branch 5 not taken.
 2320 isCoupledFFElement_[ffElemIdx] = true;
181
182
2/4
✓ Branch 1 taken 1160 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1160 times.
✗ Branch 5 not taken.
 2320 std::get<CouplingManager::porousMediumIndex>(scvfInfo_)[pmScvf.index()] = ScvfInfoPM{ffElemIdx, ffScvf.index()};
183
2/6
✓ Branch 1 taken 1160 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1160 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 2320 std::get<CouplingManager::freeFlowMomentumIndex>(scvfInfo_)[ffScvf.index()] = ScvfInfoFF{pmElemIdx, pmScvf.index(), ffScv.dofIndex()};
184 }
185 }
186 }
187
188
2/4
✓ Branch 3 taken 14 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 14 times.
✗ Branch 8 not taken.
 42 std::cout << "took " << watch.elapsed() << " seconds." << std::endl;
189 14 }
190
191 /*!
192 * \brief returns an iterable container of all indices of degrees of freedom of domain j
193 * that couple with / influence the element residual of the given element of domain i
194 *
195 * \param domainI the domain index of domain i
196 * \param eIdxI the index of the coupled element of domain í
197 * \param domainJ the domain index of domain j
198 *
199 * \note The element residual definition depends on the discretization scheme of domain i
200 * box: a container of the residuals of all sub control volumes
201 * cc : the residual of the (sub) control volume
202 * fem: the residual of the element
203 * \note This function has to be implemented by all coupling managers for all combinations of i and j
204 */
205 ✗ const std::vector<std::size_t>& couplingStencil(Dune::index_constant<CouplingManager::porousMediumIndex> domainI,
206 const std::size_t eIdxI,
207 Dune::index_constant<CouplingManager::freeFlowMomentumIndex> domainJ) const
208 {
209 ✗ if (isCoupledElement(domainI, eIdxI))
210 ✗ return stencils_[CouplingManager::porousMediumIndex].at(eIdxI);
211 else
212 ✗ return emptyStencil_;
213 }
214
215 /*!
216 * \brief returns an iterable container of all indices of degrees of freedom of domain j
217 * that couple with / influence the element residual of the given element of domain i
218 *
219 * \param domainI the domain index of domain i
220 * \param elementI the coupled element of domain i
221 * \param scvI the coupled sub control volume of domain i
222 * \param domainJ the domain index of domain j
223 *
224 * \note The element residual definition depends on the discretization scheme of domain i
225 * box: a container of the residuals of all sub control volumes
226 * cc : the residual of the (sub) control volume
227 * fem: the residual of the element
228 * \note This function has to be implemented by all coupling managers for all combinations of i and j
229 */
230 ✗ const std::vector<std::size_t>& couplingStencil(Dune::index_constant<CouplingManager::freeFlowMomentumIndex> domainI,
231 const Element<CouplingManager::freeFlowMomentumIndex>& elementI,
232 const SubControlVolume<CouplingManager::freeFlowMomentumIndex>& scvI,
233 Dune::index_constant<CouplingManager::porousMediumIndex> domainJ) const
234 {
235 ✗ if (isCoupled(domainI, scvI))
236 ✗ return stencils_[CouplingManager::freeFlowMomentumIndex].at(scvI.dofIndex());
237 else
238 ✗ return emptyStencil_;
239 }
240
241 /*!
242 * \brief Return if an element residual with index eIdx of domain i is coupled to domain j
243 */
244 template<std::size_t i>
245 ✗ bool isCoupledElement(Dune::index_constant<i>, std::size_t eIdx) const
246 {
247 if constexpr (i == CouplingManager::porousMediumIndex)
248
2/4
✓ Branch 4 taken 3160 times.
✓ Branch 5 taken 368040 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 742400 return static_cast<bool>(stencils_[i].count(eIdx));
249 else
250
4/4
✓ Branch 0 taken 181550 times.
✓ Branch 1 taken 368040 times.
✓ Branch 2 taken 181550 times.
✓ Branch 3 taken 368040 times.
 1099180 return isCoupledFFElement_[eIdx];
251 }
252
253 /*!
254 * \brief If the boundary entity is on a coupling boundary
255 * \param domainI the domain index for which to compute the flux
256 * \param scvf the sub control volume face
257 */
258 template<std::size_t i>
259 ✗ bool isCoupled(Dune::index_constant<i> domainI,
260 const SubControlVolumeFace<i>& scvf) const
261 {
262 4750704 return isCoupledScvf_[i].at(scvf.index());
263 }
264
265 /*!
266 * \brief If the boundary entity is on a coupling boundary
267 * \param domainI the domain index for which to compute the flux
268 * \param scvf the sub control volume face
269 */
270 ✗ bool isCoupledLateralScvf(Dune::index_constant<CouplingManager::freeFlowMomentumIndex> domainI,
271 const SubControlVolumeFace<CouplingManager::freeFlowMomentumIndex>& scvf) const
272 ✗ { return isCoupledLateralScvf_.count(scvf.index()); }
273
274 /*!
275 * \brief If the boundary entity is on a coupling boundary
276 * \param domainI the domain index for which to compute the flux
277 * \param scv the sub control volume
278 */
279 ✗ bool isCoupled(Dune::index_constant<CouplingManager::freeFlowMomentumIndex> domainI,
280 const SubControlVolume<CouplingManager::freeFlowMomentumIndex>& scv) const
281 ✗ { return isCoupledFFDof_[scv.dofIndex()]; }
282
283 /*!
284 * \brief Return the scvf index of the flipped scvf in the other domain
285 * \param domainI the domain index for which to compute the flux
286 * \param scvf the sub control volume face
287 */
288 template<std::size_t i>
289 ✗ std::size_t flipScvfIndex(Dune::index_constant<i> domainI,
290 const SubControlVolumeFace<i>& scvf) const
291 {
292
4/10
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 3160 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 181550 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 181550 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 181550 times.
 369420 return std::get<i>(scvfInfo_).at(scvf.index()).flipScvfIdx;
293 }
294
295 /*!
296 * \brief Return the outside element index (the element index of the other domain)
297 * \param domainI the domain index for which to compute the flux
298 * \param scvf the sub control volume face
299 */
300 template<std::size_t i>
301 ✗ std::size_t outsideElementIndex(Dune::index_constant<i> domainI,
302 const SubControlVolumeFace<i>& scvf) const
303 {
304
2/7
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 3160 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 181550 times.
 378900 return std::get<i>(scvfInfo_).at(scvf.index()).eIdxOutside;
305 }
306
307 /*!
308 * \brief Return the outside element index (the element index of the other domain)
309 * \param domainI the domain index for which to compute the flux
310 * \param scvf the sub control volume face
311 */
312 template<std::size_t i>
313 ✗ std::size_t outsideDofIndex(Dune::index_constant<i> domainI,
314 const SubControlVolumeFace<i>& scvf) const
315 {
316 if constexpr (i == CouplingManager::porousMediumIndex)
317 9480 return outsideElementIndex(domainI, scvf);
318 else
319 return std::get<i>(scvfInfo_).at(scvf.index()).dofIdxOutside;
320 }
321
322 private:
323 std::array<MapType, numSD> stencils_;
324 std::vector<std::size_t> emptyStencil_;
325 std::array<std::vector<bool>, numSD> isCoupledScvf_;
326 std::unordered_map<std::size_t, bool> isCoupledLateralScvf_;
327 std::vector<bool> isCoupledFFDof_;
328 std::vector<bool> isCoupledFFElement_;
329 std::tuple<FlipScvfMapTypeFF, FlipScvfMapTypePM> scvfInfo_;
330
331 };
332
333 } // end namespace Dumux
334
335 #endif
336