ncpnewtonmethod.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 This file is part of the Open Porous Media project (OPM).
5
6 OPM 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 2 of the License, or
9 (at your option) any later version.
10
11 OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
18
19 Consult the COPYING file in the top-level source directory of this
20 module for the precise wording of the license and the list of
21 copyright holders.
22*/
28#ifndef EWOMS_NCP_NEWTON_METHOD_HH
29#define EWOMS_NCP_NEWTON_METHOD_HH
30
31#include "ncpproperties.hh"
32
33#include <opm/common/Exceptions.hpp>
34
36
37#include <algorithm>
38
39namespace Opm::Properties {
40
41template <class TypeTag, class MyTypeTag>
42struct DiscNewtonMethod;
43
44} // namespace Opm::Properties
45
46namespace Opm {
47
53template <class TypeTag>
54class NcpNewtonMethod : public GetPropType<TypeTag, Properties::DiscNewtonMethod>
55{
57
65
66 enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
67 enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
68 enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
69 enum { fugacity0Idx = Indices::fugacity0Idx };
70 enum { saturation0Idx = Indices::saturation0Idx };
71 enum { pressure0Idx = Indices::pressure0Idx };
72 enum { ncp0EqIdx = Indices::ncp0EqIdx };
73
74public:
78 NcpNewtonMethod(Simulator& simulator) : ParentType(simulator)
79 {}
80
81protected:
82 friend ParentType;
84
85 void preSolve_(const SolutionVector&,
86 const GlobalEqVector& currentResidual)
87 {
88 const auto& constraintsMap = this->model().linearizer().constraintsMap();
89 this->lastError_ = this->error_;
90
91 // calculate the error as the maximum weighted tolerance of
92 // the solution's residual
93 this->error_ = 0;
94 for (unsigned dofIdx = 0; dofIdx < currentResidual.size(); ++dofIdx) {
95 // do not consider auxiliary DOFs for the error
96 if (dofIdx >= this->model().numGridDof() || this->model().dofTotalVolume(dofIdx) <= 0.0)
97 continue;
98
99 // also do not consider DOFs which are constraint
100 if (this->enableConstraints_()) {
101 if (constraintsMap.count(dofIdx) > 0)
102 continue;
103 }
104
105 const auto& r = currentResidual[dofIdx];
106 for (unsigned eqIdx = 0; eqIdx < r.size(); ++eqIdx) {
107 if (ncp0EqIdx <= eqIdx && eqIdx < Indices::ncp0EqIdx + numPhases)
108 continue;
109 this->error_ =
110 std::max(std::abs(r[eqIdx]*this->model().eqWeight(dofIdx, eqIdx)),
111 this->error_);
112 }
113 }
114
115 // take the other processes into account
116 this->error_ = this->comm_.max(this->error_);
117
118 // make sure that the error never grows beyond the maximum
119 // allowed one
121 throw Opm::NumericalProblem("Newton: Error "+std::to_string(double(this->error_))+
122 + " is larger than maximum allowed error of "
124 }
125
129 void updatePrimaryVariables_(unsigned globalDofIdx,
130 PrimaryVariables& nextValue,
131 const PrimaryVariables& currentValue,
132 const EqVector& update,
133 const EqVector&)
134 {
135 // normal Newton-Raphson update
136 nextValue = currentValue;
137 nextValue -= update;
138
140 // put crash barriers along the update path
142
143 // saturations: limit the change of any saturation to at most 20%
144 Scalar sumSatDelta = 0.0;
145 Scalar maxSatDelta = 0.0;
146 for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
147 maxSatDelta = std::max(std::abs(update[saturation0Idx + phaseIdx]),
148 maxSatDelta);
149 sumSatDelta += update[saturation0Idx + phaseIdx];
150 }
151 maxSatDelta = std::max(std::abs(- sumSatDelta), maxSatDelta);
152
153 if (maxSatDelta > 0.2) {
154 Scalar alpha = 0.2/maxSatDelta;
155 for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
156 nextValue[saturation0Idx + phaseIdx] =
157 currentValue[saturation0Idx + phaseIdx]
158 - alpha*update[saturation0Idx + phaseIdx];
159 }
160 }
161
162 // limit pressure reference change to 20% of the total value per iteration
163 clampValue_(nextValue[pressure0Idx],
164 currentValue[pressure0Idx]*0.8,
165 currentValue[pressure0Idx]*1.2);
166
167 // fugacities
168 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
169 Scalar& val = nextValue[fugacity0Idx + compIdx];
170 Scalar oldVal = currentValue[fugacity0Idx + compIdx];
171
172 // get the minimum activity coefficient for the component (i.e., the activity
173 // coefficient of the phase for which the component has the highest affinity)
174 Scalar minPhi = this->problem().model().minActivityCoeff(globalDofIdx, compIdx);
175 // Make sure that the activity coefficient does not get too small.
176 minPhi = std::max(0.001*currentValue[pressure0Idx], minPhi);
177
178 // allow the mole fraction of the component to change at most 70% in any
179 // phase (assuming composition independent fugacity coefficients).
180 Scalar maxDelta = 0.7 * minPhi;
181 clampValue_(val, oldVal - maxDelta, oldVal + maxDelta);
182
183 // make sure that fugacities do not become negative
184 val = std::max(val, 0.0);
185 }
186
187 // do not become grossly unphysical in a single iteration for the first few
188 // iterations of a time step
189 if (this->numIterations_ < 3) {
190 // fugacities
191 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
192 Scalar& val = nextValue[fugacity0Idx + compIdx];
193 Scalar oldVal = currentValue[fugacity0Idx + compIdx];
194 Scalar minPhi = this->problem().model().minActivityCoeff(globalDofIdx, compIdx);
195 if (oldVal < 1.0*minPhi && val > 1.0*minPhi)
196 val = 1.0*minPhi;
197 else if (oldVal > 0.0 && val < 0.0)
198 val = 0.0;
199 }
200
201 // saturations
202 for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
203 Scalar& val = nextValue[saturation0Idx + phaseIdx];
204 Scalar oldVal = currentValue[saturation0Idx + phaseIdx];
205 if (oldVal < 1.0 && val > 1.0)
206 val = 1.0;
207 else if (oldVal > 0.0 && val < 0.0)
208 val = 0.0;
209 }
210 }
211 }
212
213private:
214 void clampValue_(Scalar& val, Scalar minVal, Scalar maxVal) const
215 { val = std::max(minVal, std::min(val, maxVal)); }
216};
217} // namespace Opm
218
219#endif
A Newton solver specific to the NCP model.
Definition: ncpnewtonmethod.hh:55
friend ParentType
Definition: ncpnewtonmethod.hh:82
NcpNewtonMethod(Simulator &simulator)
Definition: ncpnewtonmethod.hh:78
void updatePrimaryVariables_(unsigned globalDofIdx, PrimaryVariables &nextValue, const PrimaryVariables &currentValue, const EqVector &update, const EqVector &)
Update a single primary variables object.
Definition: ncpnewtonmethod.hh:129
void preSolve_(const SolutionVector &, const GlobalEqVector &currentResidual)
Definition: ncpnewtonmethod.hh:85
The multi-dimensional Newton method.
Definition: newtonmethod.hh:92
auto Get(bool errorIfNotRegistered=true)
Retrieve a runtime parameter.
Definition: parametersystem.hh:840
Definition: blackoilmodel.hh:72
Definition: blackoilboundaryratevector.hh:37
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:235
Declares the properties required for the NCP compositional multi-phase model.
Definition: newtonmethodparameters.hh:31