GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/io/grid/porenetwork/parametersforgeneratedgrid.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 239 269 88.8%
Functions: 54 60 90.0%
Branches: 391 2000 19.6%
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 PoreNetworkModels
10 * \brief Helper class to assign parameters to a generated grid
11 */
12 #ifndef DUMUX_IO_PARAMETERS_FOR_GENERATED_GRID
13 #define DUMUX_IO_PARAMETERS_FOR_GENERATED_GRID
14
15 #include <algorithm>
16 #include <array>
17 #include <numeric>
18 #include <random>
19 #include <string>
20 #include <vector>
21
22 #include <dune/common/exceptions.hh>
23 #include <dune/grid/common/exceptions.hh>
24 #include <dune/geometry/axisalignedcubegeometry.hh>
25
26 #include <dumux/common/parameters.hh>
27 #include <dumux/common/random.hh>
28 #include <dumux/common/stringutilities.hh>
29 #include <dumux/geometry/intersectspointgeometry.hh>
30 #include <dumux/porenetwork/common/throatproperties.hh>
31 #include <dumux/porenetwork/common/poreproperties.hh>
32
33 namespace Dumux::PoreNetwork {
34
35 /*!
36 * \ingroup PoreNetworkModels
37 * \brief Helper class to assign parameters to a generated grid
38 */
39 template <class Grid, class Scalar>
40 25 class ParametersForGeneratedGrid
41 {
42 using GridView = typename Grid::LeafGridView;
43 using Element = typename Grid::template Codim<0>::Entity;
44 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
45 using Vertex = typename Grid::template Codim<Grid::dimension>::Entity;
46
47 static constexpr auto dim = Grid::dimension;
48 static constexpr auto dimWorld = Grid::dimensionworld;
49 using BoundaryList = std::array<int, 2*dimWorld>;
50
51 public:
52
53 37 ParametersForGeneratedGrid(const GridView& gridView, const std::string& paramGroup)
54 37 : gridView_(gridView)
55 37 , paramGroup_(paramGroup)
56 37 , priorityList_(getPriorityList_())
57 {
58
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 37 computeBoundingBox_();
59
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 37 boundaryFaceIndex_ = getBoundaryFacemarkerInput_();
60 37 }
61
62 /*!
63 * \brief Returns the boundary face marker index at given position
64 *
65 * \param pos The current position
66 */
67 17291 int boundaryFaceMarkerAtPos(const GlobalPosition& pos) const
68 {
69 // set the priority which decides the order the vertices on the boundary are indexed
70 // by default, vertices on min/max faces in x direction have the highest priority, followed by y and z
71
2/4
✓ Branch 0 taken 89274 times.
✓ Branch 1 taken 12767 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 103291 for (auto boundaryIdx : priorityList_)
72 {
73
2/4
✓ Branch 0 taken 4192 times.
✓ Branch 1 taken 85082 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 90402 if (onBoundary_(pos, boundaryIdx))
74 4402 return boundaryFaceIndex_[boundaryIdx];
75 }
76 return -1;
77 }
78
79 /*!
80 * \brief Computes and returns the label of a given throat
81 *
82 * \param poreLabels The labels of the pores adjacent to the throat.
83 */
84 44293 int throatLabel(const std::array<int, 2>& poreLabels) const
85 {
86
2/2
✓ Branch 0 taken 34035 times.
✓ Branch 1 taken 9804 times.
 44293 if (poreLabels[0] == poreLabels[1]) // both vertices are inside the domain or on the same boundary face
87 return poreLabels[0];
88
2/2
✓ Branch 0 taken 5394 times.
✓ Branch 1 taken 4410 times.
 10030 if (poreLabels[0] == -1) // vertex1 is inside the domain, vertex2 is on a boundary face
89 return poreLabels[1];
90
2/2
✓ Branch 0 taken 4487 times.
✓ Branch 1 taken 907 times.
 5536 if (poreLabels[1] == -1) // vertex2 is inside the domain, vertex1 is on a boundary face
91 return poreLabels[0];
92
93 // use the priority list to find out which pore label is favored
94
1/2
✓ Branch 0 taken 1193 times.
✗ Branch 1 not taken.
 1247 for(const auto i : priorityList_)
95 {
96
2/2
✓ Branch 0 taken 489 times.
✓ Branch 1 taken 704 times.
 1247 if (poreLabels[0] == boundaryFaceIndex_[i])
97 return poreLabels[0];
98
2/2
✓ Branch 0 taken 286 times.
✓ Branch 1 taken 418 times.
 736 if (poreLabels[1] == boundaryFaceIndex_[i])
99 return poreLabels[1];
100 }
101
102 ✗ DUNE_THROW(Dune::InvalidStateException, "Something went wrong with the throat labels");
103 }
104
105 /*!
106 * \brief Assign parameters for generically created grids
107 */
108 template <class SetParameter, class GetParameter>
109 37 void assignParameters(const SetParameter& setParameter,
110 const GetParameter& getParameter,
111 const std::size_t numSubregions)
112 {
113 using cytpe = typename GlobalPosition::value_type;
114 using InternalBoundingBox = Dune::AxisAlignedCubeGeometry<cytpe, dimWorld, dimWorld>;
115 37 std::vector<InternalBoundingBox> internalBoundingBoxes;
116
117 // divide the network into subregions, if specified
118
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 24 times.
 37 if (numSubregions > 0)
119 {
120 // get bounding boxes of subregions
121
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 for (int i = 0; i < numSubregions; ++i)
122 {
123
3/6
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
 4 auto lowerLeft = getParamFromGroup<GlobalPosition>(paramGroup_, "Grid.Subregion" + std::to_string(i) + ".LowerLeft");
124
3/6
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
 4 auto upperRight = getParamFromGroup<GlobalPosition>(paramGroup_, "Grid.Subregion" + std::to_string(i) + ".UpperRight");
125
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 internalBoundingBoxes.emplace_back(std::move(lowerLeft), std::move(upperRight));
126 }
127 }
128
129 // get the maximum possible pore body radii such that pore bodies do not intersect
130 // (requires ThroatRegionId, if specified)
131
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 37 const std::vector<Scalar> maxPoreRadius = getMaxPoreRadii_(numSubregions, getParameter);
132
2/6
✓ Branch 2 taken 25 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 25 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 37 std::vector<bool> poreRadiusLimited(gridView_.size(dim), false);
133
134 // get a helper function for getting the pore radius of a pore body not belonging to a subregion
135
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 37 auto defaultPoreRadius = poreRadiusGenerator_(-1);
136
137 // get helper functions for pore body radii on subregions
138 37 std::vector<decltype(defaultPoreRadius)> subregionPoreRadius;
139
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 25 times.
 38 for (int i = 0; i < numSubregions; ++i)
140
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 2 subregionPoreRadius.emplace_back(poreRadiusGenerator_(i));
141
142 // get helper function for pore volume
143
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 37 const auto poreVolume = poreVolumeGenerator_(getParameter);
144
145 // treat the pore body parameters (label, radius and maybe regionId)
146
3/4
✓ Branch 2 taken 16959 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 16959 times.
✓ Branch 5 taken 25 times.
 17368 for (const auto& vertex : vertices(gridView_))
147 {
148
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 16959 times.
 17331 const auto& pos = vertex.geometry().center();
149
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 16959 times.
✓ Branch 3 taken 16919 times.
✗ Branch 4 not taken.
 17331 const auto vIdxGlobal = gridView_.indexSet().index(vertex);
150 17331 const auto poreLabel = boundaryFaceMarkerAtPos(pos);
151
1/2
✓ Branch 1 taken 16959 times.
✗ Branch 2 not taken.
 34662 setParameter(vertex, "PoreLabel", poreLabel);
152
153 // sets the minimum of the given value and the maximum possible pore body radius
154 // and keeps track of capped radii
155 17828 auto setRadiusAndLogIfCapped = [&](const Scalar value)
156 {
157
4/6
✓ Branch 0 taken 5188 times.
✓ Branch 1 taken 11564 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 138 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 194 times.
 17084 if (value > maxPoreRadius[vIdxGlobal])
158 {
159 5188 poreRadiusLimited[vIdxGlobal] = true;
160
1/6
✓ Branch 2 taken 5188 times.
✗ Branch 3 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
 10376 setParameter(vertex, "PoreInscribedRadius", maxPoreRadius[vIdxGlobal]);
161 }
162 else
163
3/6
✓ Branch 3 taken 11564 times.
✗ Branch 4 not taken.
✓ Branch 9 taken 138 times.
✗ Branch 10 not taken.
✓ Branch 15 taken 194 times.
✗ Branch 16 not taken.
 35688 setParameter(vertex, "PoreInscribedRadius", value);
164 };
165
166
2/2
✓ Branch 0 taken 15959 times.
✓ Branch 1 taken 1000 times.
 17331 if (numSubregions == 0) // assign radius if no subregions are specified
167
1/2
✓ Branch 1 taken 15959 times.
✗ Branch 2 not taken.
 16331 setRadiusAndLogIfCapped(defaultPoreRadius(vertex, poreLabel));
168 else // assign region ids and radii to vertices if they are within a subregion
169 {
170 // default value for vertices not belonging to a subregion
171
2/4
✓ Branch 1 taken 1000 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1000 times.
 2000 setParameter(vertex, "PoreRegionId", -1);
172
1/2
✓ Branch 1 taken 1000 times.
✗ Branch 2 not taken.
 1000 setRadiusAndLogIfCapped(defaultPoreRadius(vertex, poreLabel));
173
174
2/2
✓ Branch 0 taken 1000 times.
✓ Branch 1 taken 1000 times.
 2000 for (int id = 0; id < numSubregions; ++id)
175 {
176
1/2
✓ Branch 1 taken 1000 times.
✗ Branch 2 not taken.
 1000 const auto& subregion = internalBoundingBoxes[id];
177
2/3
✗ Branch 0 not taken.
✓ Branch 1 taken 125 times.
✓ Branch 2 taken 875 times.
 1000 if (intersectsPointGeometry(vertex.geometry().center(), subregion))
178 {
179
1/2
✓ Branch 1 taken 125 times.
✗ Branch 2 not taken.
 250 setParameter(vertex, "PoreRegionId", id);
180
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 125 times.
✓ Branch 3 taken 125 times.
✗ Branch 4 not taken.
 250 setRadiusAndLogIfCapped(subregionPoreRadius[id](vertex, poreLabel));
181 }
182 }
183 }
184
185
2/4
✓ Branch 1 taken 16959 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16959 times.
✗ Branch 5 not taken.
 34662 setParameter(vertex, "PoreVolume", poreVolume(vertex, poreLabel));
186 }
187
188 // treat throat parameters
189
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 37 auto defaultThroatInscribedRadius = throatInscribedRadiusGenerator_(-1, getParameter);
190
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 37 auto defaultThroatLength = throatLengthGenerator_(-1, getParameter);
191
192 // get helper functions for throat radii and lengths on subregions
193 37 std::vector<decltype(defaultThroatInscribedRadius)> subregionThroatInscribedRadius;
194 37 std::vector<decltype(defaultThroatLength)> subregionThroatLength;
195
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 25 times.
 38 for (int i = 0; i < numSubregions; ++i)
196 {
197
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 subregionThroatInscribedRadius.emplace_back(throatInscribedRadiusGenerator_(i, getParameter));
198
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 subregionThroatLength.emplace_back(throatLengthGenerator_(i, getParameter));
199 }
200
201 // set throat parameters
202
2/2
✓ Branch 1 taken 43839 times.
✓ Branch 2 taken 25 times.
 44330 for (const auto& element : elements(gridView_))
203 {
204
2/2
✓ Branch 0 taken 33363 times.
✓ Branch 1 taken 10476 times.
 44293 if (numSubregions == 0) // assign values if no subregions are specified
205 {
206
2/4
✓ Branch 1 taken 33363 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 33363 times.
✗ Branch 5 not taken.
 67634 setParameter(element, "ThroatInscribedRadius", defaultThroatInscribedRadius(element));
207
2/4
✓ Branch 1 taken 33363 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 33363 times.
✗ Branch 5 not taken.
 67634 setParameter(element, "ThroatLength", defaultThroatLength(element));
208 }
209 else // assign values to throats if they are within a subregion
210 {
211 // default value for elements not belonging to a subregion
212
1/2
✓ Branch 1 taken 10476 times.
✗ Branch 2 not taken.
 20952 setParameter(element, "ThroatRegionId", -1);
213
2/4
✓ Branch 1 taken 10476 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10476 times.
✗ Branch 5 not taken.
 20952 setParameter(element, "ThroatInscribedRadius", defaultThroatInscribedRadius(element));
214
2/4
✓ Branch 1 taken 10476 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10476 times.
✗ Branch 5 not taken.
 20952 setParameter(element, "ThroatLength", defaultThroatLength(element));
215
216
2/2
✓ Branch 0 taken 10476 times.
✓ Branch 1 taken 10476 times.
 20952 for (int id = 0; id < numSubregions; ++id)
217 {
218
1/2
✓ Branch 1 taken 10476 times.
✗ Branch 2 not taken.
 10476 const auto& subregion = internalBoundingBoxes[id];
219
2/3
✗ Branch 0 not taken.
✓ Branch 1 taken 1625 times.
✓ Branch 2 taken 8851 times.
 20952 if (intersectsPointGeometry(element.geometry().center(), subregion))
220 {
221
1/2
✓ Branch 1 taken 1625 times.
✗ Branch 2 not taken.
 3250 setParameter(element, "ThroatRegionId", id);
222
2/4
✓ Branch 1 taken 1625 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1625 times.
✗ Branch 5 not taken.
 3250 setParameter(element, "ThroatInscribedRadius", subregionThroatInscribedRadius[id](element));
223
2/4
✓ Branch 1 taken 1625 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1625 times.
✗ Branch 5 not taken.
 3250 setParameter(element, "ThroatLength", subregionThroatLength[id](element));
224 }
225 }
226 }
227
228 // set the throat label
229 44293 const auto vertex0 = element.template subEntity<dim>(0);
230
1/2
✓ Branch 1 taken 43839 times.
✗ Branch 2 not taken.
 44293 const auto vertex1 = element.template subEntity<dim>(1);
231
2/4
✓ Branch 1 taken 43839 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 43839 times.
✗ Branch 5 not taken.
 88586 const std::array poreLabels{static_cast<int>(getParameter(vertex0, "PoreLabel")),
232
1/2
✓ Branch 1 taken 43839 times.
✗ Branch 2 not taken.
 88586 static_cast<int>(getParameter(vertex1, "PoreLabel"))};
233
2/4
✓ Branch 1 taken 43839 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 43839 times.
✗ Branch 5 not taken.
 88586 setParameter(element, "ThroatLabel", throatLabel(poreLabels));
234 }
235
236
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 23 times.
 37 const auto numPoreRadiusLimited = std::count(poreRadiusLimited.begin(), poreRadiusLimited.end(), true);
237
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 23 times.
 37 if (numPoreRadiusLimited > 0)
238
5/10
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
 2 std::cout << "*******\nWarning! " << numPoreRadiusLimited << " out of " << poreRadiusLimited.size()
239
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 << " pore body radii have been capped automatically in order to prevent intersecting pores\n*******" << std::endl;
240 148 }
241
242 private:
243
244 /*!
245 * \brief Returns a list of boundary face priorities from user specified input or default values if no input is given
246 *
247 * This essentially determines the index of a node on an edge or corner. For instance in a 2D case, a list of {0,1,2,3}
248 * will give highest priority to the "x-direction" and lowest to the "diagonal-direction".
249 */
250 25 BoundaryList getPriorityList_() const
251 {
252 46 const auto list = [&]()
253 {
254 BoundaryList priorityList;
255 50 std::iota(priorityList.begin(), priorityList.end(), 0);
256
257
7/12
✓ Branch 2 taken 15 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 9 times.
✓ Branch 6 taken 6 times.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✓ Branch 14 taken 2 times.
✓ Branch 18 taken 8 times.
✗ Branch 19 not taken.
✗ Branch 21 not taken.
✓ Branch 22 taken 8 times.
 50 if (hasParamInGroup(paramGroup_, "Grid.PriorityList"))
258 {
259 try {
260 // priorities can also be set in the input file
261
1/6
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 9 priorityList = getParamFromGroup<BoundaryList>(paramGroup_, "Grid.PriorityList");
262 }
263 // make sure that a priority for each direction is set
264 ✗ catch(Dune::RangeError& e) {
265 ✗ DUNE_THROW(Dumux::ParameterException, "You must specify priorities for all directions (" << dimWorld << ") \n" << e.what());
266 }
267 // make sure each direction is only set once
268
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 9 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 9 if (!isUnique_(priorityList))
269
1/72
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ 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 28 not taken.
✗ Branch 29 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 35 not taken.
✓ Branch 36 taken 9 times.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
✗ Branch 47 not taken.
✗ Branch 48 not taken.
✗ Branch 50 not taken.
✗ Branch 51 not taken.
✗ Branch 53 not taken.
✗ Branch 54 not taken.
✗ Branch 56 not taken.
✗ Branch 57 not taken.
✗ Branch 59 not taken.
✗ Branch 60 not taken.
✗ Branch 62 not taken.
✗ Branch 63 not taken.
✗ Branch 65 not taken.
✗ Branch 66 not taken.
✗ Branch 68 not taken.
✗ Branch 69 not taken.
✗ Branch 72 not taken.
✗ Branch 73 not taken.
✗ Branch 75 not taken.
✗ Branch 76 not taken.
✗ Branch 81 not taken.
✗ Branch 82 not taken.
✗ Branch 84 not taken.
✗ Branch 85 not taken.
✗ Branch 87 not taken.
✗ Branch 88 not taken.
✗ Branch 90 not taken.
✗ Branch 91 not taken.
✗ Branch 93 not taken.
✗ Branch 94 not taken.
✗ Branch 96 not taken.
✗ Branch 97 not taken.
✗ Branch 99 not taken.
✗ Branch 100 not taken.
✗ Branch 102 not taken.
✗ Branch 103 not taken.
✗ Branch 105 not taken.
✗ Branch 106 not taken.
✗ Branch 108 not taken.
✗ Branch 109 not taken.
✗ Branch 112 not taken.
✗ Branch 113 not taken.
✗ Branch 115 not taken.
✗ Branch 116 not taken.
 9 DUNE_THROW(Dumux::ParameterException, "You must specify priorities for all directions (duplicate directions)");
270
271 //make sure that the directions are correct (ranging from 0 to dimWorld-1)
272
7/48
✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 9 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 9 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 9 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 3 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 9 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
 24 if (std::any_of(priorityList.begin(), priorityList.end(), []( const int i ){ return (i < 0 || i >= 2*dimWorld); }))
273 ✗ DUNE_THROW(Dumux::ParameterException, "You must specify priorities for correct directions (0-" << 2*dimWorld-1 << ")");
274 }
275 25 return priorityList;
276
3/6
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 8 times.
✗ Branch 8 not taken.
 25 }();
277
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 16 return list;
278 }
279
280 37 void computeBoundingBox_()
281 {
282 // calculate the bounding box of the local partition of the grid view
283
2/2
✓ Branch 2 taken 17113 times.
✓ Branch 3 taken 25 times.
 35007 for (const auto& vertex : vertices(gridView_))
284 {
285
2/2
✓ Branch 0 taken 51046 times.
✓ Branch 1 taken 17113 times.
 69373 for (int i = 0; i < dimWorld; i++)
286 {
287 using std::min;
288 using std::max;
289
2/2
✓ Branch 0 taken 71 times.
✓ Branch 1 taken 50975 times.
 51991 bBoxMin_[i] = min(bBoxMin_[i], vertex.geometry().corner(0)[i]);
290
2/2
✓ Branch 0 taken 334 times.
✓ Branch 1 taken 50712 times.
 52348 bBoxMax_[i] = max(bBoxMax_[i], vertex.geometry().corner(0)[i]);
291 }
292 }
293 37 }
294
295 /*!
296 * \brief Returns a list of boundary face indices from user specified input or default values if no input is given
297 */
298 37 BoundaryList getBoundaryFacemarkerInput_() const
299 {
300 BoundaryList boundaryFaceMarker;
301 37 std::fill(boundaryFaceMarker.begin(), boundaryFaceMarker.end(), 0);
302 37 boundaryFaceMarker[0] = 1;
303 37 boundaryFaceMarker[1] = 1;
304
305
3/4
✓ Branch 2 taken 25 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 24 times.
✓ Branch 6 taken 1 times.
 74 if (hasParamInGroup(paramGroup_, "Grid.BoundaryPoreLabels"))
306 {
307 36 const auto input = getParamFromGroup<std::vector<std::string>>(paramGroup_, "Grid.BoundaryPoreLabels");
308
3/4
✓ Branch 1 taken 108 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 108 times.
✓ Branch 6 taken 24 times.
 192 for (const auto& entry : input)
309 {
310
1/2
✓ Branch 1 taken 108 times.
✗ Branch 2 not taken.
 156 const std::string errorMessage = "You must specify BoundaryPoreLabels in the format pos:num, where pos can be xMin, xMax, yMin, yMax, zMin, zMax and num is the corresponding label.\n"
311 "Example (2D, defaults are used for the remaining boundaries): xMin:2 yMax:3\n";
312
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 108 times.
 156 if (entry.find(':') == std::string::npos)
313 ✗ DUNE_THROW(Dumux::ParameterException, errorMessage);
314
315
6/10
✓ Branch 0 taken 24 times.
✓ Branch 1 taken 84 times.
✓ Branch 3 taken 24 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 24 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 144 times.
✓ Branch 11 taken 24 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 408 static const std::map<std::string, int> labels = {{"xMin", 0}, {"xMax", 1}, {"yMin", 2}, {"yMax", 3}, {"zMin", 4}, {"zMax", 5}};
316
2/4
✓ Branch 1 taken 108 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 108 times.
✗ Branch 5 not taken.
 156 const auto splitEntry = split(entry, ":");
317
2/4
✓ Branch 1 taken 108 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 108 times.
✗ Branch 5 not taken.
 312 const std::string location = std::string(splitEntry[0].begin(), splitEntry[0].end());
318
1/2
✓ Branch 1 taken 108 times.
✗ Branch 2 not taken.
 156 int value = 0;
319 try {
320
3/6
✓ Branch 1 taken 108 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 108 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 108 times.
 312 value = std::stoi(std::string(splitEntry[1].begin(), splitEntry[1].end()));
321 }
322 ✗ catch(...) {
323 ✗ DUNE_THROW(Dumux::ParameterException, errorMessage);
324 }
325
326
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 108 times.
 156 if (splitEntry.size() != 2)
327
1/26
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ 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 28 not taken.
✗ Branch 29 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 38 not taken.
✓ Branch 39 taken 108 times.
 156 DUNE_THROW(Dumux::ParameterException, errorMessage);
328 ✗ if (!labels.count(location))
329 ✗ DUNE_THROW(Dumux::ParameterException, errorMessage);
330 else
331
1/2
✓ Branch 1 taken 108 times.
✗ Branch 2 not taken.
 156 boundaryFaceMarker[labels.at(location)] = value;
332 }
333 36 }
334 37 return boundaryFaceMarker;
335 72 }
336
337 // returns the maximum possible pore body radii such that pore bodies do not intersect
338 template <class GetParameter>
339 37 std::vector<Scalar> getMaxPoreRadii_(std::size_t numSubregions, const GetParameter& getParameter) const
340 {
341 37 const auto numVertices = gridView_.size(dim);
342
1/2
✓ Branch 2 taken 25 times.
✗ Branch 3 not taken.
 37 std::vector<Scalar> maxPoreRadius(numVertices, std::numeric_limits<Scalar>::max());
343
344
2/4
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 25 times.
✗ Branch 4 not taken.
 37 if (!getParamFromGroup<bool>(paramGroup_, "Grid.CapPoreRadii", true))
345 return maxPoreRadius;
346
347 // check for a user-specified fixed throat length
348
1/4
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 37 const Scalar inputThroatLength = getParamFromGroup<Scalar>(paramGroup_, "Grid.ThroatLength", -1.0);
349 37 std::vector<Scalar> subregionInputThroatLengths;
350
351
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 24 times.
 37 if (numSubregions > 0)
352 {
353
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 for (int i = 0; i < numSubregions; ++i)
354 {
355 // adapt the parameter group if there are subregions
356
2/6
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 2 const std::string paramGroup = paramGroup_ + ".SubRegion" + std::to_string(i);
357
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 const Scalar subregionInputThroatLength = getParamFromGroup<Scalar>(paramGroup, "Grid.ThroatLength", -1.0);
358
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 subregionInputThroatLengths.push_back(subregionInputThroatLength);
359 }
360 }
361
362
2/2
✓ Branch 2 taken 43839 times.
✓ Branch 3 taken 25 times.
 88623 for (const auto& element : elements(gridView_))
363 {
364 // do not cap the pore radius if a user-specified throat length if given
365
2/2
✓ Branch 0 taken 10476 times.
✓ Branch 1 taken 33363 times.
 44293 if (numSubregions > 0)
366 {
367 // check subregions
368
1/2
✓ Branch 1 taken 10476 times.
✗ Branch 2 not taken.
 20952 const auto subregionId = getParameter(element, "ThroatRegionId");
369
1/2
✓ Branch 0 taken 10476 times.
✗ Branch 1 not taken.
 10476 if (subregionId >= 0)
370 {
371 // throat lies within a subregion
372
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10476 times.
 10476 if (subregionInputThroatLengths[subregionId] > 0.0)
373 246 continue;
374 }
375 ✗ else if (inputThroatLength > 0.0) // throat does not lie within a subregion
376 ✗ continue;
377 }
378
2/2
✓ Branch 0 taken 123 times.
✓ Branch 1 taken 33240 times.
 33817 else if (inputThroatLength > 0.0) // no subregions
379 246 continue;
380
381 // No fixed throat lengths given, check for max. pore radius.
382 // We define this as 1/2 of the length (minus a user specified value) of the shortest pore throat attached to the pore body.
383
2/2
✓ Branch 0 taken 22 times.
✓ Branch 1 taken 43694 times.
 44047 const Scalar delta = element.geometry().volume();
384
4/6
✓ Branch 0 taken 22 times.
✓ Branch 1 taken 43694 times.
✓ Branch 3 taken 22 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 22 times.
✗ Branch 7 not taken.
 44047 static const Scalar minThroatLength = getParamFromGroup<Scalar>(paramGroup_, "Grid.MinThroatLength", 1e-6);
385 44047 const Scalar maxRadius = (delta - minThroatLength)/2.0;
386
2/2
✓ Branch 0 taken 87432 times.
✓ Branch 1 taken 43716 times.
 132141 for (int vIdxLocal = 0; vIdxLocal < 2 ; ++vIdxLocal)
387 {
388
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 87432 times.
 88094 const int vIdxGlobal = gridView_.indexSet().subIndex(element, vIdxLocal, dim);
389
2/2
✓ Branch 0 taken 24200 times.
✓ Branch 1 taken 63232 times.
 112563 maxPoreRadius[vIdxGlobal] = std::min(maxPoreRadius[vIdxGlobal], maxRadius);
390 }
391 }
392
393 return maxPoreRadius;
394 37 }
395
396 // returns a function taking a vertex and a pore label and returning a radius
397 38 std::function<Scalar(const Vertex&, const int)> poreRadiusGenerator_(const int subregionId) const
398 {
399 // adapt the parameter name if there are subregions
400
9/16
✓ Branch 0 taken 25 times.
✓ Branch 1 taken 1 times.
✓ Branch 3 taken 25 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✓ Branch 14 taken 25 times.
✓ Branch 15 taken 1 times.
✓ Branch 16 taken 25 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
 41 const std::string prefix = subregionId < 0 ? "Grid." : "Grid.Subregion" + std::to_string(subregionId) + ".";
401
402 // prepare random number generation for lognormal parameter distribution
403 38 std::mt19937 generator;
404
405 // check if pores for certain labels should be treated in a special way
406
2/4
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26 times.
✗ Branch 5 not taken.
 76 const auto poreLabelsToSetFixedRadius = getParamFromGroup<std::vector<int>>(paramGroup_, prefix + "PoreLabelsToSetFixedRadius", std::vector<int>{});
407
2/4
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26 times.
✗ Branch 5 not taken.
 76 const auto poreLabelsToApplyFactorForRadius = getParamFromGroup<std::vector<int>>(paramGroup_, prefix + "PoreLabelsToApplyFactorForRadius", std::vector<int>{});
408
3/6
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 26 times.
✗ Branch 9 not taken.
 76 const auto poreRadiusForLabel = getParamFromGroup<std::vector<Scalar>>(paramGroup_, prefix + "FixedPoreRadiusForLabel", std::vector<Scalar>{});
409
3/8
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 26 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 76 const auto poreRadiusFactorForLabel = getParamFromGroup<std::vector<Scalar>>(paramGroup_, prefix + "PoreRadiusFactorForLabel", std::vector<Scalar>{});
410
411 78 const auto generateFunction = [&](auto& poreRadiusDist)
412 {
413
18/72
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 12 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 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 24 not taken.
✗ Branch 25 not taken.
✓ Branch 30 taken 4 times.
✗ Branch 31 not taken.
✓ Branch 33 taken 4 times.
✗ Branch 34 not taken.
✓ Branch 36 taken 4 times.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✓ Branch 44 taken 1 times.
✗ Branch 45 not taken.
✓ Branch 47 taken 1 times.
✗ Branch 48 not taken.
✓ Branch 50 taken 1 times.
✗ Branch 51 not taken.
✗ Branch 52 not taken.
✗ Branch 53 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 66 not taken.
✗ Branch 67 not taken.
✓ Branch 72 taken 1 times.
✗ Branch 73 not taken.
✓ Branch 75 taken 1 times.
✗ Branch 76 not taken.
✓ Branch 78 taken 1 times.
✗ Branch 79 not taken.
✗ Branch 80 not taken.
✗ Branch 81 not taken.
✓ Branch 86 taken 7 times.
✗ Branch 87 not taken.
✓ Branch 89 taken 7 times.
✗ Branch 90 not taken.
✓ Branch 92 taken 7 times.
✗ Branch 93 not taken.
✗ Branch 94 not taken.
✗ Branch 95 not taken.
✗ Branch 100 not taken.
✗ Branch 101 not taken.
✗ Branch 103 not taken.
✗ Branch 104 not taken.
✗ Branch 106 not taken.
✗ Branch 107 not taken.
✗ Branch 108 not taken.
✗ Branch 109 not taken.
✓ Branch 114 taken 1 times.
✗ Branch 115 not taken.
✓ Branch 117 taken 1 times.
✗ Branch 118 not taken.
✓ Branch 120 taken 1 times.
✗ Branch 121 not taken.
✗ Branch 122 not taken.
✗ Branch 123 not taken.
 33822 return [=](const auto& vertex, const int poreLabel) mutable
414 {
415 15489 const auto radius = poreRadiusDist(generator);
416
417 // check if pores for certain labels should be treated in a special way
418
9/30
✓ Branch 0 taken 1356 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 15396 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 15396 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 69 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 69 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 69 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 170 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✓ Branch 26 taken 24 times.
✗ Branch 27 not taken.
✓ Branch 28 taken 24 times.
✗ Branch 29 not taken.
 17084 if (poreLabelsToSetFixedRadius.empty() && poreLabelsToApplyFactorForRadius.empty())
419 return radius; // nothing special to be done
420
421 // set a fixed radius for a given label
422 ✗ else if (!poreLabelsToSetFixedRadius.empty() || !poreRadiusForLabel.empty())
423 {
424 ✗ if (poreLabelsToSetFixedRadius.size() != poreRadiusForLabel.size())
425 ✗ DUNE_THROW(Dumux::ParameterException, "PoreLabelsToSetFixedRadius must be of same size as FixedPoreRadiusForLabel");
426
427 ✗ if (const auto it = std::find(poreLabelsToSetFixedRadius.begin(), poreLabelsToSetFixedRadius.end(), poreLabel); it != poreLabelsToSetFixedRadius.end())
428 ✗ return poreRadiusForLabel[std::distance(poreLabelsToSetFixedRadius.begin(), it)];
429 }
430
431 // multiply the pore radius by a given value for a given label
432 ✗ else if (!poreLabelsToApplyFactorForRadius.empty() || !poreRadiusFactorForLabel.empty())
433 {
434 ✗ if (poreLabelsToApplyFactorForRadius.size() != poreRadiusFactorForLabel.size())
435 ✗ DUNE_THROW(Dumux::ParameterException, "PoreLabelsToApplyFactorForRadius must be of same size as PoreRadiusFactorForLabel");
436
437 ✗ if (const auto it = std::find(poreLabelsToApplyFactorForRadius.begin(), poreLabelsToApplyFactorForRadius.end(), poreLabel); it != poreLabelsToApplyFactorForRadius.end())
438 ✗ return poreRadiusFactorForLabel[std::distance(poreLabelsToApplyFactorForRadius.begin(), it)] * radius;
439 }
440
441 // default
442 return radius;
443 26 };
444 };
445
446
2/4
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26 times.
✗ Branch 5 not taken.
 76 const Scalar fixedPoreRadius = getParamFromGroup<Scalar>(paramGroup_, prefix + "PoreInscribedRadius", -1.0);
447 // return random radius according to a user-specified distribution
448
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 20 times.
 38 if (fixedPoreRadius <= 0.0)
449 {
450 // allow to specify a seed to get reproducible results
451
3/6
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6 times.
✗ Branch 8 not taken.
 27 const auto seed = getParamFromGroup<unsigned int>(paramGroup_, prefix + "ParameterRandomNumberSeed", std::random_device{}());
452 9 generator.seed(seed);
453
454
2/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
 9 const auto type = getParamFromGroup<std::string>(paramGroup_, prefix + "ParameterType", "lognormal");
455
1/2
✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
 9 if (type == "lognormal")
456 {
457 // if we use a lognormal distribution, get the mean and standard deviation from input file
458
5/10
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 6 times.
✗ Branch 11 not taken.
✓ Branch 15 taken 6 times.
✗ Branch 16 not taken.
 18 const auto [meanPoreRadius, stddevPoreRadius] = getDistributionInputParams_("Lognormal", prefix,
459 "MeanPoreInscribedRadius",
460 "StandardDeviationPoreInscribedRadius");
461 9 const Scalar variance = stddevPoreRadius*stddevPoreRadius;
462
463 using std::log;
464 using std::sqrt;
465 9 const Scalar mu = log(meanPoreRadius/sqrt(1.0 + variance/(meanPoreRadius*meanPoreRadius)));
466
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 9 const Scalar sigma = sqrt(log(1.0 + variance/(meanPoreRadius*meanPoreRadius)));
467
468 9 Dumux::SimpleLogNormalDistribution<> poreRadiusDist(mu, sigma);
469
2/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
 9 return generateFunction(poreRadiusDist);
470 }
471 ✗ else if (type == "uniform")
472 {
473 // if we use a uniform distribution, get the min and max from input file
474 ✗ const auto [minPoreRadius, maxPoreRadius] = getDistributionInputParams_("Uniform", prefix,
475 "MinPoreInscribedRadius",
476 "MaxPoreInscribedRadius");
477 ✗ Dumux::SimpleUniformDistribution<> poreRadiusDist(minPoreRadius, maxPoreRadius);
478 ✗ return generateFunction(poreRadiusDist);
479 }
480 else
481
0/24
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ 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 28 not taken.
✗ Branch 29 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
 9 DUNE_THROW(Dune::InvalidStateException, "Unknown parameter type " << type);
482 9 }
483
484 // always return a fixed constant radius
485 else
486 {
487
3/6
✓ Branch 0 taken 1356 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 69 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 170 times.
✗ Branch 5 not taken.
 1613 auto poreRadiusDist = [fixedPoreRadius](auto& gen){ return fixedPoreRadius; };
488
2/6
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 29 return generateFunction(poreRadiusDist);
489 }
490 38 }
491
492 // print helpful error message if params are not properly provided
493
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 9 std::array<Scalar, 2> getDistributionInputParams_(const std::string& distributionName,
494 const std::string& prefix,
495 const std::string& paramName0,
496 const std::string& paramName1) const
497 {
498 try
499 {
500
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 9 return std::array{getParamFromGroup<Scalar>(paramGroup_, prefix + paramName0),
501
2/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
 18 getParamFromGroup<Scalar>(paramGroup_, prefix + paramName1)};
502 }
503 ✗ catch (const Dumux::ParameterException& e)
504 {
505 std::cout << "\n" << distributionName << " pore-size distribution needs input parameters "
506 << prefix + paramName0 << " and " << prefix + paramName1 << ".\n"
507 ✗ << "Alternatively, use " << prefix << "PoreInscribedRadius to set a fixed inscribed pore radius." << std::endl;
508 ✗ DUNE_THROW(Dumux::ParameterException, e.what());
509 }
510 }
511
512 template <class GetParameter>
513 37 auto poreVolumeGenerator_(const GetParameter& getParameter) const
514 {
515
1/2
✓ Branch 2 taken 25 times.
✗ Branch 3 not taken.
 37 const auto geometry = Pore::shapeFromString(getParamFromGroup<std::string>(paramGroup_, "Grid.PoreGeometry"));
516
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 37 const auto capPoresOnBoundaries = getParamFromGroup<std::vector<int>>(paramGroup_, "Grid.CapPoresOnBoundaries", std::vector<int>{});
517
2/4
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 25 times.
 37 if (!isUnique_(capPoresOnBoundaries))
518 ✗ DUNE_THROW(Dune::InvalidStateException, "CapPoresOnBoundaries must not contain duplicates");
519
520 // automatically determine the pore volume if not provided by the grid file
521
2/4
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13 times.
✗ Branch 5 not taken.
 434 return [=,this] (const auto& vertex, const auto vIdx)
522 {
523 33918 const Scalar r = getParameter(vertex, "PoreInscribedRadius");
524 17331 const Scalar volume = [&]
525 {
526
4/6
✓ Branch 0 taken 40 times.
✓ Branch 1 taken 16587 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 138 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 194 times.
 16959 if (geometry == Pore::Shape::cylinder)
527 {
528
4/18
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 38 times.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✗ 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 20 not taken.
✗ Branch 21 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
 40 static const Scalar fixedHeight = getParamFromGroup<Scalar>(paramGroup_, "Grid.PoreHeight", -1.0);
529
2/18
✗ Branch 0 not taken.
✓ Branch 1 taken 40 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 40 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
 40 const Scalar h = fixedHeight > 0.0 ? fixedHeight : getParameter(vertex, "PoreHeight");
530 40 return Pore::volume(geometry, r, h);
531 }
532 else
533 16919 return Pore::volume(geometry, r);
534 }();
535
536
5/6
✓ Branch 0 taken 14778 times.
✓ Branch 1 taken 1849 times.
✓ Branch 2 taken 69 times.
✓ Branch 3 taken 69 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 194 times.
 16959 if (capPoresOnBoundaries.empty())
537 return volume;
538 else
539 {
540 14847 std::size_t numCaps = 0;
541 14847 Scalar factor = 1.0;
542 14847 const auto& pos = vertex.geometry().center();
543
4/6
✓ Branch 0 taken 88668 times.
✓ Branch 1 taken 14778 times.
✓ Branch 2 taken 345 times.
✓ Branch 3 taken 69 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 103860 for (auto boundaryIdx : capPoresOnBoundaries)
544 {
545
4/6
✓ Branch 0 taken 3218 times.
✓ Branch 1 taken 85450 times.
✓ Branch 2 taken 48 times.
✓ Branch 3 taken 297 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 89013 if (onBoundary_(pos, boundaryIdx))
546 {
547 3266 factor *= 0.5;
548 3266 ++numCaps;
549 }
550 }
551
552
2/6
✗ Branch 0 not taken.
✓ Branch 1 taken 14778 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 69 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 14847 if (numCaps > dimWorld)
553 ✗ DUNE_THROW(Dune::InvalidStateException, "Pore " << vIdx << " at " << pos << " capped " << numCaps << " times. Capping should not happen more than " << dimWorld << " times");
554
555 14847 return volume * factor;
556 }
557 37 };
558 37 }
559
560 // returns a lambda taking a element and returning a radius
561 template <class GetParameter>
562 38 auto throatInscribedRadiusGenerator_(const int subregionId, const GetParameter& getParameter) const
563 {
564 // adapt the parameter name if there are subregions
565
10/18
✓ Branch 0 taken 25 times.
✓ Branch 1 taken 1 times.
✓ Branch 3 taken 25 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✓ Branch 14 taken 25 times.
✓ Branch 15 taken 1 times.
✓ Branch 16 taken 25 times.
✓ Branch 18 taken 26 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
 41 const std::string prefix = subregionId < 0 ? "Grid." : "Grid.Subregion" + std::to_string(subregionId) + ".";
566
567 // check for a user-specified fixed throat radius
568
2/4
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26 times.
✗ Branch 5 not taken.
 76 const Scalar inputThroatInscribedRadius = getParamFromGroup<Scalar>(paramGroup_, prefix + "ThroatInscribedRadius", -1.0);
569
570 // shape parameter for calculation of throat radius
571
2/4
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26 times.
✗ Branch 5 not taken.
 76 const Scalar throatN = getParamFromGroup<Scalar>(paramGroup_, prefix + "ThroatInscribedRadiusN", 0.1);
572
573 45464 return [=](const Element& element)
574 {
575
6/6
✓ Branch 0 taken 43615 times.
✓ Branch 1 taken 1433 times.
✓ Branch 2 taken 69 times.
✓ Branch 3 taken 69 times.
✓ Branch 4 taken 243 times.
✓ Branch 5 taken 35 times.
 45464 const Scalar delta = element.geometry().volume();
576 45464 const std::array<Vertex, 2> vertices = {element.template subEntity<dim>(0), element.template subEntity<dim>(1)};
577
578 // the element parameters (throat radius and length)
579
6/6
✓ Branch 0 taken 43615 times.
✓ Branch 1 taken 1433 times.
✓ Branch 2 taken 69 times.
✓ Branch 3 taken 69 times.
✓ Branch 4 taken 243 times.
✓ Branch 5 taken 35 times.
 45464 if (inputThroatInscribedRadius > 0.0)
580 return inputThroatInscribedRadius;
581 else
582 {
583
3/6
✓ Branch 2 taken 43615 times.
✗ Branch 3 not taken.
✓ Branch 8 taken 69 times.
✗ Branch 9 not taken.
✓ Branch 14 taken 243 times.
✗ Branch 15 not taken.
 43927 const Scalar poreRadius0 = getParameter(vertices[0], "PoreInscribedRadius");
584
3/6
✓ Branch 2 taken 43615 times.
✗ Branch 3 not taken.
✓ Branch 7 taken 69 times.
✗ Branch 8 not taken.
✓ Branch 12 taken 243 times.
✗ Branch 13 not taken.
 43927 const Scalar poreRadius1 = getParameter(vertices[1], "PoreInscribedRadius");
585 43927 return Throat::averagedRadius(poreRadius0, poreRadius1, delta, throatN);
586 }
587 38 };
588 38 }
589
590 // returns a lambda taking a element and returning a length
591 template <class GetParameter>
592 38 auto throatLengthGenerator_(int subregionId, const GetParameter& getParameter) const
593 {
594 // adapt the parameter name if there are subregions
595
10/18
✓ Branch 0 taken 25 times.
✓ Branch 1 taken 1 times.
✓ Branch 3 taken 25 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✓ Branch 14 taken 25 times.
✓ Branch 15 taken 1 times.
✓ Branch 16 taken 25 times.
✓ Branch 18 taken 26 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
 41 const std::string prefix = subregionId < 0 ? "Grid." : "Grid.Subregion" + std::to_string(subregionId) + ".";
596
2/4
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26 times.
✗ Branch 5 not taken.
 76 const Scalar inputThroatLength = getParamFromGroup<Scalar>(paramGroup_, prefix + "ThroatLength", -1.0);
597 // decide whether to subtract the pore radii from the throat length or not
598
2/4
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26 times.
✗ Branch 5 not taken.
 76 const bool subtractRadiiFromThroatLength = getParamFromGroup<bool>(paramGroup_, prefix + "SubtractPoreInscribedRadiiFromThroatLength", true);
599
600 45464 return [=,this](const Element& element)
601 {
602
6/6
✓ Branch 0 taken 45029 times.
✓ Branch 1 taken 19 times.
✓ Branch 2 taken 69 times.
✓ Branch 3 taken 69 times.
✓ Branch 4 taken 243 times.
✓ Branch 5 taken 35 times.
 45464 if (inputThroatLength > 0.0)
603 return inputThroatLength;
604
605
3/6
✓ Branch 0 taken 45029 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 69 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 243 times.
✗ Branch 5 not taken.
 45341 const Scalar delta = element.geometry().volume();
606
3/6
✓ Branch 0 taken 45029 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 69 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 243 times.
✗ Branch 5 not taken.
 45341 if (subtractRadiiFromThroatLength)
607 {
608 45341 const std::array<Vertex, 2> vertices = {element.template subEntity<dim>(0), element.template subEntity<dim>(1)};
609
9/18
✓ Branch 2 taken 45029 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 45029 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 45029 times.
✗ Branch 9 not taken.
✓ Branch 16 taken 69 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 69 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 69 times.
✗ Branch 23 not taken.
✓ Branch 30 taken 243 times.
✗ Branch 31 not taken.
✓ Branch 33 taken 243 times.
✗ Branch 34 not taken.
✓ Branch 36 taken 243 times.
✗ Branch 37 not taken.
 45341 const Scalar result = delta - getParameter(vertices[0], "PoreInscribedRadius") - getParameter(vertices[1], "PoreInscribedRadius");
610
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 45029 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 69 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 243 times.
 45341 if (result <= 0.0)
611
0/78
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ 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 27 not taken.
✗ Branch 28 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 36 not taken.
✗ Branch 37 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 69 not taken.
✗ Branch 70 not taken.
✗ Branch 72 not taken.
✗ Branch 73 not taken.
✗ Branch 75 not taken.
✗ Branch 76 not taken.
✗ Branch 78 not taken.
✗ Branch 79 not taken.
✗ Branch 85 not taken.
✗ Branch 86 not taken.
✗ Branch 88 not taken.
✗ Branch 89 not taken.
✗ Branch 91 not taken.
✗ Branch 92 not taken.
✗ Branch 94 not taken.
✗ Branch 95 not taken.
✗ Branch 97 not taken.
✗ Branch 98 not taken.
✗ Branch 100 not taken.
✗ Branch 101 not taken.
✗ Branch 103 not taken.
✗ Branch 104 not taken.
✗ Branch 106 not taken.
✗ Branch 107 not taken.
✗ Branch 109 not taken.
✗ Branch 110 not taken.
✗ Branch 111 not taken.
✗ Branch 112 not taken.
✗ Branch 114 not taken.
✗ Branch 115 not taken.
✗ Branch 117 not taken.
✗ Branch 118 not taken.
✗ Branch 120 not taken.
✗ Branch 121 not taken.
 45341 DUNE_THROW(Dune::GridError, "Pore radii are so large they intersect! Something went wrong at element " << gridView_.indexSet().index(element));
612 else
613 return result;
614 }
615 else
616 return delta;
617 38 };
618 38 }
619
620 178287 bool onBoundary_(const GlobalPosition& pos, std::size_t boundaryIdx) const
621 {
622
6/8
✓ Branch 0 taken 43296 times.
✓ Branch 1 taken 45332 times.
✓ Branch 2 taken 44472 times.
✓ Branch 3 taken 44541 times.
✓ Branch 4 taken 254 times.
✓ Branch 5 taken 314 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 178209 constexpr auto eps = 1e-8; //TODO
623
624 static constexpr std::array boundaryIdxToCoordinateIdx{0, 0, 1, 1, 2, 2};
625 178287 const auto coordinateIdx = boundaryIdxToCoordinateIdx[boundaryIdx];
626 178287 auto isMaxBoundary = [] (int n) { return (n % 2 != 0); };
627
628
10/14
✓ Branch 0 taken 43076 times.
✓ Branch 1 taken 45070 times.
✓ Branch 2 taken 44334 times.
✓ Branch 3 taken 44334 times.
✓ Branch 4 taken 258 times.
✓ Branch 5 taken 302 times.
✓ Branch 6 taken 138 times.
✓ Branch 7 taken 207 times.
✓ Branch 8 taken 254 times.
✓ Branch 9 taken 314 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 178287 if (isMaxBoundary(boundaryIdx))
629 88060 return pos[coordinateIdx] > bBoxMax_[coordinateIdx] - eps;
630 else
631 90227 return pos[coordinateIdx] < bBoxMin_[coordinateIdx] + eps;
632 }
633
634 // check if a container is unique
635 // copy the container such that is does not get altered by std::sort
636 template <class T>
637 34 bool isUnique_(T t) const
638 {
639 34 std::sort(t.begin(), t.end());
640 34 return (std::unique(t.begin(), t.end()) == t.end());
641 }
642
643 const GridView gridView_;
644
645 const std::string paramGroup_;
646 BoundaryList priorityList_;
647 BoundaryList boundaryFaceIndex_;
648
649 //! the bounding box of the whole domain
650 GlobalPosition bBoxMin_ = GlobalPosition(std::numeric_limits<double>::max());
651 GlobalPosition bBoxMax_ = GlobalPosition(std::numeric_limits<double>::min());
652 };
653
654 } // namespace Dumux::PoreNetwork
655
656 #endif
657