28#ifndef EWOMS_BLACK_OIL_INTENSIVE_QUANTITIES_HH
29#define EWOMS_BLACK_OIL_INTENSIVE_QUANTITIES_HH
42#include <opm/common/TimingMacros.hpp>
43#include <opm/common/OpmLog/OpmLog.hpp>
45#include <opm/input/eclipse/EclipseState/Grid/FaceDir.hpp>
47#include <opm/material/fluidstates/BlackOilFluidState.hpp>
48#include <opm/material/common/Valgrind.hpp>
52#include <opm/utility/CopyablePtr.hpp>
54#include <dune/common/fmatrix.hh>
59#include <fmt/format.h>
69template <
class TypeTag>
71 :
public GetPropType<TypeTag, Properties::DiscIntensiveQuantities>
72 ,
public GetPropType<TypeTag, Properties::FluxModule>::FluxIntensiveQuantities
96 enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
97 enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
98 enum { enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>() };
99 enum { enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>() };
100 enum { enableFoam = getPropValue<TypeTag, Properties::EnableFoam>() };
101 enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
102 enum { enableVapwat = getPropValue<TypeTag, Properties::EnableVapwat>() };
103 enum { has_disgas_in_water = getPropValue<TypeTag, Properties::EnableDisgasInWater>() };
104 enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
105 enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
106 enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
107 enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
108 enum { enableDispersion = getPropValue<TypeTag, Properties::EnableDispersion>() };
109 enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
110 enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
111 enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
112 enum { waterCompIdx = FluidSystem::waterCompIdx };
113 enum { oilCompIdx = FluidSystem::oilCompIdx };
114 enum { gasCompIdx = FluidSystem::gasCompIdx };
115 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
116 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
117 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
118 enum { dimWorld = GridView::dimensionworld };
119 enum { compositionSwitchIdx = Indices::compositionSwitchIdx };
121 static constexpr bool compositionSwitchEnabled = Indices::compositionSwitchIdx >= 0;
122 static constexpr bool waterEnabled = Indices::waterEnabled;
123 static constexpr bool gasEnabled = Indices::gasEnabled;
124 static constexpr bool oilEnabled = Indices::oilEnabled;
126 using Toolbox = MathToolbox<Evaluation>;
127 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
128 using FluxIntensiveQuantities =
typename FluxModule::FluxIntensiveQuantities;
132 using DirectionalMobilityPtr = Opm::Utility::CopyablePtr<DirectionalMobility<TypeTag, Evaluation>>;
141 compositionSwitchEnabled,
144 enableSaltPrecipitation,
151 compositionSwitchEnabled,
154 enableSaltPrecipitation,
161 if (compositionSwitchEnabled) {
162 fluidState_.setRs(0.0);
163 fluidState_.setRv(0.0);
166 fluidState_.setRvw(0.0);
168 if (has_disgas_in_water) {
169 fluidState_.setRsw(0.0);
179 void update(
const ElementContext& elemCtx,
unsigned dofIdx,
unsigned timeIdx)
181 ParentType::update(elemCtx, dofIdx, timeIdx);
182 OPM_TIMEBLOCK_LOCAL(blackoilIntensiveQuanititiesUpdate);
183 const auto& problem = elemCtx.problem();
184 const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
185 const auto& linearizationType = problem.model().linearizer().getLinearizationType();
186 unsigned globalSpaceIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
187 Scalar RvMax = FluidSystem::enableVaporizedOil()
188 ? problem.maxOilVaporizationFactor(timeIdx, globalSpaceIdx)
190 Scalar RsMax = FluidSystem::enableDissolvedGas()
191 ? problem.maxGasDissolutionFactor(timeIdx, globalSpaceIdx)
194 asImp_().updateTemperature_(elemCtx, dofIdx, timeIdx);
196 unsigned pvtRegionIdx = priVars.pvtRegionIndex();
197 fluidState_.setPvtRegionIndex(pvtRegionIdx);
199 asImp_().updateSaltConcentration_(elemCtx, dofIdx, timeIdx);
203 if constexpr (waterEnabled) {
204 if (priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
205 Sw = priVars.makeEvaluation(Indices::waterSwitchIdx, timeIdx);
206 }
else if(priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Rsw ||
207 priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Disabled) {
214 if constexpr (gasEnabled) {
215 if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg) {
216 Sg = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
217 }
else if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rv) {
219 }
else if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Disabled) {
220 if constexpr (waterEnabled) {
228 Valgrind::CheckDefined(Sg);
229 Valgrind::CheckDefined(Sw);
231 Evaluation So = 1.0 - Sw - Sg;
234 if constexpr (enableSolvent) {
235 if(priVars.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
236 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
237 So -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx);
238 }
else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
239 Sg -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx);
244 if (FluidSystem::phaseIsActive(waterPhaseIdx))
245 fluidState_.setSaturation(waterPhaseIdx, Sw);
247 if (FluidSystem::phaseIsActive(gasPhaseIdx))
248 fluidState_.setSaturation(gasPhaseIdx, Sg);
250 if (FluidSystem::phaseIsActive(oilPhaseIdx))
251 fluidState_.setSaturation(oilPhaseIdx, So);
253 asImp_().solventPreSatFuncUpdate_(elemCtx, dofIdx, timeIdx);
256 std::array<Evaluation, numPhases> pC;
257 const auto& materialParams = problem.materialLawParams(globalSpaceIdx);
258 MaterialLaw::capillaryPressures(pC, materialParams, fluidState_);
259 problem.updateRelperms(mobility_, dirMob_, fluidState_, globalSpaceIdx);
263 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, dofIdx, timeIdx);
264 const Evaluation Sp = priVars.makeEvaluation(Indices::saltConcentrationIdx, timeIdx);
265 const Evaluation porosityFactor = min(1.0 - Sp, 1.0);
267 const Evaluation pcFactor = pcfactTable.eval(porosityFactor,
true);
268 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
269 if (FluidSystem::phaseIsActive(phaseIdx)) {
270 pC[phaseIdx] *= pcFactor;
275 if (priVars.primaryVarsMeaningPressure() == PrimaryVariables::PressureMeaning::Pg) {
276 const Evaluation& pg = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
277 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
278 if (FluidSystem::phaseIsActive(phaseIdx))
279 fluidState_.setPressure(phaseIdx, pg + (pC[phaseIdx] - pC[gasPhaseIdx]));
280 }
else if (priVars.primaryVarsMeaningPressure() == PrimaryVariables::PressureMeaning::Pw) {
281 const Evaluation& pw = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
282 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
283 if (FluidSystem::phaseIsActive(phaseIdx))
284 fluidState_.setPressure(phaseIdx, pw + (pC[phaseIdx] - pC[waterPhaseIdx]));
286 assert(FluidSystem::phaseIsActive(oilPhaseIdx));
287 const Evaluation& po = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
288 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
289 if (FluidSystem::phaseIsActive(phaseIdx))
290 fluidState_.setPressure(phaseIdx, po + (pC[phaseIdx] - pC[oilPhaseIdx]));
294 asImp_().solventPostSatFuncUpdate_(elemCtx, dofIdx, timeIdx);
297 asImp_().zFractionUpdate_(elemCtx, dofIdx, timeIdx);
299 Evaluation SoMax = 0.0;
300 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
301 SoMax = max(fluidState_.saturation(oilPhaseIdx),
302 problem.maxOilSaturation(globalSpaceIdx));
307 if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rs) {
308 const auto& Rs = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
309 fluidState_.setRs(Rs);
311 if (FluidSystem::enableDissolvedGas()) {
312 const Evaluation& RsSat = enableExtbo ? asImp_().rs() :
313 FluidSystem::saturatedDissolutionFactor(fluidState_,
317 fluidState_.setRs(min(RsMax, RsSat));
319 else if constexpr (compositionSwitchEnabled)
320 fluidState_.setRs(0.0);
323 if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rv) {
324 const auto& Rv = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
325 fluidState_.setRv(Rv);
327 if (FluidSystem::enableVaporizedOil() ) {
328 const Evaluation& RvSat = enableExtbo ? asImp_().rv() :
329 FluidSystem::saturatedDissolutionFactor(fluidState_,
333 fluidState_.setRv(min(RvMax, RvSat));
335 else if constexpr (compositionSwitchEnabled)
336 fluidState_.setRv(0.0);
339 if (priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Rvw) {
340 const auto& Rvw = priVars.makeEvaluation(Indices::waterSwitchIdx, timeIdx);
341 fluidState_.setRvw(Rvw);
343 if (FluidSystem::enableVaporizedWater()) {
344 const Evaluation& RvwSat = FluidSystem::saturatedVaporizationFactor(fluidState_,
347 fluidState_.setRvw(RvwSat);
351 if (priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Rsw) {
352 const auto& Rsw = priVars.makeEvaluation(Indices::waterSwitchIdx, timeIdx);
353 fluidState_.setRsw(Rsw);
355 if (FluidSystem::enableDissolvedGasInWater()) {
356 const Evaluation& RswSat = FluidSystem::saturatedDissolutionFactor(fluidState_,
359 fluidState_.setRsw(RswSat);
363 typename FluidSystem::template ParameterCache<Evaluation> paramCache;
364 paramCache.setRegionIndex(pvtRegionIdx);
365 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
366 paramCache.setMaxOilSat(SoMax);
368 paramCache.updateAll(fluidState_);
372 std::vector<std::array<Evaluation,numPhases>*> mobilities = {&mobility_};
374 for (
int i=0; i<3; i++) {
376 mobilities.push_back(&(dirMob_->getArray(i)));
379 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
380 if (!FluidSystem::phaseIsActive(phaseIdx))
382 const auto& b = FluidSystem::inverseFormationVolumeFactor(fluidState_, phaseIdx, pvtRegionIdx);
383 fluidState_.setInvB(phaseIdx, b);
384 const auto& mu = FluidSystem::viscosity(fluidState_, paramCache, phaseIdx);
385 for (
int i = 0; i<nmobilities; i++) {
386 if (enableExtbo && phaseIdx == oilPhaseIdx) {
387 (*mobilities[i])[phaseIdx] /= asImp_().oilViscosity();
389 else if (enableExtbo && phaseIdx == gasPhaseIdx) {
390 (*mobilities[i])[phaseIdx] /= asImp_().gasViscosity();
393 (*mobilities[i])[phaseIdx] /= mu;
397 Valgrind::CheckDefined(mobility_);
401 if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
402 rho = fluidState_.invB(waterPhaseIdx);
403 rho *= FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
404 if (FluidSystem::enableDissolvedGasInWater()) {
406 fluidState_.invB(waterPhaseIdx) *
408 FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
410 fluidState_.setDensity(waterPhaseIdx, rho);
413 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
414 rho = fluidState_.invB(gasPhaseIdx);
415 rho *= FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
416 if (FluidSystem::enableVaporizedOil()) {
418 fluidState_.invB(gasPhaseIdx) *
420 FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
422 if (FluidSystem::enableVaporizedWater()) {
424 fluidState_.invB(gasPhaseIdx) *
426 FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
428 fluidState_.setDensity(gasPhaseIdx, rho);
431 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
432 rho = fluidState_.invB(oilPhaseIdx);
433 rho *= FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
434 if (FluidSystem::enableDissolvedGas()) {
436 fluidState_.invB(oilPhaseIdx) *
438 FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
440 fluidState_.setDensity(oilPhaseIdx, rho);
444 referencePorosity_ = problem.porosity(elemCtx, dofIdx, timeIdx);
445 porosity_ = referencePorosity_;
449 Scalar rockCompressibility = problem.rockCompressibility(globalSpaceIdx);
450 if (rockCompressibility > 0.0) {
451 Scalar rockRefPressure = problem.rockReferencePressure(globalSpaceIdx);
453 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
454 x = rockCompressibility*(fluidState_.pressure(oilPhaseIdx) - rockRefPressure);
455 }
else if (FluidSystem::phaseIsActive(waterPhaseIdx)){
456 x = rockCompressibility*(fluidState_.pressure(waterPhaseIdx) - rockRefPressure);
458 x = rockCompressibility*(fluidState_.pressure(gasPhaseIdx) - rockRefPressure);
460 porosity_ *= 1.0 + x + 0.5*x*x;
464 porosity_ *= problem.template rockCompPoroMultiplier<Evaluation>(*
this, globalSpaceIdx);
467 if constexpr (enableMICP){
468 Evaluation biofilm_ = priVars.makeEvaluation(Indices::biofilmConcentrationIdx, timeIdx, linearizationType);
469 Evaluation calcite_ = priVars.makeEvaluation(Indices::calciteConcentrationIdx, timeIdx, linearizationType);
470 porosity_ += - biofilm_ - calcite_;
474 if (enableSaltPrecipitation && priVars.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp) {
475 Evaluation Sp = priVars.makeEvaluation(Indices::saltConcentrationIdx, timeIdx);
476 porosity_ *= (1.0 - Sp);
479 rockCompTransMultiplier_ = problem.template rockCompTransMultiplier<Evaluation>(*
this, globalSpaceIdx);
481 asImp_().solventPvtUpdate_(elemCtx, dofIdx, timeIdx);
482 asImp_().zPvtUpdate_();
483 asImp_().polymerPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
484 asImp_().updateEnergyQuantities_(elemCtx, dofIdx, timeIdx, paramCache);
485 asImp_().foamPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
486 asImp_().MICPPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
487 asImp_().saltPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
491 FluxIntensiveQuantities::update_(elemCtx, dofIdx, timeIdx);
494 DiffusionIntensiveQuantities::update_(fluidState_, paramCache, elemCtx, dofIdx, timeIdx);
497 DispersionIntensiveQuantities::update_(elemCtx, dofIdx, timeIdx);
501 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
502 if (!FluidSystem::phaseIsActive(phaseIdx))
505 assert(isfinite(fluidState_.density(phaseIdx)));
506 assert(isfinite(fluidState_.saturation(phaseIdx)));
507 assert(isfinite(fluidState_.temperature(phaseIdx)));
508 assert(isfinite(fluidState_.pressure(phaseIdx)));
509 assert(isfinite(fluidState_.invB(phaseIdx)));
511 assert(isfinite(fluidState_.Rs()));
512 assert(isfinite(fluidState_.Rv()));
521 {
return fluidState_; }
526 const Evaluation&
mobility(
unsigned phaseIdx)
const
527 {
return mobility_[phaseIdx]; }
529 const Evaluation&
mobility(
unsigned phaseIdx, FaceDir::DirEnum facedir)
const
531 using Dir = FaceDir::DirEnum;
536 return dirMob_->mobilityX_[phaseIdx];
539 return dirMob_->mobilityY_[phaseIdx];
542 return dirMob_->mobilityZ_[phaseIdx];
544 throw std::runtime_error(
"Unexpected face direction");
548 return mobility_[phaseIdx];
557 {
return porosity_; }
563 {
return rockCompTransMultiplier_; }
574 {
return fluidState_.pvtRegionIndex(); }
582 return fluidState_.viscosity(phaseIdx)*
mobility(phaseIdx);
592 {
return referencePorosity_; }
603 Implementation& asImp_()
604 {
return *
static_cast<Implementation*
>(
this); }
607 Scalar referencePorosity_;
608 Evaluation porosity_;
609 Evaluation rockCompTransMultiplier_;
610 std::array<Evaluation,numPhases> mobility_;
627 DirectionalMobilityPtr dirMob_;
Contains the classes required to extend the black-oil model by brine.
Classes required for molecular diffusion.
Classes required for mechanical dispersion.
Contains the classes required to extend the black-oil model by energy.
Contains the classes required to extend the black-oil model by solvent component. For details,...
Contains the classes required to extend the black-oil model to include the effects of foam.
Contains the classes required to extend the black-oil model by MICP.
Contains the classes required to extend the black-oil model by polymer.
Declares the properties required by the black oil model.
Contains the classes required to extend the black-oil model by solvents.
Definition: blackoilbrinemodules.hh:416
Contains the high level supplements required to extend the black oil model by brine.
Definition: blackoilbrinemodules.hh:58
static const TabulatedFunction & pcfactTable(unsigned satnumRegionIdx)
Definition: blackoilbrinemodules.hh:342
static bool hasPcfactTables()
Definition: blackoilbrinemodules.hh:386
Provides the volumetric quantities required for the calculation of molecular diffusive fluxes.
Definition: blackoildiffusionmodule.hh:274
Provides the volumetric quantities required for the calculation of dispersive fluxes.
Definition: blackoildispersionmodule.hh:278
Provides the volumetric quantities required for the equations needed by the energys extension of the ...
Definition: blackoilenergymodules.hh:333
Provides the volumetric quantities required for the equations needed by the solvents extension of the...
Definition: blackoilextbomodules.hh:562
Provides the volumetric quantities required for the equations needed by the polymers extension of the...
Definition: blackoilfoammodules.hh:461
Contains the quantities which are are constant within a finite volume in the black-oil model.
Definition: blackoilintensivequantities.hh:82
const Evaluation & porosity() const
Returns the average porosity within the control volume.
Definition: blackoilintensivequantities.hh:556
const Evaluation & mobility(unsigned phaseIdx) const
Returns the effective mobility of a given phase within the control volume.
Definition: blackoilintensivequantities.hh:526
void update(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Definition: blackoilintensivequantities.hh:179
Evaluation relativePermeability(unsigned phaseIdx) const
Returns the relative permeability of a given phase within the control volume.
Definition: blackoilintensivequantities.hh:579
auto pvtRegionIndex() const -> decltype(std::declval< FluidState >().pvtRegionIndex())
Returns the index of the PVT region used to calculate the thermodynamic quantities.
Definition: blackoilintensivequantities.hh:572
const FluidState & fluidState() const
Returns the phase state for the control-volume.
Definition: blackoilintensivequantities.hh:520
BlackOilIntensiveQuantities & operator=(const BlackOilIntensiveQuantities &other)=default
const Evaluation & rockCompTransMultiplier() const
Definition: blackoilintensivequantities.hh:562
BlackOilIntensiveQuantities(const BlackOilIntensiveQuantities &other)=default
BlackOilIntensiveQuantities()
Definition: blackoilintensivequantities.hh:159
BlackOilFluidState< Scalar, FluidSystem, enableTemperature, enableEnergy, compositionSwitchEnabled, enableVapwat, enableBrine, enableSaltPrecipitation, has_disgas_in_water, Indices::numPhases > ScalarFluidState
Definition: blackoilintensivequantities.hh:156
const Evaluation & mobility(unsigned phaseIdx, FaceDir::DirEnum facedir) const
Definition: blackoilintensivequantities.hh:529
Scalar referencePorosity() const
Returns the porosity of the rock at reference conditions.
Definition: blackoilintensivequantities.hh:591
BlackOilFluidState< Evaluation, FluidSystem, enableTemperature, enableEnergy, compositionSwitchEnabled, enableVapwat, enableBrine, enableSaltPrecipitation, has_disgas_in_water, Indices::numPhases > FluidState
Definition: blackoilintensivequantities.hh:146
Provides the volumetric quantities required for the equations needed by the MICP extension of the bla...
Definition: blackoilmicpmodules.hh:459
Provides the volumetric quantities required for the equations needed by the polymers extension of the...
Definition: blackoilpolymermodules.hh:807
Provides the volumetric quantities required for the equations needed by the solvents extension of the...
Definition: blackoilsolventmodules.hh:779
This file contains definitions related to directional mobilities.
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:242