26 #ifndef EWOMS_RICHARDS_NEWTON_METHOD_HH
27 #define EWOMS_RICHARDS_NEWTON_METHOD_HH
31 #include <opm/material/fluidstates/ImmiscibleFluidState.hpp>
33 #include <dune/common/fvector.hh>
42 template <
class TypeTag>
45 typedef typename GET_PROP_TYPE(TypeTag, DiscNewtonMethod) ParentType;
47 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
48 typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
49 typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
50 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
51 typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
52 typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
54 typedef typename GET_PROP_TYPE(TypeTag, Linearizer) Linearizer;
56 typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
57 enum { pressureWIdx = Indices::pressureWIdx };
58 enum { numPhases = FluidSystem::numPhases };
59 enum { liquidPhaseIdx =
GET_PROP_VALUE(TypeTag, LiquidPhaseIndex) };
62 typedef Dune::FieldVector<Scalar, numPhases> PhaseVector;
79 PrimaryVariables& nextValue,
80 const PrimaryVariables& currentValue,
81 const EqVector& update,
82 const EqVector& currentResidual)
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;
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:468
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:485
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:73
Definition: baseauxiliarymodule.hh:35
void updatePrimaryVariables_(int globalDofIdx, PrimaryVariables &nextValue, const PrimaryVariables ¤tValue, const EqVector &update, const EqVector ¤tResidual)
Update a single primary variables object.
Definition: richardsnewtonmethod.hh:78
RichardsNewtonMethod(Simulator &simulator)
Definition: richardsnewtonmethod.hh:65
A Richards model specific Newton method.
Definition: richardsnewtonmethod.hh:43
Contains the property declarations for the Richards model.