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 pull request #4372 from totto82/fix_gp

Browse Source 
output gas pressure for gas-water case
This commit is contained in:
Bård Skaflestad 2023-01-16 11:11:39 +01:00 committed by GitHub
commit 2691e6328f
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
3 changed files with 213 additions and 261 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
ebos/eclgenericoutputblackoilmodule.cc
View File

@ -682,7 +682,7 @@ assignToSolution(data::Solution& sol)
{"PORV_RC", UnitSystem::measure::identity, data::TargetType::RESTART_AUXILIARY, rockCompPorvMultiplier_}, {"PORV_RC", UnitSystem::measure::identity, data::TargetType::RESTART_AUXILIARY, rockCompPorvMultiplier_},
{"PPCW", UnitSystem::measure::pressure, data::TargetType::RESTART_SOLUTION, ppcw_}, {"PPCW", UnitSystem::measure::pressure, data::TargetType::RESTART_SOLUTION, ppcw_},
{"PRESROCC", UnitSystem::measure::pressure, data::TargetType::RESTART_SOLUTION, minimumOilPressure_}, {"PRESROCC", UnitSystem::measure::pressure, data::TargetType::RESTART_SOLUTION, minimumOilPressure_},
{"PRESSURE", UnitSystem::measure::pressure, data::TargetType::RESTART_SOLUTION, oilPressure_}, {"PRESSURE", UnitSystem::measure::pressure, data::TargetType::RESTART_SOLUTION, fluidPressure_},
{"PCOW", UnitSystem::measure::pressure, data::TargetType::RESTART_SOLUTION, pcow_}, {"PCOW", UnitSystem::measure::pressure, data::TargetType::RESTART_SOLUTION, pcow_},
{"PCOG", UnitSystem::measure::pressure, data::TargetType::RESTART_SOLUTION, pcog_}, {"PCOG", UnitSystem::measure::pressure, data::TargetType::RESTART_SOLUTION, pcog_},
{"PRES_OVB", UnitSystem::measure::pressure, data::TargetType::RESTART_SOLUTION, overburdenPressure_}, {"PRES_OVB", UnitSystem::measure::pressure, data::TargetType::RESTART_SOLUTION, overburdenPressure_},
@ -788,8 +788,8 @@ setRestart(const data::Solution& sol,
assert(!saturation_[oilPhaseIdx].empty()); assert(!saturation_[oilPhaseIdx].empty());
saturation_[oilPhaseIdx][elemIdx] = so; saturation_[oilPhaseIdx][elemIdx] = so;
if (!oilPressure_.empty() && sol.has("PRESSURE")) if (!fluidPressure_.empty() && sol.has("PRESSURE"))
oilPressure_[elemIdx] = sol.data("PRESSURE")[globalDofIndex]; fluidPressure_[elemIdx] = sol.data("PRESSURE")[globalDofIndex];
if (!temperature_.empty() && sol.has("TEMP")) if (!temperature_.empty() && sol.has("TEMP"))
temperature_[elemIdx] = sol.data("TEMP")[globalDofIndex]; temperature_[elemIdx] = sol.data("TEMP")[globalDofIndex];
if (!rs_.empty() && sol.has("RS")) if (!rs_.empty() && sol.has("RS"))
@ -969,7 +969,7 @@ doAllocBuffers(unsigned bufferSize,
saturation_[phaseIdx].resize(bufferSize, 0.0); saturation_[phaseIdx].resize(bufferSize, 0.0);
} }
// and oil pressure // and oil pressure
oilPressure_.resize(bufferSize, 0.0); fluidPressure_.resize(bufferSize, 0.0);
rstKeywords["PRES"] = 0; rstKeywords["PRES"] = 0;
rstKeywords["PRESSURE"] = 0; rstKeywords["PRESSURE"] = 0;

6
ebos/eclgenericoutputblackoilmodule.hh
View File

@ -48,11 +48,11 @@ template<class FluidSystem, class Scalar>
class EclGenericOutputBlackoilModule { class EclGenericOutputBlackoilModule {
public: public:
Scalar* getPRESSURE_ptr(void) { Scalar* getPRESSURE_ptr(void) {
return (this->oilPressure_.data()) ; return (this->fluidPressure_.data()) ;
}; };
int getPRESSURE_size( void ) { int getPRESSURE_size( void ) {
return (this->oilPressure_.size()) ; return (this->fluidPressure_.size()) ;
}; };
// write cumulative production and injection reports to output // write cumulative production and injection reports to output
@ -411,7 +411,7 @@ protected:
ScalarBuffer pressureTimesPoreVolume_; ScalarBuffer pressureTimesPoreVolume_;
ScalarBuffer pressureTimesHydrocarbonVolume_; ScalarBuffer pressureTimesHydrocarbonVolume_;
ScalarBuffer dynamicPoreVolume_; ScalarBuffer dynamicPoreVolume_;
ScalarBuffer oilPressure_; ScalarBuffer fluidPressure_;
ScalarBuffer temperature_; ScalarBuffer temperature_;
ScalarBuffer rs_; ScalarBuffer rs_;
ScalarBuffer rsw_; ScalarBuffer rsw_;

350
ebos/ecloutputblackoilmodule.hh
View File

@ -32,8 +32,8 @@
#include <opm/models/blackoil/blackoilproperties.hh> #include <opm/models/blackoil/blackoilproperties.hh>
#include <opm/models/utils/propertysystem.hh>
#include <opm/models/utils/parametersystem.hh> #include <opm/models/utils/parametersystem.hh>
#include <opm/models/utils/propertysystem.hh>
#include <opm/material/common/Valgrind.hpp> #include <opm/material/common/Valgrind.hpp>
@ -59,12 +59,15 @@
#include <type_traits> #include <type_traits>
#include <utility> #include <utility>
namespace Opm::Properties { namespace Opm::Properties
{
// create new type tag for the Ecl-output // create new type tag for the Ecl-output
namespace TTag { namespace TTag
struct EclOutputBlackOil {}; {
} struct EclOutputBlackOil {
};
} // namespace TTag
template <class TypeTag, class MyTypeTag> template <class TypeTag, class MyTypeTag>
struct ForceDisableFluidInPlaceOutput { struct ForceDisableFluidInPlaceOutput {
@ -89,7 +92,8 @@ struct ForceDisableResvFluidInPlaceOutput<TypeTag, TTag::EclOutputBlackOil> {
} // namespace Opm::Properties } // namespace Opm::Properties
namespace Opm { namespace Opm
{
// forward declaration // forward declaration
template <class TypeTag> template <class TypeTag>
@ -151,12 +155,10 @@ public:
} }
for (const auto& node : this->simulator_.vanguard().summaryConfig()) { for (const auto& node : this->simulator_.vanguard().summaryConfig()) {
if ((node.category() == SummaryConfigNode::Category::Block) && if ((node.category() == SummaryConfigNode::Category::Block)
collectToIORank.isCartIdxOnThisRank(node.number() - 1)) && collectToIORank.isCartIdxOnThisRank(node.number() - 1)) {
{
this->blockData_.emplace(std::piecewise_construct, this->blockData_.emplace(std::piecewise_construct,
std::forward_as_tuple(node.keyword(), std::forward_as_tuple(node.keyword(), node.number()),
node.number()),
std::forward_as_tuple(0.0)); std::forward_as_tuple(0.0));
} }
} }
@ -167,11 +169,9 @@ public:
} }
} }
this->forceDisableFipOutput_ = this->forceDisableFipOutput_ = EWOMS_GET_PARAM(TypeTag, bool, ForceDisableFluidInPlaceOutput);
EWOMS_GET_PARAM(TypeTag, bool, ForceDisableFluidInPlaceOutput);
this->forceDisableFipresvOutput_ = this->forceDisableFipresvOutput_ = EWOMS_GET_PARAM(TypeTag, bool, ForceDisableResvFluidInPlaceOutput);
EWOMS_GET_PARAM(TypeTag, bool, ForceDisableResvFluidInPlaceOutput);
} }
/*! /*!
@ -179,9 +179,15 @@ public:
*/ */
static void registerParameters() static void registerParameters()
{ {
EWOMS_REGISTER_PARAM(TypeTag, bool, ForceDisableFluidInPlaceOutput, EWOMS_REGISTER_PARAM(
TypeTag,
bool,
ForceDisableFluidInPlaceOutput,
"Do not print fluid-in-place values after each report step even if requested by the deck."); "Do not print fluid-in-place values after each report step even if requested by the deck.");
EWOMS_REGISTER_PARAM(TypeTag, bool, ForceDisableResvFluidInPlaceOutput, EWOMS_REGISTER_PARAM(
TypeTag,
bool,
ForceDisableResvFluidInPlaceOutput,
"Do not print reservoir volumes values after each report step even if requested by the deck."); "Do not print reservoir volumes values after each report step even if requested by the deck.");
} }
@ -189,7 +195,8 @@ public:
* \brief Allocate memory for the scalar fields we would like to * \brief Allocate memory for the scalar fields we would like to
* write to ECL output files * write to ECL output files
*/ */
void allocBuffers(unsigned bufferSize, unsigned reportStepNum, const bool substep, const bool log, const bool isRestart) void
allocBuffers(unsigned bufferSize, unsigned reportStepNum, const bool substep, const bool log, const bool isRestart)
{ {
if (!std::is_same<Discretization, EcfvDiscretization<TypeTag>>::value) if (!std::is_same<Discretization, EcfvDiscretization<TypeTag>>::value)
return; return;
@ -230,18 +237,18 @@ public:
Valgrind::CheckDefined(this->saturation_[phaseIdx][globalDofIdx]); Valgrind::CheckDefined(this->saturation_[phaseIdx][globalDofIdx]);
} }
if (!this->oilPressure_.empty()) { if (!this->fluidPressure_.empty()) {
if (FluidSystem::phaseIsActive(oilPhaseIdx)) { if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
this->oilPressure_[globalDofIdx] = getValue(fs.pressure(oilPhaseIdx)); // Output oil pressure as default
this->fluidPressure_[globalDofIdx] = getValue(fs.pressure(oilPhaseIdx));
} else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
// Output gas if oil is not present
this->fluidPressure_[globalDofIdx] = getValue(fs.pressure(gasPhaseIdx));
} else { } else {
// put pressure in oil pressure for output // Output water if neither oil nor gas is present
if (FluidSystem::phaseIsActive(waterPhaseIdx)) { this->fluidPressure_[globalDofIdx] = getValue(fs.pressure(waterPhaseIdx));
this->oilPressure_[globalDofIdx] = getValue(fs.pressure(waterPhaseIdx));
} else {
this->oilPressure_[globalDofIdx] = getValue(fs.pressure(gasPhaseIdx));
} }
} Valgrind::CheckDefined(this->fluidPressure_[globalDofIdx]);
Valgrind::CheckDefined(this->oilPressure_[globalDofIdx]);
} }
if (!this->temperature_.empty()) { if (!this->temperature_.empty()) {
@ -250,35 +257,34 @@ public:
} }
if (!this->gasDissolutionFactor_.empty()) { if (!this->gasDissolutionFactor_.empty()) {
Scalar SoMax = elemCtx.problem().maxOilSaturation(globalDofIdx); Scalar SoMax = elemCtx.problem().maxOilSaturation(globalDofIdx);
this->gasDissolutionFactor_[globalDofIdx] = this->gasDissolutionFactor_[globalDofIdx]
FluidSystem::template saturatedDissolutionFactor<FluidState, Scalar>(fs, oilPhaseIdx, pvtRegionIdx, SoMax); = FluidSystem::template saturatedDissolutionFactor<FluidState, Scalar>(
fs, oilPhaseIdx, pvtRegionIdx, SoMax);
Valgrind::CheckDefined(this->gasDissolutionFactor_[globalDofIdx]); Valgrind::CheckDefined(this->gasDissolutionFactor_[globalDofIdx]);
} }
if (!this->oilVaporizationFactor_.empty()) { if (!this->oilVaporizationFactor_.empty()) {
Scalar SoMax = elemCtx.problem().maxOilSaturation(globalDofIdx); Scalar SoMax = elemCtx.problem().maxOilSaturation(globalDofIdx);
this->oilVaporizationFactor_[globalDofIdx] = this->oilVaporizationFactor_[globalDofIdx]
FluidSystem::template saturatedDissolutionFactor<FluidState, Scalar>(fs, gasPhaseIdx, pvtRegionIdx, SoMax); = FluidSystem::template saturatedDissolutionFactor<FluidState, Scalar>(
fs, gasPhaseIdx, pvtRegionIdx, SoMax);
Valgrind::CheckDefined(this->oilVaporizationFactor_[globalDofIdx]); Valgrind::CheckDefined(this->oilVaporizationFactor_[globalDofIdx]);
} }
if (!this->gasFormationVolumeFactor_.empty()) { if (!this->gasFormationVolumeFactor_.empty()) {
this->gasFormationVolumeFactor_[globalDofIdx] = this->gasFormationVolumeFactor_[globalDofIdx] = 1.0
1.0/FluidSystem::template inverseFormationVolumeFactor<FluidState, Scalar>(fs, gasPhaseIdx, pvtRegionIdx); / FluidSystem::template inverseFormationVolumeFactor<FluidState, Scalar>(
fs, gasPhaseIdx, pvtRegionIdx);
Valgrind::CheckDefined(this->gasFormationVolumeFactor_[globalDofIdx]); Valgrind::CheckDefined(this->gasFormationVolumeFactor_[globalDofIdx]);
} }
if (!this->saturatedOilFormationVolumeFactor_.empty()) { if (!this->saturatedOilFormationVolumeFactor_.empty()) {
this->saturatedOilFormationVolumeFactor_[globalDofIdx] = this->saturatedOilFormationVolumeFactor_[globalDofIdx] = 1.0
1.0/FluidSystem::template saturatedInverseFormationVolumeFactor<FluidState, Scalar>(fs, oilPhaseIdx, pvtRegionIdx); / FluidSystem::template saturatedInverseFormationVolumeFactor<FluidState, Scalar>(
fs, oilPhaseIdx, pvtRegionIdx);
Valgrind::CheckDefined(this->saturatedOilFormationVolumeFactor_[globalDofIdx]); Valgrind::CheckDefined(this->saturatedOilFormationVolumeFactor_[globalDofIdx]);
} }
if (!this->oilSaturationPressure_.empty()) { if (!this->oilSaturationPressure_.empty()) {
this->oilSaturationPressure_[globalDofIdx] = this->oilSaturationPressure_[globalDofIdx]
FluidSystem::template saturationPressure<FluidState, Scalar>(fs, oilPhaseIdx, pvtRegionIdx); = FluidSystem::template saturationPressure<FluidState, Scalar>(fs, oilPhaseIdx, pvtRegionIdx);
Valgrind::CheckDefined(this->oilSaturationPressure_[globalDofIdx]); Valgrind::CheckDefined(this->oilSaturationPressure_[globalDofIdx]);
} }
if (!this->rs_.empty()) { if (!this->rs_.empty()) {
@ -341,7 +347,8 @@ public:
if (this->relativePermeability_[phaseIdx].empty()) if (this->relativePermeability_[phaseIdx].empty())
continue; continue;
this->relativePermeability_[phaseIdx][globalDofIdx] = getValue(intQuants.relativePermeability(phaseIdx)); this->relativePermeability_[phaseIdx][globalDofIdx]
= getValue(intQuants.relativePermeability(phaseIdx));
Valgrind::CheckDefined(this->relativePermeability_[phaseIdx][globalDofIdx]); Valgrind::CheckDefined(this->relativePermeability_[phaseIdx][globalDofIdx]);
} }
@ -388,10 +395,14 @@ public:
if (!this->mFracCo2_.empty()) { if (!this->mFracCo2_.empty()) {
const Scalar stdVolOil = getValue(fs.saturation(oilPhaseIdx)) * getValue(fs.invB(oilPhaseIdx)) const Scalar stdVolOil = getValue(fs.saturation(oilPhaseIdx)) * getValue(fs.invB(oilPhaseIdx))
+ getValue(fs.saturation(gasPhaseIdx)) * getValue(fs.invB(gasPhaseIdx)) * getValue(fs.Rv()); + getValue(fs.saturation(gasPhaseIdx)) * getValue(fs.invB(gasPhaseIdx)) * getValue(fs.Rv());
const Scalar stdVolGas = getValue(fs.saturation(gasPhaseIdx))*getValue(fs.invB(gasPhaseIdx))*(1.0-intQuants.yVolume().value()) const Scalar stdVolGas = getValue(fs.saturation(gasPhaseIdx)) * getValue(fs.invB(gasPhaseIdx))
+ getValue(fs.saturation(oilPhaseIdx))*getValue(fs.invB(oilPhaseIdx))*getValue(fs.Rs())*(1.0-intQuants.xVolume().value()); * (1.0 - intQuants.yVolume().value())
const Scalar stdVolCo2 = getValue(fs.saturation(gasPhaseIdx))*getValue(fs.invB(gasPhaseIdx))*intQuants.yVolume().value() + getValue(fs.saturation(oilPhaseIdx)) * getValue(fs.invB(oilPhaseIdx)) * getValue(fs.Rs())
+ getValue(fs.saturation(oilPhaseIdx))*getValue(fs.invB(oilPhaseIdx))*getValue(fs.Rs())*intQuants.xVolume().value(); * (1.0 - intQuants.xVolume().value());
const Scalar stdVolCo2 = getValue(fs.saturation(gasPhaseIdx)) * getValue(fs.invB(gasPhaseIdx))
* intQuants.yVolume().value()
+ getValue(fs.saturation(oilPhaseIdx)) * getValue(fs.invB(oilPhaseIdx)) * getValue(fs.Rs())
* intQuants.xVolume().value();
const Scalar rhoO = FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx); const Scalar rhoO = FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
const Scalar rhoG = FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx); const Scalar rhoG = FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
const Scalar rhoCO2 = intQuants.zRefDensity(); const Scalar rhoCO2 = intQuants.zRefDensity();
@ -410,7 +421,9 @@ public:
} }
if (!this->cUrea_.empty()) { if (!this->cUrea_.empty()) {
this->cUrea_[globalDofIdx] = 10 * intQuants.ureaConcentration().value(); //Reescaling back the urea concentration (see WellInterface_impl.hpp) this->cUrea_[globalDofIdx] = 10
* intQuants.ureaConcentration()
.value(); // Reescaling back the urea concentration (see WellInterface_impl.hpp)
} }
if (!this->cBiofilm_.empty()) { if (!this->cBiofilm_.empty()) {
@ -423,60 +436,55 @@ public:
if (!this->bubblePointPressure_.empty()) { if (!this->bubblePointPressure_.empty()) {
try { try {
this->bubblePointPressure_[globalDofIdx] = getValue(FluidSystem::bubblePointPressure(fs, intQuants.pvtRegionIndex())); this->bubblePointPressure_[globalDofIdx]
} = getValue(FluidSystem::bubblePointPressure(fs, intQuants.pvtRegionIndex()));
catch (const NumericalProblem&) { } catch (const NumericalProblem&) {
const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx); const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx);
this->failedCellsPb_.push_back(cartesianIdx); this->failedCellsPb_.push_back(cartesianIdx);
} }
} }
if (!this->dewPointPressure_.empty()) { if (!this->dewPointPressure_.empty()) {
try { try {
this->dewPointPressure_[globalDofIdx] = getValue(FluidSystem::dewPointPressure(fs, intQuants.pvtRegionIndex())); this->dewPointPressure_[globalDofIdx]
} = getValue(FluidSystem::dewPointPressure(fs, intQuants.pvtRegionIndex()));
catch (const NumericalProblem&) { } catch (const NumericalProblem&) {
const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx); const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx);
this->failedCellsPd_.push_back(cartesianIdx); this->failedCellsPd_.push_back(cartesianIdx);
} }
} }
if (!this->soMax_.empty()) if (!this->soMax_.empty())
this->soMax_[globalDofIdx] = this->soMax_[globalDofIdx]
std::max(getValue(fs.saturation(oilPhaseIdx)), = std::max(getValue(fs.saturation(oilPhaseIdx)), problem.maxOilSaturation(globalDofIdx));
problem.maxOilSaturation(globalDofIdx));
if (!this->swMax_.empty()) if (!this->swMax_.empty())
this->swMax_[globalDofIdx] = this->swMax_[globalDofIdx]
std::max(getValue(fs.saturation(waterPhaseIdx)), = std::max(getValue(fs.saturation(waterPhaseIdx)), problem.maxWaterSaturation(globalDofIdx));
problem.maxWaterSaturation(globalDofIdx));
if (!this->minimumOilPressure_.empty()) if (!this->minimumOilPressure_.empty())
this->minimumOilPressure_[globalDofIdx] = this->minimumOilPressure_[globalDofIdx]
std::min(getValue(fs.pressure(oilPhaseIdx)), = std::min(getValue(fs.pressure(oilPhaseIdx)), problem.minOilPressure(globalDofIdx));
problem.minOilPressure(globalDofIdx));
if (!this->overburdenPressure_.empty()) if (!this->overburdenPressure_.empty())
this->overburdenPressure_[globalDofIdx] = problem.overburdenPressure(globalDofIdx); this->overburdenPressure_[globalDofIdx] = problem.overburdenPressure(globalDofIdx);
if (!this->rockCompPorvMultiplier_.empty()) if (!this->rockCompPorvMultiplier_.empty())
this->rockCompPorvMultiplier_[globalDofIdx] = problem.template rockCompPoroMultiplier<Scalar>(intQuants, globalDofIdx); this->rockCompPorvMultiplier_[globalDofIdx]
= problem.template rockCompPoroMultiplier<Scalar>(intQuants, globalDofIdx);
if (!this->rockCompTransMultiplier_.empty()) if (!this->rockCompTransMultiplier_.empty())
this->rockCompTransMultiplier_[globalDofIdx] = problem.template rockCompTransMultiplier<Scalar>(intQuants, globalDofIdx); this->rockCompTransMultiplier_[globalDofIdx]
= problem.template rockCompTransMultiplier<Scalar>(intQuants, globalDofIdx);
const auto& matLawManager = problem.materialLawManager(); const auto& matLawManager = problem.materialLawManager();
if (matLawManager->enableHysteresis()) { if (matLawManager->enableHysteresis()) {
if (!this->pcSwMdcOw_.empty() && !this->krnSwMdcOw_.empty()) { if (!this->pcSwMdcOw_.empty() && !this->krnSwMdcOw_.empty()) {
matLawManager->oilWaterHysteresisParams( matLawManager->oilWaterHysteresisParams(
this->pcSwMdcOw_[globalDofIdx], this->pcSwMdcOw_[globalDofIdx], this->krnSwMdcOw_[globalDofIdx], globalDofIdx);
this->krnSwMdcOw_[globalDofIdx],
globalDofIdx);
} }
if (!this->pcSwMdcGo_.empty() && !this->krnSwMdcGo_.empty()) { if (!this->pcSwMdcGo_.empty() && !this->krnSwMdcGo_.empty()) {
matLawManager->gasOilHysteresisParams( matLawManager->gasOilHysteresisParams(
this->pcSwMdcGo_[globalDofIdx], this->pcSwMdcGo_[globalDofIdx], this->krnSwMdcGo_[globalDofIdx], globalDofIdx);
this->krnSwMdcGo_[globalDofIdx],
globalDofIdx);
} }
} }
@ -491,7 +499,8 @@ public:
// rs and rv with the values computed in the initially. // rs and rv with the values computed in the initially.
// Volume factors, densities and viscosities need to be recalculated with the updated rs and rv values. // Volume factors, densities and viscosities need to be recalculated with the updated rs and rv values.
// This can be removed when ebos has 100% controll over output // This can be removed when ebos has 100% controll over output
if (elemCtx.simulator().episodeIndex() < 0 && FluidSystem::phaseIsActive(oilPhaseIdx) && FluidSystem::phaseIsActive(gasPhaseIdx)) { if (elemCtx.simulator().episodeIndex() < 0 && FluidSystem::phaseIsActive(oilPhaseIdx)
&& FluidSystem::phaseIsActive(gasPhaseIdx)) {
const auto& fsInitial = problem.initialFluidState(globalDofIdx); const auto& fsInitial = problem.initialFluidState(globalDofIdx);
@ -510,30 +519,24 @@ public:
// re-compute the volume factors, viscosities and densities if asked for // re-compute the volume factors, viscosities and densities if asked for
if (!this->density_[oilPhaseIdx].empty()) if (!this->density_[oilPhaseIdx].empty())
this->density_[oilPhaseIdx][globalDofIdx] = FluidSystem::density(fsInitial, this->density_[oilPhaseIdx][globalDofIdx]
oilPhaseIdx, = FluidSystem::density(fsInitial, oilPhaseIdx, intQuants.pvtRegionIndex());
intQuants.pvtRegionIndex());
if (!this->density_[gasPhaseIdx].empty()) if (!this->density_[gasPhaseIdx].empty())
this->density_[gasPhaseIdx][globalDofIdx] = FluidSystem::density(fsInitial, this->density_[gasPhaseIdx][globalDofIdx]
gasPhaseIdx, = FluidSystem::density(fsInitial, gasPhaseIdx, intQuants.pvtRegionIndex());
intQuants.pvtRegionIndex());
if (!this->invB_[oilPhaseIdx].empty()) if (!this->invB_[oilPhaseIdx].empty())
this->invB_[oilPhaseIdx][globalDofIdx] = FluidSystem::inverseFormationVolumeFactor(fsInitial, this->invB_[oilPhaseIdx][globalDofIdx]
oilPhaseIdx, = FluidSystem::inverseFormationVolumeFactor(fsInitial, oilPhaseIdx, intQuants.pvtRegionIndex());
intQuants.pvtRegionIndex());
if (!this->invB_[gasPhaseIdx].empty()) if (!this->invB_[gasPhaseIdx].empty())
this->invB_[gasPhaseIdx][globalDofIdx] = FluidSystem::inverseFormationVolumeFactor(fsInitial, this->invB_[gasPhaseIdx][globalDofIdx]
gasPhaseIdx, = FluidSystem::inverseFormationVolumeFactor(fsInitial, gasPhaseIdx, intQuants.pvtRegionIndex());
intQuants.pvtRegionIndex());
if (!this->viscosity_[oilPhaseIdx].empty()) if (!this->viscosity_[oilPhaseIdx].empty())
this->viscosity_[oilPhaseIdx][globalDofIdx] = FluidSystem::viscosity(fsInitial, this->viscosity_[oilPhaseIdx][globalDofIdx]
oilPhaseIdx, = FluidSystem::viscosity(fsInitial, oilPhaseIdx, intQuants.pvtRegionIndex());
intQuants.pvtRegionIndex());
if (!this->viscosity_[gasPhaseIdx].empty()) if (!this->viscosity_[gasPhaseIdx].empty())
this->viscosity_[gasPhaseIdx][globalDofIdx] = FluidSystem::viscosity(fsInitial, this->viscosity_[gasPhaseIdx][globalDofIdx]
gasPhaseIdx, = FluidSystem::viscosity(fsInitial, gasPhaseIdx, intQuants.pvtRegionIndex());
intQuants.pvtRegionIndex());
} }
// Add fluid in Place values // Add fluid in Place values
@ -561,16 +564,14 @@ public:
val.second = getValue(fs.pressure(gasPhaseIdx)); val.second = getValue(fs.pressure(gasPhaseIdx));
else if (FluidSystem::phaseIsActive(waterPhaseIdx)) else if (FluidSystem::phaseIsActive(waterPhaseIdx))
val.second = getValue(fs.pressure(waterPhaseIdx)); val.second = getValue(fs.pressure(waterPhaseIdx));
} } else if ((key.first == "BTCNFHEA") || (key.first == "BTEMP")) {
else if ((key.first == "BTCNFHEA") || (key.first == "BTEMP")){
if (FluidSystem::phaseIsActive(oilPhaseIdx)) if (FluidSystem::phaseIsActive(oilPhaseIdx))
val.second = getValue(fs.temperature(oilPhaseIdx)); val.second = getValue(fs.temperature(oilPhaseIdx));
else if (FluidSystem::phaseIsActive(gasPhaseIdx)) else if (FluidSystem::phaseIsActive(gasPhaseIdx))
val.second = getValue(fs.temperature(gasPhaseIdx)); val.second = getValue(fs.temperature(gasPhaseIdx));
else if (FluidSystem::phaseIsActive(waterPhaseIdx)) else if (FluidSystem::phaseIsActive(waterPhaseIdx))
val.second = getValue(fs.temperature(waterPhaseIdx)); val.second = getValue(fs.temperature(waterPhaseIdx));
} } else if (key.first == "BWKR" || key.first == "BKRW")
else if (key.first == "BWKR" || key.first == "BKRW")
val.second = getValue(intQuants.relativePermeability(waterPhaseIdx)); val.second = getValue(intQuants.relativePermeability(waterPhaseIdx));
else if (key.first == "BGKR" || key.first == "BKRG") else if (key.first == "BGKR" || key.first == "BKRG")
val.second = getValue(intQuants.relativePermeability(gasPhaseIdx)); val.second = getValue(intQuants.relativePermeability(gasPhaseIdx));
@ -578,15 +579,15 @@ public:
val.second = getValue(intQuants.relativePermeability(oilPhaseIdx)); val.second = getValue(intQuants.relativePermeability(oilPhaseIdx));
else if (key.first == "BKROG") { else if (key.first == "BKROG") {
const auto& materialParams = problem.materialLawParams(elemCtx, dofIdx, /* timeIdx = */ 0); const auto& materialParams = problem.materialLawParams(elemCtx, dofIdx, /* timeIdx = */ 0);
const auto krog = MaterialLaw::template relpermOilInOilGasSystem<Evaluation>(materialParams, fs); const auto krog
= MaterialLaw::template relpermOilInOilGasSystem<Evaluation>(materialParams, fs);
val.second = getValue(krog); val.second = getValue(krog);
} } else if (key.first == "BKROW") {
else if (key.first == "BKROW") {
const auto& materialParams = problem.materialLawParams(elemCtx, dofIdx, /* timeIdx = */ 0); const auto& materialParams = problem.materialLawParams(elemCtx, dofIdx, /* timeIdx = */ 0);
const auto krow = MaterialLaw::template relpermOilInOilWaterSystem<Evaluation>(materialParams, fs); const auto krow
= MaterialLaw::template relpermOilInOilWaterSystem<Evaluation>(materialParams, fs);
val.second = getValue(krow); val.second = getValue(krow);
} } else if (key.first == "BWPC")
else if (key.first == "BWPC")
val.second = getValue(fs.pressure(oilPhaseIdx)) - getValue(fs.pressure(waterPhaseIdx)); val.second = getValue(fs.pressure(oilPhaseIdx)) - getValue(fs.pressure(waterPhaseIdx));
else if (key.first == "BGPC") else if (key.first == "BGPC")
val.second = getValue(fs.pressure(gasPhaseIdx)) - getValue(fs.pressure(oilPhaseIdx)); val.second = getValue(fs.pressure(gasPhaseIdx)) - getValue(fs.pressure(oilPhaseIdx));
@ -600,76 +601,56 @@ public:
val.second = getValue(fs.viscosity(gasPhaseIdx)); val.second = getValue(fs.viscosity(gasPhaseIdx));
else if (key.first == "BVOIL" || key.first == "BOVIS") else if (key.first == "BVOIL" || key.first == "BOVIS")
val.second = getValue(fs.viscosity(oilPhaseIdx)); val.second = getValue(fs.viscosity(oilPhaseIdx));
else if ((key.first == "BRPV") || else if ((key.first == "BRPV") || (key.first == "BOPV") || (key.first == "BWPV")
(key.first == "BOPV") || || (key.first == "BGPV")) {
(key.first == "BWPV") ||
(key.first == "BGPV"))
{
if (key.first == "BRPV") { if (key.first == "BRPV") {
val.second = 1.0; val.second = 1.0;
} } else if (key.first == "BOPV") {
else if (key.first == "BOPV") {
val.second = getValue(fs.saturation(oilPhaseIdx)); val.second = getValue(fs.saturation(oilPhaseIdx));
} } else if (key.first == "BWPV") {
else if (key.first == "BWPV") {
val.second = getValue(fs.saturation(waterPhaseIdx)); val.second = getValue(fs.saturation(waterPhaseIdx));
} } else {
else {
val.second = getValue(fs.saturation(gasPhaseIdx)); val.second = getValue(fs.saturation(gasPhaseIdx));
} }
// Include active pore-volume. // Include active pore-volume.
val.second *= elemCtx.simulator().model().dofTotalVolume(globalDofIdx) val.second *= elemCtx.simulator().model().dofTotalVolume(globalDofIdx)
* getValue(intQuants.porosity()); * getValue(intQuants.porosity());
} } else if (key.first == "BRS")
else if (key.first == "BRS")
val.second = getValue(fs.Rs()); val.second = getValue(fs.Rs());
else if (key.first == "BRV") else if (key.first == "BRV")
val.second = getValue(fs.Rv()); val.second = getValue(fs.Rv());
else if ((key.first == "BOIP") || else if ((key.first == "BOIP") || (key.first == "BOIPL") || (key.first == "BOIPG")
(key.first == "BOIPL") || || (key.first == "BGIP") || (key.first == "BGIPL") || (key.first == "BGIPG")
(key.first == "BOIPG") || || (key.first == "BWIP")) {
(key.first == "BGIP") ||
(key.first == "BGIPL") ||
(key.first == "BGIPG") ||
(key.first == "BWIP"))
{
if ((key.first == "BOIP") || (key.first == "BOIPL")) { if ((key.first == "BOIP") || (key.first == "BOIPL")) {
val.second = getValue(fs.invB(oilPhaseIdx)) val.second = getValue(fs.invB(oilPhaseIdx)) * getValue(fs.saturation(oilPhaseIdx));
* getValue(fs.saturation(oilPhaseIdx));
if (key.first == "BOIP") { if (key.first == "BOIP") {
val.second += getValue(fs.Rv()) * getValue(fs.invB(gasPhaseIdx)) val.second += getValue(fs.Rv()) * getValue(fs.invB(gasPhaseIdx))
* getValue(fs.saturation(gasPhaseIdx)); * getValue(fs.saturation(gasPhaseIdx));
} }
} } else if (key.first == "BOIPG") {
else if (key.first == "BOIPG") {
val.second = getValue(fs.Rv()) * getValue(fs.invB(gasPhaseIdx)) val.second = getValue(fs.Rv()) * getValue(fs.invB(gasPhaseIdx))
* getValue(fs.saturation(gasPhaseIdx)); * getValue(fs.saturation(gasPhaseIdx));
} } else if ((key.first == "BGIP") || (key.first == "BGIPG")) {
else if ((key.first == "BGIP") || (key.first == "BGIPG")) { val.second = getValue(fs.invB(gasPhaseIdx)) * getValue(fs.saturation(gasPhaseIdx));
val.second = getValue(fs.invB(gasPhaseIdx))
* getValue(fs.saturation(gasPhaseIdx));
if (key.first == "BGIP") { if (key.first == "BGIP") {
val.second += getValue(fs.Rs()) * getValue(fs.invB(oilPhaseIdx)) val.second += getValue(fs.Rs()) * getValue(fs.invB(oilPhaseIdx))
* getValue(fs.saturation(oilPhaseIdx)); * getValue(fs.saturation(oilPhaseIdx));
} }
} } else if (key.first == "BGIPL") {
else if (key.first == "BGIPL") {
val.second = getValue(fs.Rs()) * getValue(fs.invB(oilPhaseIdx)) val.second = getValue(fs.Rs()) * getValue(fs.invB(oilPhaseIdx))
* getValue(fs.saturation(oilPhaseIdx)); * getValue(fs.saturation(oilPhaseIdx));
} } else { // BWIP
else { // BWIP val.second = getValue(fs.invB(waterPhaseIdx)) * getValue(fs.saturation(waterPhaseIdx));
val.second = getValue(fs.invB(waterPhaseIdx))
* getValue(fs.saturation(waterPhaseIdx));
} }
// Include active pore-volume. // Include active pore-volume.
val.second *= elemCtx.simulator().model().dofTotalVolume(globalDofIdx) val.second *= elemCtx.simulator().model().dofTotalVolume(globalDofIdx)
* getValue(intQuants.porosity()); * getValue(intQuants.porosity());
} } else {
else {
std::string logstring = "Keyword '"; std::string logstring = "Keyword '";
logstring.append(key.first); logstring.append(key.first);
logstring.append("' is unhandled for output to file."); logstring.append("' is unhandled for output to file.");
@ -705,7 +686,8 @@ public:
if (this->tracerConcentrations_[tracerIdx].empty()) if (this->tracerConcentrations_[tracerIdx].empty())
continue; continue;
this->tracerConcentrations_[tracerIdx][globalDofIdx] = tracerModel.tracerConcentration(tracerIdx, globalDofIdx); this->tracerConcentrations_[tracerIdx][globalDofIdx]
= tracerModel.tracerConcentration(tracerIdx, globalDofIdx);
} }
} }
} }
@ -740,20 +722,13 @@ public:
* to globally unique Cartesian cell/element index. * to globally unique Cartesian cell/element index.
*/ */
template <class ActiveIndex, class CartesianIndex> template <class ActiveIndex, class CartesianIndex>
void processFluxes(const ElementContext& elemCtx, void processFluxes(const ElementContext& elemCtx, ActiveIndex&& activeIndex, CartesianIndex&& cartesianIndex)
ActiveIndex&& activeIndex,
CartesianIndex&& cartesianIndex)
{
const auto identifyCell = [&activeIndex, &cartesianIndex](const Element& elem)
-> EclInterRegFlowMap::Cell
{ {
const auto identifyCell = [&activeIndex, &cartesianIndex](const Element& elem) -> EclInterRegFlowMap::Cell {
const auto cellIndex = activeIndex(elem); const auto cellIndex = activeIndex(elem);
return { return {
static_cast<int>(cellIndex), static_cast<int>(cellIndex), cartesianIndex(cellIndex), elem.partitionType() == Dune::InteriorEntity};
cartesianIndex(cellIndex),
elem.partitionType() == Dune::InteriorEntity
};
}; };
const auto timeIdx = 0u; const auto timeIdx = 0u;
@ -765,8 +740,7 @@ public:
const auto left = identifyCell(stencil.element(face.interiorIndex())); const auto left = identifyCell(stencil.element(face.interiorIndex()));
const auto right = identifyCell(stencil.element(face.exteriorIndex())); const auto right = identifyCell(stencil.element(face.exteriorIndex()));
const auto rates = this-> const auto rates = this->getComponentSurfaceRates(elemCtx, face.area(), scvfIdx, timeIdx);
getComponentSurfaceRates(elemCtx, face.area(), scvfIdx, timeIdx);
this->interRegionFlows_.addConnection(left, right, rates); this->interRegionFlows_.addConnection(left, right, rates);
} }
@ -809,20 +783,20 @@ public:
fs.setSaturation(phaseIdx, this->saturation_[phaseIdx][elemIdx]); fs.setSaturation(phaseIdx, this->saturation_[phaseIdx][elemIdx]);
} }
if (!this->oilPressure_.empty()) { if (!this->fluidPressure_.empty()) {
// this assumes that capillary pressures only depend on the phase saturations // this assumes that capillary pressures only depend on the phase saturations
// and possibly on temperature. (this is always the case for ECL problems.) // and possibly on temperature. (this is always the case for ECL problems.)
std::array<Scalar, numPhases> pc = {0}; std::array<Scalar, numPhases> pc = {0};
const MaterialLawParams& matParams = simulator_.problem().materialLawParams(elemIdx); const MaterialLawParams& matParams = simulator_.problem().materialLawParams(elemIdx);
MaterialLaw::capillaryPressures(pc, matParams, fs); MaterialLaw::capillaryPressures(pc, matParams, fs);
Valgrind::CheckDefined(this->oilPressure_[elemIdx]); Valgrind::CheckDefined(this->fluidPressure_[elemIdx]);
Valgrind::CheckDefined(pc); Valgrind::CheckDefined(pc);
assert(FluidSystem::phaseIsActive(oilPhaseIdx)); assert(FluidSystem::phaseIsActive(oilPhaseIdx));
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) if (!FluidSystem::phaseIsActive(phaseIdx))
continue; continue;
fs.setPressure(phaseIdx, this->oilPressure_[elemIdx] + (pc[phaseIdx] - pc[oilPhaseIdx])); fs.setPressure(phaseIdx, this->fluidPressure_[elemIdx] + (pc[phaseIdx] - pc[oilPhaseIdx]));
} }
} }
@ -848,24 +822,19 @@ public:
if (!this->pcSwMdcOw_.empty() && !this->krnSwMdcOw_.empty()) { if (!this->pcSwMdcOw_.empty() && !this->krnSwMdcOw_.empty()) {
matLawManager->setOilWaterHysteresisParams( matLawManager->setOilWaterHysteresisParams(
this->pcSwMdcOw_[elemIdx], this->pcSwMdcOw_[elemIdx], this->krnSwMdcOw_[elemIdx], elemIdx);
this->krnSwMdcOw_[elemIdx],
elemIdx);
} }
if (!this->pcSwMdcGo_.empty() && !this->krnSwMdcGo_.empty()) { if (!this->pcSwMdcGo_.empty() && !this->krnSwMdcGo_.empty()) {
matLawManager->setGasOilHysteresisParams( matLawManager->setGasOilHysteresisParams(
this->pcSwMdcGo_[elemIdx], this->pcSwMdcGo_[elemIdx], this->krnSwMdcGo_[elemIdx], elemIdx);
this->krnSwMdcGo_[elemIdx],
elemIdx);
} }
} }
if (simulator_.vanguard().eclState().fieldProps().has_double("SWATINIT")) { if (simulator_.vanguard().eclState().fieldProps().has_double("SWATINIT")) {
const auto& oilWaterScaledEpsInfoDrainage = const auto& oilWaterScaledEpsInfoDrainage
simulator.problem().materialLawManager()->oilWaterScaledEpsInfoDrainage(elemIdx); = simulator.problem().materialLawManager()->oilWaterScaledEpsInfoDrainage(elemIdx);
const_cast<EclEpsScalingPointsInfo<Scalar>&>(oilWaterScaledEpsInfoDrainage).maxPcow = this->ppcw_[elemIdx]; const_cast<EclEpsScalingPointsInfo<Scalar>&>(oilWaterScaledEpsInfoDrainage).maxPcow = this->ppcw_[elemIdx];
} }
} }
private: private:
@ -875,8 +844,7 @@ private:
return false; return false;
const auto& parallelWells = simulator_.vanguard().parallelWells(); const auto& parallelWells = simulator_.vanguard().parallelWells();
std::pair<std::string, bool> value {wname, true}; std::pair<std::string, bool> value {wname, true};
auto candidate = std::lower_bound(parallelWells.begin(), parallelWells.end(), auto candidate = std::lower_bound(parallelWells.begin(), parallelWells.end(), value);
value);
return candidate == parallelWells.end() || *candidate != value; return candidate == parallelWells.end() || *candidate != value;
} }
@ -896,16 +864,14 @@ private:
// PORV, MINPV and friends). Also note that because of this, the porosity // PORV, MINPV and friends). Also note that because of this, the porosity
// returned by the intensive quantities can be outside of the physical // returned by the intensive quantities can be outside of the physical
// range [0, 1] in pathetic cases. // range [0, 1] in pathetic cases.
const auto totVolume = const auto totVolume = elemCtx.simulator().model().dofTotalVolume(globalDofIdx);
elemCtx.simulator().model().dofTotalVolume(globalDofIdx);
const double pv = totVolume * intQuants.porosity().value(); const double pv = totVolume * intQuants.porosity().value();
if (!this->pressureTimesHydrocarbonVolume_.empty() && !this->pressureTimesPoreVolume_.empty()) { if (!this->pressureTimesHydrocarbonVolume_.empty() && !this->pressureTimesPoreVolume_.empty()) {
assert(this->hydrocarbonPoreVolume_.size() == this->pressureTimesHydrocarbonVolume_.size()); assert(this->hydrocarbonPoreVolume_.size() == this->pressureTimesHydrocarbonVolume_.size());
assert(this->fip_[Inplace::Phase::PoreVolume].size() == this->pressureTimesPoreVolume_.size()); assert(this->fip_[Inplace::Phase::PoreVolume].size() == this->pressureTimesPoreVolume_.size());
this->fip_[Inplace::Phase::PoreVolume][globalDofIdx] = this->fip_[Inplace::Phase::PoreVolume][globalDofIdx] = totVolume * intQuants.referencePorosity();
totVolume * intQuants.referencePorosity();
this->dynamicPoreVolume_[globalDofIdx] = pv; this->dynamicPoreVolume_[globalDofIdx] = pv;
@ -919,10 +885,12 @@ private:
if (FluidSystem::phaseIsActive(oilPhaseIdx)) { if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
this->pressureTimesPoreVolume_[globalDofIdx] = getValue(fs.pressure(oilPhaseIdx)) * pv; this->pressureTimesPoreVolume_[globalDofIdx] = getValue(fs.pressure(oilPhaseIdx)) * pv;
this->pressureTimesHydrocarbonVolume_[globalDofIdx] = this->pressureTimesPoreVolume_[globalDofIdx] * hydrocarbon; this->pressureTimesHydrocarbonVolume_[globalDofIdx]
= this->pressureTimesPoreVolume_[globalDofIdx] * hydrocarbon;
} else if (FluidSystem::phaseIsActive(gasPhaseIdx)) { } else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
this->pressureTimesPoreVolume_[globalDofIdx] = getValue(fs.pressure(gasPhaseIdx)) * pv; this->pressureTimesPoreVolume_[globalDofIdx] = getValue(fs.pressure(gasPhaseIdx)) * pv;
this->pressureTimesHydrocarbonVolume_[globalDofIdx] = this->pressureTimesPoreVolume_[globalDofIdx] * hydrocarbon; this->pressureTimesHydrocarbonVolume_[globalDofIdx]
= this->pressureTimesPoreVolume_[globalDofIdx] * hydrocarbon;
} else if (FluidSystem::phaseIsActive(waterPhaseIdx)) { } else if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
this->pressureTimesPoreVolume_[globalDofIdx] = getValue(fs.pressure(waterPhaseIdx)) * pv; this->pressureTimesPoreVolume_[globalDofIdx] = getValue(fs.pressure(waterPhaseIdx)) * pv;
} }
@ -983,7 +951,6 @@ private:
this->fip_[Inplace::Phase::GAS][globalDofIdx] += gasInPlaceLiquid; this->fip_[Inplace::Phase::GAS][globalDofIdx] += gasInPlaceLiquid;
} }
} }
} }
void createLocalRegion_(std::vector<int>& region) void createLocalRegion_(std::vector<int>& region)
@ -1011,8 +978,7 @@ private:
* *
* \return Surface level component flow rates. * \return Surface level component flow rates.
*/ */
data::InterRegFlowMap::FlowRates data::InterRegFlowMap::FlowRates getComponentSurfaceRates(const ElementContext& elemCtx,
getComponentSurfaceRates(const ElementContext& elemCtx,
const Scalar faceArea, const Scalar faceArea,
const std::size_t scvfIdx, const std::size_t scvfIdx,
const std::size_t timeIdx) const const std::size_t timeIdx) const
@ -1026,25 +992,20 @@ private:
const auto alpha = getValue(extQuant.extrusionFactor()) * faceArea; const auto alpha = getValue(extQuant.extrusionFactor()) * faceArea;
if (FluidSystem::phaseIsActive(oilPhaseIdx)) { if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
const auto& up = elemCtx const auto& up = elemCtx.intensiveQuantities(extQuant.upstreamIndex(oilPhaseIdx), timeIdx);
.intensiveQuantities(extQuant.upstreamIndex(oilPhaseIdx), timeIdx);
using FluidState = std::remove_cv_t<std::remove_reference_t< using FluidState = std::remove_cv_t<std::remove_reference_t<decltype(up.fluidState())>>;
decltype(up.fluidState())>>;
const auto pvtReg = up.pvtRegionIndex(); const auto pvtReg = up.pvtRegionIndex();
const auto bO = getValue(getInvB_<FluidSystem, FluidState, Scalar> const auto bO = getValue(getInvB_<FluidSystem, FluidState, Scalar>(up.fluidState(), oilPhaseIdx, pvtReg));
(up.fluidState(), oilPhaseIdx, pvtReg));
const auto qO = alpha * bO * getValue(extQuant.volumeFlux(oilPhaseIdx)); const auto qO = alpha * bO * getValue(extQuant.volumeFlux(oilPhaseIdx));
rates[Component::Oil] += qO; rates[Component::Oil] += qO;
if (FluidSystem::phaseIsActive(gasPhaseIdx)) { if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
const auto Rs = getValue( const auto Rs = getValue(BlackOil::getRs_<FluidSystem, FluidState, Scalar>(up.fluidState(), pvtReg));
BlackOil::getRs_<FluidSystem, FluidState, Scalar>
(up.fluidState(), pvtReg));
rates[Component::Gas] += qO * Rs; rates[Component::Gas] += qO * Rs;
rates[Component::Disgas] += qO * Rs; rates[Component::Disgas] += qO * Rs;
@ -1052,25 +1013,20 @@ private:
} }
if (FluidSystem::phaseIsActive(gasPhaseIdx)) { if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
const auto& up = elemCtx const auto& up = elemCtx.intensiveQuantities(extQuant.upstreamIndex(gasPhaseIdx), timeIdx);
.intensiveQuantities(extQuant.upstreamIndex(gasPhaseIdx), timeIdx);
using FluidState = std::remove_cv_t<std::remove_reference_t< using FluidState = std::remove_cv_t<std::remove_reference_t<decltype(up.fluidState())>>;
decltype(up.fluidState())>>;
const auto pvtReg = up.pvtRegionIndex(); const auto pvtReg = up.pvtRegionIndex();
const auto bG = getValue(getInvB_<FluidSystem, FluidState, Scalar> const auto bG = getValue(getInvB_<FluidSystem, FluidState, Scalar>(up.fluidState(), gasPhaseIdx, pvtReg));
(up.fluidState(), gasPhaseIdx, pvtReg));
const auto qG = alpha * bG * getValue(extQuant.volumeFlux(gasPhaseIdx)); const auto qG = alpha * bG * getValue(extQuant.volumeFlux(gasPhaseIdx));
rates[Component::Gas] += qG; rates[Component::Gas] += qG;
if (FluidSystem::phaseIsActive(oilPhaseIdx)) { if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
const auto Rv = getValue( const auto Rv = getValue(BlackOil::getRv_<FluidSystem, FluidState, Scalar>(up.fluidState(), pvtReg));
BlackOil::getRv_<FluidSystem, FluidState, Scalar>
(up.fluidState(), pvtReg));
rates[Component::Oil] += qG * Rv; rates[Component::Oil] += qG * Rv;
rates[Component::Vapoil] += qG * Rv; rates[Component::Vapoil] += qG * Rv;
@ -1078,19 +1034,15 @@ private:
} }
if (FluidSystem::phaseIsActive(waterPhaseIdx)) { if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
const auto& up = elemCtx const auto& up = elemCtx.intensiveQuantities(extQuant.upstreamIndex(waterPhaseIdx), timeIdx);
.intensiveQuantities(extQuant.upstreamIndex(waterPhaseIdx), timeIdx);
using FluidState = std::remove_cv_t<std::remove_reference_t< using FluidState = std::remove_cv_t<std::remove_reference_t<decltype(up.fluidState())>>;
decltype(up.fluidState())>>;
const auto pvtReg = up.pvtRegionIndex(); const auto pvtReg = up.pvtRegionIndex();
const auto bW = getValue(getInvB_<FluidSystem, FluidState, Scalar> const auto bW = getValue(getInvB_<FluidSystem, FluidState, Scalar>(up.fluidState(), waterPhaseIdx, pvtReg));
(up.fluidState(), waterPhaseIdx, pvtReg));
rates[Component::Water] += rates[Component::Water] += alpha * bW * getValue(extQuant.volumeFlux(waterPhaseIdx));
alpha * bW * getValue(extQuant.volumeFlux(waterPhaseIdx));
} }
return rates; return rates;