3.5-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
gridadaptindicator.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 *****************************************************************************/
25#ifndef DUMUX_TWOP_ADAPTION_INDICATOR_HH
26#define DUMUX_TWOP_ADAPTION_INDICATOR_HH
27
28#include <memory>
29#include <dune/common/exceptions.hh>
30#include <dune/grid/common/partitionset.hh>
31
36
37namespace Dumux {
38
43template<class TypeTag>
45{
48 using Element = typename GridView::template Codim<0>::Entity;
52
53 enum { saturationIdx = Indices::saturationIdx };
54
55public:
66 TwoPGridAdaptIndicator(std::shared_ptr<const GridGeometry> gridGeometry, const std::string& paramGroup = "")
67 : gridGeometry_(gridGeometry)
68 , refineBound_(std::numeric_limits<Scalar>::max())
69 , coarsenBound_(std::numeric_limits<Scalar>::lowest())
70 , maxSaturationDelta_(gridGeometry_->gridView().size(0), 0.0)
71 , minLevel_(getParamFromGroup<std::size_t>(paramGroup, "Adaptive.MinLevel", 0))
72 , maxLevel_(getParamFromGroup<std::size_t>(paramGroup, "Adaptive.MaxLevel", 0))
73 {}
74
80 void setMinLevel(std::size_t minLevel)
81 {
82 minLevel_ = minLevel;
83 }
84
90 void setMaxLevel(std::size_t maxLevel)
91 {
92 maxLevel_ = maxLevel;
93 }
94
101 void setLevels(std::size_t minLevel, std::size_t maxLevel)
102 {
103 minLevel_ = minLevel;
104 maxLevel_ = maxLevel;
105 }
106
116 void calculate(const SolutionVector& sol,
117 Scalar refineTol = 0.05,
118 Scalar coarsenTol = 0.001)
119 {
121 refineBound_ = std::numeric_limits<Scalar>::max();
122 coarsenBound_ = std::numeric_limits<Scalar>::lowest();
123 maxSaturationDelta_.assign(gridGeometry_->gridView().size(0), 0.0);
124
126 if (minLevel_ >= maxLevel_)
127 return;
128
130 if (coarsenTol > refineTol)
131 DUNE_THROW(Dune::InvalidStateException, "Refine tolerance must be higher than coarsen tolerance");
132
134 Scalar globalMax = std::numeric_limits<Scalar>::lowest();
135 Scalar globalMin = std::numeric_limits<Scalar>::max();
136
138 for (const auto& element : elements(gridGeometry_->gridView()))
139 {
141 const auto globalIdxI = gridGeometry_->elementMapper().index(element);
142
144 const auto geometry = element.geometry();
145 const auto elemSol = elementSolution(element, sol, *gridGeometry_);
146 const Scalar satI = evalSolution(element, geometry, *gridGeometry_, elemSol, geometry.center())[saturationIdx];
147
149 using std::min;
150 using std::max;
151 globalMin = min(satI, globalMin);
152 globalMax = max(satI, globalMax);
153
155 for (const auto& intersection : intersections(gridGeometry_->gridView(), element))
156 {
158 if (intersection.neighbor())
159 {
161 const auto outside = intersection.outside();
162 const auto globalIdxJ = gridGeometry_->elementMapper().index(outside);
163
165 if (element.level() > outside.level() || (element.level() == outside.level() && globalIdxI < globalIdxJ))
166 {
168 const auto outsideGeometry = outside.geometry();
169 const auto elemSolJ = elementSolution(outside, sol, *gridGeometry_);
170 const Scalar satJ = evalSolution(outside, outsideGeometry, *gridGeometry_, elemSolJ, outsideGeometry.center())[saturationIdx];
171
172 using std::abs;
173 Scalar localdelta = abs(satI - satJ);
174 maxSaturationDelta_[globalIdxI] = max(maxSaturationDelta_[globalIdxI], localdelta);
175 maxSaturationDelta_[globalIdxJ] = max(maxSaturationDelta_[globalIdxJ], localdelta);
176 }
177 }
178 }
179 }
180
182 const auto globalDelta = globalMax - globalMin;
183
185 refineBound_ = refineTol*globalDelta;
186 coarsenBound_ = coarsenTol*globalDelta;
187
188// TODO: fix adaptive simulations in parallel
189//#if HAVE_MPI
190// // communicate updated values
191// using DataHandle = VectorExchange<ElementMapper, ScalarSolutionType>;
192// DataHandle dataHandle(problem_.elementMapper(), maxSaturationDelta_);
193// problem_.gridView().template communicate<DataHandle>(dataHandle,
194// Dune::InteriorBorder_All_Interface,
195// Dune::ForwardCommunication);
196//
197// using std::max;
198// refineBound_ = problem_.gridView().comm().max(refineBound_);
199// coarsenBound_ = problem_.gridView().comm().max(coarsenBound_);
200//
201//#endif
202
204 for (const auto& element : elements(gridGeometry_->gridView(), Dune::Partitions::interior))
205 if (this->operator()(element) > 0)
206 checkNeighborsRefine_(element);
207 }
208
218 int operator() (const Element& element) const
219 {
220 if (element.hasFather()
221 && maxSaturationDelta_[gridGeometry_->elementMapper().index(element)] < coarsenBound_)
222 {
223 return -1;
224 }
225 else if (element.level() < maxLevel_
226 && maxSaturationDelta_[gridGeometry_->elementMapper().index(element)] > refineBound_)
227 {
228 return 1;
229 }
230 else
231 return 0;
232 }
233
234private:
246 bool checkNeighborsRefine_(const Element &element, std::size_t level = 1)
247 {
248 for(const auto& intersection : intersections(gridGeometry_->gridView(), element))
249 {
250 if(!intersection.neighbor())
251 continue;
252
253 // obtain outside element
254 const auto outside = intersection.outside();
255
256 // only mark non-ghost elements
257 if (outside.partitionType() == Dune::GhostEntity)
258 continue;
259
260 if (outside.level() < maxLevel_ && outside.level() < element.level())
261 {
262 // ensure refinement for outside element
263 maxSaturationDelta_[gridGeometry_->elementMapper().index(outside)] = std::numeric_limits<Scalar>::max();
264 if(level < maxLevel_)
265 checkNeighborsRefine_(outside, ++level);
266 }
267 }
268
269 return true;
270 }
271
272 std::shared_ptr<const GridGeometry> gridGeometry_;
273
274 Scalar refineBound_;
275 Scalar coarsenBound_;
276 std::vector< Scalar > maxSaturationDelta_;
277 std::size_t minLevel_;
278 std::size_t maxLevel_;
279};
280
281} // end namespace Dumux
282
283#endif
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.
Element solution classes and factory functions.
free functions for the evaluation of primary variables inside elements.
PrimaryVariables evalSolution(const Element &element, const typename Element::Geometry &geometry, const typename FVElementGeometry::GridGeometry &gridGeometry, const BoxElementSolution< FVElementGeometry, PrimaryVariables > &elemSol, const typename Element::Geometry::GlobalCoordinate &globalPos, bool ignoreState=false)
Interpolates a given box element solution at a given global position. Uses the finite element cache o...
Definition: evalsolution.hh:95
auto elementSolution(const Element &element, const SolutionVector &sol, const GridGeometry &gg) -> std::enable_if_t< GridGeometry::discMethod==DiscretizationMethods::box, BoxElementSolution< typename GridGeometry::LocalView, std::decay_t< decltype(std::declval< SolutionVector >()[0])> > >
Make an element solution for box schemes.
Definition: box/elementsolution.hh:118
T getParamFromGroup(Args &&... args)
A free function to get a parameter from the parameter tree singleton with a model group.
Definition: parameters.hh:161
Definition: adapt.hh:29
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property
Definition: propertysystem.hh:150
Class defining a standard, saturation dependent indicator for grid adaptation.
Definition: gridadaptindicator.hh:45
void setMinLevel(std::size_t minLevel)
Function to set the minimum allowed level.
Definition: gridadaptindicator.hh:80
int operator()(const Element &element) const
function call operator to return mark
Definition: gridadaptindicator.hh:218
void setMaxLevel(std::size_t maxLevel)
Function to set the maximum allowed level.
Definition: gridadaptindicator.hh:90
TwoPGridAdaptIndicator(std::shared_ptr< const GridGeometry > gridGeometry, const std::string &paramGroup="")
The Constructor.
Definition: gridadaptindicator.hh:66
void calculate(const SolutionVector &sol, Scalar refineTol=0.05, Scalar coarsenTol=0.001)
Calculates the indicator used for refinement/coarsening for each grid cell.
Definition: gridadaptindicator.hh:116
void setLevels(std::size_t minLevel, std::size_t maxLevel)
Function to set the minumum/maximum allowed levels.
Definition: gridadaptindicator.hh:101
Declares all properties used in Dumux.