26 #ifndef EWOMS_PVS_INTENSIVE_QUANTITIES_HH
27 #define EWOMS_PVS_INTENSIVE_QUANTITIES_HH
34 #include <opm/material/constraintsolvers/ComputeFromReferencePhase.hpp>
35 #include <opm/material/constraintsolvers/MiscibleMultiPhaseComposition.hpp>
36 #include <opm/material/fluidstates/CompositionalFluidState.hpp>
38 #include <dune/common/fvector.hh>
39 #include <dune/common/fmatrix.hh>
52 template <
class TypeTag>
57 ,
public GET_PROP_TYPE(TypeTag, FluxModule)::FluxIntensiveQuantities
59 typedef typename GET_PROP_TYPE(TypeTag, DiscIntensiveQuantities) ParentType;
61 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
62 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
63 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
64 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
65 typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
66 typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
67 typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
68 typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
69 typedef typename GET_PROP_TYPE(TypeTag, FluxModule) FluxModule;
71 enum { switch0Idx = Indices::switch0Idx };
72 enum { pressure0Idx = Indices::pressure0Idx };
75 enum { enableDiffusion =
GET_PROP_VALUE(TypeTag, EnableDiffusion) };
77 enum { dimWorld = GridView::dimensionworld };
79 typedef Opm::MathToolbox<Evaluation> Toolbox;
80 typedef Opm::MiscibleMultiPhaseComposition<Scalar, FluidSystem, Evaluation> MiscibleMultiPhaseComposition;
81 typedef Opm::ComputeFromReferencePhase<Scalar, FluidSystem, Evaluation> ComputeFromReferencePhase;
83 typedef Dune::FieldVector<Scalar, numPhases> PhaseVector;
84 typedef Dune::FieldVector<Evaluation, numPhases> EvalPhaseVector;
85 typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
87 typedef typename FluxModule::FluxIntensiveQuantities FluxIntensiveQuantities;
93 typedef Opm::CompositionalFluidState<Evaluation, FluidSystem>
FluidState;
98 void update(
const ElementContext &elemCtx,
int dofIdx,
int timeIdx)
100 ParentType::update(elemCtx, dofIdx, timeIdx);
101 EnergyIntensiveQuantities::updateTemperatures_(fluidState_, elemCtx, dofIdx, timeIdx);
103 const auto &priVars = elemCtx.primaryVars(dofIdx, timeIdx);
104 const auto &problem = elemCtx.problem();
109 Evaluation sumSat = Toolbox::createConstant(0.0);
110 for (
int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
111 fluidState_.setSaturation(phaseIdx, priVars.explicitSaturationValue(phaseIdx, timeIdx));
112 Valgrind::CheckDefined(fluidState_.saturation(phaseIdx));
113 sumSat += fluidState_.saturation(phaseIdx);
115 Valgrind::CheckDefined(priVars.implicitSaturationIdx());
116 Valgrind::CheckDefined(sumSat);
117 fluidState_.setSaturation(priVars.implicitSaturationIdx(), 1.0 - sumSat);
124 const MaterialLawParams &materialParams =
125 problem.materialLawParams(elemCtx, dofIdx, timeIdx);
127 MaterialLaw::capillaryPressures(pC, materialParams, fluidState_);
130 const Evaluation& p0 = priVars.makeEvaluation(pressure0Idx, timeIdx);
131 for (
int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
132 fluidState_.setPressure(phaseIdx, p0 + (pC[phaseIdx] - pC[0]));
138 typename FluidSystem::ParameterCache paramCache;
139 int lowestPresentPhaseIdx = priVars.lowestPresentPhaseIdx();
140 int numNonPresentPhases = 0;
141 for (
int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
142 if (!priVars.phaseIsPresent(phaseIdx))
143 ++numNonPresentPhases;
147 if (numNonPresentPhases == numPhases - 1) {
150 Evaluation sumx = Toolbox::createConstant(0.0);
151 for (
int compIdx = 1; compIdx < numComponents; ++compIdx) {
152 const Evaluation& x = priVars.makeEvaluation(switch0Idx + compIdx - 1, timeIdx);
153 fluidState_.setMoleFraction(lowestPresentPhaseIdx, compIdx, x);
158 fluidState_.setMoleFraction(lowestPresentPhaseIdx, 0, 1 - sumx);
163 ComputeFromReferencePhase::solve(fluidState_, paramCache,
164 lowestPresentPhaseIdx,
170 int numAuxConstraints = numComponents + numNonPresentPhases - numPhases;
171 Opm::MMPCAuxConstraint<Evaluation> auxConstraints[numComponents];
175 for (; switchIdx < numPhases - 1; ++switchIdx) {
176 int compIdx = switchIdx + 1;
177 int switchPhaseIdx = switchIdx;
178 if (switchIdx >= lowestPresentPhaseIdx)
181 if (!priVars.phaseIsPresent(switchPhaseIdx)) {
182 auxConstraints[auxIdx].set(lowestPresentPhaseIdx, compIdx,
183 priVars.makeEvaluation(switch0Idx + switchIdx, timeIdx));
188 for (; auxIdx < numAuxConstraints; ++auxIdx, ++switchIdx) {
189 int compIdx = numPhases - numNonPresentPhases + auxIdx;
190 auxConstraints[auxIdx].set(lowestPresentPhaseIdx, compIdx,
191 priVars.makeEvaluation(switch0Idx + switchIdx, timeIdx));
197 MiscibleMultiPhaseComposition::solve(fluidState_, paramCache,
198 priVars.phasePresence(),
208 Scalar myNan = std::numeric_limits<Scalar>::quiet_NaN();
209 for (
int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
210 fluidState_.setEnthalpy(phaseIdx, myNan);
219 MaterialLaw::relativePermeabilities(relativePermeability_,
220 materialParams, fluidState_);
221 Valgrind::CheckDefined(relativePermeability_);
224 for (
int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
225 mobility_[phaseIdx] =
226 relativePermeability_[phaseIdx] /
fluidState().viscosity(phaseIdx);
229 porosity_ = problem.porosity(elemCtx, dofIdx, timeIdx);
230 Valgrind::CheckDefined(porosity_);
233 intrinsicPerm_ = problem.intrinsicPermeability(elemCtx, dofIdx, timeIdx);
236 FluxIntensiveQuantities::update_(elemCtx, dofIdx, timeIdx);
239 EnergyIntensiveQuantities::update_(fluidState_, paramCache, elemCtx, dofIdx, timeIdx);
242 DiffusionIntensiveQuantities::update_(fluidState_, paramCache, elemCtx, dofIdx, timeIdx);
244 fluidState_.checkDefined();
251 {
return fluidState_; }
257 {
return intrinsicPerm_; }
263 {
return relativePermeability_[phaseIdx]; }
269 {
return mobility_[phaseIdx]; }
275 {
return porosity_; }
278 FluidState fluidState_;
279 Evaluation porosity_;
280 DimMatrix intrinsicPerm_;
281 Evaluation relativePermeability_[numPhases];
282 Evaluation mobility_[numPhases];
Provides the volumetric quantities required for the energy equation.
Definition: energymodule.hh:530
Classes required for molecular diffusion.
Contains the quantities which are are constant within a finite volume in the compositional multi-phas...
Definition: pvsintensivequantities.hh:53
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:468
void update(const ElementContext &elemCtx, int dofIdx, int timeIdx)
Definition: pvsintensivequantities.hh:98
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:485
const FluidState & fluidState() const
Returns the phase state for the control-volume.
Definition: pvsintensivequantities.hh:250
const Evaluation & porosity() const
Returns the average porosity within the control volume.
Definition: pvsintensivequantities.hh:274
Definition: baseauxiliarymodule.hh:35
const DimMatrix & intrinsicPermeability() const
Returns the intrinsic permeability tensor a degree of freedom.
Definition: pvsintensivequantities.hh:256
const Evaluation & mobility(int phaseIdx) const
Returns the effective mobility of a given phase within the control volume.
Definition: pvsintensivequantities.hh:268
Provides the volumetric quantities required for the calculation of molecular diffusive fluxes...
Definition: diffusionmodule.hh:138
const Evaluation & relativePermeability(int phaseIdx) const
Returns the relative permeability of a given phase within the control volume.
Definition: pvsintensivequantities.hh:262
Opm::CompositionalFluidState< Evaluation, FluidSystem > FluidState
The type of the object returned by the fluidState() method.
Definition: pvsintensivequantities.hh:93
Declares the properties required for the compositional multi-phase primary variable switching model...
Contains the classes required to consider energy as a conservation quantity in a multi-phase module...