From 0a1210a39212d9e51d5c37fd6c9205f09f207fba Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Fri, 13 Jan 2023 10:47:23 +0100 Subject: [PATCH 1/2] output gas pressure for gas-water case --- ebos/ecloutputblackoilmodule.hh | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/ebos/ecloutputblackoilmodule.hh b/ebos/ecloutputblackoilmodule.hh index 7a3430f04..f11ce611a 100644 --- a/ebos/ecloutputblackoilmodule.hh +++ b/ebos/ecloutputblackoilmodule.hh @@ -234,11 +234,11 @@ public: if (FluidSystem::phaseIsActive(oilPhaseIdx)) { this->oilPressure_[globalDofIdx] = getValue(fs.pressure(oilPhaseIdx)); }else{ - // put pressure in oil pressure for output - if (FluidSystem::phaseIsActive(waterPhaseIdx)) { - this->oilPressure_[globalDofIdx] = getValue(fs.pressure(waterPhaseIdx)); - } else { + // put pressure in oil pressure for output. Use gas if oil is not present + if (FluidSystem::phaseIsActive(gasPhaseIdx)) { this->oilPressure_[globalDofIdx] = getValue(fs.pressure(gasPhaseIdx)); + } else { // use water if neither oil nor gas is present + this->oilPressure_[globalDofIdx] = getValue(fs.pressure(waterPhaseIdx)); } } Valgrind::CheckDefined(this->oilPressure_[globalDofIdx]); From ec2983df43ed86477bc964381c1a1d759f3e825e Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Mon, 16 Jan 2023 08:53:57 +0100 Subject: [PATCH 2/2] come code cleaning in ecloutputblackoilmodule --- ebos/eclgenericoutputblackoilmodule.cc | 8 +- ebos/eclgenericoutputblackoilmodule.hh | 6 +- ebos/ecloutputblackoilmodule.hh | 460 +++++++++++-------------- 3 files changed, 213 insertions(+), 261 deletions(-) diff --git a/ebos/eclgenericoutputblackoilmodule.cc b/ebos/eclgenericoutputblackoilmodule.cc index 08b1ed4c8..18a818fb2 100644 --- a/ebos/eclgenericoutputblackoilmodule.cc +++ b/ebos/eclgenericoutputblackoilmodule.cc @@ -682,7 +682,7 @@ assignToSolution(data::Solution& sol) {"PORV_RC", UnitSystem::measure::identity, data::TargetType::RESTART_AUXILIARY, rockCompPorvMultiplier_}, {"PPCW", UnitSystem::measure::pressure, data::TargetType::RESTART_SOLUTION, ppcw_}, {"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_}, {"PCOG", UnitSystem::measure::pressure, data::TargetType::RESTART_SOLUTION, pcog_}, {"PRES_OVB", UnitSystem::measure::pressure, data::TargetType::RESTART_SOLUTION, overburdenPressure_}, @@ -788,8 +788,8 @@ setRestart(const data::Solution& sol, assert(!saturation_[oilPhaseIdx].empty()); saturation_[oilPhaseIdx][elemIdx] = so; - if (!oilPressure_.empty() && sol.has("PRESSURE")) - oilPressure_[elemIdx] = sol.data("PRESSURE")[globalDofIndex]; + if (!fluidPressure_.empty() && sol.has("PRESSURE")) + fluidPressure_[elemIdx] = sol.data("PRESSURE")[globalDofIndex]; if (!temperature_.empty() && sol.has("TEMP")) temperature_[elemIdx] = sol.data("TEMP")[globalDofIndex]; if (!rs_.empty() && sol.has("RS")) @@ -969,7 +969,7 @@ doAllocBuffers(unsigned bufferSize, saturation_[phaseIdx].resize(bufferSize, 0.0); } // and oil pressure - oilPressure_.resize(bufferSize, 0.0); + fluidPressure_.resize(bufferSize, 0.0); rstKeywords["PRES"] = 0; rstKeywords["PRESSURE"] = 0; diff --git a/ebos/eclgenericoutputblackoilmodule.hh b/ebos/eclgenericoutputblackoilmodule.hh index 6016a9ae6..160d7d6c8 100644 --- a/ebos/eclgenericoutputblackoilmodule.hh +++ b/ebos/eclgenericoutputblackoilmodule.hh @@ -48,11 +48,11 @@ template class EclGenericOutputBlackoilModule { public: Scalar* getPRESSURE_ptr(void) { - return (this->oilPressure_.data()) ; + return (this->fluidPressure_.data()) ; }; int getPRESSURE_size( void ) { - return (this->oilPressure_.size()) ; + return (this->fluidPressure_.size()) ; }; // write cumulative production and injection reports to output @@ -411,7 +411,7 @@ protected: ScalarBuffer pressureTimesPoreVolume_; ScalarBuffer pressureTimesHydrocarbonVolume_; ScalarBuffer dynamicPoreVolume_; - ScalarBuffer oilPressure_; + ScalarBuffer fluidPressure_; ScalarBuffer temperature_; ScalarBuffer rs_; ScalarBuffer rsw_; diff --git a/ebos/ecloutputblackoilmodule.hh b/ebos/ecloutputblackoilmodule.hh index f11ce611a..e19b99091 100644 --- a/ebos/ecloutputblackoilmodule.hh +++ b/ebos/ecloutputblackoilmodule.hh @@ -32,8 +32,8 @@ #include -#include #include +#include #include @@ -59,37 +59,41 @@ #include #include -namespace Opm::Properties { +namespace Opm::Properties +{ // create new type tag for the Ecl-output -namespace TTag { -struct EclOutputBlackOil {}; -} +namespace TTag +{ + struct EclOutputBlackOil { + }; +} // namespace TTag -template +template struct ForceDisableFluidInPlaceOutput { using type = UndefinedProperty; }; -template +template struct ForceDisableFluidInPlaceOutput { static constexpr bool value = false; }; -template +template struct ForceDisableResvFluidInPlaceOutput { using type = UndefinedProperty; }; -template +template struct ForceDisableResvFluidInPlaceOutput { static constexpr bool value = false; }; } // namespace Opm::Properties -namespace Opm { +namespace Opm +{ // forward declaration template @@ -128,9 +132,9 @@ class EclOutputBlackOilModule : public EclGenericOutputBlackoilModule - EclOutputBlackOilModule(const Simulator& simulator, + EclOutputBlackOilModule(const Simulator& simulator, const std::vector& wbp_index_list, - const CollectDataToIORankType& collectToIORank) + const CollectDataToIORankType& collectToIORank) : BaseType(simulator.vanguard().eclState(), simulator.vanguard().schedule(), simulator.vanguard().summaryConfig(), @@ -144,19 +148,17 @@ public: getPropValue(), getPropValue(), getPropValue()) - , simulator_(simulator) + , simulator_(simulator) { for (auto& region_pair : this->regions_) { this->createLocalRegion_(region_pair.second); } for (const auto& node : this->simulator_.vanguard().summaryConfig()) { - if ((node.category() == SummaryConfigNode::Category::Block) && - collectToIORank.isCartIdxOnThisRank(node.number() - 1)) - { + if ((node.category() == SummaryConfigNode::Category::Block) + && collectToIORank.isCartIdxOnThisRank(node.number() - 1)) { this->blockData_.emplace(std::piecewise_construct, - std::forward_as_tuple(node.keyword(), - node.number()), + std::forward_as_tuple(node.keyword(), node.number()), std::forward_as_tuple(0.0)); } } @@ -167,11 +169,9 @@ public: } } - this->forceDisableFipOutput_ = - EWOMS_GET_PARAM(TypeTag, bool, ForceDisableFluidInPlaceOutput); + this->forceDisableFipOutput_ = EWOMS_GET_PARAM(TypeTag, bool, ForceDisableFluidInPlaceOutput); - this->forceDisableFipresvOutput_ = - EWOMS_GET_PARAM(TypeTag, bool, ForceDisableResvFluidInPlaceOutput); + this->forceDisableFipresvOutput_ = EWOMS_GET_PARAM(TypeTag, bool, ForceDisableResvFluidInPlaceOutput); } /*! @@ -179,19 +179,26 @@ public: */ static void registerParameters() { - EWOMS_REGISTER_PARAM(TypeTag, bool, ForceDisableFluidInPlaceOutput, - "Do not print fluid-in-place values after each report step even if requested by the deck."); - EWOMS_REGISTER_PARAM(TypeTag, bool, ForceDisableResvFluidInPlaceOutput, - "Do not print reservoir volumes values after each report step even if requested by the deck."); + EWOMS_REGISTER_PARAM( + TypeTag, + bool, + ForceDisableFluidInPlaceOutput, + "Do not print fluid-in-place values after each report step even if requested by the deck."); + EWOMS_REGISTER_PARAM( + TypeTag, + bool, + ForceDisableResvFluidInPlaceOutput, + "Do not print reservoir volumes values after each report step even if requested by the deck."); } /*! * \brief Allocate memory for the scalar fields we would like to * 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 >::value) + if (!std::is_same>::value) return; this->doAllocBuffers(bufferSize, @@ -210,7 +217,7 @@ public: */ void processElement(const ElementContext& elemCtx) { - if (!std::is_same >::value) + if (!std::is_same>::value) return; const auto& problem = elemCtx.simulator().problem(); @@ -222,7 +229,7 @@ public: unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0); unsigned pvtRegionIdx = elemCtx.primaryVars(dofIdx, /*timeIdx=*/0).pvtRegionIndex(); - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { if (this->saturation_[phaseIdx].empty()) continue; @@ -230,18 +237,18 @@ public: Valgrind::CheckDefined(this->saturation_[phaseIdx][globalDofIdx]); } - if (!this->oilPressure_.empty()) { + if (!this->fluidPressure_.empty()) { if (FluidSystem::phaseIsActive(oilPhaseIdx)) { - this->oilPressure_[globalDofIdx] = getValue(fs.pressure(oilPhaseIdx)); - }else{ - // put pressure in oil pressure for output. Use gas if oil is not present - if (FluidSystem::phaseIsActive(gasPhaseIdx)) { - this->oilPressure_[globalDofIdx] = getValue(fs.pressure(gasPhaseIdx)); - } else { // use water if neither oil nor gas is present - this->oilPressure_[globalDofIdx] = getValue(fs.pressure(waterPhaseIdx)); - } + // 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 { + // Output water if neither oil nor gas is present + this->fluidPressure_[globalDofIdx] = getValue(fs.pressure(waterPhaseIdx)); } - Valgrind::CheckDefined(this->oilPressure_[globalDofIdx]); + Valgrind::CheckDefined(this->fluidPressure_[globalDofIdx]); } if (!this->temperature_.empty()) { @@ -250,35 +257,34 @@ public: } if (!this->gasDissolutionFactor_.empty()) { Scalar SoMax = elemCtx.problem().maxOilSaturation(globalDofIdx); - this->gasDissolutionFactor_[globalDofIdx] = - FluidSystem::template saturatedDissolutionFactor(fs, oilPhaseIdx, pvtRegionIdx, SoMax); + this->gasDissolutionFactor_[globalDofIdx] + = FluidSystem::template saturatedDissolutionFactor( + fs, oilPhaseIdx, pvtRegionIdx, SoMax); Valgrind::CheckDefined(this->gasDissolutionFactor_[globalDofIdx]); - } if (!this->oilVaporizationFactor_.empty()) { Scalar SoMax = elemCtx.problem().maxOilSaturation(globalDofIdx); - this->oilVaporizationFactor_[globalDofIdx] = - FluidSystem::template saturatedDissolutionFactor(fs, gasPhaseIdx, pvtRegionIdx, SoMax); + this->oilVaporizationFactor_[globalDofIdx] + = FluidSystem::template saturatedDissolutionFactor( + fs, gasPhaseIdx, pvtRegionIdx, SoMax); Valgrind::CheckDefined(this->oilVaporizationFactor_[globalDofIdx]); - } if (!this->gasFormationVolumeFactor_.empty()) { - this->gasFormationVolumeFactor_[globalDofIdx] = - 1.0/FluidSystem::template inverseFormationVolumeFactor(fs, gasPhaseIdx, pvtRegionIdx); + this->gasFormationVolumeFactor_[globalDofIdx] = 1.0 + / FluidSystem::template inverseFormationVolumeFactor( + fs, gasPhaseIdx, pvtRegionIdx); Valgrind::CheckDefined(this->gasFormationVolumeFactor_[globalDofIdx]); - } if (!this->saturatedOilFormationVolumeFactor_.empty()) { - this->saturatedOilFormationVolumeFactor_[globalDofIdx] = - 1.0/FluidSystem::template saturatedInverseFormationVolumeFactor(fs, oilPhaseIdx, pvtRegionIdx); + this->saturatedOilFormationVolumeFactor_[globalDofIdx] = 1.0 + / FluidSystem::template saturatedInverseFormationVolumeFactor( + fs, oilPhaseIdx, pvtRegionIdx); Valgrind::CheckDefined(this->saturatedOilFormationVolumeFactor_[globalDofIdx]); - } if (!this->oilSaturationPressure_.empty()) { - this->oilSaturationPressure_[globalDofIdx] = - FluidSystem::template saturationPressure(fs, oilPhaseIdx, pvtRegionIdx); + this->oilSaturationPressure_[globalDofIdx] + = FluidSystem::template saturationPressure(fs, oilPhaseIdx, pvtRegionIdx); Valgrind::CheckDefined(this->oilSaturationPressure_[globalDofIdx]); - } if (!this->rs_.empty()) { @@ -308,7 +314,7 @@ public: Valgrind::CheckDefined(this->rvw_[globalDofIdx]); } - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { if (this->invB_[phaseIdx].empty()) continue; @@ -316,7 +322,7 @@ public: Valgrind::CheckDefined(this->invB_[phaseIdx][globalDofIdx]); } - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { if (this->density_[phaseIdx].empty()) continue; @@ -324,24 +330,25 @@ public: Valgrind::CheckDefined(this->density_[phaseIdx][globalDofIdx]); } - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { if (this->viscosity_[phaseIdx].empty()) continue; - if (!this->extboX_.empty() && phaseIdx==oilPhaseIdx) + if (!this->extboX_.empty() && phaseIdx == oilPhaseIdx) this->viscosity_[phaseIdx][globalDofIdx] = getValue(intQuants.oilViscosity()); - else if (!this->extboX_.empty() && phaseIdx==gasPhaseIdx) + else if (!this->extboX_.empty() && phaseIdx == gasPhaseIdx) this->viscosity_[phaseIdx][globalDofIdx] = getValue(intQuants.gasViscosity()); else this->viscosity_[phaseIdx][globalDofIdx] = getValue(fs.viscosity(phaseIdx)); Valgrind::CheckDefined(this->viscosity_[phaseIdx][globalDofIdx]); } - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { if (this->relativePermeability_[phaseIdx].empty()) continue; - this->relativePermeability_[phaseIdx][globalDofIdx] = getValue(intQuants.relativePermeability(phaseIdx)); + this->relativePermeability_[phaseIdx][globalDofIdx] + = getValue(intQuants.relativePermeability(phaseIdx)); Valgrind::CheckDefined(this->relativePermeability_[phaseIdx][globalDofIdx]); } @@ -372,7 +379,7 @@ public: if (!this->permFact_.empty()) { this->permFact_[globalDofIdx] = intQuants.permFactor().value(); } - + if (!this->extboX_.empty()) { this->extboX_[globalDofIdx] = intQuants.xVolume().value(); } @@ -386,19 +393,23 @@ public: } if (!this->mFracCo2_.empty()) { - const Scalar stdVolOil = getValue(fs.saturation(oilPhaseIdx))*getValue(fs.invB(oilPhaseIdx)) - + 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()) - + getValue(fs.saturation(oilPhaseIdx))*getValue(fs.invB(oilPhaseIdx))*getValue(fs.Rs())*(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 rhoG= FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx); - const Scalar rhoCO2= intQuants.zRefDensity(); - const Scalar stdMassTotal= 1.0e-10 + stdVolOil*rhoO + stdVolGas*rhoG + stdVolCo2*rhoCO2; - this->mFracOil_[globalDofIdx] = stdVolOil*rhoO/stdMassTotal; - this->mFracGas_[globalDofIdx] = stdVolGas*rhoG/stdMassTotal; - this->mFracCo2_[globalDofIdx] = stdVolCo2*rhoCO2/stdMassTotal; + const Scalar stdVolOil = getValue(fs.saturation(oilPhaseIdx)) * getValue(fs.invB(oilPhaseIdx)) + + 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()) + + getValue(fs.saturation(oilPhaseIdx)) * getValue(fs.invB(oilPhaseIdx)) * getValue(fs.Rs()) + * (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 rhoG = FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx); + const Scalar rhoCO2 = intQuants.zRefDensity(); + const Scalar stdMassTotal = 1.0e-10 + stdVolOil * rhoO + stdVolGas * rhoG + stdVolCo2 * rhoCO2; + this->mFracOil_[globalDofIdx] = stdVolOil * rhoO / stdMassTotal; + this->mFracGas_[globalDofIdx] = stdVolGas * rhoG / stdMassTotal; + this->mFracCo2_[globalDofIdx] = stdVolCo2 * rhoCO2 / stdMassTotal; } if (!this->cMicrobes_.empty()) { @@ -410,7 +421,9 @@ public: } 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()) { @@ -423,67 +436,62 @@ public: if (!this->bubblePointPressure_.empty()) { try { - this->bubblePointPressure_[globalDofIdx] = getValue(FluidSystem::bubblePointPressure(fs, intQuants.pvtRegionIndex())); - } - catch (const NumericalProblem&) { + this->bubblePointPressure_[globalDofIdx] + = getValue(FluidSystem::bubblePointPressure(fs, intQuants.pvtRegionIndex())); + } catch (const NumericalProblem&) { const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx); this->failedCellsPb_.push_back(cartesianIdx); } } if (!this->dewPointPressure_.empty()) { try { - this->dewPointPressure_[globalDofIdx] = getValue(FluidSystem::dewPointPressure(fs, intQuants.pvtRegionIndex())); - } - catch (const NumericalProblem&) { + this->dewPointPressure_[globalDofIdx] + = getValue(FluidSystem::dewPointPressure(fs, intQuants.pvtRegionIndex())); + } catch (const NumericalProblem&) { const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx); this->failedCellsPd_.push_back(cartesianIdx); } } if (!this->soMax_.empty()) - this->soMax_[globalDofIdx] = - std::max(getValue(fs.saturation(oilPhaseIdx)), - problem.maxOilSaturation(globalDofIdx)); + this->soMax_[globalDofIdx] + = std::max(getValue(fs.saturation(oilPhaseIdx)), problem.maxOilSaturation(globalDofIdx)); if (!this->swMax_.empty()) - this->swMax_[globalDofIdx] = - std::max(getValue(fs.saturation(waterPhaseIdx)), - problem.maxWaterSaturation(globalDofIdx)); + this->swMax_[globalDofIdx] + = std::max(getValue(fs.saturation(waterPhaseIdx)), problem.maxWaterSaturation(globalDofIdx)); if (!this->minimumOilPressure_.empty()) - this->minimumOilPressure_[globalDofIdx] = - std::min(getValue(fs.pressure(oilPhaseIdx)), - problem.minOilPressure(globalDofIdx)); + this->minimumOilPressure_[globalDofIdx] + = std::min(getValue(fs.pressure(oilPhaseIdx)), problem.minOilPressure(globalDofIdx)); if (!this->overburdenPressure_.empty()) this->overburdenPressure_[globalDofIdx] = problem.overburdenPressure(globalDofIdx); if (!this->rockCompPorvMultiplier_.empty()) - this->rockCompPorvMultiplier_[globalDofIdx] = problem.template rockCompPoroMultiplier(intQuants, globalDofIdx); + this->rockCompPorvMultiplier_[globalDofIdx] + = problem.template rockCompPoroMultiplier(intQuants, globalDofIdx); if (!this->rockCompTransMultiplier_.empty()) - this->rockCompTransMultiplier_[globalDofIdx] = problem.template rockCompTransMultiplier(intQuants, globalDofIdx); + this->rockCompTransMultiplier_[globalDofIdx] + = problem.template rockCompTransMultiplier(intQuants, globalDofIdx); const auto& matLawManager = problem.materialLawManager(); if (matLawManager->enableHysteresis()) { if (!this->pcSwMdcOw_.empty() && !this->krnSwMdcOw_.empty()) { matLawManager->oilWaterHysteresisParams( - this->pcSwMdcOw_[globalDofIdx], - this->krnSwMdcOw_[globalDofIdx], - globalDofIdx); + this->pcSwMdcOw_[globalDofIdx], this->krnSwMdcOw_[globalDofIdx], globalDofIdx); } if (!this->pcSwMdcGo_.empty() && !this->krnSwMdcGo_.empty()) { matLawManager->gasOilHysteresisParams( - this->pcSwMdcGo_[globalDofIdx], - this->krnSwMdcGo_[globalDofIdx], - globalDofIdx); + this->pcSwMdcGo_[globalDofIdx], this->krnSwMdcGo_[globalDofIdx], globalDofIdx); } } if (!this->ppcw_.empty()) { this->ppcw_[globalDofIdx] = matLawManager->oilWaterScaledEpsInfoDrainage(globalDofIdx).maxPcow; - //printf("ppcw_[%d] = %lg\n", globalDofIdx, ppcw_[globalDofIdx]); + // printf("ppcw_[%d] = %lg\n", globalDofIdx, ppcw_[globalDofIdx]); } // hack to make the intial output of rs and rv Ecl compatible. // For cells with swat == 1 Ecl outputs; rs = rsSat and rv=rvSat, in all but the initial step @@ -491,7 +499,8 @@ public: // 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. // 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); @@ -510,30 +519,24 @@ public: // re-compute the volume factors, viscosities and densities if asked for if (!this->density_[oilPhaseIdx].empty()) - this->density_[oilPhaseIdx][globalDofIdx] = FluidSystem::density(fsInitial, - oilPhaseIdx, - intQuants.pvtRegionIndex()); + this->density_[oilPhaseIdx][globalDofIdx] + = FluidSystem::density(fsInitial, oilPhaseIdx, intQuants.pvtRegionIndex()); if (!this->density_[gasPhaseIdx].empty()) - this->density_[gasPhaseIdx][globalDofIdx] = FluidSystem::density(fsInitial, - gasPhaseIdx, - intQuants.pvtRegionIndex()); + this->density_[gasPhaseIdx][globalDofIdx] + = FluidSystem::density(fsInitial, gasPhaseIdx, intQuants.pvtRegionIndex()); if (!this->invB_[oilPhaseIdx].empty()) - this->invB_[oilPhaseIdx][globalDofIdx] = FluidSystem::inverseFormationVolumeFactor(fsInitial, - oilPhaseIdx, - intQuants.pvtRegionIndex()); + this->invB_[oilPhaseIdx][globalDofIdx] + = FluidSystem::inverseFormationVolumeFactor(fsInitial, oilPhaseIdx, intQuants.pvtRegionIndex()); if (!this->invB_[gasPhaseIdx].empty()) - this->invB_[gasPhaseIdx][globalDofIdx] = FluidSystem::inverseFormationVolumeFactor(fsInitial, - gasPhaseIdx, - intQuants.pvtRegionIndex()); + this->invB_[gasPhaseIdx][globalDofIdx] + = FluidSystem::inverseFormationVolumeFactor(fsInitial, gasPhaseIdx, intQuants.pvtRegionIndex()); if (!this->viscosity_[oilPhaseIdx].empty()) - this->viscosity_[oilPhaseIdx][globalDofIdx] = FluidSystem::viscosity(fsInitial, - oilPhaseIdx, - intQuants.pvtRegionIndex()); + this->viscosity_[oilPhaseIdx][globalDofIdx] + = FluidSystem::viscosity(fsInitial, oilPhaseIdx, intQuants.pvtRegionIndex()); if (!this->viscosity_[gasPhaseIdx].empty()) - this->viscosity_[gasPhaseIdx][globalDofIdx] = FluidSystem::viscosity(fsInitial, - gasPhaseIdx, - intQuants.pvtRegionIndex()); + this->viscosity_[gasPhaseIdx][globalDofIdx] + = FluidSystem::viscosity(fsInitial, gasPhaseIdx, intQuants.pvtRegionIndex()); } // Add fluid in Place values @@ -541,7 +544,7 @@ public: // Adding block data const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx); - for (auto& val: this->blockData_) { + for (auto& val : this->blockData_) { const auto& key = val.first; assert(key.second > 0); unsigned int cartesianIdxBlock = key.second - 1; @@ -554,23 +557,21 @@ public: val.second = getValue(fs.saturation(oilPhaseIdx)); else if (key.first == "BNSAT") val.second = intQuants.solventSaturation().value(); - else if ((key.first == "BPR") || (key.first == "BPRESSUR")){ - if (FluidSystem::phaseIsActive(oilPhaseIdx)) + else if ((key.first == "BPR") || (key.first == "BPRESSUR")) { + if (FluidSystem::phaseIsActive(oilPhaseIdx)) val.second = getValue(fs.pressure(oilPhaseIdx)); - else if (FluidSystem::phaseIsActive(gasPhaseIdx)) + else if (FluidSystem::phaseIsActive(gasPhaseIdx)) val.second = getValue(fs.pressure(gasPhaseIdx)); - else if (FluidSystem::phaseIsActive(waterPhaseIdx)) + else if (FluidSystem::phaseIsActive(waterPhaseIdx)) val.second = getValue(fs.pressure(waterPhaseIdx)); - } - else if ((key.first == "BTCNFHEA") || (key.first == "BTEMP")){ - if (FluidSystem::phaseIsActive(oilPhaseIdx)) + } else if ((key.first == "BTCNFHEA") || (key.first == "BTEMP")) { + if (FluidSystem::phaseIsActive(oilPhaseIdx)) val.second = getValue(fs.temperature(oilPhaseIdx)); - else if (FluidSystem::phaseIsActive(gasPhaseIdx)) + else if (FluidSystem::phaseIsActive(gasPhaseIdx)) val.second = getValue(fs.temperature(gasPhaseIdx)); - else if (FluidSystem::phaseIsActive(waterPhaseIdx)) + else if (FluidSystem::phaseIsActive(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)); else if (key.first == "BGKR" || key.first == "BKRG") val.second = getValue(intQuants.relativePermeability(gasPhaseIdx)); @@ -578,15 +579,15 @@ public: val.second = getValue(intQuants.relativePermeability(oilPhaseIdx)); else if (key.first == "BKROG") { const auto& materialParams = problem.materialLawParams(elemCtx, dofIdx, /* timeIdx = */ 0); - const auto krog = MaterialLaw::template relpermOilInOilGasSystem(materialParams, fs); + const auto krog + = MaterialLaw::template relpermOilInOilGasSystem(materialParams, fs); 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 krow = MaterialLaw::template relpermOilInOilWaterSystem(materialParams, fs); + const auto krow + = MaterialLaw::template relpermOilInOilWaterSystem(materialParams, fs); val.second = getValue(krow); - } - else if (key.first == "BWPC") + } else if (key.first == "BWPC") val.second = getValue(fs.pressure(oilPhaseIdx)) - getValue(fs.pressure(waterPhaseIdx)); else if (key.first == "BGPC") val.second = getValue(fs.pressure(gasPhaseIdx)) - getValue(fs.pressure(oilPhaseIdx)); @@ -600,76 +601,56 @@ public: val.second = getValue(fs.viscosity(gasPhaseIdx)); else if (key.first == "BVOIL" || key.first == "BOVIS") val.second = getValue(fs.viscosity(oilPhaseIdx)); - else if ((key.first == "BRPV") || - (key.first == "BOPV") || - (key.first == "BWPV") || - (key.first == "BGPV")) - { + else if ((key.first == "BRPV") || (key.first == "BOPV") || (key.first == "BWPV") + || (key.first == "BGPV")) { if (key.first == "BRPV") { val.second = 1.0; - } - else if (key.first == "BOPV") { + } else if (key.first == "BOPV") { val.second = getValue(fs.saturation(oilPhaseIdx)); - } - else if (key.first == "BWPV") { + } else if (key.first == "BWPV") { val.second = getValue(fs.saturation(waterPhaseIdx)); - } - else { + } else { val.second = getValue(fs.saturation(gasPhaseIdx)); } // Include active pore-volume. val.second *= elemCtx.simulator().model().dofTotalVolume(globalDofIdx) * getValue(intQuants.porosity()); - } - else if (key.first == "BRS") + } else if (key.first == "BRS") val.second = getValue(fs.Rs()); else if (key.first == "BRV") val.second = getValue(fs.Rv()); - else if ((key.first == "BOIP") || - (key.first == "BOIPL") || - (key.first == "BOIPG") || - (key.first == "BGIP") || - (key.first == "BGIPL") || - (key.first == "BGIPG") || - (key.first == "BWIP")) - { + else if ((key.first == "BOIP") || (key.first == "BOIPL") || (key.first == "BOIPG") + || (key.first == "BGIP") || (key.first == "BGIPL") || (key.first == "BGIPG") + || (key.first == "BWIP")) { if ((key.first == "BOIP") || (key.first == "BOIPL")) { - val.second = getValue(fs.invB(oilPhaseIdx)) - * getValue(fs.saturation(oilPhaseIdx)); + val.second = getValue(fs.invB(oilPhaseIdx)) * getValue(fs.saturation(oilPhaseIdx)); if (key.first == "BOIP") { val.second += getValue(fs.Rv()) * getValue(fs.invB(gasPhaseIdx)) * getValue(fs.saturation(gasPhaseIdx)); } - } - else if (key.first == "BOIPG") { + } else if (key.first == "BOIPG") { val.second = getValue(fs.Rv()) * getValue(fs.invB(gasPhaseIdx)) * getValue(fs.saturation(gasPhaseIdx)); - } - else if ((key.first == "BGIP") || (key.first == "BGIPG")) { - val.second = getValue(fs.invB(gasPhaseIdx)) - * getValue(fs.saturation(gasPhaseIdx)); + } else if ((key.first == "BGIP") || (key.first == "BGIPG")) { + val.second = getValue(fs.invB(gasPhaseIdx)) * getValue(fs.saturation(gasPhaseIdx)); if (key.first == "BGIP") { val.second += getValue(fs.Rs()) * getValue(fs.invB(oilPhaseIdx)) * getValue(fs.saturation(oilPhaseIdx)); } - } - else if (key.first == "BGIPL") { + } else if (key.first == "BGIPL") { val.second = getValue(fs.Rs()) * getValue(fs.invB(oilPhaseIdx)) * getValue(fs.saturation(oilPhaseIdx)); - } - else { // BWIP - val.second = getValue(fs.invB(waterPhaseIdx)) - * getValue(fs.saturation(waterPhaseIdx)); + } else { // BWIP + val.second = getValue(fs.invB(waterPhaseIdx)) * getValue(fs.saturation(waterPhaseIdx)); } // Include active pore-volume. val.second *= elemCtx.simulator().model().dofTotalVolume(globalDofIdx) * getValue(intQuants.porosity()); - } - else { + } else { std::string logstring = "Keyword '"; logstring.append(key.first); logstring.append("' is unhandled for output to file."); @@ -701,11 +682,12 @@ public: // tracers const auto& tracerModel = simulator_.problem().tracerModel(); if (!this->tracerConcentrations_.empty()) { - for (int tracerIdx = 0; tracerIdx < tracerModel.numTracers(); tracerIdx++){ + for (int tracerIdx = 0; tracerIdx < tracerModel.numTracers(); tracerIdx++) { if (this->tracerConcentrations_[tracerIdx].empty()) continue; - this->tracerConcentrations_[tracerIdx][globalDofIdx] = tracerModel.tracerConcentration(tracerIdx, globalDofIdx); + this->tracerConcentrations_[tracerIdx][globalDofIdx] + = tracerModel.tracerConcentration(tracerIdx, globalDofIdx); } } } @@ -740,33 +722,25 @@ public: * to globally unique Cartesian cell/element index. */ template - void processFluxes(const ElementContext& elemCtx, - ActiveIndex&& activeIndex, - CartesianIndex&& cartesianIndex) + void processFluxes(const ElementContext& elemCtx, 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); return { - static_cast(cellIndex), - cartesianIndex(cellIndex), - elem.partitionType() == Dune::InteriorEntity - }; + static_cast(cellIndex), cartesianIndex(cellIndex), elem.partitionType() == Dune::InteriorEntity}; }; const auto timeIdx = 0u; const auto& stencil = elemCtx.stencil(timeIdx); const auto numInteriorFaces = elemCtx.numInteriorFaces(timeIdx); - for (auto scvfIdx = 0*numInteriorFaces; scvfIdx < numInteriorFaces; ++scvfIdx) { + for (auto scvfIdx = 0 * numInteriorFaces; scvfIdx < numInteriorFaces; ++scvfIdx) { const auto& face = stencil.interiorFace(scvfIdx); - 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 rates = this-> - getComponentSurfaceRates(elemCtx, face.area(), scvfIdx, timeIdx); + const auto rates = this->getComponentSurfaceRates(elemCtx, face.area(), scvfIdx, timeIdx); this->interRegionFlows_.addConnection(left, right, rates); } @@ -802,40 +776,40 @@ public: template void assignToFluidState(FluidState& fs, unsigned elemIdx) const { - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { if (this->saturation_[phaseIdx].empty()) continue; 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 // and possibly on temperature. (this is always the case for ECL problems.) std::array pc = {0}; const MaterialLawParams& matParams = simulator_.problem().materialLawParams(elemIdx); MaterialLaw::capillaryPressures(pc, matParams, fs); - Valgrind::CheckDefined(this->oilPressure_[elemIdx]); + Valgrind::CheckDefined(this->fluidPressure_[elemIdx]); Valgrind::CheckDefined(pc); assert(FluidSystem::phaseIsActive(oilPhaseIdx)); for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { if (!FluidSystem::phaseIsActive(phaseIdx)) continue; - fs.setPressure(phaseIdx, this->oilPressure_[elemIdx] + (pc[phaseIdx] - pc[oilPhaseIdx])); + fs.setPressure(phaseIdx, this->fluidPressure_[elemIdx] + (pc[phaseIdx] - pc[oilPhaseIdx])); } } if (!this->temperature_.empty()) fs.setTemperature(this->temperature_[elemIdx]); if (!this->rs_.empty()) - fs.setRs(this->rs_[elemIdx]); + fs.setRs(this->rs_[elemIdx]); if (!this->rsw_.empty()) - fs.setRsw(this->rsw_[elemIdx]); + fs.setRsw(this->rsw_[elemIdx]); if (!this->rv_.empty()) - fs.setRv(this->rv_[elemIdx]); + fs.setRv(this->rv_[elemIdx]); if (!this->rvw_.empty()) - fs.setRvw(this->rvw_[elemIdx]); + fs.setRvw(this->rvw_[elemIdx]); } void initHysteresisParams(Simulator& simulator, unsigned elemIdx) const @@ -848,35 +822,29 @@ public: if (!this->pcSwMdcOw_.empty() && !this->krnSwMdcOw_.empty()) { matLawManager->setOilWaterHysteresisParams( - this->pcSwMdcOw_[elemIdx], - this->krnSwMdcOw_[elemIdx], - elemIdx); + this->pcSwMdcOw_[elemIdx], this->krnSwMdcOw_[elemIdx], elemIdx); } if (!this->pcSwMdcGo_.empty() && !this->krnSwMdcGo_.empty()) { matLawManager->setGasOilHysteresisParams( - this->pcSwMdcGo_[elemIdx], - this->krnSwMdcGo_[elemIdx], - elemIdx); + this->pcSwMdcGo_[elemIdx], this->krnSwMdcGo_[elemIdx], elemIdx); } } if (simulator_.vanguard().eclState().fieldProps().has_double("SWATINIT")) { - const auto& oilWaterScaledEpsInfoDrainage = - simulator.problem().materialLawManager()->oilWaterScaledEpsInfoDrainage(elemIdx); + const auto& oilWaterScaledEpsInfoDrainage + = simulator.problem().materialLawManager()->oilWaterScaledEpsInfoDrainage(elemIdx); const_cast&>(oilWaterScaledEpsInfoDrainage).maxPcow = this->ppcw_[elemIdx]; } - } private: bool isDefunctParallelWell(std::string wname) const override { - if (simulator_.gridView().comm().size()==1) + if (simulator_.gridView().comm().size() == 1) return false; const auto& parallelWells = simulator_.vanguard().parallelWells(); - std::pair value{wname, true}; - auto candidate = std::lower_bound(parallelWells.begin(), parallelWells.end(), - value); + std::pair value {wname, true}; + auto candidate = std::lower_bound(parallelWells.begin(), parallelWells.end(), value); return candidate == parallelWells.end() || *candidate != value; } @@ -896,16 +864,14 @@ private: // PORV, MINPV and friends). Also note that because of this, the porosity // returned by the intensive quantities can be outside of the physical // range [0, 1] in pathetic cases. - const auto totVolume = - elemCtx.simulator().model().dofTotalVolume(globalDofIdx); + const auto totVolume = elemCtx.simulator().model().dofTotalVolume(globalDofIdx); const double pv = totVolume * intQuants.porosity().value(); 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()); - this->fip_[Inplace::Phase::PoreVolume][globalDofIdx] = - totVolume * intQuants.referencePorosity(); + this->fip_[Inplace::Phase::PoreVolume][globalDofIdx] = totVolume * intQuants.referencePorosity(); this->dynamicPoreVolume_[globalDofIdx] = pv; @@ -919,10 +885,12 @@ private: if (FluidSystem::phaseIsActive(oilPhaseIdx)) { 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)) { 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)) { this->pressureTimesPoreVolume_[globalDofIdx] = getValue(fs.pressure(waterPhaseIdx)) * pv; } @@ -983,7 +951,6 @@ private: this->fip_[Inplace::Phase::GAS][globalDofIdx] += gasInPlaceLiquid; } } - } void createLocalRegion_(std::vector& region) @@ -1011,86 +978,71 @@ private: * * \return Surface level component flow rates. */ - data::InterRegFlowMap::FlowRates - getComponentSurfaceRates(const ElementContext& elemCtx, - const Scalar faceArea, - const std::size_t scvfIdx, - const std::size_t timeIdx) const + data::InterRegFlowMap::FlowRates getComponentSurfaceRates(const ElementContext& elemCtx, + const Scalar faceArea, + const std::size_t scvfIdx, + const std::size_t timeIdx) const { using Component = data::InterRegFlowMap::Component; - auto rates = data::InterRegFlowMap::FlowRates{}; + auto rates = data::InterRegFlowMap::FlowRates {}; const auto& extQuant = elemCtx.extensiveQuantities(scvfIdx, timeIdx); const auto alpha = getValue(extQuant.extrusionFactor()) * faceArea; if (FluidSystem::phaseIsActive(oilPhaseIdx)) { - const auto& up = elemCtx - .intensiveQuantities(extQuant.upstreamIndex(oilPhaseIdx), timeIdx); + const auto& up = elemCtx.intensiveQuantities(extQuant.upstreamIndex(oilPhaseIdx), timeIdx); - using FluidState = std::remove_cv_t>; + using FluidState = std::remove_cv_t>; const auto pvtReg = up.pvtRegionIndex(); - const auto bO = getValue(getInvB_ - (up.fluidState(), oilPhaseIdx, pvtReg)); + const auto bO = getValue(getInvB_(up.fluidState(), oilPhaseIdx, pvtReg)); const auto qO = alpha * bO * getValue(extQuant.volumeFlux(oilPhaseIdx)); rates[Component::Oil] += qO; if (FluidSystem::phaseIsActive(gasPhaseIdx)) { - const auto Rs = getValue( - BlackOil::getRs_ - (up.fluidState(), pvtReg)); + const auto Rs = getValue(BlackOil::getRs_(up.fluidState(), pvtReg)); - rates[Component::Gas] += qO * Rs; + rates[Component::Gas] += qO * Rs; rates[Component::Disgas] += qO * Rs; } } if (FluidSystem::phaseIsActive(gasPhaseIdx)) { - const auto& up = elemCtx - .intensiveQuantities(extQuant.upstreamIndex(gasPhaseIdx), timeIdx); + const auto& up = elemCtx.intensiveQuantities(extQuant.upstreamIndex(gasPhaseIdx), timeIdx); - using FluidState = std::remove_cv_t>; + using FluidState = std::remove_cv_t>; const auto pvtReg = up.pvtRegionIndex(); - const auto bG = getValue(getInvB_ - (up.fluidState(), gasPhaseIdx, pvtReg)); + const auto bG = getValue(getInvB_(up.fluidState(), gasPhaseIdx, pvtReg)); const auto qG = alpha * bG * getValue(extQuant.volumeFlux(gasPhaseIdx)); rates[Component::Gas] += qG; if (FluidSystem::phaseIsActive(oilPhaseIdx)) { - const auto Rv = getValue( - BlackOil::getRv_ - (up.fluidState(), pvtReg)); + const auto Rv = getValue(BlackOil::getRv_(up.fluidState(), pvtReg)); - rates[Component::Oil] += qG * Rv; + rates[Component::Oil] += qG * Rv; rates[Component::Vapoil] += qG * Rv; } } if (FluidSystem::phaseIsActive(waterPhaseIdx)) { - const auto& up = elemCtx - .intensiveQuantities(extQuant.upstreamIndex(waterPhaseIdx), timeIdx); + const auto& up = elemCtx.intensiveQuantities(extQuant.upstreamIndex(waterPhaseIdx), timeIdx); - using FluidState = std::remove_cv_t>; + using FluidState = std::remove_cv_t>; const auto pvtReg = up.pvtRegionIndex(); - const auto bW = getValue(getInvB_ - (up.fluidState(), waterPhaseIdx, pvtReg)); + const auto bW = getValue(getInvB_(up.fluidState(), waterPhaseIdx, pvtReg)); - rates[Component::Water] += - alpha * bW * getValue(extQuant.volumeFlux(waterPhaseIdx)); + rates[Component::Water] += alpha * bW * getValue(extQuant.volumeFlux(waterPhaseIdx)); } return rates;