27#ifndef OPM_OUTPUT_BLACK_OIL_MODULE_HPP
28#define OPM_OUTPUT_BLACK_OIL_MODULE_HPP
30#include <dune/common/fvector.hh>
34#include <opm/common/Exceptions.hpp>
35#include <opm/common/TimingMacros.hpp>
36#include <opm/common/OpmLog/OpmLog.hpp>
37#include <opm/common/utility/Visitor.hpp>
39#include <opm/input/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
41#include <opm/material/common/Valgrind.hpp>
42#include <opm/material/fluidmatrixinteractions/EclEpsScalingPoints.hpp>
43#include <opm/material/fluidstates/BlackOilFluidState.hpp>
44#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
51#include <opm/output/data/Cells.hpp>
52#include <opm/output/eclipse/EclipseIO.hpp>
53#include <opm/output/eclipse/Inplace.hpp>
74template <
class TypeTag>
75class EcfvDiscretization;
83template <
class TypeTag>
95 using FluidState =
typename IntensiveQuantities::FluidState;
97 using Element =
typename GridView::template Codim<0>::Entity;
98 using ElementIterator =
typename GridView::template Codim<0>::Iterator;
101 using Dir = FaceDir::DirEnum;
108 static constexpr int conti0EqIdx = Indices::conti0EqIdx;
109 static constexpr int numPhases = FluidSystem::numPhases;
110 static constexpr int oilPhaseIdx = FluidSystem::oilPhaseIdx;
111 static constexpr int gasPhaseIdx = FluidSystem::gasPhaseIdx;
112 static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
113 static constexpr int gasCompIdx = FluidSystem::gasCompIdx;
114 static constexpr int oilCompIdx = FluidSystem::oilCompIdx;
115 static constexpr int waterCompIdx = FluidSystem::waterCompIdx;
116 enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
117 enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
118 enum { enableVapwat = getPropValue<TypeTag, Properties::EnableVapwat>() };
119 enum { enableDisgasInWater = getPropValue<TypeTag, Properties::EnableDisgasInWater>() };
120 enum { enableDissolvedGas = Indices::compositionSwitchIdx >= 0 };
122 template<
int idx,
class VectorType>
123 static Scalar value_or_zero(
const VectorType& v)
125 if constexpr (idx == -1) {
128 return v.empty() ? 0.0 : v[idx];
134 const SummaryConfig& smryCfg,
136 :
BaseType(simulator.vanguard().eclState(),
137 simulator.vanguard().schedule(),
139 simulator.vanguard().summaryState(),
141 [this](const int idx)
142 {
return simulator_.problem().eclWriter().collectOnIORank().localIdxToGlobalIdx(idx); },
143 simulator.vanguard().grid().comm(),
144 getPropValue<TypeTag, Properties::EnableEnergy>(),
145 getPropValue<TypeTag, Properties::EnableTemperature>(),
146 getPropValue<TypeTag, Properties::EnableMech>(),
147 getPropValue<TypeTag, Properties::EnableSolvent>(),
148 getPropValue<TypeTag, Properties::EnablePolymer>(),
149 getPropValue<TypeTag, Properties::EnableFoam>(),
150 getPropValue<TypeTag, Properties::EnableBrine>(),
151 getPropValue<TypeTag, Properties::EnableSaltPrecipitation>(),
152 getPropValue<TypeTag, Properties::EnableExtbo>(),
153 getPropValue<TypeTag, Properties::EnableMICP>())
154 , simulator_(simulator)
155 , collectOnIORank_(collectOnIORank)
157 for (
auto& region_pair : this->
regions_) {
158 this->createLocalRegion_(region_pair.second);
161 auto isCartIdxOnThisRank = [&collectOnIORank](
const int idx) {
162 return collectOnIORank.isCartIdxOnThisRank(idx);
167 if (! Parameters::Get<Parameters::OwnerCellsFirst>()) {
168 const std::string msg =
"The output code does not support --owner-cells-first=false.";
169 if (collectOnIORank.isIORank()) {
172 OPM_THROW_NOLOG(std::runtime_error, msg);
175 if (smryCfg.match(
"[FB]PP[OGW]") || smryCfg.match(
"RPP[OGW]*")) {
176 auto rset = this->
eclState_.fieldProps().fip_regions();
177 rset.push_back(
"PVTNUM");
183 .emplace(this->simulator_.gridView().comm(),
184 FluidSystem::numPhases, rset,
185 [fp = std::cref(this->eclState_.fieldProps())]
186 (
const std::string& rsetName) ->
decltype(
auto)
187 { return fp.get().get_int(rsetName); });
197 const unsigned reportStepNum,
200 const bool isRestart)
206 const auto& problem = this->simulator_.problem();
213 &problem.materialLawManager()->hysteresisConfig(),
214 problem.eclWriter().getOutputNnc().size());
219 const int reportStepNum)
221 this->setupElementExtractors_();
222 this->setupBlockExtractors_(isSubStep, reportStepNum);
228 this->extractors_.clear();
229 this->blockExtractors_.clear();
230 this->extraBlockExtractors_.clear();
244 if (this->extractors_.empty()) {
248 const auto& matLawManager = simulator_.problem().materialLawManager();
251 for (
unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(0); ++dofIdx) {
252 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, 0);
253 const auto& fs = intQuants.fluidState();
256 elemCtx.globalSpaceIndex(dofIdx, 0),
257 elemCtx.primaryVars(dofIdx, 0).pvtRegionIndex(),
264 if (matLawManager->enableHysteresis()) {
265 if (FluidSystem::phaseIsActive(oilPhaseIdx) && FluidSystem::phaseIsActive(waterPhaseIdx)) {
266 matLawManager->oilWaterHysteresisParams(hysterParams.
somax,
271 if (FluidSystem::phaseIsActive(oilPhaseIdx) && FluidSystem::phaseIsActive(gasPhaseIdx)) {
272 matLawManager->gasOilHysteresisParams(hysterParams.
sgmax,
290 if (this->blockExtractors_.empty() && this->extraBlockExtractors_.empty()) {
294 for (
unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(0); ++dofIdx) {
296 const auto globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, 0);
297 const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx);
299 const auto be_it = this->blockExtractors_.find(cartesianIdx);
300 const auto bee_it = this->extraBlockExtractors_.find(cartesianIdx);
301 if (be_it == this->blockExtractors_.end() &&
302 bee_it == this->extraBlockExtractors_.end())
307 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, 0);
308 const auto& fs = intQuants.fluidState();
318 if (be_it != this->blockExtractors_.end()) {
321 if (bee_it != this->extraBlockExtractors_.end()) {
328 const std::size_t reportStepNum,
330 boost::posix_time::ptime currentDate,
335 if (comm.rank() != 0) {
340 std::unique_ptr<FIPConfig> fipSched;
341 if (reportStepNum > 0) {
342 const auto& rpt = this->
schedule_[reportStepNum-1].rpt_config.get();
343 fipSched = std::make_unique<FIPConfig>(rpt);
345 const FIPConfig& fipc = reportStepNum == 0 ? this->
eclState_.getEclipseConfig().fip()
350 this->
logOutput_.timeStamp(
"BALANCE", elapsed, reportStepNum, currentDate);
353 this->
logOutput_.fip(inplace, initial_inplace,
"");
355 if (fipc.output(FIPConfig::OutputField::FIPNUM)) {
356 this->
logOutput_.fip(inplace, initial_inplace,
"FIPNUM");
358 if (fipc.output(FIPConfig::OutputField::RESV))
362 if (fipc.output(FIPConfig::OutputField::FIP)) {
363 for (
const auto& reg : this->regions_) {
364 if (reg.first !=
"FIPNUM") {
365 std::ostringstream ss;
366 ss <<
"BAL" << reg.first.substr(3);
367 this->
logOutput_.timeStamp(ss.str(), elapsed, reportStepNum, currentDate);
368 this->
logOutput_.fip(inplace, initial_inplace, reg.first);
370 if (fipc.output(FIPConfig::OutputField::RESV))
382 if (comm.rank() != 0) {
386 if ((reportStepNum == 0) && (!substep) &&
387 (this->
schedule_.initialReportConfiguration().has_value()) &&
388 (this->schedule_.initialReportConfiguration()->contains(
"CSVFIP"))) {
390 std::ostringstream csv_stream;
396 this->
logOutput_.fip_csv(csv_stream, initial_inplace,
"FIPNUM");
398 for (
const auto& reg : this->regions_) {
399 if (reg.first !=
"FIPNUM") {
400 this->
logOutput_.fip_csv(csv_stream, initial_inplace, reg.first);
404 const IOConfig& io = this->
eclState_.getIOConfig();
405 auto csv_fname = io.getOutputDir() +
"/" + io.getBaseName() +
".CSV";
407 std::ofstream outputFile(csv_fname);
409 outputFile << csv_stream.str();
443 template <
class ActiveIndex,
class CartesianIndex>
445 ActiveIndex&& activeIndex,
446 CartesianIndex&& cartesianIndex)
449 const auto identifyCell = [&activeIndex, &cartesianIndex](
const Element& elem)
452 const auto cellIndex = activeIndex(elem);
455 static_cast<int>(cellIndex),
456 cartesianIndex(cellIndex),
457 elem.partitionType() == Dune::InteriorEntity
461 const auto timeIdx = 0u;
462 const auto& stencil = elemCtx.stencil(timeIdx);
463 const auto numInteriorFaces = elemCtx.numInteriorFaces(timeIdx);
465 for (
auto scvfIdx = 0 * numInteriorFaces; scvfIdx < numInteriorFaces; ++scvfIdx) {
466 const auto& face = stencil.interiorFace(scvfIdx);
467 const auto left = identifyCell(stencil.element(face.interiorIndex()));
468 const auto right = identifyCell(stencil.element(face.exteriorIndex()));
470 const auto rates = this->
471 getComponentSurfaceRates(elemCtx, face.area(), scvfIdx, timeIdx);
504 template <
class Flu
idState>
507 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
511 fs.setSaturation(phaseIdx, this->
saturation_[phaseIdx][elemIdx]);
517 std::array<Scalar, numPhases> pc = {0};
518 const MaterialLawParams& matParams = simulator_.problem().materialLawParams(elemIdx);
519 MaterialLaw::capillaryPressures(pc, matParams, fs);
521 Valgrind::CheckDefined(pc);
523 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
524 if (!FluidSystem::phaseIsActive(phaseIdx))
527 if (Indices::oilEnabled)
528 fs.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx]));
529 else if (Indices::gasEnabled)
530 fs.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx]));
531 else if (Indices::waterEnabled)
533 fs.setPressure(phaseIdx, pressure);
539 if constexpr (enableDissolvedGas) {
540 if (!this->
rs_.empty())
541 fs.setRs(this->rs_[elemIdx]);
542 if (!this->
rv_.empty())
543 fs.setRv(this->rv_[elemIdx]);
545 if constexpr (enableDisgasInWater) {
546 if (!this->
rsw_.empty())
547 fs.setRsw(this->rsw_[elemIdx]);
549 if constexpr (enableVapwat) {
550 if (!this->
rvw_.empty())
551 fs.setRvw(this->rvw_[elemIdx]);
557 if (!this->
soMax_.empty())
558 simulator.problem().setMaxOilSaturation(elemIdx, this->
soMax_[elemIdx]);
560 if (simulator.problem().materialLawManager()->enableHysteresis()) {
561 auto matLawManager = simulator.problem().materialLawManager();
563 if (FluidSystem::phaseIsActive(oilPhaseIdx)
564 && FluidSystem::phaseIsActive(waterPhaseIdx)) {
569 if (matLawManager->enableNonWettingHysteresis()) {
570 if (!this->
soMax_.empty()) {
571 somax = this->
soMax_[elemIdx];
574 if (matLawManager->enableWettingHysteresis()) {
575 if (!this->
swMax_.empty()) {
576 swmax = this->
swMax_[elemIdx];
579 if (matLawManager->enablePCHysteresis()) {
580 if (!this->
swmin_.empty()) {
581 swmin = this->
swmin_[elemIdx];
584 matLawManager->setOilWaterHysteresisParams(
585 somax, swmax, swmin, elemIdx);
587 if (FluidSystem::phaseIsActive(oilPhaseIdx)
588 && FluidSystem::phaseIsActive(gasPhaseIdx)) {
593 if (matLawManager->enableNonWettingHysteresis()) {
594 if (!this->
sgmax_.empty()) {
595 sgmax = this->
sgmax_[elemIdx];
598 if (matLawManager->enableWettingHysteresis()) {
599 if (!this->
shmax_.empty()) {
600 shmax = this->
shmax_[elemIdx];
603 if (matLawManager->enablePCHysteresis()) {
604 if (!this->
somin_.empty()) {
605 somin = this->
somin_[elemIdx];
608 matLawManager->setGasOilHysteresisParams(
609 sgmax, shmax, somin, elemIdx);
614 if (simulator_.vanguard().eclState().fieldProps().has_double(
"SWATINIT")) {
615 simulator.problem().materialLawManager()
616 ->applyRestartSwatInit(elemIdx, this->
ppcw_[elemIdx]);
622 for (
unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(0); ++dofIdx) {
623 updateFluidInPlace_(elemCtx, dofIdx);
628 const IntensiveQuantities& intQuants,
629 const double totVolume)
631 this->updateFluidInPlace_(globalDofIdx, intQuants, totVolume);
635 template <
typename T>
636 using RemoveCVR = std::remove_cv_t<std::remove_reference_t<T>>;
638 template <
typename,
class =
void>
639 struct HasGeoMech :
public std::false_type {};
641 template <
typename Problem>
643 Problem, std::void_t<decltype(std::declval<Problem>().geoMechModel())>
644 > :
public std::true_type {};
646 bool isDefunctParallelWell(
const std::string& wname)
const override
648 if (simulator_.gridView().comm().size() == 1)
650 const auto& parallelWells = simulator_.vanguard().parallelWells();
651 std::pair<std::string, bool> value {wname,
true};
652 auto candidate = std::lower_bound(parallelWells.begin(), parallelWells.end(), value);
653 return candidate == parallelWells.end() || *candidate != value;
656 bool isOwnedByCurrentRank(
const std::string& wname)
const override
658 return this->simulator_.problem().wellModel().isOwner(wname);
661 void updateFluidInPlace_(
const ElementContext& elemCtx,
const unsigned dofIdx)
663 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, 0);
664 const unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, 0);
665 const auto totVolume = elemCtx.simulator().model().dofTotalVolume(globalDofIdx);
667 this->updateFluidInPlace_(globalDofIdx, intQuants, totVolume);
670 void updateFluidInPlace_(
const unsigned globalDofIdx,
671 const IntensiveQuantities& intQuants,
672 const double totVolume)
676 this->updateTotalVolumesAndPressures_(globalDofIdx, intQuants, totVolume);
679 this->updatePhaseInplaceVolumes_(globalDofIdx, intQuants, totVolume);
683 void createLocalRegion_(std::vector<int>& region)
689 region.resize(simulator_.gridView().size(0));
690 std::size_t elemIdx = 0;
691 for (
const auto& elem : elements(simulator_.gridView())) {
692 if (elem.partitionType() != Dune::InteriorEntity) {
700 template <
typename Flu
idState>
701 void aggregateAverageDensityContributions_(
const FluidState& fs,
702 const unsigned int globalDofIdx,
705 auto pvCellValue = RegionPhasePoreVolAverage::CellValue{};
706 pvCellValue.porv = porv;
708 for (
auto phaseIdx = 0*FluidSystem::numPhases;
709 phaseIdx < FluidSystem::numPhases; ++phaseIdx)
711 if (! FluidSystem::phaseIsActive(phaseIdx)) {
715 pvCellValue.value = getValue(fs.density(phaseIdx));
716 pvCellValue.sat = getValue(fs.saturation(phaseIdx));
719 ->addCell(globalDofIdx,
720 RegionPhasePoreVolAverage::Phase { phaseIdx },
740 data::InterRegFlowMap::FlowRates
741 getComponentSurfaceRates(
const ElementContext& elemCtx,
742 const Scalar faceArea,
743 const std::size_t scvfIdx,
744 const std::size_t timeIdx)
const
746 using Component = data::InterRegFlowMap::Component;
748 auto rates = data::InterRegFlowMap::FlowRates {};
750 const auto& extQuant = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
752 const auto alpha = getValue(extQuant.extrusionFactor()) * faceArea;
754 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
755 const auto& up = elemCtx
756 .intensiveQuantities(extQuant.upstreamIndex(oilPhaseIdx), timeIdx);
758 const auto pvtReg = up.pvtRegionIndex();
760 const auto bO = getValue(getInvB_<FluidSystem, FluidState, Scalar>
761 (up.fluidState(), oilPhaseIdx, pvtReg));
763 const auto qO = alpha * bO * getValue(extQuant.volumeFlux(oilPhaseIdx));
765 rates[Component::Oil] += qO;
767 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
768 const auto Rs = getValue(
769 BlackOil::getRs_<FluidSystem, FluidState, Scalar>
770 (up.fluidState(), pvtReg));
772 rates[Component::Gas] += qO * Rs;
773 rates[Component::Disgas] += qO * Rs;
777 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
778 const auto& up = elemCtx
779 .intensiveQuantities(extQuant.upstreamIndex(gasPhaseIdx), timeIdx);
781 const auto pvtReg = up.pvtRegionIndex();
783 const auto bG = getValue(getInvB_<FluidSystem, FluidState, Scalar>
784 (up.fluidState(), gasPhaseIdx, pvtReg));
786 const auto qG = alpha * bG * getValue(extQuant.volumeFlux(gasPhaseIdx));
788 rates[Component::Gas] += qG;
790 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
791 const auto Rv = getValue(
792 BlackOil::getRv_<FluidSystem, FluidState, Scalar>
793 (up.fluidState(), pvtReg));
795 rates[Component::Oil] += qG * Rv;
796 rates[Component::Vapoil] += qG * Rv;
800 if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
801 const auto& up = elemCtx
802 .intensiveQuantities(extQuant.upstreamIndex(waterPhaseIdx), timeIdx);
804 const auto pvtReg = up.pvtRegionIndex();
806 const auto bW = getValue(getInvB_<FluidSystem, FluidState, Scalar>
807 (up.fluidState(), waterPhaseIdx, pvtReg));
809 rates[Component::Water] +=
810 alpha * bW * getValue(extQuant.volumeFlux(waterPhaseIdx));
816 template <
typename Flu
idState>
817 Scalar hydroCarbonFraction(
const FluidState& fs)
const
819 if (this->
eclState_.runspec().co2Storage()) {
826 auto hydrocarbon = Scalar {0};
827 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
828 hydrocarbon += getValue(fs.saturation(oilPhaseIdx));
831 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
832 hydrocarbon += getValue(fs.saturation(gasPhaseIdx));
838 void updateTotalVolumesAndPressures_(
const unsigned globalDofIdx,
839 const IntensiveQuantities& intQuants,
840 const double totVolume)
842 const auto& fs = intQuants.fluidState();
844 const double pv = totVolume * intQuants.porosity().value();
845 const auto hydrocarbon = this->hydroCarbonFraction(fs);
849 totVolume * intQuants.referencePorosity());
856 !this->pressureTimesPoreVolume_.empty())
859 assert(this->
fipC_.
get(Inplace::Phase::PoreVolume).size() == this->pressureTimesPoreVolume_.size());
861 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
863 getValue(fs.pressure(oilPhaseIdx)) * pv;
868 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
870 getValue(fs.pressure(gasPhaseIdx)) * pv;
875 else if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
877 getValue(fs.pressure(waterPhaseIdx)) * pv;
882 void updatePhaseInplaceVolumes_(
const unsigned globalDofIdx,
883 const IntensiveQuantities& intQuants,
884 const double totVolume)
886 std::array<Scalar, FluidSystem::numPhases> fip {};
887 std::array<Scalar, FluidSystem::numPhases> fipr{};
889 const auto& fs = intQuants.fluidState();
890 const auto pv = totVolume * intQuants.porosity().value();
892 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
893 if (!FluidSystem::phaseIsActive(phaseIdx)) {
897 const auto b = getValue(fs.invB(phaseIdx));
898 const auto s = getValue(fs.saturation(phaseIdx));
900 fipr[phaseIdx] = s * pv;
901 fip [phaseIdx] = b * fipr[phaseIdx];
906 fs.saltConcentration().value(),
909 if (FluidSystem::phaseIsActive(oilPhaseIdx) &&
910 FluidSystem::phaseIsActive(gasPhaseIdx))
912 this->updateOilGasDistribution(globalDofIdx, fs, fip);
915 if (FluidSystem::phaseIsActive(waterPhaseIdx) &&
916 FluidSystem::phaseIsActive(gasPhaseIdx))
918 this->updateGasWaterDistribution(globalDofIdx, fs, fip);
921 if (FluidSystem::phaseIsActive(gasPhaseIdx) &&
924 this->updateCO2InGas(globalDofIdx, pv, intQuants);
928 (FluidSystem::phaseIsActive(waterPhaseIdx) ||
929 FluidSystem::phaseIsActive(oilPhaseIdx)))
931 this->updateCO2InWater(globalDofIdx, pv, fs);
934 if constexpr(enableMICP) {
935 const auto surfVolWat = pv * getValue(fs.invB(waterPhaseIdx));
937 this->updateMicrobialMass(globalDofIdx, intQuants, surfVolWat);
940 this->updateOxygenMass(globalDofIdx, intQuants, surfVolWat);
943 this->updateUreaMass(globalDofIdx, intQuants, surfVolWat);
946 this->updateBiofilmMass(globalDofIdx, intQuants, totVolume);
949 this->updateCalciteMass(globalDofIdx, intQuants, totVolume);
955 this->updateWaterMass(globalDofIdx, fs, fip);
959 template <
typename Flu
idState,
typename FIPArray>
960 void updateOilGasDistribution(
const unsigned globalDofIdx,
961 const FluidState& fs,
965 const auto gasInPlaceLiquid = getValue(fs.Rs()) * fip[oilPhaseIdx];
966 const auto oilInPlaceGas = getValue(fs.Rv()) * fip[gasPhaseIdx];
971 template <
typename Flu
idState,
typename FIPArray>
972 void updateGasWaterDistribution(
const unsigned globalDofIdx,
973 const FluidState& fs,
977 const auto gasInPlaceWater = getValue(fs.Rsw()) * fip[waterPhaseIdx];
978 const auto waterInPlaceGas = getValue(fs.Rvw()) * fip[gasPhaseIdx];
983 template <
typename IntensiveQuantities>
984 void updateCO2InGas(
const unsigned globalDofIdx,
986 const IntensiveQuantities& intQuants)
988 const auto& scaledDrainageInfo = this->simulator_.problem().materialLawManager()
989 ->oilWaterScaledEpsInfoDrainage(globalDofIdx);
991 const auto& fs = intQuants.fluidState();
992 Scalar sgcr = scaledDrainageInfo.Sgcr;
993 if (this->simulator_.problem().materialLawManager()->enableHysteresis()) {
994 const auto& matParams = simulator_.problem().materialLawParams(globalDofIdx);
995 sgcr = MaterialLaw::trappedGasSaturation(matParams,
false);
998 Scalar trappedGasSaturation = scaledDrainageInfo.Sgcr;
999 if (this->
fipC_.
has(Inplace::Phase::CO2MassInGasPhaseMaximumTrapped) ||
1000 this->fipC_.has(Inplace::Phase::CO2MassInGasPhaseMaximumUnTrapped))
1002 if (this->simulator_.problem().materialLawManager()->enableHysteresis()) {
1003 const auto& matParams = simulator_.problem().materialLawParams(globalDofIdx);
1005 trappedGasSaturation = MaterialLaw::trappedGasSaturation(matParams,
true);
1009 const Scalar sg = getValue(fs.saturation(gasPhaseIdx));
1010 Scalar strandedGasSaturation = scaledDrainageInfo.Sgcr;
1011 if (this->
fipC_.
has(Inplace::Phase::CO2MassInGasPhaseEffectiveTrapped) ||
1012 this->fipC_.has(Inplace::Phase::CO2MassInGasPhaseEffectiveUnTrapped))
1014 if (this->simulator_.problem().materialLawManager()->enableHysteresis()) {
1015 const auto& matParams = simulator_.problem().materialLawParams(globalDofIdx);
1016 const double krg = getValue(intQuants.relativePermeability(gasPhaseIdx));
1017 strandedGasSaturation = MaterialLaw::strandedGasSaturation(matParams, sg, krg);
1021 const typename FIPContainer<FluidSystem>::Co2InGasInput v{
1025 getValue(fs.density(gasPhaseIdx)),
1026 FluidSystem::phaseIsActive(waterPhaseIdx)
1027 ? FluidSystem::convertRvwToXgW(getValue(fs.Rvw()), fs.pvtRegionIndex())
1028 : FluidSystem::convertRvToXgO(getValue(fs.Rv()), fs.pvtRegionIndex()),
1029 FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex()),
1030 trappedGasSaturation,
1031 strandedGasSaturation,
1037 template <
typename Flu
idState>
1038 void updateCO2InWater(
const unsigned globalDofIdx,
1040 const FluidState& fs)
1042 const auto co2InWater = FluidSystem::phaseIsActive(oilPhaseIdx)
1043 ? this->co2InWaterFromOil(fs, pv)
1044 : this->co2InWaterFromWater(fs, pv);
1046 const Scalar mM = FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex());
1051 template <
typename Flu
idState>
1052 Scalar co2InWaterFromWater(
const FluidState& fs,
const double pv)
const
1054 const double rhow = getValue(fs.density(waterPhaseIdx));
1055 const double sw = getValue(fs.saturation(waterPhaseIdx));
1056 const double xwG = FluidSystem::convertRswToXwG(getValue(fs.Rsw()), fs.pvtRegionIndex());
1058 const Scalar mM = FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex());
1060 return xwG * pv * rhow * sw / mM;
1063 template <
typename Flu
idState>
1064 Scalar co2InWaterFromOil(
const FluidState& fs,
const double pv)
const
1066 const double rhoo = getValue(fs.density(oilPhaseIdx));
1067 const double so = getValue(fs.saturation(oilPhaseIdx));
1068 const double xoG = FluidSystem::convertRsToXoG(getValue(fs.Rs()), fs.pvtRegionIndex());
1070 const Scalar mM = FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex());
1072 return xoG * pv * rhoo * so / mM;
1075 template <
typename Flu
idState,
typename FIPArray>
1076 void updateWaterMass(
const unsigned globalDofIdx,
1077 const FluidState& fs,
1081 const Scalar rhoW = FluidSystem::referenceDensity(waterPhaseIdx, fs.pvtRegionIndex());
1086 template <
typename IntensiveQuantities>
1087 void updateMicrobialMass(
const unsigned globalDofIdx,
1088 const IntensiveQuantities& intQuants,
1089 const double surfVolWat)
1091 const Scalar mass = surfVolWat * intQuants.microbialConcentration().value();
1096 template <
typename IntensiveQuantities>
1097 void updateOxygenMass(
const unsigned globalDofIdx,
1098 const IntensiveQuantities& intQuants,
1099 const double surfVolWat)
1101 const Scalar mass = surfVolWat * intQuants.oxygenConcentration().value();
1106 template <
typename IntensiveQuantities>
1107 void updateUreaMass(
const unsigned globalDofIdx,
1108 const IntensiveQuantities& intQuants,
1109 const double surfVolWat)
1111 const Scalar mass = surfVolWat * intQuants.ureaConcentration().value();
1116 template <
typename IntensiveQuantities>
1117 void updateBiofilmMass(
const unsigned globalDofIdx,
1118 const IntensiveQuantities& intQuants,
1119 const double totVolume)
1121 const Scalar mass = totVolume * intQuants.biofilmMass().value();
1126 template <
typename IntensiveQuantities>
1127 void updateCalciteMass(
const unsigned globalDofIdx,
1128 const IntensiveQuantities& intQuants,
1129 const double totVolume)
1131 const Scalar mass = totVolume * intQuants.calciteMass().value();
1137 void setupElementExtractors_()
1139 using Entry =
typename Extractor::Entry;
1140 using Context =
typename Extractor::Context;
1141 using ScalarEntry =
typename Extractor::ScalarEntry;
1142 using PhaseEntry =
typename Extractor::PhaseEntry;
1144 const bool hasResidual = simulator_.model().linearizer().residual().size() > 0;
1145 const auto& hysteresisConfig = simulator_.problem().materialLawManager()->hysteresisConfig();
1147 auto extractors = std::array{
1149 [](
const unsigned phase,
const Context& ectx)
1150 {
return getValue(ectx.fs.saturation(phase)); }
1153 Entry{PhaseEntry{&this->
invB_,
1154 [](
const unsigned phase,
const Context& ectx)
1155 {
return getValue(ectx.fs.invB(phase)); }
1159 [](
const unsigned phase,
const Context& ectx)
1160 {
return getValue(ectx.fs.density(phase)); }
1164 [](
const unsigned phase,
const Context& ectx)
1165 {
return getValue(ectx.intQuants.relativePermeability(phase)); }
1169 [
this](
const unsigned phaseIdx,
const Context& ectx)
1171 if (this->
extboC_.allocated() && phaseIdx == oilPhaseIdx) {
1172 return getValue(ectx.intQuants.oilViscosity());
1174 else if (this->
extboC_.allocated() && phaseIdx == gasPhaseIdx) {
1175 return getValue(ectx.intQuants.gasViscosity());
1178 return getValue(ectx.fs.viscosity(phaseIdx));
1184 [&modelResid = this->simulator_.model().linearizer().residual()]
1185 (
const unsigned phaseIdx,
const Context& ectx)
1187 const unsigned sIdx = FluidSystem::solventComponentIndex(phaseIdx);
1188 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(sIdx);
1189 return modelResid[ectx.globalDofIdx][activeCompIdx];
1195 [&problem = this->simulator_.problem()](
const Context& ectx)
1197 return problem.template
1198 rockCompPoroMultiplier<Scalar>(ectx.intQuants,
1204 [&problem = this->simulator_.problem()](
const Context& ectx)
1207 template rockCompTransMultiplier<Scalar>(ectx.intQuants,
1212 [&problem = this->simulator_.problem()](
const Context& ectx)
1214 return std::min(getValue(ectx.fs.pressure(oilPhaseIdx)),
1215 problem.minOilPressure(ectx.globalDofIdx));
1221 &vanguard = this->simulator_.vanguard()](
const Context& ectx)
1225 FluidSystem::bubblePointPressure(ectx.fs,
1226 ectx.intQuants.pvtRegionIndex())
1228 }
catch (
const NumericalProblem&) {
1229 const auto cartesianIdx = vanguard.cartesianIndex(ectx.globalDofIdx);
1230 failedCells.push_back(cartesianIdx);
1238 &vanguard = this->simulator_.vanguard()](
const Context& ectx)
1242 FluidSystem::dewPointPressure(ectx.fs,
1243 ectx.intQuants.pvtRegionIndex())
1245 }
catch (
const NumericalProblem&) {
1246 const auto cartesianIdx = vanguard.cartesianIndex(ectx.globalDofIdx);
1247 failedCells.push_back(cartesianIdx);
1254 [&problem = simulator_.problem()](
const Context& ectx)
1255 {
return problem.overburdenPressure(ectx.globalDofIdx); }
1259 [](
const Context& ectx)
1260 {
return getValue(ectx.fs.temperature(oilPhaseIdx)); }
1263 Entry{ScalarEntry{&this->
sSol_,
1264 [](
const Context& ectx)
1265 {
return getValue(ectx.intQuants.solventSaturation()); }
1268 Entry{ScalarEntry{&this->
rswSol_,
1269 [](
const Context& ectx)
1270 {
return getValue(ectx.intQuants.rsSolw()); }
1274 [](
const Context& ectx)
1275 {
return getValue(ectx.intQuants.polymerConcentration()); }
1278 Entry{ScalarEntry{&this->
cFoam_,
1279 [](
const Context& ectx)
1280 {
return getValue(ectx.intQuants.foamConcentration()); }
1283 Entry{ScalarEntry{&this->
cSalt_,
1284 [](
const Context& ectx)
1285 {
return getValue(ectx.fs.saltConcentration()); }
1288 Entry{ScalarEntry{&this->
pSalt_,
1289 [](
const Context& ectx)
1290 {
return getValue(ectx.fs.saltSaturation()); }
1294 [](
const Context& ectx)
1295 {
return getValue(ectx.intQuants.permFactor()); }
1298 Entry{ScalarEntry{&this->
rPorV_,
1299 [&model = this->simulator_.model()](
const Context& ectx)
1301 const auto totVolume = model.dofTotalVolume(ectx.globalDofIdx);
1302 return totVolume * getValue(ectx.intQuants.porosity());
1306 Entry{ScalarEntry{&this->
rs_,
1307 [](
const Context& ectx)
1308 {
return getValue(ectx.fs.Rs()); }
1311 Entry{ScalarEntry{&this->
rv_,
1312 [](
const Context& ectx)
1313 {
return getValue(ectx.fs.Rv()); }
1316 Entry{ScalarEntry{&this->
rsw_,
1317 [](
const Context& ectx)
1318 {
return getValue(ectx.fs.Rsw()); }
1321 Entry{ScalarEntry{&this->
rvw_,
1322 [](
const Context& ectx)
1323 {
return getValue(ectx.fs.Rvw()); }
1326 Entry{ScalarEntry{&this->
ppcw_,
1327 [&matLawManager = *this->simulator_.problem().materialLawManager()]
1328 (
const Context& ectx)
1330 return matLawManager.
1331 oilWaterScaledEpsInfoDrainage(ectx.globalDofIdx).maxPcow;
1336 [&problem = this->simulator_.problem()](
const Context& ectx)
1338 return problem.drsdtcon(ectx.globalDofIdx,
1343 Entry{ScalarEntry{&this->
pcgw_,
1344 [](
const Context& ectx)
1346 return getValue(ectx.fs.pressure(gasPhaseIdx)) -
1347 getValue(ectx.fs.pressure(waterPhaseIdx));
1351 Entry{ScalarEntry{&this->
pcow_,
1352 [](
const Context& ectx)
1354 return getValue(ectx.fs.pressure(oilPhaseIdx)) -
1355 getValue(ectx.fs.pressure(waterPhaseIdx));
1359 Entry{ScalarEntry{&this->
pcog_,
1360 [](
const Context& ectx)
1362 return getValue(ectx.fs.pressure(gasPhaseIdx)) -
1363 getValue(ectx.fs.pressure(oilPhaseIdx));
1368 [](
const Context& ectx)
1370 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1372 return getValue(ectx.fs.pressure(oilPhaseIdx));
1374 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1376 return getValue(ectx.fs.pressure(gasPhaseIdx));
1380 return getValue(ectx.fs.pressure(waterPhaseIdx));
1386 [&problem = this->simulator_.problem()](
const Context& ectx)
1388 const Scalar SoMax = problem.maxOilSaturation(ectx.globalDofIdx);
1389 return FluidSystem::template
1390 saturatedDissolutionFactor<FluidState, Scalar>(ectx.fs,
1398 [&problem = this->simulator_.problem()](
const Context& ectx)
1400 const Scalar SoMax = problem.maxOilSaturation(ectx.globalDofIdx);
1401 return FluidSystem::template
1402 saturatedDissolutionFactor<FluidState, Scalar>(ectx.fs,
1410 [&problem = this->simulator_.problem()](
const Context& ectx)
1412 const Scalar SwMax = problem.maxWaterSaturation(ectx.globalDofIdx);
1413 return FluidSystem::template
1414 saturatedDissolutionFactor<FluidState, Scalar>(ectx.fs,
1422 [](
const Context& ectx)
1424 return FluidSystem::template
1425 saturatedVaporizationFactor<FluidState, Scalar>(ectx.fs,
1432 [](
const Context& ectx)
1434 return 1.0 / FluidSystem::template
1435 inverseFormationVolumeFactor<FluidState, Scalar>(ectx.fs,
1442 [](
const Context& ectx)
1444 return 1.0 / FluidSystem::template
1445 saturatedInverseFormationVolumeFactor<FluidState, Scalar>(ectx.fs,
1452 [](
const Context& ectx)
1454 return FluidSystem::template
1455 saturationPressure<FluidState, Scalar>(ectx.fs,
1461 Entry{ScalarEntry{&this->
soMax_,
1462 [&problem = this->simulator_.problem()](
const Context& ectx)
1464 return std::max(getValue(ectx.fs.saturation(oilPhaseIdx)),
1465 problem.maxOilSaturation(ectx.globalDofIdx));
1468 !hysteresisConfig.enableHysteresis()
1470 Entry{ScalarEntry{&this->
swMax_,
1471 [&problem = this->simulator_.problem()](
const Context& ectx)
1473 return std::max(getValue(ectx.fs.saturation(waterPhaseIdx)),
1474 problem.maxWaterSaturation(ectx.globalDofIdx));
1477 !hysteresisConfig.enableHysteresis()
1479 Entry{ScalarEntry{&this->
soMax_,
1480 [](
const Context& ectx)
1481 {
return ectx.hParams.somax; }
1483 hysteresisConfig.enableHysteresis() &&
1484 hysteresisConfig.enableNonWettingHysteresis() &&
1485 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1486 FluidSystem::phaseIsActive(waterPhaseIdx)
1488 Entry{ScalarEntry{&this->
swMax_,
1489 [](
const Context& ectx)
1490 {
return ectx.hParams.swmax; }
1492 hysteresisConfig.enableHysteresis() &&
1493 hysteresisConfig.enableWettingHysteresis() &&
1494 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1495 FluidSystem::phaseIsActive(waterPhaseIdx)
1497 Entry{ScalarEntry{&this->
swmin_,
1498 [](
const Context& ectx)
1499 {
return ectx.hParams.swmin; }
1501 hysteresisConfig.enableHysteresis() &&
1502 hysteresisConfig.enablePCHysteresis() &&
1503 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1504 FluidSystem::phaseIsActive(waterPhaseIdx)
1506 Entry{ScalarEntry{&this->
sgmax_,
1507 [](
const Context& ectx)
1508 {
return ectx.hParams.sgmax; }
1510 hysteresisConfig.enableHysteresis() &&
1511 hysteresisConfig.enableNonWettingHysteresis() &&
1512 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1513 FluidSystem::phaseIsActive(gasPhaseIdx)
1515 Entry{ScalarEntry{&this->
shmax_,
1516 [](
const Context& ectx)
1517 {
return ectx.hParams.shmax; }
1519 hysteresisConfig.enableHysteresis() &&
1520 hysteresisConfig.enableWettingHysteresis() &&
1521 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1522 FluidSystem::phaseIsActive(gasPhaseIdx)
1524 Entry{ScalarEntry{&this->
somin_,
1525 [](
const Context& ectx)
1526 {
return ectx.hParams.somin; }
1528 hysteresisConfig.enableHysteresis() &&
1529 hysteresisConfig.enablePCHysteresis() &&
1530 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1531 FluidSystem::phaseIsActive(gasPhaseIdx)
1533 Entry{[&model = this->simulator_.model(),
this](
const Context& ectx)
1537 const auto porv = ectx.intQuants.referencePorosity()
1538 * model.dofTotalVolume(ectx.globalDofIdx);
1540 this->aggregateAverageDensityContributions_(ectx.fs, ectx.globalDofIdx,
1541 static_cast<double>(porv));
1544 Entry{[&extboC = this->
extboC_](
const Context& ectx)
1546 extboC.assignVolumes(ectx.globalDofIdx,
1547 ectx.intQuants.xVolume().value(),
1548 ectx.intQuants.yVolume().value());
1549 extboC.assignZFraction(ectx.globalDofIdx,
1550 ectx.intQuants.zFraction().value());
1552 const Scalar stdVolOil = getValue(ectx.fs.saturation(oilPhaseIdx)) *
1553 getValue(ectx.fs.invB(oilPhaseIdx)) +
1554 getValue(ectx.fs.saturation(gasPhaseIdx)) *
1555 getValue(ectx.fs.invB(gasPhaseIdx)) *
1556 getValue(ectx.fs.Rv());
1557 const Scalar stdVolGas = getValue(ectx.fs.saturation(gasPhaseIdx)) *
1558 getValue(ectx.fs.invB(gasPhaseIdx)) *
1559 (1.0 - ectx.intQuants.yVolume().value()) +
1560 getValue(ectx.fs.saturation(oilPhaseIdx)) *
1561 getValue(ectx.fs.invB(oilPhaseIdx)) *
1562 getValue(ectx.fs.Rs()) *
1563 (1.0 - ectx.intQuants.xVolume().value());
1564 const Scalar stdVolCo2 = getValue(ectx.fs.saturation(gasPhaseIdx)) *
1565 getValue(ectx.fs.invB(gasPhaseIdx)) *
1566 ectx.intQuants.yVolume().value() +
1567 getValue(ectx.fs.saturation(oilPhaseIdx)) *
1568 getValue(ectx.fs.invB(oilPhaseIdx)) *
1569 getValue(ectx.fs.Rs()) *
1570 ectx.intQuants.xVolume().value();
1571 const Scalar rhoO = FluidSystem::referenceDensity(gasPhaseIdx, ectx.pvtRegionIdx);
1572 const Scalar rhoG = FluidSystem::referenceDensity(gasPhaseIdx, ectx.pvtRegionIdx);
1573 const Scalar rhoCO2 = ectx.intQuants.zRefDensity();
1574 const Scalar stdMassTotal = 1.0e-10 + stdVolOil * rhoO + stdVolGas * rhoG + stdVolCo2 * rhoCO2;
1575 extboC.assignMassFractions(ectx.globalDofIdx,
1576 stdVolGas * rhoG / stdMassTotal,
1577 stdVolOil * rhoO / stdMassTotal,
1578 stdVolCo2 * rhoCO2 / stdMassTotal);
1581 Entry{[&micpC = this->
micpC_](
const Context& ectx)
1583 micpC.assign(ectx.globalDofIdx,
1584 ectx.intQuants.microbialConcentration().value(),
1585 ectx.intQuants.oxygenConcentration().value(),
1586 ectx.intQuants.ureaConcentration().value(),
1587 ectx.intQuants.biofilmConcentration().value(),
1588 ectx.intQuants.calciteConcentration().value());
1589 }, this->
micpC_.allocated()
1591 Entry{[&rftC = this->
rftC_,
1592 &vanguard = this->simulator_.vanguard()](
const Context& ectx)
1594 const auto cartesianIdx = vanguard.cartesianIndex(ectx.globalDofIdx);
1595 rftC.assign(cartesianIdx,
1596 [&fs = ectx.fs]() {
return getValue(fs.pressure(oilPhaseIdx)); },
1597 [&fs = ectx.fs]() {
return getValue(fs.saturation(waterPhaseIdx)); },
1598 [&fs = ectx.fs]() {
return getValue(fs.saturation(gasPhaseIdx)); });
1602 &tM = this->simulator_.problem().tracerModel()](
const Context& ectx)
1604 tC.assignFreeConcentrations(ectx.globalDofIdx,
1605 [gIdx = ectx.globalDofIdx, &tM](
const unsigned tracerIdx)
1606 {
return tM.freeTracerConcentration(tracerIdx, gIdx); });
1607 tC.assignSolConcentrations(ectx.globalDofIdx,
1608 [gIdx = ectx.globalDofIdx, &tM](
const unsigned tracerIdx)
1609 {
return tM.solTracerConcentration(tracerIdx, gIdx); });
1612 Entry{[&flowsInf = this->simulator_.problem().model().linearizer().getFlowsInfo(),
1613 &flowsC = this->
flowsC_](
const Context& ectx)
1615 constexpr auto gas_idx = Indices::gasEnabled ?
1616 conti0EqIdx + Indices::canonicalToActiveComponentIndex(gasCompIdx) : -1;
1617 constexpr auto oil_idx = Indices::oilEnabled ?
1618 conti0EqIdx + Indices::canonicalToActiveComponentIndex(oilCompIdx) : -1;
1619 constexpr auto water_idx = Indices::waterEnabled ?
1620 conti0EqIdx + Indices::canonicalToActiveComponentIndex(waterCompIdx) : -1;
1621 const auto& flowsInfos = flowsInf[ectx.globalDofIdx];
1622 for (
const auto& flowsInfo : flowsInfos) {
1623 flowsC.assignFlows(ectx.globalDofIdx,
1626 value_or_zero<gas_idx>(flowsInfo.flow),
1627 value_or_zero<oil_idx>(flowsInfo.flow),
1628 value_or_zero<water_idx>(flowsInfo.flow));
1630 }, !this->simulator_.problem().model().linearizer().getFlowsInfo().empty()
1632 Entry{[&floresInf = this->simulator_.problem().model().linearizer().getFloresInfo(),
1633 &flowsC = this->
flowsC_](
const Context& ectx)
1635 constexpr auto gas_idx = Indices::gasEnabled ?
1636 conti0EqIdx + Indices::canonicalToActiveComponentIndex(gasCompIdx) : -1;
1637 constexpr auto oil_idx = Indices::oilEnabled ?
1638 conti0EqIdx + Indices::canonicalToActiveComponentIndex(oilCompIdx) : -1;
1639 constexpr auto water_idx = Indices::waterEnabled ?
1640 conti0EqIdx + Indices::canonicalToActiveComponentIndex(waterCompIdx) : -1;
1641 const auto& floresInfos = floresInf[ectx.globalDofIdx];
1642 for (
const auto& floresInfo : floresInfos) {
1643 flowsC.assignFlores(ectx.globalDofIdx,
1646 value_or_zero<gas_idx>(floresInfo.flow),
1647 value_or_zero<oil_idx>(floresInfo.flow),
1648 value_or_zero<water_idx>(floresInfo.flow));
1650 }, !this->simulator_.problem().model().linearizer().getFloresInfo().empty()
1657 Entry{ScalarEntry{&this->
rv_,
1658 [&problem = this->simulator_.problem()](
const Context& ectx)
1659 {
return problem.initialFluidState(ectx.globalDofIdx).Rv(); }
1661 simulator_.episodeIndex() < 0 &&
1662 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1663 FluidSystem::phaseIsActive(gasPhaseIdx)
1665 Entry{ScalarEntry{&this->
rs_,
1666 [&problem = this->simulator_.problem()](
const Context& ectx)
1667 {
return problem.initialFluidState(ectx.globalDofIdx).Rs(); }
1669 simulator_.episodeIndex() < 0 &&
1670 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1671 FluidSystem::phaseIsActive(gasPhaseIdx)
1673 Entry{ScalarEntry{&this->
rsw_,
1674 [&problem = this->simulator_.problem()](
const Context& ectx)
1675 {
return problem.initialFluidState(ectx.globalDofIdx).Rsw(); }
1677 simulator_.episodeIndex() < 0 &&
1678 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1679 FluidSystem::phaseIsActive(gasPhaseIdx)
1681 Entry{ScalarEntry{&this->
rvw_,
1682 [&problem = this->simulator_.problem()](
const Context& ectx)
1683 {
return problem.initialFluidState(ectx.globalDofIdx).Rvw(); }
1685 simulator_.episodeIndex() < 0 &&
1686 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1687 FluidSystem::phaseIsActive(gasPhaseIdx)
1691 [&problem = this->simulator_.problem()](
const unsigned phase,
1692 const Context& ectx)
1694 const auto& fsInitial = problem.initialFluidState(ectx.globalDofIdx);
1695 return FluidSystem::density(fsInitial,
1697 ectx.intQuants.pvtRegionIndex());
1700 simulator_.episodeIndex() < 0 &&
1701 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1702 FluidSystem::phaseIsActive(gasPhaseIdx)
1704 Entry{PhaseEntry{&this->
invB_,
1705 [&problem = this->simulator_.problem()](
const unsigned phase,
1706 const Context& ectx)
1708 const auto& fsInitial = problem.initialFluidState(ectx.globalDofIdx);
1709 return FluidSystem::inverseFormationVolumeFactor(fsInitial,
1711 ectx.intQuants.pvtRegionIndex());
1714 simulator_.episodeIndex() < 0 &&
1715 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1716 FluidSystem::phaseIsActive(gasPhaseIdx)
1719 [&problem = this->simulator_.problem()](
const unsigned phase,
1720 const Context& ectx)
1722 const auto& fsInitial = problem.initialFluidState(ectx.globalDofIdx);
1723 return FluidSystem::viscosity(fsInitial,
1725 ectx.intQuants.pvtRegionIndex());
1728 simulator_.episodeIndex() < 0 &&
1729 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1730 FluidSystem::phaseIsActive(gasPhaseIdx)
1737 if constexpr (getPropValue<TypeTag, Properties::EnableMech>()) {
1738 if (this->
mech_.allocated()) {
1739 this->extractors_.push_back(
1740 Entry{[&mech = this->
mech_,
1741 &model = simulator_.problem().geoMechModel()](
const Context& ectx)
1743 mech.assignDelStress(ectx.globalDofIdx,
1744 model.delstress(ectx.globalDofIdx));
1746 mech.assignDisplacement(ectx.globalDofIdx,
1747 model.disp(ectx.globalDofIdx,
true));
1750 mech.assignFracStress(ectx.globalDofIdx,
1751 model.fractureStress(ectx.globalDofIdx));
1753 mech.assignLinStress(ectx.globalDofIdx,
1754 model.linstress(ectx.globalDofIdx));
1756 mech.assignPotentialForces(ectx.globalDofIdx,
1757 model.mechPotentialForce(ectx.globalDofIdx),
1758 model.mechPotentialPressForce(ectx.globalDofIdx),
1759 model.mechPotentialTempForce(ectx.globalDofIdx));
1761 mech.assignStrain(ectx.globalDofIdx,
1762 model.strain(ectx.globalDofIdx,
true));
1765 mech.assignStress(ectx.globalDofIdx,
1766 model.stress(ectx.globalDofIdx,
true));
1775 void setupBlockExtractors_(
const bool isSubStep,
1776 const int reportStepNum)
1779 using Context =
typename BlockExtractor::Context;
1780 using PhaseEntry =
typename BlockExtractor::PhaseEntry;
1781 using ScalarEntry =
typename BlockExtractor::ScalarEntry;
1783 using namespace std::string_view_literals;
1785 const auto pressure_handler =
1786 Entry{ScalarEntry{std::vector{
"BPR"sv,
"BPRESSUR"sv},
1787 [](
const Context& ectx)
1789 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1790 return getValue(ectx.fs.pressure(oilPhaseIdx));
1792 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1793 return getValue(ectx.fs.pressure(gasPhaseIdx));
1796 return getValue(ectx.fs.pressure(waterPhaseIdx));
1802 const auto handlers = std::array{
1804 Entry{PhaseEntry{std::array{
1805 std::array{
"BWSAT"sv,
"BOSAT"sv,
"BGSAT"sv},
1806 std::array{
"BSWAT"sv,
"BSOIL"sv,
"BSGAS"sv}
1808 [](
const unsigned phaseIdx,
const Context& ectx)
1810 return getValue(ectx.fs.saturation(phaseIdx));
1814 Entry{ScalarEntry{
"BNSAT",
1815 [](
const Context& ectx)
1817 return ectx.intQuants.solventSaturation().value();
1821 Entry{ScalarEntry{std::vector{
"BTCNFHEA"sv,
"BTEMP"sv},
1822 [](
const Context& ectx)
1824 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1825 return getValue(ectx.fs.temperature(oilPhaseIdx));
1827 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1828 return getValue(ectx.fs.temperature(gasPhaseIdx));
1831 return getValue(ectx.fs.temperature(waterPhaseIdx));
1836 Entry{PhaseEntry{std::array{
1837 std::array{
"BWKR"sv,
"BOKR"sv,
"BGKR"sv},
1838 std::array{
"BKRW"sv,
"BKRO"sv,
"BKRG"sv}
1840 [](
const unsigned phaseIdx,
const Context& ectx)
1842 return getValue(ectx.intQuants.relativePermeability(phaseIdx));
1846 Entry{ScalarEntry{
"BKROG",
1847 [&problem = this->simulator_.problem()](
const Context& ectx)
1849 const auto& materialParams =
1850 problem.materialLawParams(ectx.elemCtx,
1853 return getValue(MaterialLaw::template
1854 relpermOilInOilGasSystem<Evaluation>(materialParams,
1859 Entry{ScalarEntry{
"BKROW",
1860 [&problem = this->simulator_.problem()](
const Context& ectx)
1862 const auto& materialParams = problem.materialLawParams(ectx.elemCtx,
1865 return getValue(MaterialLaw::template
1866 relpermOilInOilWaterSystem<Evaluation>(materialParams,
1871 Entry{ScalarEntry{
"BWPC",
1872 [](
const Context& ectx)
1874 return getValue(ectx.fs.pressure(oilPhaseIdx)) -
1875 getValue(ectx.fs.pressure(waterPhaseIdx));
1879 Entry{ScalarEntry{
"BGPC",
1880 [](
const Context& ectx)
1882 return getValue(ectx.fs.pressure(gasPhaseIdx)) -
1883 getValue(ectx.fs.pressure(oilPhaseIdx));
1887 Entry{ScalarEntry{
"BWPR",
1888 [](
const Context& ectx)
1890 return getValue(ectx.fs.pressure(waterPhaseIdx));
1894 Entry{ScalarEntry{
"BGPR",
1895 [](
const Context& ectx)
1897 return getValue(ectx.fs.pressure(gasPhaseIdx));
1901 Entry{PhaseEntry{std::array{
1902 std::array{
"BVWAT"sv,
"BVOIL"sv,
"BVGAS"sv},
1903 std::array{
"BWVIS"sv,
"BOVIS"sv,
"BGVIS"sv}
1905 [](
const unsigned phaseIdx,
const Context& ectx)
1907 return getValue(ectx.fs.viscosity(phaseIdx));
1911 Entry{PhaseEntry{std::array{
1912 std::array{
"BWDEN"sv,
"BODEN"sv,
"BGDEN"sv},
1913 std::array{
"BDENW"sv,
"BDENO"sv,
"BDENG"sv}
1915 [](
const unsigned phaseIdx,
const Context& ectx)
1917 return getValue(ectx.fs.density(phaseIdx));
1921 Entry{ScalarEntry{
"BFLOWI",
1922 [&flowsC = this->
flowsC_](
const Context& ectx)
1924 return flowsC.getFlow(ectx.globalDofIdx, Dir::XPlus, waterCompIdx);
1928 Entry{ScalarEntry{
"BFLOWJ",
1929 [&flowsC = this->
flowsC_](
const Context& ectx)
1931 return flowsC.getFlow(ectx.globalDofIdx, Dir::YPlus, waterCompIdx);
1935 Entry{ScalarEntry{
"BFLOWK",
1936 [&flowsC = this->
flowsC_](
const Context& ectx)
1938 return flowsC.getFlow(ectx.globalDofIdx, Dir::ZPlus, waterCompIdx);
1942 Entry{ScalarEntry{
"BRPV",
1943 [&model = this->simulator_.model()](
const Context& ectx)
1945 return getValue(ectx.intQuants.porosity()) *
1946 model.dofTotalVolume(ectx.globalDofIdx);
1950 Entry{PhaseEntry{std::array{
"BWPV"sv,
"BOPV"sv,
"BGPV"sv},
1951 [&model = this->simulator_.model()](
const unsigned phaseIdx,
1952 const Context& ectx)
1954 return getValue(ectx.fs.saturation(phaseIdx)) *
1955 getValue(ectx.intQuants.porosity()) *
1956 model.dofTotalVolume(ectx.globalDofIdx);
1960 Entry{ScalarEntry{
"BRS",
1961 [](
const Context& ectx)
1963 return getValue(ectx.fs.Rs());
1967 Entry{ScalarEntry{
"BRV",
1968 [](
const Context& ectx)
1970 return getValue(ectx.fs.Rv());
1974 Entry{ScalarEntry{
"BOIP",
1975 [&model = this->simulator_.model()](
const Context& ectx)
1977 return (getValue(ectx.fs.invB(oilPhaseIdx)) *
1978 getValue(ectx.fs.saturation(oilPhaseIdx)) +
1979 getValue(ectx.fs.Rv()) *
1980 getValue(ectx.fs.invB(gasPhaseIdx)) *
1981 getValue(ectx.fs.saturation(gasPhaseIdx))) *
1982 model.dofTotalVolume(ectx.globalDofIdx) *
1983 getValue(ectx.intQuants.porosity());
1987 Entry{ScalarEntry{
"BGIP",
1988 [&model = this->simulator_.model()](
const Context& ectx)
1990 Scalar result = getValue(ectx.fs.invB(gasPhaseIdx)) *
1991 getValue(ectx.fs.saturation(gasPhaseIdx));
1993 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1994 result += getValue(ectx.fs.Rs()) *
1995 getValue(ectx.fs.invB(oilPhaseIdx)) *
1996 getValue(ectx.fs.saturation(oilPhaseIdx));
1999 result += getValue(ectx.fs.Rsw()) *
2000 getValue(ectx.fs.invB(waterPhaseIdx)) *
2001 getValue(ectx.fs.saturation(waterPhaseIdx));
2005 model.dofTotalVolume(ectx.globalDofIdx) *
2006 getValue(ectx.intQuants.porosity());
2010 Entry{ScalarEntry{
"BWIP",
2011 [&model = this->simulator_.model()](
const Context& ectx)
2013 return getValue(ectx.fs.invB(waterPhaseIdx)) *
2014 getValue(ectx.fs.saturation(waterPhaseIdx)) *
2015 model.dofTotalVolume(ectx.globalDofIdx) *
2016 getValue(ectx.intQuants.porosity());
2020 Entry{ScalarEntry{
"BOIPL",
2021 [&model = this->simulator_.model()](
const Context& ectx)
2023 return getValue(ectx.fs.invB(oilPhaseIdx)) *
2024 getValue(ectx.fs.saturation(oilPhaseIdx)) *
2025 model.dofTotalVolume(ectx.globalDofIdx) *
2026 getValue(ectx.intQuants.porosity());
2030 Entry{ScalarEntry{
"BGIPL",
2031 [&model = this->simulator_.model()](
const Context& ectx)
2034 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
2035 result = getValue(ectx.fs.Rs()) *
2036 getValue(ectx.fs.invB(oilPhaseIdx)) *
2037 getValue(ectx.fs.saturation(oilPhaseIdx));
2040 result = getValue(ectx.fs.Rsw()) *
2041 getValue(ectx.fs.invB(waterPhaseIdx)) *
2042 getValue(ectx.fs.saturation(waterPhaseIdx));
2045 model.dofTotalVolume(ectx.globalDofIdx) *
2046 getValue(ectx.intQuants.porosity());
2050 Entry{ScalarEntry{
"BGIPG",
2051 [&model = this->simulator_.model()](
const Context& ectx)
2053 return getValue(ectx.fs.invB(gasPhaseIdx)) *
2054 getValue(ectx.fs.saturation(gasPhaseIdx)) *
2055 model.dofTotalVolume(ectx.globalDofIdx) *
2056 getValue(ectx.intQuants.porosity());
2060 Entry{ScalarEntry{
"BOIPG",
2061 [&model = this->simulator_.model()](
const Context& ectx)
2063 return getValue(ectx.fs.Rv()) *
2064 getValue(ectx.fs.invB(gasPhaseIdx)) *
2065 getValue(ectx.fs.saturation(gasPhaseIdx)) *
2066 model.dofTotalVolume(ectx.globalDofIdx) *
2067 getValue(ectx.intQuants.porosity());
2071 Entry{PhaseEntry{std::array{
"BPPW"sv,
"BPPO"sv,
"BPPG"sv},
2072 [&simConfig = this->
eclState_.getSimulationConfig(),
2073 &grav = this->simulator_.problem().gravity(),
2075 &problem = this->simulator_.problem(),
2076 ®ions = this->
regions_](
const unsigned phaseIdx,
const Context& ectx)
2078 auto phase = RegionPhasePoreVolAverage::Phase{};
2079 phase.ix = phaseIdx;
2088 const auto datum = simConfig.datumDepths()(regions[
"FIPNUM"][ectx.dofIdx] - 1);
2091 const auto region = RegionPhasePoreVolAverage::Region {
2092 ectx.elemCtx.primaryVars(ectx.dofIdx, 0).pvtRegionIndex() + 1
2095 const auto density = regionAvgDensity->value(
"PVTNUM", phase, region);
2097 const auto press = getValue(ectx.fs.pressure(phase.ix));
2098 const auto dz = problem.dofCenterDepth(ectx.globalDofIdx) - datum;
2099 return press - density*dz*grav[GridView::dimensionworld - 1];
2103 Entry{ScalarEntry{
"BAMIP",
2104 [&model = this->simulator_.model()](
const Context& ectx)
2106 const Scalar rhoW = FluidSystem::referenceDensity(waterPhaseIdx,
2107 ectx.intQuants.pvtRegionIndex());
2108 return getValue(ectx.fs.invB(waterPhaseIdx)) *
2109 getValue(ectx.fs.saturation(waterPhaseIdx)) *
2111 model.dofTotalVolume(ectx.globalDofIdx) *
2112 getValue(ectx.intQuants.porosity());
2116 Entry{ScalarEntry{
"BMMIP",
2117 [&model = this->simulator_.model()](
const Context& ectx)
2119 return getValue(ectx.intQuants.microbialConcentration()) *
2120 getValue(ectx.intQuants.porosity()) *
2121 model.dofTotalVolume(ectx.globalDofIdx);
2125 Entry{ScalarEntry{
"BMOIP",
2126 [&model = this->simulator_.model()](
const Context& ectx)
2128 return getValue(ectx.intQuants.oxygenConcentration()) *
2129 getValue(ectx.intQuants.porosity()) *
2130 model.dofTotalVolume(ectx.globalDofIdx);
2134 Entry{ScalarEntry{
"BMUIP",
2135 [&model = this->simulator_.model()](
const Context& ectx)
2137 return getValue(ectx.intQuants.ureaConcentration()) *
2138 getValue(ectx.intQuants.porosity()) *
2139 model.dofTotalVolume(ectx.globalDofIdx) * 1;
2143 Entry{ScalarEntry{
"BMBIP",
2144 [&model = this->simulator_.model()](
const Context& ectx)
2146 return model.dofTotalVolume(ectx.globalDofIdx) *
2147 getValue(ectx.intQuants.biofilmMass());
2151 Entry{ScalarEntry{
"BMCIP",
2152 [&model = this->simulator_.model()](
const Context& ectx)
2154 return model.dofTotalVolume(ectx.globalDofIdx) *
2155 getValue(ectx.intQuants.calciteMass());
2159 Entry{ScalarEntry{
"BGMIP",
2160 [&model = this->simulator_.model()](
const Context& ectx)
2162 Scalar result = getValue(ectx.fs.invB(gasPhaseIdx)) *
2163 getValue(ectx.fs.saturation(gasPhaseIdx));
2165 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
2166 result += getValue(ectx.fs.Rs()) *
2167 getValue(ectx.fs.invB(oilPhaseIdx)) *
2168 getValue(ectx.fs.saturation(oilPhaseIdx));
2171 result += getValue(ectx.fs.Rsw()) *
2172 getValue(ectx.fs.invB(waterPhaseIdx)) *
2173 getValue(ectx.fs.saturation(waterPhaseIdx));
2175 const Scalar rhoG = FluidSystem::referenceDensity(gasPhaseIdx,
2176 ectx.intQuants.pvtRegionIndex());
2178 model.dofTotalVolume(ectx.globalDofIdx) *
2179 getValue(ectx.intQuants.porosity()) *
2184 Entry{ScalarEntry{
"BGMGP",
2185 [&model = this->simulator_.model()](
const Context& ectx)
2187 const Scalar rhoG = FluidSystem::referenceDensity(gasPhaseIdx,
2188 ectx.intQuants.pvtRegionIndex());
2189 return getValue(ectx.fs.invB(gasPhaseIdx)) *
2190 getValue(ectx.fs.saturation(gasPhaseIdx)) *
2191 model.dofTotalVolume(ectx.globalDofIdx) *
2192 getValue(ectx.intQuants.porosity()) *
2197 Entry{ScalarEntry{
"BGMDS",
2198 [&model = this->simulator_.model()](
const Context& ectx)
2201 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
2202 result = getValue(ectx.fs.Rs()) *
2203 getValue(ectx.fs.invB(oilPhaseIdx)) *
2204 getValue(ectx.fs.saturation(oilPhaseIdx));
2207 result = getValue(ectx.fs.Rsw()) *
2208 getValue(ectx.fs.invB(waterPhaseIdx)) *
2209 getValue(ectx.fs.saturation(waterPhaseIdx));
2211 const Scalar rhoG = FluidSystem::referenceDensity(gasPhaseIdx,
2212 ectx.intQuants.pvtRegionIndex());
2214 model.dofTotalVolume(ectx.globalDofIdx) *
2215 getValue(ectx.intQuants.porosity()) *
2220 Entry{ScalarEntry{
"BGMST",
2221 [&model = this->simulator_.model(),
2222 &problem = this->simulator_.problem()](
const Context& ectx)
2224 const auto& scaledDrainageInfo = problem.materialLawManager()
2225 ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2226 const Scalar sg = getValue(ectx.fs.saturation(gasPhaseIdx));
2227 Scalar strandedGas = scaledDrainageInfo.Sgcr;
2228 if (problem.materialLawManager()->enableHysteresis()) {
2229 const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2230 const Scalar krg = getValue(ectx.intQuants.relativePermeability(gasPhaseIdx));
2231 strandedGas = MaterialLaw::strandedGasSaturation(matParams, sg, krg);
2233 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2234 FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2235 : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2236 return (1.0 - xgW) *
2237 model.dofTotalVolume(ectx.globalDofIdx) *
2238 getValue(ectx.intQuants.porosity()) *
2239 getValue(ectx.fs.density(gasPhaseIdx)) *
2240 std::min(strandedGas, sg);
2244 Entry{ScalarEntry{
"BGMUS",
2245 [&model = this->simulator_.model(),
2246 &problem = this->simulator_.problem()](
const Context& ectx)
2248 const auto& scaledDrainageInfo = problem.materialLawManager()
2249 ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2250 const Scalar sg = getValue(ectx.fs.saturation(gasPhaseIdx));
2251 Scalar strandedGas = scaledDrainageInfo.Sgcr;
2252 if (problem.materialLawManager()->enableHysteresis()) {
2253 const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2254 const Scalar krg = getValue(ectx.intQuants.relativePermeability(gasPhaseIdx));
2255 strandedGas = MaterialLaw::strandedGasSaturation(matParams, sg, krg);
2257 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2258 FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2259 : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2260 return (1.0 - xgW) *
2261 model.dofTotalVolume(ectx.globalDofIdx) *
2262 getValue(ectx.intQuants.porosity()) *
2263 getValue(ectx.fs.density(gasPhaseIdx)) *
2264 std::max(Scalar{0.0}, sg - strandedGas);
2268 Entry{ScalarEntry{
"BGMTR",
2269 [&model = this->simulator_.model(),
2270 &problem = this->simulator_.problem()](
const Context& ectx)
2272 const auto& scaledDrainageInfo = problem.materialLawManager()
2273 ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2274 Scalar trappedGas = scaledDrainageInfo.Sgcr;
2275 if (problem.materialLawManager()->enableHysteresis()) {
2276 const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2277 trappedGas = MaterialLaw::trappedGasSaturation(matParams,
true);
2279 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2280 FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2281 : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2282 return (1.0 - xgW) *
2283 model.dofTotalVolume(ectx.globalDofIdx) *
2284 getValue(ectx.intQuants.porosity()) *
2285 getValue(ectx.fs.density(gasPhaseIdx)) *
2286 std::min(trappedGas, getValue(ectx.fs.saturation(gasPhaseIdx)));
2290 Entry{ScalarEntry{
"BGMMO",
2291 [&model = this->simulator_.model(),
2292 &problem = this->simulator_.problem()](
const Context& ectx)
2294 const auto& scaledDrainageInfo = problem.materialLawManager()
2295 ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2296 Scalar trappedGas = scaledDrainageInfo.Sgcr;
2297 if (problem.materialLawManager()->enableHysteresis()) {
2298 const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2299 trappedGas = MaterialLaw::trappedGasSaturation(matParams,
true);
2301 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2302 FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2303 : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2304 return (1.0 - xgW) *
2305 model.dofTotalVolume(ectx.globalDofIdx) *
2306 getValue(ectx.intQuants.porosity()) *
2307 getValue(ectx.fs.density(gasPhaseIdx)) *
2308 std::max(Scalar{0.0}, getValue(ectx.fs.saturation(gasPhaseIdx)) - trappedGas);
2312 Entry{ScalarEntry{
"BGKTR",
2313 [&model = this->simulator_.model(),
2314 &problem = this->simulator_.problem()](
const Context& ectx)
2316 const auto& scaledDrainageInfo = problem.materialLawManager()
2317 ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2318 const Scalar sg = getValue(ectx.fs.saturation(gasPhaseIdx));
2319 Scalar sgcr = scaledDrainageInfo.Sgcr;
2320 if (problem.materialLawManager()->enableHysteresis()) {
2321 const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2322 sgcr = MaterialLaw::trappedGasSaturation(matParams,
false);
2328 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2329 FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2330 : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2331 return (1.0 - xgW) *
2332 model.dofTotalVolume(ectx.globalDofIdx) *
2333 getValue(ectx.intQuants.porosity()) *
2334 getValue(ectx.fs.density(gasPhaseIdx)) *
2335 getValue(ectx.fs.saturation(gasPhaseIdx));
2340 Entry{ScalarEntry{
"BGKMO",
2341 [&model = this->simulator_.model(),
2342 &problem = this->simulator_.problem()](
const Context& ectx)
2344 const auto& scaledDrainageInfo = problem.materialLawManager()
2345 ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2346 const Scalar sg = getValue(ectx.fs.saturation(gasPhaseIdx));
2347 Scalar sgcr = scaledDrainageInfo.Sgcr;
2348 if (problem.materialLawManager()->enableHysteresis()) {
2349 const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2350 sgcr = MaterialLaw::trappedGasSaturation(matParams,
false);
2356 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2357 FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2358 : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2359 return (1.0 - xgW) *
2360 model.dofTotalVolume(ectx.globalDofIdx) *
2361 getValue(ectx.intQuants.porosity()) *
2362 getValue(ectx.fs.density(gasPhaseIdx)) *
2363 getValue(ectx.fs.saturation(gasPhaseIdx));
2368 Entry{ScalarEntry{
"BGCDI",
2369 [&model = this->simulator_.model(),
2370 &problem = this->simulator_.problem()](
const Context& ectx)
2372 const auto& scaledDrainageInfo = problem.materialLawManager()
2373 ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2374 Scalar sgcr = scaledDrainageInfo.Sgcr;
2375 if (problem.materialLawManager()->enableHysteresis()) {
2376 const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2377 sgcr = MaterialLaw::trappedGasSaturation(matParams,
false);
2379 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2380 FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2381 : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2382 return (1.0 - xgW) *
2383 model.dofTotalVolume(ectx.globalDofIdx) *
2384 getValue(ectx.intQuants.porosity()) *
2385 getValue(ectx.fs.density(gasPhaseIdx)) *
2386 std::min(sgcr, getValue(ectx.fs.saturation(gasPhaseIdx))) /
2387 FluidSystem::molarMass(gasCompIdx, ectx.intQuants.pvtRegionIndex());
2391 Entry{ScalarEntry{
"BGCDM",
2392 [&model = this->simulator_.model(),
2393 &problem = this->simulator_.problem()](
const Context& ectx)
2395 const auto& scaledDrainageInfo = problem.materialLawManager()
2396 ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2397 Scalar sgcr = scaledDrainageInfo.Sgcr;
2398 if (problem.materialLawManager()->enableHysteresis()) {
2399 const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2400 sgcr = MaterialLaw::trappedGasSaturation(matParams,
false);
2402 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2403 FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2404 : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2405 return (1.0 - xgW) *
2406 model.dofTotalVolume(ectx.globalDofIdx) *
2407 getValue(ectx.intQuants.porosity()) *
2408 getValue(ectx.fs.density(gasPhaseIdx)) *
2409 std::max(Scalar{0.0}, getValue(ectx.fs.saturation(gasPhaseIdx)) - sgcr) /
2410 FluidSystem::molarMass(gasCompIdx, ectx.intQuants.pvtRegionIndex());
2414 Entry{ScalarEntry{
"BGKDI",
2415 [&model = this->simulator_.model(),
2416 &problem = this->simulator_.problem()](
const Context& ectx)
2418 const auto& scaledDrainageInfo = problem.materialLawManager()
2419 ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2420 const Scalar sg = getValue(ectx.fs.saturation(gasPhaseIdx));
2421 Scalar sgcr = scaledDrainageInfo.Sgcr;
2422 if (problem.materialLawManager()->enableHysteresis()) {
2423 const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2424 sgcr = MaterialLaw::trappedGasSaturation(matParams,
false);
2430 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2431 FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2432 : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2433 return (1.0 - xgW) *
2434 model.dofTotalVolume(ectx.globalDofIdx) *
2435 getValue(ectx.intQuants.porosity()) *
2436 getValue(ectx.fs.density(gasPhaseIdx)) *
2437 getValue(ectx.fs.saturation(gasPhaseIdx)) /
2438 FluidSystem::molarMass(gasCompIdx, ectx.intQuants.pvtRegionIndex());
2443 Entry{ScalarEntry{
"BGKDM",
2444 [&model = this->simulator_.model(),
2445 &problem = this->simulator_.problem()](
const Context& ectx)
2447 const auto& scaledDrainageInfo = problem.materialLawManager()
2448 ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2449 const Scalar sg = getValue(ectx.fs.saturation(gasPhaseIdx));
2450 Scalar sgcr = scaledDrainageInfo.Sgcr;
2451 if (problem.materialLawManager()->enableHysteresis()) {
2452 const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2453 sgcr = MaterialLaw::trappedGasSaturation(matParams,
false);
2459 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2460 FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2461 : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2462 return (1.0 - xgW) *
2463 model.dofTotalVolume(ectx.globalDofIdx) *
2464 getValue(ectx.intQuants.porosity()) *
2465 getValue(ectx.fs.density(gasPhaseIdx)) *
2466 getValue(ectx.fs.saturation(gasPhaseIdx)) /
2467 FluidSystem::molarMass(gasCompIdx, ectx.intQuants.pvtRegionIndex());
2472 Entry{ScalarEntry{
"BWCD",
2473 [&model = this->simulator_.model()](
const Context& ectx)
2476 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
2477 result = getValue(ectx.fs.Rs()) *
2478 getValue(ectx.fs.invB(oilPhaseIdx)) *
2479 getValue(ectx.fs.saturation(oilPhaseIdx));
2482 result = getValue(ectx.fs.Rsw()) *
2483 getValue(ectx.fs.invB(waterPhaseIdx)) *
2484 getValue(ectx.fs.saturation(waterPhaseIdx));
2486 const Scalar rhoG = FluidSystem::referenceDensity(gasPhaseIdx,
2487 ectx.intQuants.pvtRegionIndex());
2489 model.dofTotalVolume(ectx.globalDofIdx) *
2490 getValue(ectx.intQuants.porosity()) *
2492 FluidSystem::molarMass(gasCompIdx, ectx.intQuants.pvtRegionIndex());
2496 Entry{ScalarEntry{
"BWIPG",
2497 [&model = this->simulator_.model()](
const Context& ectx)
2499 Scalar result = 0.0;
2500 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
2501 result = getValue(ectx.fs.Rvw()) *
2502 getValue(ectx.fs.invB(gasPhaseIdx)) *
2503 getValue(ectx.fs.saturation(gasPhaseIdx));
2506 model.dofTotalVolume(ectx.globalDofIdx) *
2507 getValue(ectx.intQuants.porosity());
2511 Entry{ScalarEntry{
"BWIPL",
2512 [&model = this->simulator_.model()](
const Context& ectx)
2514 return getValue(ectx.fs.invB(waterPhaseIdx)) *
2515 getValue(ectx.fs.saturation(waterPhaseIdx)) *
2516 model.dofTotalVolume(ectx.globalDofIdx) *
2517 getValue(ectx.intQuants.porosity());
2526 if (reportStepNum > 0 && !isSubStep) {
2528 const auto& rpt = this->
schedule_[reportStepNum - 1].rpt_config.get();
2529 if (rpt.contains(
"WELLS") && rpt.at(
"WELLS") > 1) {
2531 [&c = this->collectOnIORank_](
const int idx)
2532 {
return c.isCartIdxOnThisRank(idx); });
2534 const auto extraHandlers = std::array{
2543 const Simulator& simulator_;
2544 const CollectDataOnIORankType& collectOnIORank_;
2545 std::vector<typename Extractor::Entry> extractors_;
Declares the properties required by the black oil model.
Definition: CollectDataOnIORank.hpp:56
The base class for the element-centered finite-volume discretization scheme.
Definition: ecfvdiscretization.hh:160
void assignMicrobialMass(const unsigned globalDofIdx, const Scalar microbialMass)
void assignCalciteMass(const unsigned globalDofIdx, const Scalar calciteMass)
void assignCo2InWater(const unsigned globalDofIdx, const Scalar co2InWater, const Scalar mM)
void assignVolumesSurface(const unsigned globalDofIdx, const std::array< Scalar, numPhases > &fip)
bool has(const Inplace::Phase phase) const
bool hasMicrobialMass() const
void assignWaterMass(const unsigned globalDofIdx, const std::array< Scalar, numPhases > &fip, const Scalar rhoW)
void assignCo2InGas(const unsigned globalDofIdx, const Co2InGasInput &v)
bool hasOxygenMass() const
void assignVolumesReservoir(const unsigned globalDofIdx, const Scalar saltConcentration, const std::array< Scalar, numPhases > &fipr)
void assignPoreVolume(const unsigned globalDofIdx, const Scalar value)
void assignOxygenMass(const unsigned globalDofIdx, const Scalar oxygenMass)
void assignOilGasDistribution(const unsigned globalDofIdx, const Scalar gasInPlaceLiquid, const Scalar oilInPlaceGas)
void assignBiofilmMass(const unsigned globalDofIdx, const Scalar biofilmMass)
bool hasWaterMass() const
bool hasCo2InWater() const
void assignUreaMass(const unsigned globalDofIdx, const Scalar ureaMass)
bool hasCalciteMass() const
bool hasBiofilmMass() const
const std::vector< Scalar > & get(const Inplace::Phase phase) const
void assignGasWater(const unsigned globalDofIdx, const std::array< Scalar, numPhases > &fip, const Scalar gasInPlaceWater, const Scalar waterInPlaceGas)
Definition: GenericOutputBlackoilModule.hpp:76
std::map< std::pair< std::string, int >, double > blockData_
Definition: GenericOutputBlackoilModule.hpp:434
std::array< ScalarBuffer, numPhases > relativePermeability_
Definition: GenericOutputBlackoilModule.hpp:423
ScalarBuffer fluidPressure_
Definition: GenericOutputBlackoilModule.hpp:377
std::array< ScalarBuffer, numPhases > density_
Definition: GenericOutputBlackoilModule.hpp:421
ScalarBuffer saturatedOilFormationVolumeFactor_
Definition: GenericOutputBlackoilModule.hpp:409
ScalarBuffer overburdenPressure_
Definition: GenericOutputBlackoilModule.hpp:383
ScalarBuffer gasDissolutionFactorInWater_
Definition: GenericOutputBlackoilModule.hpp:403
const EclipseState & eclState_
Definition: GenericOutputBlackoilModule.hpp:336
ScalarBuffer swmin_
Definition: GenericOutputBlackoilModule.hpp:399
ScalarBuffer rockCompPorvMultiplier_
Definition: GenericOutputBlackoilModule.hpp:407
MICPContainer< Scalar > micpC_
Definition: GenericOutputBlackoilModule.hpp:411
RFTContainer< GetPropType< TypeTag, Properties::FluidSystem > > rftC_
Definition: GenericOutputBlackoilModule.hpp:431
bool computeFip_
Definition: GenericOutputBlackoilModule.hpp:359
ScalarBuffer dewPointPressure_
Definition: GenericOutputBlackoilModule.hpp:406
LogOutputHelper< Scalar > logOutput_
Definition: GenericOutputBlackoilModule.hpp:343
std::vector< int > failedCellsPb_
Definition: GenericOutputBlackoilModule.hpp:368
ScalarBuffer permFact_
Definition: GenericOutputBlackoilModule.hpp:392
ScalarBuffer rsw_
Definition: GenericOutputBlackoilModule.hpp:380
ScalarBuffer pcog_
Definition: GenericOutputBlackoilModule.hpp:414
std::optional< RegionPhasePoreVolAverage > regionAvgDensity_
Definition: GenericOutputBlackoilModule.hpp:442
std::array< ScalarBuffer, numPhases > invB_
Definition: GenericOutputBlackoilModule.hpp:420
ScalarBuffer pSalt_
Definition: GenericOutputBlackoilModule.hpp:391
ScalarBuffer cFoam_
Definition: GenericOutputBlackoilModule.hpp:389
ScalarBuffer bubblePointPressure_
Definition: GenericOutputBlackoilModule.hpp:405
ScalarBuffer temperature_
Definition: GenericOutputBlackoilModule.hpp:378
ScalarBuffer ppcw_
Definition: GenericOutputBlackoilModule.hpp:400
FIPContainer< GetPropType< TypeTag, Properties::FluidSystem > > fipC_
Definition: GenericOutputBlackoilModule.hpp:361
ScalarBuffer rockCompTransMultiplier_
Definition: GenericOutputBlackoilModule.hpp:410
MechContainer< Scalar > mech_
Definition: GenericOutputBlackoilModule.hpp:417
ScalarBuffer dynamicPoreVolume_
Definition: GenericOutputBlackoilModule.hpp:375
ScalarBuffer minimumOilPressure_
Definition: GenericOutputBlackoilModule.hpp:408
ScalarBuffer gasFormationVolumeFactor_
Definition: GenericOutputBlackoilModule.hpp:371
std::array< ScalarBuffer, numPhases > residual_
Definition: GenericOutputBlackoilModule.hpp:427
void doAllocBuffers(unsigned bufferSize, unsigned reportStepNum, const bool substep, const bool log, const bool isRestart, const EclHysteresisConfig *hysteresisConfig, unsigned numOutputNnc=0, std::map< std::string, int > rstKeywords={})
ScalarBuffer shmax_
Definition: GenericOutputBlackoilModule.hpp:397
const Schedule & schedule_
Definition: GenericOutputBlackoilModule.hpp:337
FlowsContainer< GetPropType< TypeTag, Properties::FluidSystem > > flowsC_
Definition: GenericOutputBlackoilModule.hpp:429
ExtboContainer< Scalar > extboC_
Definition: GenericOutputBlackoilModule.hpp:393
void setupExtraBlockData(const std::size_t reportStepNum, std::function< bool(int)> isCartIdxOnThisRank)
ScalarBuffer oilSaturationPressure_
Definition: GenericOutputBlackoilModule.hpp:384
InterRegFlowMap interRegionFlows_
Definition: GenericOutputBlackoilModule.hpp:342
ScalarBuffer pcgw_
Definition: GenericOutputBlackoilModule.hpp:412
ScalarBuffer cPolymer_
Definition: GenericOutputBlackoilModule.hpp:388
void setupBlockData(std::function< bool(int)> isCartIdxOnThisRank)
ScalarBuffer rvw_
Definition: GenericOutputBlackoilModule.hpp:382
std::array< ScalarBuffer, numPhases > saturation_
Definition: GenericOutputBlackoilModule.hpp:419
std::unordered_map< std::string, std::vector< int > > regions_
Definition: GenericOutputBlackoilModule.hpp:362
ScalarBuffer rPorV_
Definition: GenericOutputBlackoilModule.hpp:376
ScalarBuffer oilVaporizationFactor_
Definition: GenericOutputBlackoilModule.hpp:402
std::vector< int > failedCellsPd_
Definition: GenericOutputBlackoilModule.hpp:369
ScalarBuffer rs_
Definition: GenericOutputBlackoilModule.hpp:379
ScalarBuffer drsdtcon_
Definition: GenericOutputBlackoilModule.hpp:385
ScalarBuffer sSol_
Definition: GenericOutputBlackoilModule.hpp:386
std::map< std::pair< std::string, int >, double > extraBlockData_
Definition: GenericOutputBlackoilModule.hpp:437
ScalarBuffer pressureTimesPoreVolume_
Definition: GenericOutputBlackoilModule.hpp:373
ScalarBuffer gasDissolutionFactor_
Definition: GenericOutputBlackoilModule.hpp:401
std::array< ScalarBuffer, numPhases > viscosity_
Definition: GenericOutputBlackoilModule.hpp:422
bool forceDisableFipOutput_
Definition: GenericOutputBlackoilModule.hpp:357
ScalarBuffer soMax_
Definition: GenericOutputBlackoilModule.hpp:394
ScalarBuffer sgmax_
Definition: GenericOutputBlackoilModule.hpp:396
ScalarBuffer somin_
Definition: GenericOutputBlackoilModule.hpp:398
ScalarBuffer hydrocarbonPoreVolume_
Definition: GenericOutputBlackoilModule.hpp:372
const std::optional< Inplace > & initialInplace() const
Definition: GenericOutputBlackoilModule.hpp:225
ScalarBuffer waterVaporizationFactor_
Definition: GenericOutputBlackoilModule.hpp:404
ScalarBuffer cSalt_
Definition: GenericOutputBlackoilModule.hpp:390
TracerContainer< GetPropType< TypeTag, Properties::FluidSystem > > tracerC_
Definition: GenericOutputBlackoilModule.hpp:425
ScalarBuffer rv_
Definition: GenericOutputBlackoilModule.hpp:381
ScalarBuffer pcow_
Definition: GenericOutputBlackoilModule.hpp:413
ScalarBuffer swMax_
Definition: GenericOutputBlackoilModule.hpp:395
ScalarBuffer pressureTimesHydrocarbonVolume_
Definition: GenericOutputBlackoilModule.hpp:374
ScalarBuffer rswSol_
Definition: GenericOutputBlackoilModule.hpp:387
Inter-region flow accumulation maps for all region definition arrays.
Definition: InterRegFlows.hpp:179
void addConnection(const Cell &source, const Cell &destination, const data::InterRegFlowMap::FlowRates &rates)
void clear()
Clear all internal buffers, but preserve allocated capacity.
Output module for the results black oil model writing in ECL binary format.
Definition: OutputBlackoilModule.hpp:85
void processElement(const ElementContext &elemCtx)
Modify the internal buffers according to the intensive quanties relevant for an element.
Definition: OutputBlackoilModule.hpp:237
void initializeFluxData()
Prepare for capturing connection fluxes, particularly to account for inter-region flows.
Definition: OutputBlackoilModule.hpp:481
void setupExtractors(const bool isSubStep, const int reportStepNum)
Setup list of active element-level data extractors.
Definition: OutputBlackoilModule.hpp:218
void allocBuffers(const unsigned bufferSize, const unsigned reportStepNum, const bool substep, const bool log, const bool isRestart)
Allocate memory for the scalar fields we would like to write to ECL output files.
Definition: OutputBlackoilModule.hpp:196
void processFluxes(const ElementContext &elemCtx, ActiveIndex &&activeIndex, CartesianIndex &&cartesianIndex)
Capture connection fluxes, particularly to account for inter-region flows.
Definition: OutputBlackoilModule.hpp:444
void clearExtractors()
Clear list of active element-level data extractors.
Definition: OutputBlackoilModule.hpp:226
void outputFipAndResvLogToCSV(const std::size_t reportStepNum, const bool substep, const Parallel::Communication &comm)
Definition: OutputBlackoilModule.hpp:378
void assignToFluidState(FluidState &fs, unsigned elemIdx) const
Definition: OutputBlackoilModule.hpp:505
void initHysteresisParams(Simulator &simulator, unsigned elemIdx) const
Definition: OutputBlackoilModule.hpp:555
void updateFluidInPlace(const ElementContext &elemCtx)
Definition: OutputBlackoilModule.hpp:620
OutputBlackOilModule(const Simulator &simulator, const SummaryConfig &smryCfg, const CollectDataOnIORankType &collectOnIORank)
Definition: OutputBlackoilModule.hpp:133
void outputFipAndResvLog(const Inplace &inplace, const std::size_t reportStepNum, double elapsed, boost::posix_time::ptime currentDate, const bool substep, const Parallel::Communication &comm)
Definition: OutputBlackoilModule.hpp:327
const InterRegFlowMap & getInterRegFlows() const
Get read-only access to collection of inter-region flows.
Definition: OutputBlackoilModule.hpp:499
void processElementBlockData(const ElementContext &elemCtx)
Definition: OutputBlackoilModule.hpp:283
void finalizeFluxData()
Finalize capturing connection fluxes.
Definition: OutputBlackoilModule.hpp:491
void updateFluidInPlace(const unsigned globalDofIdx, const IntensiveQuantities &intQuants, const double totVolume)
Definition: OutputBlackoilModule.hpp:627
Declare the properties used by the infrastructure code of the finite volume discretizations.
Dune::Communication< MPIComm > Communication
Definition: ParallelCommunication.hpp:30
Definition: blackoilboundaryratevector.hh:39
std::string moduleVersionName()
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:233
This file provides the infrastructure to retrieve run-time parameters.
The Opm property system, traits with inheritance.
Minimal characteristics of a cell from a simulation grid.
Definition: InterRegFlows.hpp:50