3.5-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
newtonconsistencychecks.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_NEWTON_CONSISTENCY_CHECKS
25#define DUMUX_PNM_NEWTON_CONSISTENCY_CHECKS
26
27#include <vector>
28#include <iostream>
29#include <dune/common/exceptions.hh>
33
34namespace Dumux::PoreNetwork {
35
36#ifndef DOXYGEN
37namespace Detail {
38
39// Primary template
40template <typename T, typename U = int>
41struct SaturationIndex {};
42
43// Specialization for 2p
44template <typename T>
45struct SaturationIndex <T, decltype((void) T::saturationIdx, 0)>
46{
47 static constexpr auto value = T::saturationIdx;
48};
49
50// Specialization for 2pnc
51template <typename T>
52struct SaturationIndex <T, decltype((void) T::switchIdx, 0)>
53{
54 static constexpr auto value = T::switchIdx;
55};
56
57// Specialization for MPNC
58template <typename T>
59struct SaturationIndex <T, decltype((void) T::s0Idx, 0)>
60{
61 static constexpr auto value = T::s0Idx;
62};
63} // end namespace Detail
64#endif
65
70template<class GridVariables, class SolutionVector>
72{
73public:
74
78 void performChecks(const GridVariables& gridVariables, const SolutionVector& uCurrentIter, const SolutionVector& prevSol) const
79 {
80 checkRelativeSaturationShift(gridVariables, uCurrentIter, prevSol);
81 checkIfValuesArePhysical(gridVariables, uCurrentIter);
82 checkIfCapillaryPressureIsCloseToEntryPressure(gridVariables, uCurrentIter);
83 }
84
88 void checkRelativeSaturationShift(const GridVariables& gridVariables, const SolutionVector& uCurrentIter, const SolutionVector& prevSol) const
89 {
90 using Scalar = typename SolutionVector::block_type::value_type;
91 const auto& problem = gridVariables.curGridVolVars().problem();
92
93 static const Scalar allowedSaturationChange = getParamFromGroup<Scalar>(problem.paramGroup(), "Newton.AllowedSaturationChange", -1.0);
94 if (allowedSaturationChange < 0.0)
95 return;
96
97 const auto& curGridVolVars = gridVariables.curGridVolVars();
98 const auto& prevGridVolVars = gridVariables.prevGridVolVars();
99
100 std::vector<bool> dofVisited(uCurrentIter.size(), false);
101
102 auto fvGeometry = localView(problem.gridGeometry());
103 auto curElemVolVars = localView(curGridVolVars);
104 auto prevElemVolVars = localView(prevGridVolVars);
105
106 for (const auto& element : elements(problem.gridGeometry().gridView()))
107 {
108 fvGeometry.bindElement(element);
109 curElemVolVars.bindElement(element, fvGeometry, uCurrentIter);
110 prevElemVolVars.bindElement(element, fvGeometry, prevSol);
111
112 for (const auto& scv : scvs(fvGeometry))
113 {
114 if (dofVisited[scv.dofIndex()])
115 continue;
116 dofVisited[scv.dofIndex()] = true;
117
118 const Scalar satNew = curElemVolVars[scv].saturation(0);
119 const Scalar satOld = prevElemVolVars[scv].saturation(0);
120 using std::abs;
121 const Scalar deltaS = abs(satNew - satOld);
122 static const bool considerRelativeShift = getParamFromGroup<bool>(problem.paramGroup(), "Newton.SaturationChangeIsRelative", false);
123
124 if (considerRelativeShift)
125 {
126 // satOld mgiht be (close to) zero, so have to take care of this
127 if (satOld > 1e-3 && deltaS / satOld > allowedSaturationChange)
128 DUNE_THROW(NumericalProblem, "Saturation change too high at dof " << scv.dofIndex() << ", old sat. " << satOld << ", new sat. " << satNew << std::endl);
129 }
130 else if (deltaS > allowedSaturationChange)
131 DUNE_THROW(NumericalProblem, "Saturation change too high at dof " << scv.dofIndex() << ", old sat. " << satOld << ", new sat. " << satNew << std::endl);
132 }
133 }
134 }
135
139 void checkIfValuesArePhysical(const GridVariables& gridVariables, const SolutionVector& uCurrentIter) const
140 {
141 const auto& problem = gridVariables.curGridVolVars().problem();
142
143 static const bool doPlausibilityCheck = getParamFromGroup<bool>(problem.paramGroup(), "Newton.PlausibilityCheck", false);
144 if (!doPlausibilityCheck)
145 return;
146
147 for (std::size_t i = 0; i < uCurrentIter.size(); ++i)
148 {
149 const auto& priVars = uCurrentIter[i];
150 using Indices = typename GridVariables::VolumeVariables::Indices;
151 static constexpr auto saturationOrMoleFractionIndex = Detail::SaturationIndex<Indices>::value;
152 if (priVars[saturationOrMoleFractionIndex] < 0.0 || priVars[saturationOrMoleFractionIndex] > 1.0)
153 {
154 std::cout << "at dof " << i << ": saturation " << priVars[saturationOrMoleFractionIndex] << std::endl;
155 DUNE_THROW(NumericalProblem, "Saturation is below 0 or above 1");
156 }
157 }
158 }
159
163 void checkIfCapillaryPressureIsCloseToEntryPressure(const GridVariables& gridVariables, const SolutionVector& uCurrentIter) const
164 {
165 // this check is implemented in the invasion state itself
166 const auto& invasionState = gridVariables.gridFluxVarsCache().invasionState();
167 invasionState.checkIfCapillaryPressureIsCloseToEntryPressure(uCurrentIter, gridVariables.curGridVolVars(), gridVariables.gridFluxVarsCache());
168 }
169};
170} // end namespace Dumux
171#endif
Some exceptions thrown in DuMux
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.
Free function to get the local view of a grid cache object.
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
Exception thrown if a fixable numerical problem occurs.
Definition: exceptions.hh:39
Consistency checks for the PNM two-phase model.
Definition: newtonconsistencychecks.hh:72
void checkRelativeSaturationShift(const GridVariables &gridVariables, const SolutionVector &uCurrentIter, const SolutionVector &prevSol) const
Checks if the relative shift of saturation between to consecutive time steps is below a given treshol...
Definition: newtonconsistencychecks.hh:88
void checkIfValuesArePhysical(const GridVariables &gridVariables, const SolutionVector &uCurrentIter) const
Checks if the saturation is between zero and one.
Definition: newtonconsistencychecks.hh:139
void performChecks(const GridVariables &gridVariables, const SolutionVector &uCurrentIter, const SolutionVector &prevSol) const
Perform all checks.
Definition: newtonconsistencychecks.hh:78
void checkIfCapillaryPressureIsCloseToEntryPressure(const GridVariables &gridVariables, const SolutionVector &uCurrentIter) const
Checks if the capillary pressure at pore from which a throat was invaded is sufficiently close to the...
Definition: newtonconsistencychecks.hh:163