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
35#include <dune/common/fvector.hh>
36
37namespace Opm {
38
44template <class TypeTag>
45class RichardsNewtonMethod : public GetPropType<TypeTag, Properties::DiscNewtonMethod>
46{
48
57
59 enum { pressureWIdx = Indices::pressureWIdx };
60 enum { numPhases = FluidSystem::numPhases };
61 enum { liquidPhaseIdx = getPropValue<TypeTag, Properties::LiquidPhaseIndex>() };
62 enum { gasPhaseIdx = getPropValue<TypeTag, Properties::GasPhaseIndex>() };
63
64 using PhaseVector = Dune::FieldVector<Scalar, numPhases>;
65
66public:
67 RichardsNewtonMethod(Simulator& simulator) : ParentType(simulator)
68 {}
69
70protected:
72 friend ParentType;
73
77 void updatePrimaryVariables_(unsigned globalDofIdx,
78 PrimaryVariables& nextValue,
79 const PrimaryVariables& currentValue,
80 const EqVector& update,
81 const EqVector&)
82 {
83 // normal Newton-Raphson update
84 nextValue = currentValue;
85 nextValue -= update;
86
87 // do not clamp anything after 4 iterations
88 if (this->numIterations_ > 4)
89 return;
90
91 const auto& problem = this->simulator_.problem();
92
93 // calculate the old wetting phase saturation
94 const MaterialLawParams& matParams =
95 problem.materialLawParams(globalDofIdx, /*timeIdx=*/0);
96
97 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
98
99 // set the temperature
100 Scalar T = problem.temperature(globalDofIdx, /*timeIdx=*/0);
101 fs.setTemperature(T);
102
104 // calculate the phase pressures of the previous iteration
106
107 // first, we have to find the minimum capillary pressure
108 // (i.e. Sw = 0)
109 fs.setSaturation(liquidPhaseIdx, 1.0);
110 fs.setSaturation(gasPhaseIdx, 0.0);
111 PhaseVector pC;
112 MaterialLaw::capillaryPressures(pC, matParams, fs);
113
114 // non-wetting pressure can be larger than the
115 // reference pressure if the medium is fully
116 // saturated by the wetting phase
117 Scalar pWOld = currentValue[pressureWIdx];
118 Scalar pNOld =
119 std::max(problem.referencePressure(globalDofIdx, /*timeIdx=*/0),
120 pWOld + (pC[gasPhaseIdx] - pC[liquidPhaseIdx]));
121
123 // find the saturations of the previous iteration
125 fs.setPressure(liquidPhaseIdx, pWOld);
126 fs.setPressure(gasPhaseIdx, pNOld);
127
128 PhaseVector satOld;
129 MaterialLaw::saturations(satOld, matParams, fs);
130 satOld[liquidPhaseIdx] = std::max<Scalar>(0.0, satOld[liquidPhaseIdx]);
131
133 // find the wetting phase pressures which
134 // corrospond to a 20% increase and a 20% decrease
135 // of the wetting saturation
137 fs.setSaturation(liquidPhaseIdx, satOld[liquidPhaseIdx] - 0.2);
138 fs.setSaturation(gasPhaseIdx, 1.0 - (satOld[liquidPhaseIdx] - 0.2));
139 MaterialLaw::capillaryPressures(pC, matParams, fs);
140 Scalar pwMin = pNOld - (pC[gasPhaseIdx] - pC[liquidPhaseIdx]);
141
142 fs.setSaturation(liquidPhaseIdx, satOld[liquidPhaseIdx] + 0.2);
143 fs.setSaturation(gasPhaseIdx, 1.0 - (satOld[liquidPhaseIdx] + 0.2));
144 MaterialLaw::capillaryPressures(pC, matParams, fs);
145 Scalar pwMax = pNOld - (pC[gasPhaseIdx] - pC[liquidPhaseIdx]);
146
148 // clamp the result to the minimum and the maximum
149 // pressures we just calculated
151 Scalar pW = nextValue[pressureWIdx];
152 pW = std::max(pwMin, std::min(pW, pwMax));
153 nextValue[pressureWIdx] = pW;
154 }
155};
156} // namespace Opm
157
158#endif
The multi-dimensional Newton method.
Definition: newtonmethod.hh:92
A Richards model specific Newton method.
Definition: richardsnewtonmethod.hh:46
RichardsNewtonMethod(Simulator &simulator)
Definition: richardsnewtonmethod.hh:67
void updatePrimaryVariables_(unsigned globalDofIdx, PrimaryVariables &nextValue, const PrimaryVariables &currentValue, const EqVector &update, const EqVector &)
Update a single primary variables object.
Definition: richardsnewtonmethod.hh:77
friend ParentType
Definition: richardsnewtonmethod.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
Contains the property declarations for the Richards model.