GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/porousmediumflow/2p/gridadaptindicator.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 53 56 94.6%
Functions: 20 20 100.0%
Branches: 161 237 67.9%
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 TwoPModel
10 * \brief Class defining a standard, saturation dependent indicator for grid adaptation.
11 */
12
13 #ifndef DUMUX_TWOP_ADAPTION_INDICATOR_HH
14 #define DUMUX_TWOP_ADAPTION_INDICATOR_HH
15
16 #include <memory>
17 #include <dune/common/exceptions.hh>
18 #include <dune/grid/common/partitionset.hh>
19
20 #include <dumux/common/properties.hh>
21 #include <dumux/common/parameters.hh>
22 #include <dumux/discretization/elementsolution.hh>
23 #include <dumux/discretization/evalsolution.hh>
24
25 namespace Dumux {
26
27 /*!
28 * \ingroup TwoPModel
29 * \brief Class defining a standard, saturation dependent indicator for grid adaptation.
30 */
31 template<class TypeTag>
32 class TwoPGridAdaptIndicator
33 {
34 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
35 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
36 using Element = typename GridView::template Codim<0>::Entity;
37 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
38 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
39 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
40
41 enum { saturationIdx = Indices::saturationIdx };
42
43 public:
44 /*!
45 * \brief The Constructor
46 *
47 * \param gridGeometry The finite volume grid geometry
48 * \param paramGroup The parameter group in which to look for runtime parameters first (default is "")
49 *
50 * Note: refineBound_, coarsenBound_ & maxSaturationDelta_ are chosen
51 * in a way such that the indicator returns false for all elements
52 * before having been calculated.
53 */
54 5 TwoPGridAdaptIndicator(std::shared_ptr<const GridGeometry> gridGeometry, const std::string& paramGroup = "")
55 : gridGeometry_(gridGeometry)
56 , refineBound_(std::numeric_limits<Scalar>::max())
57 , coarsenBound_(std::numeric_limits<Scalar>::lowest())
58
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 5 , maxSaturationDelta_(gridGeometry_->gridView().size(0), 0.0)
59 10 , minLevel_(getParamFromGroup<std::size_t>(paramGroup, "Adaptive.MinLevel", 0))
60
7/20
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 5 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 5 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 5 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 5 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 4 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✓ Branch 20 taken 4 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
 5 , maxLevel_(getParamFromGroup<std::size_t>(paramGroup, "Adaptive.MaxLevel", 0))
61 5 {}
62
63 /*!
64 * \brief Function to set the minimum allowed level.
65 *
66 * \param minLevel The minimum level
67 */
68 void setMinLevel(std::size_t minLevel)
69 {
70 minLevel_ = minLevel;
71 }
72
73 /*!
74 * \brief Function to set the maximum allowed level.
75 *
76 *\param maxLevel The maximum level
77 */
78 void setMaxLevel(std::size_t maxLevel)
79 {
80 maxLevel_ = maxLevel;
81 }
82
83 /*!
84 * \brief Function to set the minimum/maximum allowed levels.
85 *
86 * \param minLevel The minimum level
87 * \param maxLevel The maximum level
88 */
89 void setLevels(std::size_t minLevel, std::size_t maxLevel)
90 {
91 minLevel_ = minLevel;
92 maxLevel_ = maxLevel;
93 }
94
95 /*!
96 * \brief Calculates the indicator used for refinement/coarsening for each grid cell.
97 *
98 * \param sol The solution vector
99 * \param refineTol The refinement tolerance
100 * \param coarsenTol The coarsening tolerance
101 *
102 * This standard two-phase indicator is based on the saturation gradient.
103 */
104 10 void calculate(const SolutionVector& sol,
105 Scalar refineTol = 0.05,
106 Scalar coarsenTol = 0.001)
107 {
108 //! Reset the indicator to a state that returns false for all elements
109 10 refineBound_ = std::numeric_limits<Scalar>::max();
110 10 coarsenBound_ = std::numeric_limits<Scalar>::lowest();
111 30 maxSaturationDelta_.assign(gridGeometry_->gridView().size(0), 0.0);
112
113 //! maxLevel_ must be higher than minLevel_ to allow for refinement
114
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 10 if (minLevel_ >= maxLevel_)
115 ✗ return;
116
117 //! Check for inadmissible tolerance combination
118
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 10 if (coarsenTol > refineTol)
119 ✗ DUNE_THROW(Dune::InvalidStateException, "Refine tolerance must be higher than coarsen tolerance");
120
121 //! Variables to hold the max/mon saturation values on the leaf
122 10 Scalar globalMax = std::numeric_limits<Scalar>::lowest();
123 10 Scalar globalMin = std::numeric_limits<Scalar>::max();
124
125 //! Calculate minimum and maximum saturation
126
10/12
✓ Branch 3 taken 6228 times.
✓ Branch 4 taken 10 times.
✓ Branch 5 taken 6228 times.
✓ Branch 6 taken 1 times.
✓ Branch 7 taken 9 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 53079 times.
✓ Branch 10 taken 9 times.
✓ Branch 11 taken 53079 times.
✓ Branch 12 taken 9 times.
✓ Branch 14 taken 53079 times.
✗ Branch 15 not taken.
 118644 for (const auto& element : elements(gridGeometry_->gridView()))
127 {
128 //! Index of the current leaf-element
129
3/6
✓ Branch 1 taken 53079 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 53079 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 53079 times.
✗ Branch 8 not taken.
 177921 const auto globalIdxI = gridGeometry_->elementMapper().index(element);
130
131 //! Obtain the saturation at the center of the element
132
2/4
✓ Branch 1 taken 53079 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 53079 times.
✗ Branch 5 not taken.
 112386 const auto geometry = element.geometry();
133
2/4
✓ Branch 1 taken 53079 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 53079 times.
✗ Branch 5 not taken.
 118614 const auto elemSol = elementSolution(element, sol, *gridGeometry_);
134
7/8
✓ Branch 1 taken 53079 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 32792 times.
✓ Branch 4 taken 26515 times.
✓ Branch 5 taken 31834 times.
✓ Branch 6 taken 21245 times.
✓ Branch 7 taken 31834 times.
✓ Branch 8 taken 21245 times.
 59307 const Scalar satI = evalSolution(element, geometry, *gridGeometry_, elemSol, geometry.center())[saturationIdx];
135
136 //! Maybe update the global minimum/maximum
137 using std::min;
138 using std::max;
139
2/2
✓ Branch 0 taken 32792 times.
✓ Branch 1 taken 26515 times.
 59307 globalMin = min(satI, globalMin);
140
2/2
✓ Branch 0 taken 47115 times.
✓ Branch 1 taken 12192 times.
 59307 globalMax = max(satI, globalMax);
141
142 //! Calculate maximum delta in saturation for this cell
143
9/17
✓ Branch 1 taken 53079 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 59307 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 59307 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 31122 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 267245 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 24558 times.
✓ Branch 16 taken 336 times.
✓ Branch 17 taken 214166 times.
✓ Branch 18 taken 24894 times.
✗ Branch 19 not taken.
 774319 for (const auto& intersection : intersections(gridGeometry_->gridView(), element))
144 {
145 //! Only consider internal intersections
146
4/4
✓ Branch 0 taken 24558 times.
✓ Branch 1 taken 212162 times.
✓ Branch 2 taken 26898 times.
✓ Branch 3 taken 336 times.
 263954 if (intersection.neighbor())
147 {
148 //! Access neighbor
149
1/2
✓ Branch 1 taken 236384 times.
✗ Branch 2 not taken.
 448210 const auto outside = intersection.outside();
150
3/6
✓ Branch 1 taken 236384 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 236384 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 236384 times.
✗ Branch 8 not taken.
 709152 const auto globalIdxJ = gridGeometry_->elementMapper().index(outside);
151
152 //! Visit intersection only once
153
18/18
✓ Branch 0 taken 24528 times.
✓ Branch 1 taken 211856 times.
✓ Branch 2 taken 24528 times.
✓ Branch 3 taken 30 times.
✓ Branch 4 taken 236354 times.
✓ Branch 5 taken 30 times.
✓ Branch 6 taken 232624 times.
✓ Branch 7 taken 3730 times.
✓ Branch 8 taken 24498 times.
✓ Branch 9 taken 208156 times.
✓ Branch 10 taken 24498 times.
✓ Branch 11 taken 30 times.
✓ Branch 12 taken 220375 times.
✓ Branch 13 taken 12249 times.
✓ Branch 14 taken 204426 times.
✓ Branch 15 taken 3700 times.
✓ Branch 16 taken 102213 times.
✓ Branch 17 taken 102213 times.
 285500 if (element.level() > outside.level() || (element.level() == outside.level() && globalIdxI < globalIdxJ))
154 {
155 //! Obtain saturation in the neighbor
156
1/2
✓ Branch 1 taken 118192 times.
✗ Branch 2 not taken.
 224105 const auto outsideGeometry = outside.geometry();
157
2/4
✓ Branch 1 taken 118192 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 118192 times.
✗ Branch 5 not taken.
 236384 const auto elemSolJ = elementSolution(outside, sol, *gridGeometry_);
158
9/10
✓ Branch 1 taken 118192 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 21216 times.
✓ Branch 4 taken 96976 times.
✓ Branch 5 taken 21216 times.
✓ Branch 6 taken 84697 times.
✓ Branch 7 taken 33495 times.
✓ Branch 8 taken 84697 times.
✓ Branch 9 taken 664 times.
✓ Branch 10 taken 11615 times.
 118192 const Scalar satJ = evalSolution(outside, outsideGeometry, *gridGeometry_, elemSolJ, outsideGeometry.center())[saturationIdx];
159
160 using std::abs;
161
2/2
✓ Branch 0 taken 21880 times.
✓ Branch 1 taken 96312 times.
 118192 Scalar localdelta = abs(satI - satJ);
162
4/4
✓ Branch 0 taken 21880 times.
✓ Branch 1 taken 96312 times.
✓ Branch 2 taken 21880 times.
✓ Branch 3 taken 96312 times.
 236384 maxSaturationDelta_[globalIdxI] = max(maxSaturationDelta_[globalIdxI], localdelta);
163
4/4
✓ Branch 0 taken 45793 times.
✓ Branch 1 taken 72399 times.
✓ Branch 2 taken 45793 times.
✓ Branch 3 taken 72399 times.
 282177 maxSaturationDelta_[globalIdxJ] = max(maxSaturationDelta_[globalIdxJ], localdelta);
164 }
165 }
166 }
167 }
168
169 //! Compute the maximum delta in saturation
170 10 const auto globalDelta = globalMax - globalMin;
171
172 //! Compute the refinement/coarsening bounds
173 10 refineBound_ = refineTol*globalDelta;
174 10 coarsenBound_ = coarsenTol*globalDelta;
175
176 // TODO: fix adaptive simulations in parallel
177 //#if HAVE_MPI
178 // // communicate updated values
179 // using DataHandle = VectorExchange<ElementMapper, ScalarSolutionType>;
180 // DataHandle dataHandle(problem_.elementMapper(), maxSaturationDelta_);
181 // problem_.gridView().template communicate<DataHandle>(dataHandle,
182 // Dune::InteriorBorder_All_Interface,
183 // Dune::ForwardCommunication);
184 //
185 // using std::max;
186 // refineBound_ = problem_.gridView().comm().max(refineBound_);
187 // coarsenBound_ = problem_.gridView().comm().max(coarsenBound_);
188 //
189 //#endif
190
191 //! check if neighbors have to be refined too
192
11/14
✓ Branch 3 taken 6228 times.
✓ Branch 4 taken 10 times.
✓ Branch 5 taken 6228 times.
✓ Branch 6 taken 1 times.
✓ Branch 7 taken 9 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 53079 times.
✓ Branch 10 taken 9 times.
✓ Branch 11 taken 53079 times.
✓ Branch 12 taken 9 times.
✓ Branch 14 taken 53079 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 53079 times.
✗ Branch 18 not taken.
 59346 for (const auto& element : elements(gridGeometry_->gridView(), Dune::Partitions::interior))
193
4/4
✓ Branch 1 taken 53877 times.
✓ Branch 2 taken 5430 times.
✓ Branch 3 taken 4039 times.
✓ Branch 4 taken 49040 times.
 59307 if (this->operator()(element) > 0)
194
1/2
✓ Branch 1 taken 4039 times.
✗ Branch 2 not taken.
 4837 checkNeighborsRefine_(element);
195 }
196
197 /*!
198 * \brief function call operator to return mark
199 *
200 * \return 1 if an element should be refined
201 * -1 if an element should be coarsened
202 * 0 otherwise
203 *
204 * \param element A grid element
205 */
206 118614 int operator() (const Element& element) const
207 {
208 if (element.hasFather()
209
8/9
✓ Branch 0 taken 12456 times.
✓ Branch 1 taken 94182 times.
✓ Branch 2 taken 24432 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 46564 times.
✓ Branch 7 taken 49414 times.
✓ Branch 8 taken 57224 times.
✓ Branch 9 taken 49414 times.
✓ Branch 10 taken 10660 times.
 131070 && maxSaturationDelta_[gridGeometry_->elementMapper().index(element)] < coarsenBound_)
210 {
211 return -1;
212 }
213
2/2
✓ Branch 0 taken 1596 times.
✓ Branch 1 taken 200 times.
 60336 else if (element.level() < maxLevel_
214
6/6
✓ Branch 0 taken 21610 times.
✓ Branch 1 taken 38726 times.
✓ Branch 5 taken 11931 times.
✓ Branch 6 taken 9679 times.
✓ Branch 7 taken 11931 times.
✓ Branch 8 taken 9679 times.
 60336 && maxSaturationDelta_[gridGeometry_->elementMapper().index(element)] > refineBound_)
215 {
216 return 1;
217 }
218 else
219 50657 return 0;
220 }
221
222 private:
223 /*!
224 * \brief Method ensuring the refinement ratio of 2:1
225 *
226 * For any given element, a loop over the neighbors checks if the
227 * entities refinement would require that any of the neighbors has
228 * to be refined, too. This is done recursively over all levels of the grid.
229 *
230 * \param element Element of interest that is to be refined
231 * \param level level of the refined element: it is at least 1
232 * \return true if everything was successful
233 */
234 4872 bool checkNeighborsRefine_(const Element &element, std::size_t level = 1)
235 {
236
8/13
✓ Branch 4 taken 798 times.
✓ Branch 5 taken 20584 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 798 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 3990 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 16510 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 3184 times.
✓ Branch 15 taken 8 times.
✓ Branch 17 taken 3192 times.
✗ Branch 18 not taken.
 58892 for(const auto& intersection : intersections(gridGeometry_->gridView(), element))
237 {
238
4/4
✓ Branch 0 taken 3184 times.
✓ Branch 1 taken 16469 times.
✓ Branch 2 taken 3233 times.
✓ Branch 3 taken 8 times.
 22894 if(!intersection.neighbor())
239 57 continue;
240
241 // obtain outside element
242
1/2
✓ Branch 1 taken 19645 times.
✗ Branch 2 not taken.
 36106 const auto outside = intersection.outside();
243
244 // only mark non-ghost elements
245
2/3
✓ Branch 0 taken 3184 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 16461 times.
 19645 if (outside.partitionType() == Dune::GhostEntity)
246 ✗ continue;
247
248
11/13
✓ Branch 0 taken 3159 times.
✓ Branch 1 taken 16486 times.
✓ Branch 2 taken 3159 times.
✓ Branch 3 taken 16232 times.
✓ Branch 4 taken 254 times.
✓ Branch 5 taken 3159 times.
✓ Branch 6 taken 16207 times.
✓ Branch 7 taken 3159 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 19366 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 41 times.
✓ Branch 12 taken 16166 times.
 22829 if (outside.level() < maxLevel_ && outside.level() < element.level())
249 {
250 // ensure refinement for outside element
251
5/8
✓ Branch 1 taken 41 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 41 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 41 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 35 times.
✓ Branch 10 taken 6 times.
 123 maxSaturationDelta_[gridGeometry_->elementMapper().index(outside)] = std::numeric_limits<Scalar>::max();
252
2/2
✓ Branch 0 taken 35 times.
✓ Branch 1 taken 6 times.
 41 if(level < maxLevel_)
253
1/2
✓ Branch 1 taken 35 times.
✗ Branch 2 not taken.
 35 checkNeighborsRefine_(outside, ++level);
254 }
255 }
256
257 4872 return true;
258 }
259
260 std::shared_ptr<const GridGeometry> gridGeometry_;
261
262 Scalar refineBound_;
263 Scalar coarsenBound_;
264 std::vector< Scalar > maxSaturationDelta_;
265 std::size_t minLevel_;
266 std::size_t maxLevel_;
267 };
268
269 } // end namespace Dumux
270
271 #endif
272