blackoilnewtonmethod.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_BLACK_OIL_NEWTON_METHOD_HH
29 #define EWOMS_BLACK_OIL_NEWTON_METHOD_HH
30 
31 #include "blackoilproperties.hh"
32 
33 #include <ewoms/common/signum.hh>
34 
35 #include <opm/common/Unused.hpp>
36 
37 namespace Ewoms {
38 
44 template <class TypeTag>
45 class BlackOilNewtonMethod : public GET_PROP_TYPE(TypeTag, DiscNewtonMethod)
46 {
47  typedef typename GET_PROP_TYPE(TypeTag, DiscNewtonMethod) ParentType;
48  typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
49  typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
50  typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) GlobalEqVector;
51  typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
52  typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
53  typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
54  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
55  typedef typename GET_PROP_TYPE(TypeTag, Linearizer) Linearizer;
56 
57  static const unsigned numEq = GET_PROP_VALUE(TypeTag, NumEq);
58 
59 public:
60  BlackOilNewtonMethod(Simulator& simulator) : ParentType(simulator)
61  { }
62 
66  static void registerParameters()
67  {
68  ParentType::registerParameters();
69  }
70 
75  unsigned numPriVarsSwitched() const
76  { return numPriVarsSwitched_; }
77 
78 protected:
79  friend NewtonMethod<TypeTag>;
80  friend ParentType;
81 
86  {
87  numPriVarsSwitched_ = 0;
88  ParentType::beginIteration_();
89  }
90 
94  void endIteration_(SolutionVector& uCurrentIter,
95  const SolutionVector& uLastIter)
96  {
97 #if HAVE_MPI
98  // in the MPI enabled case we need to add up the number of DOF
99  // for which the interpretation changed over all processes.
100  int localSwitched = numPriVarsSwitched_;
101  MPI_Allreduce(&localSwitched,
102  &numPriVarsSwitched_,
103  /*num=*/1,
104  MPI_INT,
105  MPI_SUM,
106  MPI_COMM_WORLD);
107 #endif // HAVE_MPI
108 
109  this->simulator_.model().newtonMethod().endIterMsg()
110  << ", num switched=" << numPriVarsSwitched_;
111 
112  ParentType::endIteration_(uCurrentIter, uLastIter);
113  }
114 
115  void update_(SolutionVector& nextSolution,
116  const SolutionVector& currentSolution,
117  const GlobalEqVector& solutionUpdate,
118  const GlobalEqVector& currentResidual)
119  {
120  const auto& comm = this->simulator_.gridView().comm();
121 
122  int succeeded;
123  try {
124  ParentType::update_(nextSolution,
125  currentSolution,
126  solutionUpdate,
127  currentResidual);
128  succeeded = 1;
129  }
130  catch (...) {
131  std::cout << "Newton update threw an exception on rank "
132  << comm.rank() << "\n";
133  succeeded = 0;
134  }
135  succeeded = comm.min(succeeded);
136 
137  if (!succeeded)
138  OPM_THROW(Opm::NumericalProblem,
139  "A process did not succeed in adapting the primary variables");
140 
141  numPriVarsSwitched_ = comm.sum(numPriVarsSwitched_);
142  }
143 
147  void updatePrimaryVariables_(unsigned globalDofIdx,
148  PrimaryVariables& nextValue,
149  const PrimaryVariables& currentValue,
150  const EqVector& update,
151  const EqVector& currentResidual OPM_UNUSED)
152  {
153  for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
154  // calculate the update of the current primary variable. For the
155  // black-oil model we limit the pressure and saturation updates, but do
156  // we not clamp anything after the specified number of iterations was
157  // reached
158  Scalar delta = update[eqIdx];
159 
160  // limit changes in water saturation to 20%
161  if (eqIdx == Indices::waterSaturationIdx
162  && std::abs(delta) > 0.2)
163  {
164  delta = Ewoms::signum(delta)*0.2;
165  }
166  else if (eqIdx == Indices::compositionSwitchIdx) {
167  // the switching primary variable for composition is tricky because the
168  // "reasonable" value ranges it exhibits vary widely depending on its
169  // interpretation (it can represent Sg, Rs or Rv). so far, we only limit
170  // changes in gas saturation to 20%
171  if (currentValue.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg
172  && std::abs(delta) > 0.2)
173  {
174  delta = Ewoms::signum(delta)*0.2;
175  }
176  }
177 
178  // do the actual update
179  nextValue[eqIdx] = currentValue[eqIdx] - delta;
180  }
181 
182  // switch the new primary variables to something which is physically meaningful
183  if (nextValue.adaptPrimaryVariables(this->problem(), globalDofIdx))
184  ++ numPriVarsSwitched_;
185  }
186 
187 private:
188  int numPriVarsSwitched_;
189 };
190 } // namespace Ewoms
191 
192 #endif
Definition: baseauxiliarymodule.hh:37
int signum(Scalar val)
Template function which returns the sign of a floating point value.
Definition: signum.hh:40
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
void endIteration_(SolutionVector &uCurrentIter, const SolutionVector &uLastIter)
Indicates that one Newton iteration was finished.
Definition: blackoilnewtonmethod.hh:94
A newton solver which is specific to the black oil model.
Definition: blackoilnewtonmethod.hh:45
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:486
void updatePrimaryVariables_(unsigned globalDofIdx, PrimaryVariables &nextValue, const PrimaryVariables &currentValue, const EqVector &update, const EqVector &currentResidual OPM_UNUSED)
Update a single primary variables object.
Definition: blackoilnewtonmethod.hh:147
Declares the properties required by the black oil model.
unsigned numPriVarsSwitched() const
Returns the number of degrees of freedom for which the interpretation has changed for the most recent...
Definition: blackoilnewtonmethod.hh:75
The multi-dimensional Newton method.
Definition: newtonmethod.hh:56
static void registerParameters()
Register all run-time parameters for the immiscible model.
Definition: blackoilnewtonmethod.hh:66
void beginIteration_()
Indicates the beginning of a Newton iteration.
Definition: blackoilnewtonmethod.hh:85
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:75