28 #ifndef EWOMS_RICHARDS_NEWTON_METHOD_HH
29 #define EWOMS_RICHARDS_NEWTON_METHOD_HH
33 #include <opm/material/fluidstates/ImmiscibleFluidState.hpp>
34 #include <opm/common/Unused.hpp>
36 #include <dune/common/fvector.hh>
45 template <
class TypeTag>
48 typedef typename GET_PROP_TYPE(TypeTag, DiscNewtonMethod) ParentType;
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;
57 typedef typename GET_PROP_TYPE(TypeTag, Linearizer) Linearizer;
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) };
65 typedef Dune::FieldVector<Scalar, numPhases> PhaseVector;
79 PrimaryVariables& nextValue,
80 const PrimaryVariables& currentValue,
81 const EqVector& update,
82 const EqVector& currentResidual OPM_UNUSED)
85 nextValue = currentValue;
89 if (this->numIterations_ > 4)
92 const auto& problem = this->simulator_.problem();
95 const MaterialLawParams& matParams =
96 problem.materialLawParams(globalDofIdx, 0);
98 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
101 Scalar T = problem.temperature(globalDofIdx, 0);
102 fs.setTemperature(T);
110 fs.setSaturation(liquidPhaseIdx, 1.0);
111 fs.setSaturation(gasPhaseIdx, 0.0);
113 MaterialLaw::capillaryPressures(pC, matParams, fs);
118 Scalar pWOld = currentValue[pressureWIdx];
120 std::max(problem.referencePressure(globalDofIdx, 0),
121 pWOld + (pC[gasPhaseIdx] - pC[liquidPhaseIdx]));
126 fs.setPressure(liquidPhaseIdx, pWOld);
127 fs.setPressure(gasPhaseIdx, pNOld);
130 MaterialLaw::saturations(satOld, matParams, fs);
131 satOld[liquidPhaseIdx] = std::max<Scalar>(0.0, satOld[liquidPhaseIdx]);
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]);
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]);
152 Scalar pW = nextValue[pressureWIdx];
153 pW = std::max(pwMin, std::min(pW, pwMax));
154 nextValue[pressureWIdx] = pW;
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 ¤tValue, const EqVector &update, const EqVector ¤tResidual OPM_UNUSED)
Update a single primary variables object.
Definition: richardsnewtonmethod.hh:78