GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/multidomain/embedded/pointsourcedata.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 70 71 98.6%
Functions: 14 16 87.5%
Branches: 45 60 75.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-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 EmbeddedCoupling
10 * \brief Data associated with a point source
11 */
12
13 #ifndef DUMUX_MULTIDOMAIN_EMBEDDED_POINTSOURCEDATA_HH
14 #define DUMUX_MULTIDOMAIN_EMBEDDED_POINTSOURCEDATA_HH
15
16 #include <vector>
17 #include <dune/common/fvector.hh>
18 #include <dumux/common/properties.hh>
19 #include <dumux/common/indextraits.hh>
20 #include <dumux/discretization/method.hh>
21
22 namespace Dumux {
23
24 /*!
25 * \ingroup EmbeddedCoupling
26 * \brief A point source data class used for integration in multidimensional models
27 * \note The point source and related data are connected via an identifier (id)
28 */
29 template<class MDTraits>
30
1/2
✓ Branch 1 taken 28038 times.
✗ Branch 2 not taken.
 86988 class PointSourceData
31 {
32 using Scalar = typename MDTraits::Scalar;
33 using ShapeValues = typename std::vector<Dune::FieldVector<Scalar, 1> >;
34
35 template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
36 template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
37 template<std::size_t id> using GridView = typename GridGeometry<id>::GridView;
38 template<std::size_t id> using GridIndex = typename IndexTraits<GridView<id>>::GridIndex;
39 template<std::size_t id> using SolutionVector = GetPropType<SubDomainTypeTag<id>, Properties::SolutionVector>;
40 template<std::size_t id> using PrimaryVariables = GetPropType<SubDomainTypeTag<id>, Properties::PrimaryVariables>;
41
42 static constexpr auto bulkIdx = typename MDTraits::template SubDomain<0>::Index();
43 static constexpr auto lowDimIdx = typename MDTraits::template SubDomain<1>::Index();
44
45 template<std::size_t id>
46 static constexpr bool isBox()
47 { return GridGeometry<id>::discMethod == DiscretizationMethods::box; }
48
49 public:
50 4494 void addBulkInterpolation(const ShapeValues& shapeValues,
51 const std::vector<GridIndex<bulkIdx>>& cornerIndices,
52 GridIndex<bulkIdx> eIdx)
53 {
54 static_assert(isBox<bulkIdx>(), "This interface is only available for the box method.");
55 4494 bulkElementIdx_ = eIdx;
56
1/2
✓ Branch 1 taken 4494 times.
✗ Branch 2 not taken.
 4494 bulkCornerIndices_ = cornerIndices;
57
1/2
✓ Branch 1 taken 4494 times.
✗ Branch 2 not taken.
 4494 bulkShapeValues_ = shapeValues;
58 4494 }
59
60 236 void addLowDimInterpolation(const ShapeValues& shapeValues,
61 const std::vector<GridIndex<lowDimIdx>>& cornerIndices,
62 GridIndex<lowDimIdx> eIdx)
63 {
64 static_assert(isBox<lowDimIdx>(), "This interface is only available for the box method.");
65 236 lowDimElementIdx_ = eIdx;
66
1/2
✓ Branch 1 taken 236 times.
✗ Branch 2 not taken.
 236 lowDimCornerIndices_ = cornerIndices;
67
1/2
✓ Branch 1 taken 236 times.
✗ Branch 2 not taken.
 236 lowDimShapeValues_ = shapeValues;
68 236 }
69
70 82494 void addBulkInterpolation(GridIndex<bulkIdx> eIdx)
71 {
72 static_assert(!isBox<bulkIdx>(), "This interface is not available for the box method.");
73
1/2
✓ Branch 1 taken 82494 times.
✗ Branch 2 not taken.
 82494 bulkElementIdx_ = eIdx;
74 }
75
76 86752 void addLowDimInterpolation(GridIndex<lowDimIdx> eIdx)
77 {
78 static_assert(!isBox<lowDimIdx>(), "This interface is not available for the box method.");
79
2/3
✓ Branch 1 taken 30204 times.
✗ Branch 2 not taken.
✓ Branch 0 taken 48640 times.
 86752 lowDimElementIdx_ = eIdx;
80 }
81
82 314123547 PrimaryVariables<bulkIdx> interpolateBulk(const SolutionVector<bulkIdx>& sol) const
83 {
84 314123547 PrimaryVariables<bulkIdx> bulkPriVars(0.0);
85 if (isBox<bulkIdx>())
86 {
87
2/2
✓ Branch 0 taken 55668796 times.
✓ Branch 1 taken 13917199 times.
 69585995 for (int i = 0; i < bulkCornerIndices_.size(); ++i)
88
2/2
✓ Branch 0 taken 55668796 times.
✓ Branch 1 taken 55668796 times.
 111337592 for (int priVarIdx = 0; priVarIdx < bulkPriVars.size(); ++priVarIdx)
89 55668796 bulkPriVars[priVarIdx] += sol[bulkCornerIndices_[i]][priVarIdx]*bulkShapeValues_[i];
90 }
91 else
92 {
93 300206348 bulkPriVars = sol[bulkElementIdx()];
94 }
95 13917199 return bulkPriVars;
96 }
97
98 328245343 PrimaryVariables<lowDimIdx> interpolateLowDim(const SolutionVector<lowDimIdx>& sol) const
99 {
100 328245343 PrimaryVariables<lowDimIdx> lowDimPriVars(0.0);
101 if (isBox<lowDimIdx>())
102 {
103
2/2
✓ Branch 0 taken 135360 times.
✓ Branch 1 taken 67680 times.
 203040 for (int i = 0; i < lowDimCornerIndices_.size(); ++i)
104
2/2
✓ Branch 0 taken 135360 times.
✓ Branch 1 taken 135360 times.
 270720 for (int priVarIdx = 0; priVarIdx < lowDimPriVars.size(); ++priVarIdx)
105 135360 lowDimPriVars[priVarIdx] += sol[lowDimCornerIndices_[i]][priVarIdx]*lowDimShapeValues_[i];
106 }
107 else
108 {
109 328177663 lowDimPriVars = sol[lowDimElementIdx()];
110 }
111 67680 return lowDimPriVars;
112 }
113
114
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 138475 times.
 921992113 GridIndex<lowDimIdx> lowDimElementIdx() const
115
4/4
✓ Branch 0 taken 670560 times.
✓ Branch 1 taken 9416769 times.
✓ Branch 2 taken 33754 times.
✓ Branch 3 taken 4511722 times.
 314077831 { return lowDimElementIdx_; }
116
117 300206348 GridIndex<bulkIdx> bulkElementIdx() const
118
4/4
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 1018726 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 464198 times.
 300206348 { return bulkElementIdx_; }
119
120 private:
121 ShapeValues bulkShapeValues_, lowDimShapeValues_;
122 std::vector<GridIndex<bulkIdx>> bulkCornerIndices_;
123 std::vector<GridIndex<lowDimIdx>> lowDimCornerIndices_;
124 GridIndex<bulkIdx> bulkElementIdx_;
125 GridIndex<lowDimIdx> lowDimElementIdx_;
126 };
127
128 /*!
129 * \ingroup EmbeddedCoupling
130 * \brief A point source data class used for integration in multidimensional models
131 * \note The point source and related data are connected via an identifier (id)
132 * When explicitly computing the circle average, i.e. the pressure for
133 * the source term is computed as an integral over the circle describing
134 * the surface of the one-dimensional tube. This exact determination of the bulk pressure
135 * is necessary for convergence studies.
136 */
137 template<class MDTraits>
138 class PointSourceDataCircleAverage : public PointSourceData<MDTraits>
139 {
140 using ParentType = PointSourceData<MDTraits>;
141 using Scalar = typename MDTraits::Scalar;
142 using ShapeValues = typename std::vector<Dune::FieldVector<Scalar, 1> >;
143
144 template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
145 template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
146 template<std::size_t id> using GridView = typename GridGeometry<id>::GridView;
147 template<std::size_t id> using GridIndex = typename IndexTraits<GridView<id>>::GridIndex;
148 template<std::size_t id> using SolutionVector = GetPropType<SubDomainTypeTag<id>, Properties::SolutionVector>;
149 template<std::size_t id> using PrimaryVariables = GetPropType<SubDomainTypeTag<id>, Properties::PrimaryVariables>;
150
151 static constexpr auto bulkIdx = typename MDTraits::template SubDomain<0>::Index();
152 static constexpr auto lowDimIdx = typename MDTraits::template SubDomain<1>::Index();
153
154 template<std::size_t id>
155 static constexpr bool isBox()
156 { return GridGeometry<id>::discMethod == DiscretizationMethods::box; }
157
158 public:
159
2/3
✓ Branch 1 taken 10310 times.
✗ Branch 2 not taken.
✓ Branch 0 taken 48640 times.
 58950 PointSourceDataCircleAverage() : enableBulkCircleInterpolation_(false) {}
160
161 14121796 PrimaryVariables<bulkIdx> interpolateBulk(const SolutionVector<bulkIdx>& sol) const
162 {
163 // bulk interpolation on the circle is only enabled for source in the
164 // lower dimensional domain if we use a circle distributed bulk sources
165 // (see coupling manager circle sources)
166
1/2
✓ Branch 0 taken 10597764 times.
✗ Branch 1 not taken.
 14121796 PrimaryVariables<bulkIdx> bulkPriVars(sol[0]);
167
1/2
✓ Branch 0 taken 10597764 times.
✗ Branch 1 not taken.
 14121796 bulkPriVars = 0.0;
168
1/2
✓ Branch 0 taken 14121796 times.
✗ Branch 1 not taken.
 14121796 if (enableBulkCircleInterpolation_)
169 {
170 // compute the average of the bulk privars over the circle around the integration point
171 // this computes $\bar{p} = \frac{1}{2\pi R} int_0^2*\pi p R \text{d}\theta.
172 if (isBox<bulkIdx>())
173 {
174
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4379284 times.
 4379284 assert(circleCornerIndices_.size() == circleShapeValues_.size());
175
176 Scalar weightSum = 0.0;
177
2/2
✓ Branch 0 taken 71220086 times.
✓ Branch 1 taken 4379284 times.
 75599370 for (std::size_t j = 0; j < circleStencil_.size(); ++j)
178 {
179 71220086 PrimaryVariables<bulkIdx> priVars(0.0);
180 71220086 const auto& cornerIndices = *(circleCornerIndices_[j]);
181 71220086 const auto& shapeValues = circleShapeValues_[j];
182
2/2
✓ Branch 0 taken 569760688 times.
✓ Branch 1 taken 71220086 times.
 640980774 for (int i = 0; i < cornerIndices.size(); ++i)
183 {
184 569760688 const auto& localSol = sol[cornerIndices[i]];
185 569760688 const auto& shapeValue = shapeValues[i];
186
2/2
✓ Branch 0 taken 569760688 times.
✓ Branch 1 taken 569760688 times.
 1139521376 for (int priVarIdx = 0; priVarIdx < PrimaryVariables<bulkIdx>::size(); ++priVarIdx)
187 569760688 priVars[priVarIdx] += localSol[priVarIdx]*shapeValue;
188 }
189 // multiply with weight and add
190 71220086 priVars *= circleIpWeight_[j];
191 71220086 weightSum += circleIpWeight_[j];
192 71220086 bulkPriVars += priVars;
193 }
194 4379284 bulkPriVars /= weightSum;
195 }
196 else
197 {
198 Scalar weightSum = 0.0;
199
2/2
✓ Branch 0 taken 902061139 times.
✓ Branch 1 taken 9742512 times.
 911803651 for (int j = 0; j < circleStencil_.size(); ++j)
200 {
201
2/2
✓ Branch 0 taken 937277183 times.
✓ Branch 1 taken 902061139 times.
 1839338322 for (int priVarIdx = 0; priVarIdx < bulkPriVars.size(); ++priVarIdx)
202 937277183 bulkPriVars[priVarIdx] += sol[circleStencil_[j]][priVarIdx]*circleIpWeight_[j];
203
204 902061139 weightSum += circleIpWeight_[j];
205 }
206 9742512 bulkPriVars /= weightSum;
207 }
208 14121796 return bulkPriVars;
209 }
210 else
211 {
212 ✗ return ParentType::interpolateBulk(sol);
213 }
214 }
215
216 2242 void addCircleInterpolation(const std::vector<const std::vector<GridIndex<bulkIdx>>*>& circleCornerIndices,
217 const std::vector<ShapeValues>& circleShapeValues,
218 const std::vector<Scalar>& circleIpWeight,
219 const std::vector<GridIndex<bulkIdx>>& circleStencil)
220 {
221 2242 circleCornerIndices_ = circleCornerIndices;
222 2242 circleShapeValues_ = circleShapeValues;
223 2242 circleIpWeight_ = circleIpWeight;
224 2242 circleStencil_ = circleStencil;
225 2242 enableBulkCircleInterpolation_ = true;
226
227 2242 }
228
229 56708 void addCircleInterpolation(const std::vector<Scalar>& circleIpWeight,
230 const std::vector<GridIndex<bulkIdx>>& circleStencil)
231 {
232
1/2
✓ Branch 1 taken 56708 times.
✗ Branch 2 not taken.
 56708 circleIpWeight_ = circleIpWeight;
233
1/2
✓ Branch 1 taken 56708 times.
✗ Branch 2 not taken.
 56708 circleStencil_ = circleStencil;
234
1/2
✓ Branch 1 taken 8068 times.
✗ Branch 2 not taken.
 56708 enableBulkCircleInterpolation_ = true;
235 48640 }
236
237 const std::vector<GridIndex<bulkIdx>>& circleStencil() const
238 { return circleStencil_; }
239
240 private:
241 std::vector<const std::vector<GridIndex<bulkIdx>>*> circleCornerIndices_;
242 std::vector<ShapeValues> circleShapeValues_;
243 std::vector<Scalar> circleIpWeight_;
244 std::vector<GridIndex<bulkIdx>> circleStencil_;
245 bool enableBulkCircleInterpolation_;
246 };
247
248 } // end namespace Dumux
249
250 #endif
251