GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/adaptive/initializationindicator.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 57 62 91.9%
Functions: 15 15 100.0%
Branches: 127 203 62.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-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 Adaptive
10 * \brief Class defining an initialization indicator for grid adaption
11 */
12 #ifndef DUMUX_GRIDADAPTINITIALIZATIONINDICATOR_HH
13 #define DUMUX_GRIDADAPTINITIALIZATIONINDICATOR_HH
14
15 #include <memory>
16
17 #include <dune/geometry/type.hh>
18 #include <dumux/common/properties.hh>
19 #include <dumux/common/parameters.hh>
20 #include <dumux/common/typetraits/problem.hh>
21 #include <dumux/discretization/method.hh>
22
23 namespace Dumux {
24
25 /*!
26 * \ingroup Adaptive
27 * \brief Class defining an initialization indicator for grid adaption.
28 * Refines at sources and boundaries. Use before computing on the grid.
29 *
30 * \tparam TypeTag The problem TypeTag
31 */
32 template<class TypeTag>
33 class GridAdaptInitializationIndicator
34 {
35 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
36 using Problem = GetPropType<TypeTag, Properties::Problem>;
37 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
38 using Element = typename GridView::Traits::template Codim<0>::Entity;
39
40 using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
41 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
42
43 static constexpr bool isBox = GetPropType<TypeTag, Properties::GridGeometry>::discMethod == DiscretizationMethods::box;
44
45 public:
46
47 /*!
48 * \brief Constructor
49 * \note Reads the following parameters from the parameter tree
50 * - Adaptive.MinLevel The minimum refinement level
51 * - Adaptive.MaxLevel The maximum refinement level
52 * - Adaptive.RefineAtDirichletBC If to refine at Dirichlet boundaries (default: true)
53 * - Adaptive.RefineAtFluxBC If to refine at Neumann/Robin boundaries (default: true)
54 * - Adaptive.RefineAtSource If to refine where source terms are specified (default: true)
55 * - Adaptive.BCRefinementThreshold The threshold above which fluxes are treated as non-zero (default: 1e-10)
56 * \param problem The problem object
57 * \param gridGeometry The finite volume geometry of the grid
58 * \param gridVariables The secondary variables on the grid
59 */
60 5 GridAdaptInitializationIndicator(std::shared_ptr<const Problem> problem,
61 std::shared_ptr<const GridGeometry> gridGeometry,
62 std::shared_ptr<const GridVariables> gridVariables)
63 : problem_(problem)
64 , gridGeometry_(gridGeometry)
65 , gridVariables_(gridVariables)
66
2/4
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
 10 , minLevel_(getParamFromGroup<std::size_t>(problem->paramGroup(), "Adaptive.MinLevel"))
67
2/4
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
 10 , maxLevel_(getParamFromGroup<std::size_t>(problem->paramGroup(), "Adaptive.MaxLevel"))
68
2/4
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
 10 , refineAtDirichletBC_(getParamFromGroup<bool>(problem->paramGroup(), "Adaptive.RefineAtDirichletBC", true))
69
2/4
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
 10 , refineAtFluxBC_(getParamFromGroup<bool>(problem->paramGroup(), "Adaptive.RefineAtFluxBC", true))
70
2/4
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
 10 , refineAtSource_(getParamFromGroup<bool>(problem->paramGroup(), "Adaptive.RefineAtSource", true))
71
9/24
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 5 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 5 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 5 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 5 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 5 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 5 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 5 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 5 times.
✗ Branch 29 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.
 5 , eps_(getParamFromGroup<Scalar>(problem->paramGroup(), "Adaptive.BCRefinementThreshold", 1e-10))
72 5 {}
73
74 /*!
75 * \brief Function to set the minimum allowed level.
76 */
77 void setMinLevel(std::size_t minLevel)
78 {
79 minLevel_ = minLevel;
80 }
81
82 /*!
83 * \brief Function to set the maximum allowed level.
84 */
85 void setMaxLevel(std::size_t maxLevel)
86 {
87 maxLevel_ = maxLevel;
88 }
89
90 /*!
91 * \brief Function to set the minimum/maximum allowed levels.
92 */
93 void setLevels(std::size_t minLevel, std::size_t maxLevel)
94 {
95 minLevel_ = minLevel;
96 maxLevel_ = maxLevel;
97 }
98
99 /*!
100 * \brief Function to set if refinement at Dirichlet boundaries is to be used.
101 */
102 void setRefinementAtDirichletBC(bool refine)
103 {
104 refineAtDirichletBC_ = refine;
105 }
106
107 /*!
108 * \brief Function to set if refinement at sources is to be used.
109 */
110 void setRefinementAtSources(bool refine)
111 {
112 refineAtSource_ = refine;
113 }
114
115 /*!
116 * \brief Function to set if refinement at sources is to be used.
117 */
118 void setRefinementAtFluxBoundaries(bool refine)
119 {
120 refineAtFluxBC_ = refine;
121 }
122
123 /*!
124 * \brief Calculates the indicator used for refinement/coarsening for each grid cell.
125 *
126 * \param sol The solution for which indicator is to be calculated
127 */
128 template<class SolutionVector>
129 10 void calculate(const SolutionVector& sol)
130 {
131 //! prepare an indicator for refinement
132 30 indicatorVector_.assign(gridGeometry_->gridView().size(0), false);
133
134 // get the fvGeometry and elementVolVars needed for the bc and source interfaces
135
2/4
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10 times.
✗ Branch 5 not taken.
 28 auto fvGeometry = localView(*gridGeometry_);
136
3/8
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 10 times.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 40 auto elemVolVars = localView(gridVariables_->curGridVolVars());
137 // elemFluxVarsCache for neumann interface
138
4/8
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 10 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
 32 auto elemFluxVarsCache = localView(gridVariables_->gridFluxVarsCache());
139
140
15/20
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 10 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 12372 times.
✓ Branch 10 taken 10 times.
✓ Branch 11 taken 12372 times.
✓ Branch 12 taken 2 times.
✓ Branch 13 taken 8 times.
✓ Branch 14 taken 12372 times.
✓ Branch 15 taken 49350 times.
✓ Branch 16 taken 8 times.
✓ Branch 17 taken 61722 times.
✓ Branch 18 taken 8 times.
✓ Branch 20 taken 49350 times.
✗ Branch 21 not taken.
✓ Branch 23 taken 49350 times.
✗ Branch 24 not taken.
 61760 for (const auto& element : elements(gridGeometry_->gridView()))
141 {
142
3/6
✓ Branch 1 taken 61722 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 61722 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 61722 times.
✗ Branch 8 not taken.
 185166 const auto eIdx = gridGeometry_->elementMapper().index(element);
143
144 //! refine any element being below the minimum level
145
3/5
✗ Branch 0 not taken.
✓ Branch 1 taken 61722 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 12372 times.
✓ Branch 4 taken 49350 times.
 74094 if (element.level() < minLevel_)
146 {
147 ✗ indicatorVector_[eIdx] = true;
148 ✗ continue; // proceed to the next element
149 }
150
151 // If refinement at sources/BCs etc is deactivated, skip the rest
152
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 61722 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 61722 if (!refineAtSource_ && !refineAtFluxBC_ && !refineAtDirichletBC_)
153 continue;
154
155 // if the element is already on the maximum permissive level, skip rest
156
5/5
✓ Branch 0 taken 12202 times.
✓ Branch 1 taken 49520 times.
✓ Branch 2 taken 12202 times.
✓ Branch 3 taken 49256 times.
✓ Branch 4 taken 264 times.
 74094 if (element.level() == maxLevel_)
157 continue;
158
159 // Bind all of the local views
160
1/2
✓ Branch 1 taken 61288 times.
✗ Branch 2 not taken.
 61288 fvGeometry.bind(element);
161
1/2
✓ Branch 1 taken 61288 times.
✗ Branch 2 not taken.
 61288 elemVolVars.bind(element, fvGeometry, sol);
162
1/2
✓ Branch 1 taken 61288 times.
✗ Branch 2 not taken.
 61288 elemFluxVarsCache.bind(element, fvGeometry, elemVolVars);
163
164
165 //! Check if we have to refine around a source term
166
1/2
✓ Branch 0 taken 61288 times.
✗ Branch 1 not taken.
 61288 if (refineAtSource_)
167 {
168
4/4
✓ Branch 0 taken 97894 times.
✓ Branch 1 taken 61286 times.
✓ Branch 2 taken 48808 times.
✓ Branch 3 taken 12202 times.
 171382 for (const auto& scv : scvs(fvGeometry))
169 {
170 195788 auto source = problem_->source(element, fvGeometry, elemVolVars, scv);
171
2/4
✓ Branch 1 taken 97894 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 97894 times.
✗ Branch 5 not taken.
 195788 auto pointSource = problem_->scvPointSources(element, fvGeometry, elemVolVars, scv);
172
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 97892 times.
 293682 if (source.infinity_norm() + pointSource.infinity_norm() > eps_)
173 {
174 4 indicatorVector_[eIdx] = true;
175 2 break; // element is marked, escape scv loop
176 }
177 }
178 }
179
180 //! Check if we have to refine at the boundary
181
4/4
✓ Branch 0 taken 61286 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 61286 times.
✓ Branch 3 taken 2 times.
 122576 if (!indicatorVector_[eIdx] // proceed if element is not already marked
182
5/5
✓ Branch 0 taken 610 times.
✓ Branch 1 taken 60676 times.
✓ Branch 2 taken 610 times.
✓ Branch 3 taken 14056 times.
✓ Branch 4 taken 46620 times.
 73488 && element.hasBoundaryIntersections() // proceed if element is on boundary
183
4/6
✓ Branch 0 taken 61286 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 3074 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 3074 times.
✗ Branch 5 not taken.
 64362 && (refineAtDirichletBC_ || refineAtFluxBC_)) // proceed if boundary refinement is active
184 {
185 // cell-centered schemes
186 if (!isBox)
187 {
188
4/4
✓ Branch 0 taken 12314 times.
✓ Branch 1 taken 2400 times.
✓ Branch 2 taken 12314 times.
✓ Branch 3 taken 2400 times.
 17178 for (const auto& scvf : scvfs(fvGeometry))
189 {
190 // skip non-boundary scvfs
191
2/2
✓ Branch 0 taken 9210 times.
✓ Branch 1 taken 3104 times.
 12314 if (!scvf.boundary())
192 9210 continue;
193
194 6208 const auto bcTypes = problem_->boundaryTypes(element, scvf);
195 // We assume pure BCs, mixed boundary types are not allowed anyway!
196
5/6
✓ Branch 0 taken 1280 times.
✓ Branch 1 taken 1824 times.
✓ Branch 2 taken 1280 times.
✓ Branch 3 taken 1824 times.
✓ Branch 4 taken 1280 times.
✗ Branch 5 not taken.
 6208 if(bcTypes.hasOnlyDirichlet() && refineAtDirichletBC_)
197 {
198 ✗ indicatorVector_[eIdx] = true;
199 64 break; // element is marked, escape scvf loop
200 }
201
202 // we are on a pure Neumann boundary
203
1/2
✓ Branch 0 taken 3104 times.
✗ Branch 1 not taken.
 3104 else if(refineAtFluxBC_)
204 {
205 6208 const auto fluxes = problem_->neumann(element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf);
206
2/2
✓ Branch 0 taken 64 times.
✓ Branch 1 taken 3040 times.
 6208 if (fluxes.infinity_norm() > eps_)
207 {
208 128 indicatorVector_[eIdx] = true;
209 64 break; // element is marked, escape scvf loop
210 }
211 }
212 }
213 }
214 // box-scheme
215 else
216 {
217 // container to store bcTypes
218 using BoundaryTypes = typename ProblemTraits<Problem>::BoundaryTypes;
219
4/8
✓ Branch 1 taken 610 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 610 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 610 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 610 times.
✗ Branch 10 not taken.
 3050 std::vector<BoundaryTypes> bcTypes(fvGeometry.numScv());
220
221 // Get bcTypes and maybe mark for refinement on Dirichlet boundaries
222
4/4
✓ Branch 0 taken 2440 times.
✓ Branch 1 taken 610 times.
✓ Branch 2 taken 2440 times.
✓ Branch 3 taken 610 times.
 3660 for (const auto& scv : scvs(fvGeometry))
223 {
224
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 2440 times.
 4880 bcTypes[scv.localDofIndex()] = problem_->boundaryTypes(element, scv);
225
1/8
✗ Branch 0 not taken.
✓ Branch 1 taken 2440 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 2440 if (refineAtDirichletBC_ && bcTypes[scv.localDofIndex()].hasDirichlet())
226 {
227 ✗ indicatorVector_[eIdx] = true;
228 ✗ break; // element is marked, escape scv loop
229 }
230 }
231
232 // If element hasn't been marked because of Dirichlet BCS, check Neumann BCs
233
3/6
✓ Branch 0 taken 610 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 610 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 610 times.
✗ Branch 5 not taken.
 1220 if (!indicatorVector_[eIdx] && refineAtFluxBC_)
234 {
235
4/4
✓ Branch 0 taken 3660 times.
✓ Branch 1 taken 594 times.
✓ Branch 2 taken 3660 times.
✓ Branch 3 taken 594 times.
 4864 for (const auto& scvf : scvfs(fvGeometry))
236 {
237 //! check if scvf is on Neumann boundary
238
10/10
✓ Branch 0 taken 1220 times.
✓ Branch 1 taken 2440 times.
✓ Branch 2 taken 700 times.
✓ Branch 3 taken 520 times.
✓ Branch 4 taken 700 times.
✓ Branch 5 taken 520 times.
✓ Branch 6 taken 700 times.
✓ Branch 7 taken 520 times.
✓ Branch 8 taken 700 times.
✓ Branch 9 taken 520 times.
 3660 if (scvf.boundary() && bcTypes[scvf.insideScvIdx()].hasNeumann())
239 {
240 1400 const auto fluxes = problem_->neumann(element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf);
241
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 684 times.
 1400 if (fluxes.infinity_norm() > eps_)
242 {
243 32 indicatorVector_[eIdx] = true;
244 16 break; // element is marked, escape scvf loop
245 }
246 }
247 }
248 }
249 }
250 }
251 }
252 10 }
253
254 /*! \brief function call operator to return mark
255 *
256 * \return 1 if an element should be refined
257 * -1 if an element should be coarsened
258 * 0 otherwise
259 *
260 * \note In this initialization indicator implementation
261 * element coarsening is not considered.
262 *
263 * \param element A grid element
264 */
265 61722 int operator() (const Element& element) const
266 {
267
4/4
✓ Branch 3 taken 82 times.
✓ Branch 4 taken 61640 times.
✓ Branch 5 taken 82 times.
✓ Branch 6 taken 61640 times.
 185166 if (indicatorVector_[gridGeometry_->elementMapper().index(element)])
268 82 return 1;
269 return 0;
270 }
271
272 private:
273 std::shared_ptr<const Problem> problem_; //!< The problem to be solved
274 std::shared_ptr<const GridGeometry> gridGeometry_; //!< The finite volume grid geometry
275 std::shared_ptr<const GridVariables> gridVariables_; //!< The secondary variables on the grid
276 std::vector<bool> indicatorVector_; //!< Indicator for BCs/sources
277
278 int minLevel_; //!< The minimum allowed level
279 int maxLevel_; //!< The maximum allowed level
280 bool refineAtDirichletBC_; //!< Specifies if it should be refined at Dirichlet BCs
281 bool refineAtFluxBC_; //!< Specifies if it should be refined at non-zero Neumann BCs
282 bool refineAtSource_; //!< Specifies if it should be refined at sources
283 Scalar eps_; //!< Threshold for refinement at sources/BCS
284 };
285
286 } // end namespace Dumux
287
288 #endif
289