28 #ifndef EWOMS_BLACK_OIL_PRIMARY_VARIABLES_HH 29 #define EWOMS_BLACK_OIL_PRIMARY_VARIABLES_HH 37 #include <dune/common/fvector.hh> 39 #include <opm/material/constraintsolvers/NcpFlash.hpp> 40 #include <opm/material/fluidstates/CompositionalFluidState.hpp> 41 #include <opm/material/fluidstates/SimpleModularFluidState.hpp> 42 #include <opm/material/fluidsystems/BlackOilFluidSystem.hpp> 43 #include <opm/common/Valgrind.hpp> 46 template <
class TypeTag,
bool enableSolvent>
49 template <
class TypeTag,
bool enablePolymer>
57 template <
class TypeTag>
61 typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) Implementation;
63 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
64 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
65 typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
66 typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
67 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
68 typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
69 typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
75 enum { waterSaturationIdx = Indices::waterSaturationIdx };
76 enum { pressureSwitchIdx = Indices::pressureSwitchIdx };
77 enum { compositionSwitchIdx = Indices::compositionSwitchIdx };
79 static const bool compositionSwitchEnabled = Indices::gasEnabled;
80 static const bool waterEnabled = Indices::waterEnabled;
84 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
85 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
86 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
92 enum { gasCompIdx = FluidSystem::gasCompIdx };
93 enum { waterCompIdx = FluidSystem::waterCompIdx };
94 enum { oilCompIdx = FluidSystem::oilCompIdx };
96 typedef typename Opm::MathToolbox<Evaluation> Toolbox;
97 typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
102 static_assert(numPhases == 3,
"The black-oil model assumes three phases!");
103 static_assert(numComponents == 3,
"The black-oil model assumes three components!");
106 enum PrimaryVarsMeaning {
115 Opm::Valgrind::SetUndefined(*
this);
125 Opm::Valgrind::SetUndefined(primaryVarsMeaning_);
143 { pvtRegionIdx_ =
static_cast<unsigned short>(value); }
149 {
return pvtRegionIdx_; }
156 {
return primaryVarsMeaning_; }
163 { primaryVarsMeaning_ = newMeaning; }
168 template <
class Flu
idState>
170 const MaterialLawParams& matParams,
171 bool isInEquilibrium =
false)
173 typedef typename std::remove_reference<typename FluidState::Scalar>::type ConstEvaluation;
174 typedef typename std::remove_const<ConstEvaluation>::type FsEvaluation;
175 typedef typename Opm::MathToolbox<FsEvaluation> FsToolbox;
179 for (
unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
180 Opm::Valgrind::CheckDefined(fluidState.temperature(0));
181 Opm::Valgrind::CheckDefined(fluidState.temperature(phaseIdx));
183 assert(fluidState.temperature(0) == fluidState.temperature(phaseIdx));
189 if (isInEquilibrium) {
196 typename FluidSystem::template ParameterCache<Scalar> paramCache;
197 paramCache.setRegionIndex(pvtRegionIdx_);
198 paramCache.setMaxOilSat(FsToolbox::value(fluidState.saturation(oilPhaseIdx)));
201 typedef Opm::NcpFlash<Scalar, FluidSystem> NcpFlash;
202 typedef Opm::CompositionalFluidState<Scalar, FluidSystem> FlashFluidState;
203 FlashFluidState fsFlash;
204 fsFlash.setTemperature(FsToolbox::value(fluidState.temperature(0)));
205 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
206 fsFlash.setPressure(phaseIdx, FsToolbox::value(fluidState.pressure(phaseIdx)));
207 fsFlash.setSaturation(phaseIdx, FsToolbox::value(fluidState.saturation(phaseIdx)));
208 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
209 fsFlash.setMoleFraction(phaseIdx, compIdx, FsToolbox::value(fluidState.moleFraction(phaseIdx, compIdx)));
212 paramCache.updateAll(fsFlash);
213 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
214 if (!FluidSystem::phaseIsActive(phaseIdx))
217 Scalar rho = FluidSystem::template density<FlashFluidState, Scalar>(fsFlash, paramCache, phaseIdx);
218 fsFlash.setDensity(phaseIdx, rho);
222 ComponentVector globalMolarities(0.0);
223 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
224 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
225 if (!FluidSystem::phaseIsActive(phaseIdx))
228 globalMolarities[compIdx] +=
229 fsFlash.saturation(phaseIdx) * fsFlash.molarity(phaseIdx, compIdx);
237 NcpFlash::template solve<MaterialLaw>(fsFlash, matParams, paramCache, globalMolarities);
243 template <
class Flu
idState,
class SolventContainer>
245 const MaterialLawParams& matParams,
247 bool isInEquilibrium =
false)
261 template <
class Flu
idState>
264 typedef typename std::remove_reference<typename FluidState::Scalar>::type ConstEvaluation;
265 typedef typename std::remove_const<ConstEvaluation>::type FsEvaluation;
266 typedef typename Opm::MathToolbox<FsEvaluation> FsToolbox;
268 bool gasPresent = (fluidState.saturation(gasPhaseIdx) > 0.0);
269 bool oilPresent = (fluidState.saturation(oilPhaseIdx) > 0.0);
272 if ((gasPresent && oilPresent) || (!gasPresent && !oilPresent))
275 primaryVarsMeaning_ = Sw_po_Sg;
276 else if (oilPresent) {
279 if (FluidSystem::enableDissolvedGas())
280 primaryVarsMeaning_ = Sw_po_Rs;
282 primaryVarsMeaning_ = Sw_po_Sg;
288 if (FluidSystem::enableVaporizedOil())
289 primaryVarsMeaning_ = Sw_pg_Rv;
291 primaryVarsMeaning_ = Sw_po_Sg;
297 (*this)[waterSaturationIdx] = FsToolbox::value(fluidState.saturation(waterPhaseIdx));
298 (*this)[pressureSwitchIdx] = FsToolbox::value(fluidState.pressure(oilPhaseIdx));
299 if( compositionSwitchEnabled )
300 (*this)[compositionSwitchIdx] = FsToolbox::value(fluidState.saturation(gasPhaseIdx));
303 const auto& Rs = Opm::BlackOil::getRs_<FluidSystem, Scalar, FluidState>(fluidState, pvtRegionIdx_);
306 (*this)[waterSaturationIdx] = FsToolbox::value(fluidState.saturation(waterPhaseIdx));
307 (*this)[pressureSwitchIdx] = FsToolbox::value(fluidState.pressure(oilPhaseIdx));
308 if( compositionSwitchEnabled )
309 (*this)[compositionSwitchIdx] = Rs;
314 const auto& Rv = Opm::BlackOil::getRv_<FluidSystem, Scalar, FluidState>(fluidState, pvtRegionIdx_);
316 (*this)[waterSaturationIdx] = FsToolbox::value(fluidState.saturation(waterPhaseIdx));
318 (*this)[pressureSwitchIdx] = FsToolbox::value(fluidState.pressure(gasPhaseIdx));
319 if( compositionSwitchEnabled )
320 (*this)[compositionSwitchIdx] = Rv;
339 static const Scalar thresholdWaterFilledCell = 1.0 - eps;
348 Sw = (*this)[Indices::waterSaturationIdx];
353 if (Sw >= thresholdWaterFilledCell) {
357 (*this)[Indices::waterSaturationIdx] = 1.0;
359 if (compositionSwitchEnabled)
360 (*this)[Indices::compositionSwitchIdx] = 0.0;
367 if (compositionSwitchEnabled)
368 Sg = (*this)[Indices::compositionSwitchIdx];
370 Scalar So = 1.0 - Sw - Sg - solventSaturation();
372 Scalar So2 = 1.0 - Sw - solventSaturation();
373 if (Sg < -eps && So2 > 0.0 && FluidSystem::enableDissolvedGas()) {
381 Scalar po = (*this)[Indices::pressureSwitchIdx];
382 Scalar T = asImp_().temperature_();
383 Scalar SoMax = problem.model().maxOilSaturation(globalDofIdx);
384 Scalar RsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_,
390 if (compositionSwitchEnabled)
391 (*this)[Indices::compositionSwitchIdx] = RsSat;
395 asImp_().adaptPrimaryVariables(problem, globalDofIdx, eps);
400 Scalar Sg2 = 1.0 - Sw - solventSaturation();
401 if (So < -eps && Sg2 > 0.0 && FluidSystem::enableVaporizedOil()) {
404 Scalar po = (*this)[Indices::pressureSwitchIdx];
408 Scalar pC[numPhases] = { 0.0 };
409 const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
410 computeCapillaryPressures_(pC, 0.0, Sg2 + solventSaturation(), Sw, matParams);
411 Scalar pg = po + (pC[gasPhaseIdx] - pC[oilPhaseIdx]);
415 Scalar T = asImp_().temperature_();
416 Scalar SoMax = problem.model().maxOilSaturation(globalDofIdx);
418 FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_,
425 (*this)[Indices::pressureSwitchIdx] = pg;
426 if (compositionSwitchEnabled)
427 (*this)[Indices::compositionSwitchIdx] = RvSat;
431 asImp_().adaptPrimaryVariables(problem, globalDofIdx, eps);
440 assert(compositionSwitchEnabled);
443 if (Sw >= thresholdWaterFilledCell) {
448 (*this)[Indices::waterSaturationIdx] = 1.0;
450 (*this)[Indices::compositionSwitchIdx] = 0.0;
454 asImp_().adaptPrimaryVariables(problem, globalDofIdx, eps);
462 Scalar T = asImp_().temperature_();
463 Scalar po = (*this)[Indices::pressureSwitchIdx];
464 Scalar So = 1.0 - Sw - solventSaturation();
465 Scalar SoMax = std::max(So, problem.model().maxOilSaturation(globalDofIdx));
467 FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, T, po, So, SoMax);
468 Scalar Rs = (*this)[Indices::compositionSwitchIdx];
469 if (Rs > RsSat*(1.0 + eps)) {
473 (*this)[Indices::compositionSwitchIdx] = 0.0;
477 asImp_().adaptPrimaryVariables(problem, globalDofIdx, eps);
486 assert(compositionSwitchEnabled);
488 Scalar pg = (*this)[Indices::pressureSwitchIdx];
489 Scalar Sg = 1.0 - Sw - solventSaturation();
492 if (Sw >= thresholdWaterFilledCell) {
496 Scalar pC[numPhases] = { 0.0 };
497 const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
498 computeCapillaryPressures_(pC,
500 Sg + solventSaturation(),
503 Scalar po = pg + (pC[oilPhaseIdx] - pC[gasPhaseIdx]);
507 (*this)[Indices::waterSaturationIdx] = 1.0;
509 (*this)[Indices::pressureSwitchIdx] = po;
510 (*this)[Indices::compositionSwitchIdx] = 0.0;
514 asImp_().adaptPrimaryVariables(problem, globalDofIdx, eps);
523 Scalar T = asImp_().temperature_();
524 Scalar SoMax = problem.model().maxOilSaturation(globalDofIdx);
526 FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_,
532 Scalar Rv = (*this)[Indices::compositionSwitchIdx];
533 if (Rv > RvSat*(1.0 + eps)) {
537 Scalar pC[numPhases] = { 0.0 };
538 const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
539 computeCapillaryPressures_(pC,
541 Sg + solventSaturation(),
544 Scalar po = pg + (pC[oilPhaseIdx] - pC[gasPhaseIdx]);
547 (*this)[Indices::pressureSwitchIdx] = po;
548 (*this)[Indices::compositionSwitchIdx] = Sg;
552 asImp_().adaptPrimaryVariables(problem, globalDofIdx, eps);
567 for (
unsigned i = 0; i < numEq; ++i)
574 Implementation& asImp_()
575 {
return *
static_cast<Implementation*
>(
this); }
577 const Implementation& asImp_()
const 578 {
return *
static_cast<const Implementation*
>(
this); }
580 Scalar solventSaturation()
const 585 return (*
this)[Indices::solventSaturationIdx];
588 Scalar polymerConcentration()
const 593 return (*
this)[Indices::polymerConcentrationIdx];
596 template <
class Container>
597 void computeCapillaryPressures_(Container& result,
601 const MaterialLawParams& matParams)
const 603 typedef Opm::SimpleModularFluidState<Scalar,
614 false> SatOnlyFluidState;
615 SatOnlyFluidState fluidState;
616 fluidState.setSaturation(waterPhaseIdx, Sw);
617 fluidState.setSaturation(oilPhaseIdx, So);
618 fluidState.setSaturation(gasPhaseIdx, Sg);
620 MaterialLaw::capillaryPressures(result, matParams, fluidState);
624 Scalar temperature_()
const 625 {
return FluidSystem::surfaceTemperature; }
627 PrimaryVarsMeaning primaryVarsMeaning_;
628 unsigned short pvtRegionIdx_;
Definition: baseauxiliarymodule.hh:37
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
Contains the classes required to extend the black-oil model by solvents.
static void assignPrimaryVars(PrimaryVariables &priVars, Scalar polymerConcentration)
Assign the polymer specific primary variables to a PrimaryVariables object.
Definition: blackoilpolymermodules.hh:511
unsigned pvtRegionIndex() const
Return the index of the region which should be used for PVT properties.
Definition: blackoilprimaryvariables.hh:148
Declares the properties required by the black oil model.
bool adaptPrimaryVariables(const Problem &problem, unsigned globalDofIdx, Scalar eps=0.0)
Adapt the interpretation of the switching variables to be physically meaningful.
Definition: blackoilprimaryvariables.hh:337
Contains the high level supplements required to extend the black oil model by solvents.
Definition: blackoilprimaryvariables.hh:47
Represents the primary variables used by the a model.
Definition: fvbaseprimaryvariables.hh:48
Contains the high level supplements required to extend the black oil model by polymer.
Definition: blackoilpolymermodules.hh:64
PrimaryVarsMeaning primaryVarsMeaning() const
Return the interpretation which should be applied to the switching primary variables.
Definition: blackoilprimaryvariables.hh:155
void setPvtRegionIndex(unsigned value)
Set the index of the region which should be used for PVT properties.
Definition: blackoilprimaryvariables.hh:142
void setPrimaryVarsMeaning(PrimaryVarsMeaning newMeaning)
Set the interpretation which should be applied to the switching primary variables.
Definition: blackoilprimaryvariables.hh:162
Represents the primary variables used by the a model.
void assignMassConservative(const FluidState &fluidState, const MaterialLawParams &matParams, bool isInEquilibrium=false)
Set the primary variables from an arbitrary fluid state in a mass conservative way.
Definition: blackoilprimaryvariables.hh:169
Represents the primary variables used by the black-oil model.
Definition: blackoilprimaryvariables.hh:58
void assignNaive(const FluidState &fluidState)
Directly retrieve the primary variables from an arbitrary fluid state.
Definition: blackoilprimaryvariables.hh:262
BlackOilPrimaryVariables(Scalar value)
Constructor with assignment from scalar.
Definition: blackoilprimaryvariables.hh:122
static void assignPrimaryVars(PrimaryVariables &priVars, Scalar solventSaturation)
Assign the solvent specific primary variables to a PrimaryVariables object.
Definition: blackoilsolventmodules.hh:615
Contains the classes required to extend the black-oil model by polymer.