All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
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/Unused.hpp>
34 #include <opm/common/ErrorMacros.hpp>
35 #include <opm/common/Exceptions.hpp>
36 
37 #include <algorithm>
38 
39 namespace Ewoms {
40 
46 template <class TypeTag>
47 class NcpNewtonMethod : public GET_PROP_TYPE(TypeTag, DiscNewtonMethod)
48 {
49  typedef typename GET_PROP_TYPE(TypeTag, DiscNewtonMethod) ParentType;
50 
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;
55  typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
56  typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
57  typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) GlobalEqVector;
58 
59  enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
60  enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
61  enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
62  enum { fugacity0Idx = Indices::fugacity0Idx };
63  enum { saturation0Idx = Indices::saturation0Idx };
64  enum { pressure0Idx = Indices::pressure0Idx };
65  enum { ncp0EqIdx = Indices::ncp0EqIdx };
66 
67 public:
71  NcpNewtonMethod(Simulator& simulator) : ParentType(simulator)
72  {
73  Dune::FMatrixPrecision<Scalar>::set_singular_limit(1e-35);
74  }
75 
76 protected:
77  friend ParentType;
78  friend NewtonMethod<TypeTag>;
79 
80  void preSolve_(const SolutionVector& currentSolution OPM_UNUSED,
81  const GlobalEqVector& currentResidual)
82  {
83  const auto& constraintsMap = this->model().linearizer().constraintsMap();
84  this->lastError_ = this->error_;
85 
86  // calculate the error as the maximum weighted tolerance of
87  // the solution's residual
88  this->error_ = 0;
89  for (unsigned dofIdx = 0; dofIdx < currentResidual.size(); ++dofIdx) {
90  // do not consider auxiliary DOFs for the error
91  if (dofIdx >= this->model().numGridDof() || this->model().dofTotalVolume(dofIdx) <= 0.0)
92  continue;
93 
94  // also do not consider DOFs which are constraint
95  if (this->enableConstraints_()) {
96  if (constraintsMap.count(dofIdx) > 0)
97  continue;
98  }
99 
100  const auto& r = currentResidual[dofIdx];
101  for (unsigned eqIdx = 0; eqIdx < r.size(); ++eqIdx) {
102  if (ncp0EqIdx <= eqIdx && eqIdx < Indices::ncp0EqIdx + numPhases)
103  continue;
104  this->error_ =
105  std::max(std::abs(r[eqIdx]*this->model().eqWeight(dofIdx, eqIdx)),
106  this->error_);
107  }
108  }
109 
110  // take the other processes into account
111  this->error_ = this->comm_.max(this->error_);
112 
113  // make sure that the error never grows beyond the maximum
114  // allowed one
115  if (this->error_ > EWOMS_GET_PARAM(TypeTag, Scalar, NewtonMaxError))
116  OPM_THROW(Opm::NumericalProblem,
117  "Newton: Error " << this->error_
118  << " is larger than maximum allowed error of "
119  << EWOMS_GET_PARAM(TypeTag, Scalar, NewtonMaxError));
120  }
121 
125  void updatePrimaryVariables_(unsigned globalDofIdx,
126  PrimaryVariables& nextValue,
127  const PrimaryVariables& currentValue,
128  const EqVector& update,
129  const EqVector& currentResidual OPM_UNUSED)
130  {
131  // normal Newton-Raphson update
132  nextValue = currentValue;
133  nextValue -= update;
134 
136  // put crash barriers along the update path
138 
139  // saturations: limit the change of any saturation to at most 20%
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]),
144  maxSatDelta);
145  sumSatDelta += update[saturation0Idx + phaseIdx];
146  }
147  maxSatDelta = std::max(std::abs(- sumSatDelta), maxSatDelta);
148 
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];
155  }
156  }
157 
158  // limit pressure reference change to 20% of the total value per iteration
159  clampValue_(nextValue[pressure0Idx],
160  currentValue[pressure0Idx]*0.8,
161  currentValue[pressure0Idx]*1.2);
162 
163  // fugacities
164  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
165  Scalar& val = nextValue[fugacity0Idx + compIdx];
166  Scalar oldVal = currentValue[fugacity0Idx + compIdx];
167 
168  // allow the mole fraction of the component to change
169  // at most 70% (assuming composition independent
170  // fugacity coefficients)
171  Scalar minPhi = this->problem().model().minActivityCoeff(globalDofIdx, compIdx);
172  Scalar maxDelta = 0.7 * minPhi;
173 
174  clampValue_(val, oldVal - maxDelta, oldVal + maxDelta);
175  }
176 
177  // do not become grossly unphysical in a single iteration for the first few
178  // iterations of a time step
179  if (this->numIterations_ < 3) {
180  // fugacities
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)
186  val = 1.0*minPhi;
187  else if (oldVal > 0.0 && val < 0.0)
188  val = 0.0;
189  }
190 
191  // saturations
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)
196  val = 1.0;
197  else if (oldVal > 0.0 && val < 0.0)
198  val = 0.0;
199  }
200  }
201  }
202 
203 private:
204  void clampValue_(Scalar& val, Scalar minVal, Scalar maxVal) const
205  { val = std::max(minVal, std::min(val, maxVal)); }
206 };
207 } // namespace Ewoms
208 
209 #endif
void updatePrimaryVariables_(unsigned globalDofIdx, PrimaryVariables &nextValue, const PrimaryVariables &currentValue, const EqVector &update, const EqVector &currentResidual OPM_UNUSED)
Update a single primary variables object.
Definition: ncpnewtonmethod.hh:125
#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