All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
richardsnewtonmethod.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_RICHARDS_NEWTON_METHOD_HH
29 #define EWOMS_RICHARDS_NEWTON_METHOD_HH
30 
31 #include "richardsproperties.hh"
32 
33 #include <opm/material/fluidstates/ImmiscibleFluidState.hpp>
34 #include <opm/common/Unused.hpp>
35 
36 #include <dune/common/fvector.hh>
37 
38 namespace Ewoms {
39 
45 template <class TypeTag>
46 class RichardsNewtonMethod : public GET_PROP_TYPE(TypeTag, DiscNewtonMethod)
47 {
48  typedef typename GET_PROP_TYPE(TypeTag, DiscNewtonMethod) ParentType;
49 
50  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
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, FluidSystem) FluidSystem;
54  typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
55  typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
56  typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
57  typedef typename GET_PROP_TYPE(TypeTag, Linearizer) Linearizer;
58 
59  typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
60  enum { pressureWIdx = Indices::pressureWIdx };
61  enum { numPhases = FluidSystem::numPhases };
62  enum { liquidPhaseIdx = GET_PROP_VALUE(TypeTag, LiquidPhaseIndex) };
63  enum { gasPhaseIdx = GET_PROP_VALUE(TypeTag, GasPhaseIndex) };
64 
65  typedef Dune::FieldVector<Scalar, numPhases> PhaseVector;
66 
67 public:
68  RichardsNewtonMethod(Simulator& simulator) : ParentType(simulator)
69  {}
70 
71 protected:
72  friend NewtonMethod<TypeTag>;
73  friend ParentType;
74 
78  void updatePrimaryVariables_(unsigned globalDofIdx,
79  PrimaryVariables& nextValue,
80  const PrimaryVariables& currentValue,
81  const EqVector& update,
82  const EqVector& currentResidual OPM_UNUSED)
83  {
84  // normal Newton-Raphson update
85  nextValue = currentValue;
86  nextValue -= update;
87 
88  // do not clamp anything after 4 iterations
89  if (this->numIterations_ > 4)
90  return;
91 
92  const auto& problem = this->simulator_.problem();
93 
94  // calculate the old wetting phase saturation
95  const MaterialLawParams& matParams =
96  problem.materialLawParams(globalDofIdx, /*timeIdx=*/0);
97 
98  Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
99 
100  // set the temperature
101  Scalar T = problem.temperature(globalDofIdx, /*timeIdx=*/0);
102  fs.setTemperature(T);
103 
105  // calculate the phase pressures of the previous iteration
107 
108  // first, we have to find the minimum capillary pressure
109  // (i.e. Sw = 0)
110  fs.setSaturation(liquidPhaseIdx, 1.0);
111  fs.setSaturation(gasPhaseIdx, 0.0);
112  PhaseVector pC;
113  MaterialLaw::capillaryPressures(pC, matParams, fs);
114 
115  // non-wetting pressure can be larger than the
116  // reference pressure if the medium is fully
117  // saturated by the wetting phase
118  Scalar pWOld = currentValue[pressureWIdx];
119  Scalar pNOld =
120  std::max(problem.referencePressure(globalDofIdx, /*timeIdx=*/0),
121  pWOld + (pC[gasPhaseIdx] - pC[liquidPhaseIdx]));
122 
124  // find the saturations of the previous iteration
126  fs.setPressure(liquidPhaseIdx, pWOld);
127  fs.setPressure(gasPhaseIdx, pNOld);
128 
129  PhaseVector satOld;
130  MaterialLaw::saturations(satOld, matParams, fs);
131  satOld[liquidPhaseIdx] = std::max<Scalar>(0.0, satOld[liquidPhaseIdx]);
132 
134  // find the wetting phase pressures which
135  // corrospond to a 20% increase and a 20% decrease
136  // of the wetting saturation
138  fs.setSaturation(liquidPhaseIdx, satOld[liquidPhaseIdx] - 0.2);
139  fs.setSaturation(gasPhaseIdx, 1.0 - (satOld[liquidPhaseIdx] - 0.2));
140  MaterialLaw::capillaryPressures(pC, matParams, fs);
141  Scalar pwMin = pNOld - (pC[gasPhaseIdx] - pC[liquidPhaseIdx]);
142 
143  fs.setSaturation(liquidPhaseIdx, satOld[liquidPhaseIdx] + 0.2);
144  fs.setSaturation(gasPhaseIdx, 1.0 - (satOld[liquidPhaseIdx] + 0.2));
145  MaterialLaw::capillaryPressures(pC, matParams, fs);
146  Scalar pwMax = pNOld - (pC[gasPhaseIdx] - pC[liquidPhaseIdx]);
147 
149  // clamp the result to the minimum and the maximum
150  // pressures we just calculated
152  Scalar pW = nextValue[pressureWIdx];
153  pW = std::max(pwMin, std::min(pW, pwMax));
154  nextValue[pressureWIdx] = pW;
155  }
156 };
157 } // namespace Ewoms
158 
159 #endif
A Richards model specific Newton method.
Definition: richardsnewtonmethod.hh:46
#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
Contains the property declarations for the Richards model.
The multi-dimensional Newton method.
Definition: newtonmethod.hh:56
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:75
void updatePrimaryVariables_(unsigned globalDofIdx, PrimaryVariables &nextValue, const PrimaryVariables &currentValue, const EqVector &update, const EqVector &currentResidual OPM_UNUSED)
Update a single primary variables object.
Definition: richardsnewtonmethod.hh:78