3.5-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
invasionstate.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_PNM_2P_INVASIONSTATE_HH
25#define DUMUX_PNM_2P_INVASIONSTATE_HH
26
27#include <vector>
28#include <type_traits>
29#include <dune/common/std/type_traits.hh>
32
33namespace Dumux::PoreNetwork {
34
39template<class P>
41{
42 using Problem = P;
43
44 template <class T>
45 using GlobalCapillaryPressureDetector = decltype(std::declval<T>().globalCapillaryPressure());
46
47 template<class T>
48 static constexpr bool hasGlobalCapillaryPressure()
49 { return Dune::Std::is_detected<GlobalCapillaryPressureDetector, T>::value; }
50
51 enum class EventType {invasion, snapOff, none};
52
53public:
54
55 TwoPInvasionState(const Problem& problem) : problem_(problem)
56 {
57 // initialize the invasion state
58 invadedCurrentIteration_.resize(problem.gridGeometry().gridView().size(0));
59 invadedPreviousTimeStep_.resize(problem.gridGeometry().gridView().size(0));
60
61 for (auto&& element : elements(problem.gridGeometry().gridView()))
62 {
63 const auto eIdx = problem.gridGeometry().elementMapper().index(element);
64 invadedCurrentIteration_[eIdx] = problem.initialInvasionState(element);
65 invadedPreviousTimeStep_[eIdx] = invadedCurrentIteration_[eIdx];
66 }
67
68 numThroatsInvaded_ = std::count(invadedCurrentIteration_.begin(), invadedCurrentIteration_.end(), true);
69 verbose_ = getParamFromGroup<bool>(problem.paramGroup(), "InvasionState.Verbosity", true);
70 restrictToGlobalCapillaryPressure_ = getParamFromGroup<bool>(problem.paramGroup(), "InvasionState.RestrictInvasionToGlobalCapillaryPressure", false);
71
72 if constexpr (hasGlobalCapillaryPressure<Problem>())
73 {
74 if (restrictToGlobalCapillaryPressure_)
75 std::cout << "\n *** Invasion behavior is restricted by a global capillary pressure defined in the problem! *** \n" << std::endl;
76 else
77 std::cout << "\n *** WARNING: global capillary pressure defined in the problem but InvasionState.RestrictInvasionToGlobalCapillaryPressure is set to false.\n"
78 << " Invasion behavior will NOT be restricted! ***\n" << std::endl;
79 }
80 }
81
83 template<class Element>
84 bool invaded(const Element& element) const
85 {
86 const auto eIdx = problem_.gridGeometry().elementMapper().index(element);
87 return invadedCurrentIteration_[eIdx];
88 }
89
91 std::size_t numThroatsInvaded() const
92 { return numThroatsInvaded_; }
93
95 template<class SolutionVector, class GridVolumeVariables, class GridFluxVariablesCache>
96 bool update(const SolutionVector& sol, const GridVolumeVariables& gridVolVars, GridFluxVariablesCache& gridFluxVarsCache)
97 {
98 hasChangedInCurrentIteration_ = false;
99 auto fvGeometry = localView(problem_.gridGeometry());
100 auto elemVolVars = localView(gridVolVars);
101 auto elemFluxVarsCache = localView(gridFluxVarsCache);
102 for (auto&& element : elements(problem_.gridGeometry().gridView()))
103 {
104 fvGeometry.bindElement(element);
105 elemVolVars.bind(element, fvGeometry, sol);
106 elemFluxVarsCache.bind(element, fvGeometry, elemVolVars);
107
108 for (auto&& scvf : scvfs(fvGeometry))
109 {
110 // checks if invasion or snap-off occured after Newton iteration step
111 if (const auto invasionResult = invasionSwitch_(element, elemVolVars, elemFluxVarsCache[scvf]); invasionResult)
112 {
113 hasChangedInCurrentIteration_ = true;
114 if constexpr (GridFluxVariablesCache::cachingEnabled)
115 {
116 const auto eIdx = problem_.gridGeometry().elementMapper().index(element);
117 gridFluxVarsCache.cache(eIdx, scvf.index()).update(problem_, element, fvGeometry, elemVolVars, scvf, invadedCurrentIteration_[eIdx]);
118 }
119 }
120 }
121 }
122 numThroatsInvaded_ = std::count(invadedCurrentIteration_.begin(), invadedCurrentIteration_.end(), true);
123 return hasChangedInCurrentIteration_;
124 }
125
127 void reset()
128 {
129 hasChangedInCurrentIteration_ = false;
130 invadedCurrentIteration_ = invadedPreviousTimeStep_;
131 }
132
134 bool hasChanged() const
135 { return hasChangedComparedToPreviousTimestep_; }
136
139 { return hasChangedInCurrentIteration_; }
140
142 void advance()
143 {
144 hasChangedComparedToPreviousTimestep_ = (invadedPreviousTimeStep_ != invadedCurrentIteration_);
145 invadedPreviousTimeStep_ = invadedCurrentIteration_;
146 }
147
148 template<class SolutionVector, class GridVolumeVariables, class GridFluxVariablesCache>
150 const GridVolumeVariables& gridVolVars,
151 const GridFluxVariablesCache& gridFluxVarsCache) const
152 {
153 using Scalar = typename SolutionVector::block_type::value_type;
154 static const Scalar accuracyCriterion = getParamFromGroup<Scalar>(problem_.paramGroup(), "InvasionState.AccuracyCriterion", -1.0);
155
156 if (accuracyCriterion < 0.0)
157 return;
158
159 auto fvGeometry = localView(problem_.gridGeometry());
160 auto elemVolVars = localView(gridVolVars);
161 auto elemFluxVarsCache = localView(gridFluxVarsCache);
162 for (auto&& element : elements(problem_.gridGeometry().gridView()))
163 {
164 // Only consider throats which have been invaded during the current time step
165 const auto eIdx = problem_.gridGeometry().elementMapper().index(element);
166 if (!invadedCurrentIteration_[eIdx] || invadedPreviousTimeStep_[eIdx] == invadedCurrentIteration_[eIdx])
167 continue;
168
169 fvGeometry.bindElement(element);
170 elemVolVars.bind(element, fvGeometry, sol);
171 elemFluxVarsCache.bind(element, fvGeometry, elemVolVars);
172
173 for (auto&& scvf : scvfs(fvGeometry))
174 {
175 // checks if pc is close enough to the entry pressure value
176 const auto& fluxVarsCache = elemFluxVarsCache[scvf];
177
178 using std::max;
179 const Scalar pc = max(elemVolVars[0].capillaryPressure(), elemVolVars[1].capillaryPressure());
180
181 if (pc < accuracyCriterion * fluxVarsCache.pcEntry())
182 DUNE_THROW(NumericalProblem, "At element " << eIdx << ": pc " << pc << " too far away form pcEntry " << fluxVarsCache.pcEntry());
183 }
184 }
185 }
186
187private:
188
190 template<class Element, class ElementVolumeVariables, class FluxVariablesCache>
191 auto invasionSwitch_(const Element& element,
192 const ElementVolumeVariables& elemVolVars,
193 const FluxVariablesCache& fluxVarsCache)
194
195 {
196 using Scalar = typename ElementVolumeVariables::VolumeVariables::PrimaryVariables::value_type;
197 const auto& gridGeometry = problem_.gridGeometry();
198 const auto& spatialParams = problem_.spatialParams();
199 const auto eIdx = gridGeometry.elementMapper().index(element);
200 bool invadedBeforeSwitch = invadedCurrentIteration_[eIdx];
201 bool invadedAfterSwitch = invadedBeforeSwitch;
202
203 // Result type, containing the local scv index of the pore from which the invasion/snap-off occurred
204 // Evaluates to 'false' if no invasion/snap-off occurred
205 struct Result
206 {
207 std::uint8_t localScvIdxWithCriticalPc;
208 Scalar criticalPc;
209 EventType event = EventType::none;
210
211 operator bool() const
212 { return event != EventType::none; }
213 };
214
215 // Block non-wetting phase flux out of the outlet
216 static const auto blockNonwettingPhase = getParamFromGroup<std::vector<int>>(problem_.paramGroup(), "InvasionState.BlockNonwettingPhaseAtThroatLabel", std::vector<int>{Labels::outlet});
217 if (!blockNonwettingPhase.empty() && std::find(blockNonwettingPhase.begin(), blockNonwettingPhase.end(), gridGeometry.throatLabel(eIdx)) != blockNonwettingPhase.end())
218 {
219 invadedCurrentIteration_[eIdx] = false;
220 return Result{}; // nothing happened
221 }
222
223 //Determine whether throat gets invaded or snap-off occurs
224 const std::array<Scalar, 2> pc = { elemVolVars[0].capillaryPressure(), elemVolVars[1].capillaryPressure() };
225 const auto pcMax = std::max_element(pc.begin(), pc.end());
226 const Scalar pcEntry = fluxVarsCache.pcEntry();
227 const Scalar pcSnapoff = fluxVarsCache.pcSnapoff();
228
229 // check if there is a user-specified global capillary pressure which needs to be obeyed
230 if (maybeRestrictToGlobalCapillaryPressure_(pcEntry))
231 {
232 if (*pcMax > pcEntry)
233 {
234 std::cout << "Throat " << eIdx << " would have been invaded by pc of " << *pcMax << "but a global capillary pressure restricion was set in the problem.";
235 std::cout << ". pcEntry: " << spatialParams.pcEntry(element, elemVolVars) << std::endl;
236 }
237
238 invadedCurrentIteration_[eIdx] = false;
239 return Result{}; //nothing happened
240 }
241
242 if (*pcMax > pcEntry)
243 invadedAfterSwitch = true;
244 else if (*pcMax <= pcSnapoff)
245 invadedAfterSwitch = false;
246
247 invadedCurrentIteration_[eIdx] = invadedAfterSwitch;
248
249 if (invadedBeforeSwitch == invadedAfterSwitch)
250 return Result{}; // nothing happened
251 else
252 {
253 Result result;
254 result.localScvIdxWithCriticalPc = std::distance(pc.begin(), pcMax);
255 result.criticalPc = *pcMax;
256 result.event = !invadedBeforeSwitch && invadedAfterSwitch ? EventType::invasion : EventType::snapOff;
257
258 if (verbose_)
259 {
260 const auto wPhaseIdx = spatialParams.template wettingPhase<typename ElementVolumeVariables::VolumeVariables::FluidSystem>(element, elemVolVars);
261 const std::array sw = { elemVolVars[0].saturation(wPhaseIdx), elemVolVars[1].saturation(wPhaseIdx) };
262 const auto vIdx = gridGeometry.gridView().indexSet().subIndex(element, result.localScvIdxWithCriticalPc, 1);
263 if (result.event == EventType::invasion)
264 {
265 std::cout << "Throat " << eIdx << " was invaded from pore " << vIdx << " :";
266 std::cout << " pc: " << *pcMax;
267 std::cout << ", pcEntry: " << spatialParams.pcEntry(element, elemVolVars);
268 std::cout << ", sw: " << sw[result.localScvIdxWithCriticalPc] << std::endl;
269 }
270 else
271 {
272 std::cout << "Snap-off occured at throat " << eIdx << " from pore " << vIdx << " :";
273 std::cout << " pc: " << *pcMax;
274 std::cout << ", pcSnapoff: " << spatialParams.pcSnapoff(element, elemVolVars);
275 std::cout << ", sw: " << sw[result.localScvIdxWithCriticalPc] << std::endl;
276 }
277 }
278
279 return result;
280 }
281 }
282
285 template<class Scalar>
286 bool maybeRestrictToGlobalCapillaryPressure_(const Scalar pcEntry) const
287 {
288 if constexpr (hasGlobalCapillaryPressure<Problem>())
289 return restrictToGlobalCapillaryPressure_ && (pcEntry > problem_.globalCapillaryPressure());
290 else
291 return false;
292 }
293
294 std::vector<bool> invadedCurrentIteration_;
295 std::vector<bool> invadedPreviousTimeStep_;
296 bool hasChangedInCurrentIteration_ = false;
297 bool hasChangedComparedToPreviousTimestep_ = false;
298 std::size_t numThroatsInvaded_;
299 bool verbose_;
300 bool restrictToGlobalCapillaryPressure_;
301
302 const Problem& problem_;
303};
304
305} // end namespace Dumux::PoreNetwork
306
307#endif
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.
Defines labels for pores and throats.
static ctype distance(const Dune::FieldVector< ctype, dimWorld > &a, const Dune::FieldVector< ctype, dimWorld > &b)
Compute the shortest distance between two points.
Definition: distance.hh:292
GridCache::LocalView localView(const GridCache &gridCache)
Free function to get the local view of a grid cache object.
Definition: localview.hh:38
Definition: discretization/porenetwork/fvelementgeometry.hh:34
std::string capillaryPressure() noexcept
I/O name of capillary pressure.
Definition: name.hh:135
Exception thrown if a fixable numerical problem occurs.
Definition: exceptions.hh:39
This class updates the invasion state for the two-phase PNM.
Definition: invasionstate.hh:41
void checkIfCapillaryPressureIsCloseToEntryPressure(const SolutionVector &sol, const GridVolumeVariables &gridVolVars, const GridFluxVariablesCache &gridFluxVarsCache) const
Definition: invasionstate.hh:149
bool hasChangedInCurrentIteration() const
Return whether an invasion or snap-off occurred anywhere during the current Newton iteration.
Definition: invasionstate.hh:138
bool invaded(const Element &element) const
Return whether a given throat is invaded or not.
Definition: invasionstate.hh:84
bool hasChanged() const
Return whether an invasion or snap-off occurred anywhere. Can be used, e.g., for output file writing ...
Definition: invasionstate.hh:134
bool update(const SolutionVector &sol, const GridVolumeVariables &gridVolVars, GridFluxVariablesCache &gridFluxVarsCache)
Update the invasion state of all throats. This is done after each Newton step by a call from the Newt...
Definition: invasionstate.hh:96
TwoPInvasionState(const Problem &problem)
Definition: invasionstate.hh:55
void advance()
This is called after the Newton method has successfully finished one time step.
Definition: invasionstate.hh:142
std::size_t numThroatsInvaded() const
Return the number of currently invaded throats.
Definition: invasionstate.hh:91
void reset()
Restore the old invasion state after a Newton iteration has failed.
Definition: invasionstate.hh:127
static constexpr int outlet
Label for pores/throats on an outlet.
Definition: labels.hh:41