28 #ifndef EWOMS_NCP_NEWTON_METHOD_HH 29 #define EWOMS_NCP_NEWTON_METHOD_HH 33 #include <opm/common/Unused.hpp> 34 #include <opm/common/ErrorMacros.hpp> 35 #include <opm/common/Exceptions.hpp> 46 template <
class TypeTag>
49 typedef typename GET_PROP_TYPE(TypeTag, DiscNewtonMethod) ParentType;
51 typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
52 typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
53 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
54 typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
56 typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
57 typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) GlobalEqVector;
62 enum { fugacity0Idx = Indices::fugacity0Idx };
63 enum { saturation0Idx = Indices::saturation0Idx };
64 enum { pressure0Idx = Indices::pressure0Idx };
65 enum { ncp0EqIdx = Indices::ncp0EqIdx };
73 Dune::FMatrixPrecision<Scalar>::set_singular_limit(1e-35);
80 void preSolve_(
const SolutionVector& currentSolution OPM_UNUSED,
81 const GlobalEqVector& currentResidual)
83 const auto& constraintsMap = this->model().linearizer().constraintsMap();
84 this->lastError_ = this->error_;
89 for (
unsigned dofIdx = 0; dofIdx < currentResidual.size(); ++dofIdx) {
91 if (dofIdx >= this->model().numGridDof() || this->model().dofTotalVolume(dofIdx) <= 0.0)
95 if (this->enableConstraints_()) {
96 if (constraintsMap.count(dofIdx) > 0)
100 const auto& r = currentResidual[dofIdx];
101 for (
unsigned eqIdx = 0; eqIdx < r.size(); ++eqIdx) {
102 if (ncp0EqIdx <= eqIdx && eqIdx < Indices::ncp0EqIdx + numPhases)
105 std::max(std::abs(r[eqIdx]*this->model().eqWeight(dofIdx, eqIdx)),
111 this->error_ = this->comm_.max(this->error_);
116 OPM_THROW(Opm::NumericalProblem,
117 "Newton: Error " << this->error_
118 <<
" is larger than maximum allowed error of " 126 PrimaryVariables& nextValue,
127 const PrimaryVariables& currentValue,
128 const EqVector& update,
129 const EqVector& currentResidual OPM_UNUSED)
132 nextValue = currentValue;
140 Scalar sumSatDelta = 0.0;
141 Scalar maxSatDelta = 0.0;
142 for (
unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
143 maxSatDelta = std::max(std::abs(update[saturation0Idx + phaseIdx]),
145 sumSatDelta += update[saturation0Idx + phaseIdx];
147 maxSatDelta = std::max(std::abs(- sumSatDelta), maxSatDelta);
149 if (maxSatDelta > 0.2) {
150 Scalar alpha = 0.2/maxSatDelta;
151 for (
unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
152 nextValue[saturation0Idx + phaseIdx] =
153 currentValue[saturation0Idx + phaseIdx]
154 - alpha*update[saturation0Idx + phaseIdx];
159 clampValue_(nextValue[pressure0Idx],
160 currentValue[pressure0Idx]*0.8,
161 currentValue[pressure0Idx]*1.2);
164 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
165 Scalar& val = nextValue[fugacity0Idx + compIdx];
166 Scalar oldVal = currentValue[fugacity0Idx + compIdx];
171 Scalar minPhi = this->problem().model().minActivityCoeff(globalDofIdx, compIdx);
172 Scalar maxDelta = 0.7 * minPhi;
174 clampValue_(val, oldVal - maxDelta, oldVal + maxDelta);
179 if (this->numIterations_ < 3) {
181 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
182 Scalar& val = nextValue[fugacity0Idx + compIdx];
183 Scalar oldVal = currentValue[fugacity0Idx + compIdx];
184 Scalar minPhi = this->problem().model().minActivityCoeff(globalDofIdx, compIdx);
185 if (oldVal < 1.0*minPhi && val > 1.0*minPhi)
187 else if (oldVal > 0.0 && val < 0.0)
192 for (
unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
193 Scalar& val = nextValue[saturation0Idx + phaseIdx];
194 Scalar oldVal = currentValue[saturation0Idx + phaseIdx];
195 if (oldVal < 1.0 && val > 1.0)
197 else if (oldVal > 0.0 && val < 0.0)
204 void clampValue_(Scalar& val, Scalar minVal, Scalar maxVal)
const 205 { val = std::max(minVal, std::min(val, maxVal)); }
void updatePrimaryVariables_(unsigned globalDofIdx, PrimaryVariables &nextValue, const PrimaryVariables ¤tValue, const EqVector &update, const EqVector ¤tResidual OPM_UNUSED)
Update a single primary variables object.
Definition: ncpnewtonmethod.hh:125
Definition: baseauxiliarymodule.hh:37
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:486
Declares the properties required for the NCP compositional multi-phase model.
#define EWOMS_GET_PARAM(TypeTag, ParamType, ParamName)
Retrieve a runtime parameter.
Definition: parametersystem.hh:99
NcpNewtonMethod(Simulator &simulator)
Definition: ncpnewtonmethod.hh:71
The multi-dimensional Newton method.
Definition: newtonmethod.hh:56
A Newton solver specific to the NCP model.
Definition: ncpnewtonmethod.hh:47
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:75