3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
gridadaptionindicatorlocal.hh
Go to the documentation of this file.
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 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
24#ifndef DUMUX_GRIDADAPTIONINDICATOR2PLOCAL_HH
25#define DUMUX_GRIDADAPTIONINDICATOR2PLOCAL_HH
26
29
30namespace Dumux {
31
38template<class TypeTag>
40{
41private:
45 using Element = typename GridView::Traits::template Codim<0>::Entity;
46
48 using ScalarSolutionType = typename SolutionTypes::ScalarSolution;
49
52
54
55 enum
56 {
57 sw = Indices::saturationW,
58 sn = Indices::saturationNw
59 };
60 enum
61 {
62 wPhaseIdx = Indices::wPhaseIdx,
63 nPhaseIdx = Indices::nPhaseIdx
64 };
65
66public:
73 {
74 // prepare an indicator for refinement
75 if(indicatorVector_.size() != problem_.variables().cellDataGlobal().size())
76 {
77 indicatorVector_.resize(problem_.variables().cellDataGlobal().size());
78 }
79 indicatorVector_ = -1e100;
80
81 Scalar globalMax = -1e100;
82 Scalar globalMin = 1e100;
83
84 Scalar maxLocalDelta = 0;
85
86 // 1) calculate Indicator -> min, maxvalues
87 // loop over all leaf-elements
88 for (const auto& element : elements(problem_.gridView()))
89 {
90 // calculate minimum and maximum saturation
91 // index of the current leaf-elements
92 int globalIdxI = problem_.variables().index(element);
93
94 if (refineAtSource_)
95 {
96 PrimaryVariables source(0.0);
97 problem_.source(source, element);
98 for (int i = 0; i < 2; i++)
99 {
100 using std::abs;
101 if (abs(source[i]) > 1e-10)
102 {
103 indicatorVector_[globalIdxI] = 10;
104 break;
105 }
106 }
107 }
108
109 if (indicatorVector_[globalIdxI] == 10)
110 {
111 continue;
112 }
113
114 Scalar satI = 0.0;
115 switch (saturationType_)
116 {
117 case sw:
118 satI = problem_.variables().cellData(globalIdxI).saturation(wPhaseIdx);
119 break;
120 case sn:
121 satI = problem_.variables().cellData(globalIdxI).saturation(nPhaseIdx);
122 break;
123 }
124
125 using std::min;
126 globalMin = min(satI, globalMin);
127 using std::max;
128 globalMax = max(satI, globalMax);
129
130 // calculate refinement indicator in all cells
131 for (const auto& intersection : intersections(problem_.gridView(), element))
132 {
133 if (indicatorVector_[globalIdxI] == 10)
134 {
135 break;
136 }
137
138 // exit, if it is not a neighbor
139 if (intersection.boundary())
140 {
141 BoundaryTypes bcTypes;
142 problem_.boundaryTypes(bcTypes, intersection);
143
144 for (int i = 0; i < 2; i++)
145 {
146 if (bcTypes.isNeumann(i))
147 {
148 PrimaryVariables flux(0.0);
149 problem_.neumann(flux, intersection);
150
151 bool fluxBound = false;
152 for (int j = 0; j < 2; j++)
153 {
154 using std::abs;
155 if (abs(flux[j]) > 1e-10)
156 {
157 if (refineAtFluxBC_)
158 {
159 indicatorVector_[globalIdxI] = 10;
160 fluxBound = true;
161 break;
162 }
163 }
164 }
165 if (fluxBound)
166 break;
167 }
168 else if (bcTypes.isDirichlet(i))
169 {
170 if (refineAtDirichletBC_)
171 {
172 indicatorVector_[globalIdxI] = 10;
173 break;
174 }
175 }
176 }
177 }
178 else
179 {
180 // get element pointer from neighbor
181 auto outside = intersection.outside();
182 int globalIdxJ = problem_.variables().index(outside);
183
184 // treat each intersection only from one side
185 if (element.level() > outside.level()
186 || (element.level() == outside.level() && globalIdxI < globalIdxJ))
187 {
188 Scalar satJ = 0.;
189 switch (saturationType_)
190 {
191 case sw:
192 satJ = problem_.variables().cellData(globalIdxJ).saturation(wPhaseIdx);
193 break;
194 case sn:
195 satJ = problem_.variables().cellData(globalIdxJ).saturation(nPhaseIdx);
196 break;
197 }
198
199 using std::abs;
200 Scalar localdelta = abs(satI - satJ);
201 using std::max;
202 indicatorVector_[globalIdxI][0] = max(indicatorVector_[globalIdxI][0], localdelta);
203 indicatorVector_[globalIdxJ][0] = max(indicatorVector_[globalIdxJ][0], localdelta);
204
205 maxLocalDelta = max(maxLocalDelta, localdelta);
206 }
207 }
208 }
209 }
210
211 Scalar globaldelta = globalMax - globalMin;
212
213 if (maxLocalDelta > 0.)
214 {
215 indicatorVector_ /= maxLocalDelta;
216 }
217 if (globaldelta > 0.)
218 {
219 maxLocalDelta /= globaldelta;
220 }
221
222 refineBound_ = refinetol_*maxLocalDelta;
223 coarsenBound_ = coarsentol_*maxLocalDelta;
224 }
225
233 bool refine(const Element& element)
234 {
235 return (indicatorVector_[problem_.elementMapper().index(element)] > refineBound_);
236 }
237
245 bool coarsen(const Element& element)
246 {
247 return (indicatorVector_[problem_.elementMapper().index(element)] < coarsenBound_);
248 }
249
253 void init()
254 {
255 refineBound_ = 0.;
256 coarsenBound_ = 0.;
257 };
258
265 {
266 indicatorVector_[idx] = coarsenBound_+1;
267 }
268
275 {
276 indicatorVector_[idx] = refineBound_ - 1;
277 }
278
290 problem_(problem)
291 {
292 refinetol_ = getParam<Scalar>("GridAdapt.RefineTolerance");
293 coarsentol_ = getParam<Scalar>("GridAdapt.CoarsenTolerance");
294 refineAtDirichletBC_ = getParam<bool>("GridAdapt.RefineAtDirichletBC");
295 refineAtFluxBC_ = getParam<bool>("GridAdapt.RefineAtFluxBC");
296 refineAtSource_ = getParam<bool>("GridAdapt.RefineAtSource");
297 }
298
299private:
300 Problem& problem_;
301 Scalar refinetol_;
302 Scalar coarsentol_;
303 Scalar refineBound_;
304 Scalar coarsenBound_;
305 ScalarSolutionType indicatorVector_;
306 static const int saturationType_ = getPropValue<TypeTag, Properties::SaturationFormulation>();
307 bool refineAtDirichletBC_; // switch for refinement at Dirichlet BC's
308 bool refineAtFluxBC_; // switch for refinement at Neumann BC's
309 bool refineAtSource_; // switch for refinement at sources
310};
311} // end namespace Dumux
312
313#endif
Base file for properties related to sequential IMPET algorithms.
Definition: adapt.hh:29
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type GetProp
get the type of a property (equivalent to old macro GET_PROP(...))
Definition: propertysystem.hh:140
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:149
Class defining a standard, saturation dependent indicator for grid adaption.
Definition: gridadaptionindicatorlocal.hh:40
void setIndicatorToCoarse(int idx)
Function for changing the indicatorVector values for coarsening.
Definition: gridadaptionindicatorlocal.hh:274
bool refine(const Element &element)
Indicator function for marking of grid cells for refinement.
Definition: gridadaptionindicatorlocal.hh:233
void calculateIndicator()
Calculates the indicator used for refinement/coarsening for each grid cell.
Definition: gridadaptionindicatorlocal.hh:72
void init()
Initializes the adaption indicator class.
Definition: gridadaptionindicatorlocal.hh:253
GridAdaptionIndicator2PLocal(Problem &problem)
Constructs a GridAdaptionIndicator instance.
Definition: gridadaptionindicatorlocal.hh:289
bool coarsen(const Element &element)
Indicator function for marking of grid cells for coarsening.
Definition: gridadaptionindicatorlocal.hh:245
void setIndicatorToRefine(int idx)
Function for changing the indicatorVector values for refinement.
Definition: gridadaptionindicatorlocal.hh:264
Defines the properties required for (immiscible) two-phase sequential models.