OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge 8341381c42 into 363d2eb902

Browse Source
This commit is contained in:
Tobias Meyer Andersen 2025-02-13 13:10:42 +01:00 committed by GitHub
commit b05604b256
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
11 changed files with 157 additions and 108 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

69
flowexperimental/BlackOilIntensiveQuantitiesGlobalIndex.hpp
View File

@ -154,28 +154,31 @@ public:
/*!
* \copydoc IntensiveQuantities::update
*/
void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
template <class DynamicFluidSystem = FluidSystem>
void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx, const DynamicFluidSystem& fluidSystem = DynamicFluidSystem{})
{
ParentType::update(elemCtx, dofIdx, timeIdx);
const auto& problem = elemCtx.problem();
const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
unsigned globalSpaceIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
this->update(problem,priVars,globalSpaceIdx,timeIdx);
this->update(problem,priVars,globalSpaceIdx,timeIdx, fluidSystem);
}
template <class DynamicFluidSystem = FluidSystem>
void update(const Problem& problem,
const PrimaryVariables& priVars,
unsigned globalSpaceIdx,
unsigned timeIdx)
unsigned timeIdx,
const DynamicFluidSystem& fluidSystem = DynamicFluidSystem{})
{
this->extrusionFactor_ = 1.0;// to avoid fixing parent update
OPM_TIMEBLOCK_LOCAL(UpdateIntensiveQuantities);
Scalar RvMax = FluidSystem::enableVaporizedOil()
Scalar RvMax = fluidSystem.enableVaporizedOil()
? problem.maxOilVaporizationFactor(timeIdx, globalSpaceIdx)
: 0.0;
Scalar RsMax = FluidSystem::enableDissolvedGas()
Scalar RsMax = fluidSystem.enableDissolvedGas()
? problem.maxGasDissolutionFactor(timeIdx, globalSpaceIdx)
: 0.0;
@ -221,15 +224,15 @@ public:
So -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx);
}
if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
if (fluidSystem.phaseIsActive(waterPhaseIdx)) {
fluidState_.setSaturation(waterPhaseIdx, Sw);
}
if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
if (fluidSystem.phaseIsActive(gasPhaseIdx)) {
fluidState_.setSaturation(gasPhaseIdx, Sg);
}
if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
if (fluidSystem.phaseIsActive(oilPhaseIdx)) {
fluidState_.setSaturation(oilPhaseIdx, So);
}
@ -239,21 +242,21 @@ public:
if (priVars.primaryVarsMeaningPressure() == PrimaryVariables::PressureMeaning::Pg) {
const Evaluation& pg = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
if (FluidSystem::phaseIsActive(phaseIdx)) {
if (fluidSystem.phaseIsActive(phaseIdx)) {
fluidState_.setPressure(phaseIdx, pg + (pC[phaseIdx] - pC[gasPhaseIdx]));
}
}
} else {
const Evaluation& po = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
if (FluidSystem::phaseIsActive(phaseIdx)) {
if (fluidSystem.phaseIsActive(phaseIdx)) {
fluidState_.setPressure(phaseIdx, po + (pC[phaseIdx] - pC[oilPhaseIdx]));
}
}
}
Evaluation SoMax = 0.0;
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
if (fluidSystem.phaseIsActive(fluidSystem.oilPhaseIdx)) {
SoMax = max(fluidState_.saturation(oilPhaseIdx),
problem.maxOilSaturation(globalSpaceIdx));
}
@ -264,10 +267,10 @@ public:
const auto& Rs = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
fluidState_.setRs(Rs);
} else {
if (FluidSystem::enableDissolvedGas()) { // Add So > 0? i.e. if only water set rs = 0)
if (fluidSystem.enableDissolvedGas()) { // Add So > 0? i.e. if only water set rs = 0)
OPM_TIMEBLOCK_LOCAL(UpdateSaturatedRs);
const Evaluation& RsSat = enableExtbo ? asImp_().rs() :
FluidSystem::saturatedDissolutionFactor(fluidState_,
fluidSystem.saturatedDissolutionFactor(fluidState_,
oilPhaseIdx,
pvtRegionIdx,
SoMax);
@ -282,11 +285,11 @@ public:
const auto& Rv = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
fluidState_.setRv(Rv);
} else {
if (FluidSystem::enableVaporizedOil() ) { // Add Sg > 0? i.e. if only water set rv = 0)
if (fluidSystem.enableVaporizedOil() ) { // Add Sg > 0? i.e. if only water set rv = 0)
OPM_TIMEBLOCK_LOCAL(UpdateSaturatedRv);
//NB! should save the indexing for later evalustion
const Evaluation& RvSat = enableExtbo ? asImp_().rv() :
FluidSystem::saturatedDissolutionFactor(fluidState_,
fluidSystem.saturatedDissolutionFactor(fluidState_,
gasPhaseIdx,
pvtRegionIdx,
SoMax);
@ -302,9 +305,9 @@ public:
fluidState_.setRvw(Rvw);
} else {
//NB! should save the indexing for later evaluation
if (FluidSystem::enableVaporizedWater()) { // Add Sg > 0? i.e. if only water set rv = 0)
if (fluidSystem.enableVaporizedWater()) { // Add Sg > 0? i.e. if only water set rv = 0)
OPM_TIMEBLOCK_LOCAL(UpdateSaturatedRv);
const Evaluation& RvwSat = FluidSystem::saturatedVaporizationFactor(fluidState_,
const Evaluation& RvwSat = fluidSystem.saturatedVaporizationFactor(fluidState_,
gasPhaseIdx,
pvtRegionIdx);
fluidState_.setRvw(RvwSat);
@ -312,7 +315,7 @@ public:
}
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
if (!fluidSystem.phaseIsActive(phaseIdx)) {
continue;
}
computeInverseFormationVolumeFactorAndViscosity(fluidState_, phaseIdx, pvtRegionIdx, SoMax);
@ -321,38 +324,38 @@ public:
// calculate the phase densities
Evaluation rho;
if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
if (fluidSystem.phaseIsActive(waterPhaseIdx)) {
OPM_TIMEBLOCK_LOCAL(UpdateWDensity);
rho = fluidState_.invB(waterPhaseIdx);
rho *= FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
rho *= fluidSystem.referenceDensity(waterPhaseIdx, pvtRegionIdx);
fluidState_.setDensity(waterPhaseIdx, rho);
}
if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
if (fluidSystem.phaseIsActive(gasPhaseIdx)) {
OPM_TIMEBLOCK_LOCAL(UpdateGDensity);
rho = fluidState_.invB(gasPhaseIdx);
rho *= FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
if (FluidSystem::enableVaporizedOil()) {
rho *= fluidSystem.referenceDensity(gasPhaseIdx, pvtRegionIdx);
if (fluidSystem.enableVaporizedOil()) {
rho += fluidState_.invB(gasPhaseIdx) *
fluidState_.Rv() *
FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
fluidSystem.referenceDensity(oilPhaseIdx, pvtRegionIdx);
}
if (FluidSystem::enableVaporizedWater()) {
if (fluidSystem.enableVaporizedWater()) {
rho += fluidState_.invB(gasPhaseIdx) *
fluidState_.Rvw() *
FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
fluidSystem.referenceDensity(waterPhaseIdx, pvtRegionIdx);
}
fluidState_.setDensity(gasPhaseIdx, rho);
}
if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
if (fluidSystem.phaseIsActive(oilPhaseIdx)) {
OPM_TIMEBLOCK_LOCAL(UpdateODensity);
rho = fluidState_.invB(oilPhaseIdx);
rho *= FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
if (FluidSystem::enableDissolvedGas()) {
rho *= fluidSystem.referenceDensity(oilPhaseIdx, pvtRegionIdx);
if (fluidSystem.enableDissolvedGas()) {
rho += fluidState_.invB(oilPhaseIdx) *
fluidState_.Rs() *
FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
fluidSystem.referenceDensity(gasPhaseIdx, pvtRegionIdx);
}
fluidState_.setDensity(oilPhaseIdx, rho);
}
@ -367,9 +370,9 @@ public:
OPM_TIMEBLOCK_LOCAL(UpdateRockCompressibility);
Scalar rockRefPressure = problem.rockReferencePressure(globalSpaceIdx);
Evaluation x;
if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
if (fluidSystem.phaseIsActive(oilPhaseIdx)) {
x = rockCompressibility*(fluidState_.pressure(oilPhaseIdx) - rockRefPressure);
} else if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
} else if (fluidSystem.phaseIsActive(waterPhaseIdx)) {
x = rockCompressibility*(fluidState_.pressure(waterPhaseIdx) - rockRefPressure);
} else {
x = rockCompressibility*(fluidState_.pressure(gasPhaseIdx) - rockRefPressure);
@ -385,7 +388,7 @@ public:
#ifndef NDEBUG
// some safety checks in debug mode
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
if (!fluidSystem.phaseIsActive(phaseIdx)) {
continue;
}

10
opm/models/blackoil/blackoilenergymodules.hh
View File

@ -384,10 +384,12 @@ public:
* \brief Compute the intensive quantities needed to handle energy conservation
*
*/
template <class FluidSystemT = FluidSystem>
void updateEnergyQuantities_(const ElementContext& elemCtx,
unsigned dofIdx,
unsigned timeIdx,
const typename FluidSystem::template ParameterCache<Evaluation>& paramCache)
const typename FluidSystem::template ParameterCache<Evaluation>& paramCache,
const FluidSystemT& fluidSystem = FluidSystemT{})
{
auto& fs = asImp_().fluidState_;
@ -398,7 +400,7 @@ public:
continue;
}
const auto& h = FluidSystem::enthalpy(fs, paramCache, phaseIdx);
const auto& h = fluidSystem.enthalpy(fs, paramCache, phaseIdx);
fs.setEnthalpy(phaseIdx, h);
}
@ -479,10 +481,12 @@ public:
}
}
template <class FluidSystemT = FluidSystem>
void updateEnergyQuantities_(const ElementContext&,
unsigned,
unsigned,
const typename FluidSystem::template ParameterCache<Evaluation>&)
const typename FluidSystem::template ParameterCache<Evaluation>&,
const FluidSystemT&)
{ }
const Evaluation& rockInternalEnergy() const

125
opm/models/blackoil/blackoilintensivequantities.hh
View File

@ -182,7 +182,9 @@ public:
}
}
void updateSaturations(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
template <class DynamicFluidSystem = FluidSystem>
void updateSaturations(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx, const DynamicFluidSystem& fluidSystem = DynamicFluidSystem{})
{
const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
@ -227,25 +229,26 @@ public:
// deal with solvent
if constexpr (enableSolvent) {
if(priVars.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
if (fluidSystem.phaseIsActive(oilPhaseIdx)) {
So -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx);
} else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
} else if (fluidSystem.phaseIsActive(gasPhaseIdx)) {
Sg -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx);
}
}
}
if (FluidSystem::phaseIsActive(waterPhaseIdx))
if (fluidSystem.phaseIsActive(waterPhaseIdx))
fluidState_.setSaturation(waterPhaseIdx, Sw);
if (FluidSystem::phaseIsActive(gasPhaseIdx))
if (fluidSystem.phaseIsActive(gasPhaseIdx))
fluidState_.setSaturation(gasPhaseIdx, Sg);
if (FluidSystem::phaseIsActive(oilPhaseIdx))
if (fluidSystem.phaseIsActive(oilPhaseIdx))
fluidState_.setSaturation(oilPhaseIdx, So);
}
void updateRelpermAndPressures(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
template <class DynamicFluidSystem = FluidSystem>
void updateRelpermAndPressures(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx, const DynamicFluidSystem& fluidSystem = DynamicFluidSystem{})
{
const auto& problem = elemCtx.problem();
const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
@ -276,7 +279,7 @@ public:
const auto& pcfactTable = BrineModule::pcfactTable(satnumRegionIdx);
const Evaluation pcFactor = pcfactTable.eval(porosityFactor, /*extrapolation=*/true);
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
if (FluidSystem::phaseIsActive(phaseIdx)) {
if (fluidSystem.phaseIsActive(phaseIdx)) {
pC[phaseIdx] *= pcFactor;
}
}
@ -286,18 +289,18 @@ public:
if (priVars.primaryVarsMeaningPressure() == PrimaryVariables::PressureMeaning::Pg) {
const Evaluation& pg = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
if (FluidSystem::phaseIsActive(phaseIdx))
if (fluidSystem.phaseIsActive(phaseIdx))
fluidState_.setPressure(phaseIdx, pg + (pC[phaseIdx] - pC[gasPhaseIdx]));
} else if (priVars.primaryVarsMeaningPressure() == PrimaryVariables::PressureMeaning::Pw) {
const Evaluation& pw = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
if (FluidSystem::phaseIsActive(phaseIdx))
if (fluidSystem.phaseIsActive(phaseIdx))
fluidState_.setPressure(phaseIdx, pw + (pC[phaseIdx] - pC[waterPhaseIdx]));
} else {
assert(FluidSystem::phaseIsActive(oilPhaseIdx));
assert(fluidSystem.phaseIsActive(oilPhaseIdx));
const Evaluation& po = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
if (FluidSystem::phaseIsActive(phaseIdx))
if (fluidSystem.phaseIsActive(phaseIdx))
fluidState_.setPressure(phaseIdx, po + (pC[phaseIdx] - pC[oilPhaseIdx]));
}
@ -311,25 +314,26 @@ public:
}
Evaluation updateRsRvRsw(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
template <class DynamicFluidSystem = FluidSystem>
Evaluation updateRsRvRsw(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx, const DynamicFluidSystem& fluidSystem = DynamicFluidSystem{})
{
const auto& problem = elemCtx.problem();
const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
const unsigned globalSpaceIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
const unsigned pvtRegionIdx = priVars.pvtRegionIndex();
Scalar RvMax = FluidSystem::enableVaporizedOil()
Scalar RvMax = fluidSystem.enableVaporizedOil()
? problem.maxOilVaporizationFactor(timeIdx, globalSpaceIdx)
: 0.0;
Scalar RsMax = FluidSystem::enableDissolvedGas()
Scalar RsMax = fluidSystem.enableDissolvedGas()
? problem.maxGasDissolutionFactor(timeIdx, globalSpaceIdx)
: 0.0;
Scalar RswMax = FluidSystem::enableDissolvedGasInWater()
Scalar RswMax = fluidSystem.enableDissolvedGasInWater()
? problem.maxGasDissolutionFactor(timeIdx, globalSpaceIdx)
: 0.0;
Evaluation SoMax = 0.0;
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
if (fluidSystem.phaseIsActive(fluidSystem.oilPhaseIdx)) {
SoMax = max(fluidState_.saturation(oilPhaseIdx),
problem.maxOilSaturation(globalSpaceIdx));
}
@ -339,9 +343,9 @@ public:
const auto& Rs = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
fluidState_.setRs(Rs);
} else {
if (FluidSystem::enableDissolvedGas()) { // Add So > 0? i.e. if only water set rs = 0)
if (fluidSystem.enableDissolvedGas()) { // Add So > 0? i.e. if only water set rs = 0)
const Evaluation& RsSat = enableExtbo ? asImp_().rs() :
FluidSystem::saturatedDissolutionFactor(fluidState_,
fluidSystem.saturatedDissolutionFactor(fluidState_,
oilPhaseIdx,
pvtRegionIdx,
SoMax);
@ -354,9 +358,9 @@ public:
const auto& Rv = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
fluidState_.setRv(Rv);
} else {
if (FluidSystem::enableVaporizedOil() ) { // Add Sg > 0? i.e. if only water set rv = 0)
if (fluidSystem.enableVaporizedOil() ) { // Add Sg > 0? i.e. if only water set rv = 0)
const Evaluation& RvSat = enableExtbo ? asImp_().rv() :
FluidSystem::saturatedDissolutionFactor(fluidState_,
fluidSystem.saturatedDissolutionFactor(fluidState_,
gasPhaseIdx,
pvtRegionIdx,
SoMax);
@ -370,8 +374,8 @@ public:
const auto& Rvw = priVars.makeEvaluation(Indices::waterSwitchIdx, timeIdx);
fluidState_.setRvw(Rvw);
} else {
if (FluidSystem::enableVaporizedWater()) { // Add Sg > 0? i.e. if only water set rv = 0)
const Evaluation& RvwSat = FluidSystem::saturatedVaporizationFactor(fluidState_,
if (fluidSystem.enableVaporizedWater()) { // Add Sg > 0? i.e. if only water set rv = 0)
const Evaluation& RvwSat = fluidSystem.saturatedVaporizationFactor(fluidState_,
gasPhaseIdx,
pvtRegionIdx);
fluidState_.setRvw(RvwSat);
@ -382,8 +386,8 @@ public:
const auto& Rsw = priVars.makeEvaluation(Indices::waterSwitchIdx, timeIdx);
fluidState_.setRsw(Rsw);
} else {
if (FluidSystem::enableDissolvedGasInWater()) {
const Evaluation& RswSat = FluidSystem::saturatedDissolutionFactor(fluidState_,
if (fluidSystem.enableDissolvedGasInWater()) {
const Evaluation& RswSat = fluidSystem.saturatedDissolutionFactor(fluidState_,
waterPhaseIdx,
pvtRegionIdx);
fluidState_.setRsw(min(RswMax, RswSat));
@ -393,7 +397,8 @@ public:
return SoMax;
}
void updateMobilityAndInvB()
template <class DynamicFluidSystem = FluidSystem>
void updateMobilityAndInvB(const DynamicFluidSystem& fluidSystem = DynamicFluidSystem{})
{
const unsigned pvtRegionIdx = fluidState_.pvtRegionIndex();
@ -407,11 +412,11 @@ public:
}
}
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx))
if (!fluidSystem.phaseIsActive(phaseIdx))
continue;
const auto& b = FluidSystem::inverseFormationVolumeFactor(fluidState_, phaseIdx, pvtRegionIdx);
const auto& b = fluidSystem.inverseFormationVolumeFactor(fluidState_, phaseIdx, pvtRegionIdx);
fluidState_.setInvB(phaseIdx, b);
const auto& mu = FluidSystem::viscosity(fluidState_, phaseIdx, pvtRegionIdx);
const auto& mu = fluidSystem.viscosity(fluidState_, phaseIdx, pvtRegionIdx);
for (int i = 0; i<nmobilities; i++) {
if (enableExtbo && phaseIdx == oilPhaseIdx) {
(*mobilities[i])[phaseIdx] /= asImp_().oilViscosity();
@ -427,56 +432,58 @@ public:
Valgrind::CheckDefined(mobility_);
}
void updatePhaseDensities()
template <class DynamicFluidSystem = FluidSystem>
void updatePhaseDensities(const DynamicFluidSystem& fluidSystem = DynamicFluidSystem{})
{
const unsigned pvtRegionIdx = fluidState_.pvtRegionIndex();
// calculate the phase densities
Evaluation rho;
if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
if (fluidSystem.phaseIsActive(waterPhaseIdx)) {
rho = fluidState_.invB(waterPhaseIdx);
rho *= FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
if (FluidSystem::enableDissolvedGasInWater()) {
rho *= fluidSystem.referenceDensity(waterPhaseIdx, pvtRegionIdx);
if (fluidSystem.enableDissolvedGasInWater()) {
rho +=
fluidState_.invB(waterPhaseIdx) *
fluidState_.Rsw() *
FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
fluidSystem.referenceDensity(gasPhaseIdx, pvtRegionIdx);
}
fluidState_.setDensity(waterPhaseIdx, rho);
}
if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
if (fluidSystem.phaseIsActive(gasPhaseIdx)) {
rho = fluidState_.invB(gasPhaseIdx);
rho *= FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
if (FluidSystem::enableVaporizedOil()) {
rho *= fluidSystem.referenceDensity(gasPhaseIdx, pvtRegionIdx);
if (fluidSystem.enableVaporizedOil()) {
rho +=
fluidState_.invB(gasPhaseIdx) *
fluidState_.Rv() *
FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
fluidSystem.referenceDensity(oilPhaseIdx, pvtRegionIdx);
}
if (FluidSystem::enableVaporizedWater()) {
if (fluidSystem.enableVaporizedWater()) {
rho +=
fluidState_.invB(gasPhaseIdx) *
fluidState_.Rvw() *
FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
fluidSystem.referenceDensity(waterPhaseIdx, pvtRegionIdx);
}
fluidState_.setDensity(gasPhaseIdx, rho);
}
if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
if (fluidSystem.phaseIsActive(oilPhaseIdx)) {
rho = fluidState_.invB(oilPhaseIdx);
rho *= FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
if (FluidSystem::enableDissolvedGas()) {
rho *= fluidSystem.referenceDensity(oilPhaseIdx, pvtRegionIdx);
if (fluidSystem.enableDissolvedGas()) {
rho +=
fluidState_.invB(oilPhaseIdx) *
fluidState_.Rs() *
FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
fluidSystem.referenceDensity(gasPhaseIdx, pvtRegionIdx);
}
fluidState_.setDensity(oilPhaseIdx, rho);
}
}
void updatePorosity(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
template <class DynamicFluidSystem = FluidSystem>
void updatePorosity(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx, const DynamicFluidSystem& fluidSystem = DynamicFluidSystem{})
{
const auto& problem = elemCtx.problem();
const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
@ -493,9 +500,9 @@ public:
if (rockCompressibility > 0.0) {
Scalar rockRefPressure = problem.rockReferencePressure(globalSpaceIdx);
Evaluation x;
if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
if (fluidSystem.phaseIsActive(oilPhaseIdx)) {
x = rockCompressibility*(fluidState_.pressure(oilPhaseIdx) - rockRefPressure);
} else if (FluidSystem::phaseIsActive(waterPhaseIdx)){
} else if (fluidSystem.phaseIsActive(waterPhaseIdx)){
x = rockCompressibility*(fluidState_.pressure(waterPhaseIdx) - rockRefPressure);
} else {
x = rockCompressibility*(fluidState_.pressure(gasPhaseIdx) - rockRefPressure);
@ -520,11 +527,12 @@ public:
}
}
void assertFiniteMembers()
template <class DynamicFluidSystem = FluidSystem>
void assertFiniteMembers(const DynamicFluidSystem& fluidSystem = DynamicFluidSystem{})
{
// some safety checks in debug mode
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx))
if (!fluidSystem.phaseIsActive(phaseIdx))
continue;
assert(isfinite(fluidState_.density(phaseIdx)));
@ -540,7 +548,8 @@ public:
/*!
* \copydoc IntensiveQuantities::update
*/
void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
template <class DynamicFluidSystem = FluidSystem>
void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx, const DynamicFluidSystem& fluidSystem = DynamicFluidSystem{})
{
ParentType::update(elemCtx, dofIdx, timeIdx);
@ -554,19 +563,19 @@ public:
fluidState_.setPvtRegionIndex(pvtRegionIdx);
updateTempSalt(elemCtx, dofIdx, timeIdx);
updateSaturations(elemCtx, dofIdx, timeIdx);
updateRelpermAndPressures(elemCtx, dofIdx, timeIdx);
updateSaturations(elemCtx, dofIdx, timeIdx, fluidSystem);
updateRelpermAndPressures(elemCtx, dofIdx, timeIdx, fluidSystem);
// update extBO parameters
if constexpr (enableExtbo) {
asImp_().zFractionUpdate_(elemCtx, dofIdx, timeIdx);
}
Evaluation SoMax = updateRsRvRsw(elemCtx, dofIdx, timeIdx);
Evaluation SoMax = updateRsRvRsw(elemCtx, dofIdx, timeIdx, fluidSystem);
updateMobilityAndInvB();
updatePhaseDensities();
updatePorosity(elemCtx, dofIdx, timeIdx);
updateMobilityAndInvB(fluidSystem);
updatePhaseDensities(fluidSystem);
updatePorosity(elemCtx, dofIdx, timeIdx, fluidSystem);
rockCompTransMultiplier_ = problem.template rockCompTransMultiplier<Evaluation>(*this, globalSpaceIdx);
@ -582,7 +591,7 @@ public:
typename FluidSystem::template ParameterCache<Evaluation> paramCache;
paramCache.setRegionIndex(pvtRegionIdx);
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
if (fluidSystem.phaseIsActive(fluidSystem.oilPhaseIdx)) {
paramCache.setMaxOilSat(SoMax);
}
paramCache.updateAll(fluidState_);
@ -611,7 +620,7 @@ public:
DispersionIntensiveQuantities::update_(elemCtx, dofIdx, timeIdx);
#ifndef NDEBUG
assertFiniteMembers();
assertFiniteMembers(fluidSystem);
#endif
}

32
opm/models/discretization/common/fvbaseelementcontext.hh
View File

@ -209,8 +209,16 @@ public:
*
* \param timeIdx The index of the solution vector used by the time discretization.
*/
void updatePrimaryIntensiveQuantities(unsigned timeIdx)
{ updateIntensiveQuantities_(timeIdx, numPrimaryDof(timeIdx)); }
template <class FluidSystemT = std::nullptr_t>
void updatePrimaryIntensiveQuantities(unsigned timeIdx, const FluidSystemT& Fsystem = nullptr)
{
if constexpr (std::is_same_v<FluidSystemT, std::nullptr_t>) {
updateIntensiveQuantities_(timeIdx, numPrimaryDof(timeIdx));
}
else {
updateIntensiveQuantities_(timeIdx, numPrimaryDof(timeIdx), Fsystem);
}
}
/*!
* \brief Compute the intensive quantities of a single sub-control volume of the
@ -545,7 +553,8 @@ protected:
*
* This method considers the intensive quantities cache.
*/
void updateIntensiveQuantities_(unsigned timeIdx, size_t numDof)
template <class FluidSystemT = std::nullptr_t>
void updateIntensiveQuantities_(unsigned timeIdx, size_t numDof, const FluidSystemT& Fsystem = nullptr)
{
// update the intensive quantities for the whole history
const SolutionVector& globalSol = model().solution(timeIdx);
@ -564,7 +573,13 @@ protected:
dofVars_[dofIdx].intensiveQuantities[timeIdx] = *cachedIntQuants;
}
else {
updateSingleIntQuants_(dofSol, dofIdx, timeIdx);
if constexpr (std::is_same_v<FluidSystemT, std::nullptr_t>){
updateSingleIntQuants_(dofSol, dofIdx, timeIdx);
}
else
{
updateSingleIntQuants_(dofSol, dofIdx, timeIdx, Fsystem);
}
model().updateCachedIntensiveQuantities(dofVars_[dofIdx].intensiveQuantities[timeIdx],
globalIdx,
timeIdx);
@ -572,7 +587,8 @@ protected:
}
}
void updateSingleIntQuants_(const PrimaryVariables& priVars, unsigned dofIdx, unsigned timeIdx)
template <class FluidSystemT = std::nullptr_t>
void updateSingleIntQuants_(const PrimaryVariables& priVars, unsigned dofIdx, unsigned timeIdx, const FluidSystemT& Fsystem = nullptr)
{
#ifndef NDEBUG
if (enableStorageCache_ && timeIdx != 0 && problem().recycleFirstIterationStorage())
@ -581,7 +597,11 @@ protected:
#endif
dofVars_[dofIdx].priVars[timeIdx] = &priVars;
dofVars_[dofIdx].intensiveQuantities[timeIdx].update(/*context=*/asImp_(), dofIdx, timeIdx);
if constexpr (std::is_same_v<FluidSystemT, std::nullptr_t>) {
dofVars_[dofIdx].intensiveQuantities[timeIdx].update(/*context=*/asImp_(), dofIdx, timeIdx);
} else {
dofVars_[dofIdx].intensiveQuantities[timeIdx].update(/*context=*/asImp_(), dofIdx, timeIdx, Fsystem);
}
}
IntensiveQuantities intensiveQuantitiesStashed_;

3
opm/models/flash/flashintensivequantities.hh
View File

@ -100,7 +100,8 @@ public:
/*!
* \copydoc IntensiveQuantities::update
*/
void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
template <class DynamicFluidSystem = FluidSystem>
void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx, [[maybe_unused]] const DynamicFluidSystem& fluidSystem = DynamicFluidSystem{})
{
ParentType::update(elemCtx, dofIdx, timeIdx);
EnergyIntensiveQuantities::updateTemperatures_(fluidState_, elemCtx, dofIdx, timeIdx);

3
opm/models/immiscible/immiscibleintensivequantities.hh
View File

@ -90,7 +90,8 @@ public:
/*!
* \copydoc IntensiveQuantities::update
*/
void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
template <class DynamicFluidSystem = FluidSystem>
void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx, const DynamicFluidSystem& fluidSystem = DynamicFluidSystem{})
{
ParentType::update(elemCtx, dofIdx, timeIdx);
EnergyIntensiveQuantities::updateTemperatures_(fluidState_, elemCtx, dofIdx, timeIdx);

4
opm/models/ncp/ncpintensivequantities.hh
View File

@ -97,9 +97,11 @@ public:
/*!
* \brief IntensiveQuantities::update
*/
template <class DynamicFluidSystem = FluidSystem>
void update(const ElementContext& elemCtx,
unsigned dofIdx,
unsigned timeIdx)
unsigned timeIdx,
const DynamicFluidSystem& fluidSystem = DynamicFluidSystem{})
{
ParentType::update(elemCtx, dofIdx, timeIdx);
ParentType::checkDefined();

3
opm/models/ptflash/flashintensivequantities.hh
View File

@ -105,7 +105,8 @@ public:
/*!
* \copydoc IntensiveQuantities::update
*/
void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
template <class DynamicFluidSystem = FluidSystem>
void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx, [[maybe_unused]] const DynamicFluidSystem& fluidSystem = DynamicFluidSystem{})
{
ParentType::update(elemCtx, dofIdx, timeIdx);
EnergyIntensiveQuantities::updateTemperatures_(fluidState_, elemCtx, dofIdx, timeIdx);

3
opm/models/pvs/pvsintensivequantities.hh
View File

@ -105,7 +105,8 @@ public:
/*!
* \copydoc ImmiscibleIntensiveQuantities::update
*/
void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
template <class DynamicFluidSystem = FluidSystem>
void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx, [[maybe_unused]] const DynamicFluidSystem& fluidSystem = DynamicFluidSystem{})
{
ParentType::update(elemCtx, dofIdx, timeIdx);
EnergyIntensiveQuantities::updateTemperatures_(fluidState_, elemCtx, dofIdx, timeIdx);

3
opm/models/richards/richardsintensivequantities.hh
View File

@ -84,7 +84,8 @@ public:
/*!
* \copydoc IntensiveQuantities::update
*/
void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
template <class DynamicFluidSystem = FluidSystem>
void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx, const DynamicFluidSystem& fluidSystem = DynamicFluidSystem{})
{
ParentType::update(elemCtx, dofIdx, timeIdx);

10
opm/simulators/flow/FIBlackoilModel.hpp
View File

@ -59,6 +59,7 @@ class FIBlackOilModel : public BlackOilModel<TypeTag>
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using ThreadManager = GetPropType<TypeTag, Properties::ThreadManager>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using LocalFluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using Element = typename GridView::template Codim<0>::Entity;
using ElementIterator = typename GridView::template Codim<0>::Iterator;
enum {
@ -80,6 +81,10 @@ public:
void invalidateAndUpdateIntensiveQuantities(unsigned timeIdx) const
{
using DynamicFluidSystem = std::remove_reference_t<decltype(LocalFluidSystem::getNonStatic())>;
DynamicFluidSystem* fluidSystemInstance = &LocalFluidSystem::getNonStatic();
this->invalidateIntensiveQuantitiesCache(timeIdx);
OPM_BEGIN_PARALLEL_TRY_CATCH();
if constexpr (gridIsUnchanging) {
@ -90,7 +95,8 @@ public:
ElementContext elemCtx(this->simulator_);
for (const auto& elem : chunk) {
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(timeIdx);
elemCtx.updatePrimaryIntensiveQuantities(timeIdx, *fluidSystemInstance);
}
}
} else {
@ -98,7 +104,7 @@ public:
ElementContext elemCtx(this->simulator_);
for (const auto& elem : elements(this->gridView_)) {
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(timeIdx);
elemCtx.updatePrimaryIntensiveQuantities(timeIdx, *fluidSystemInstance);
}
}
OPM_END_PARALLEL_TRY_CATCH("InvalideAndUpdateIntensiveQuantities: state error", this->simulator_.vanguard().grid().comm());