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<
class VectorType>
123 static Scalar value_or_zero(
int idx,
const VectorType& v)
128 return v.empty() ? 0.0 : v[idx];
133 const SummaryConfig& smryCfg,
135 :
BaseType(simulator.vanguard().eclState(),
136 simulator.vanguard().schedule(),
138 simulator.vanguard().summaryState(),
140 [this](const int idx)
141 {
return simulator_.problem().eclWriter().collectOnIORank().localIdxToGlobalIdx(idx); },
142 simulator.vanguard().grid().comm(),
143 getPropValue<TypeTag, Properties::EnableEnergy>(),
144 getPropValue<TypeTag, Properties::EnableTemperature>(),
145 getPropValue<TypeTag, Properties::EnableMech>(),
146 getPropValue<TypeTag, Properties::EnableSolvent>(),
147 getPropValue<TypeTag, Properties::EnablePolymer>(),
148 getPropValue<TypeTag, Properties::EnableFoam>(),
149 getPropValue<TypeTag, Properties::EnableBrine>(),
150 getPropValue<TypeTag, Properties::EnableSaltPrecipitation>(),
151 getPropValue<TypeTag, Properties::EnableExtbo>(),
152 getPropValue<TypeTag, Properties::EnableMICP>())
153 , simulator_(simulator)
154 , collectOnIORank_(collectOnIORank)
156 for (
auto& region_pair : this->
regions_) {
157 this->createLocalRegion_(region_pair.second);
160 auto isCartIdxOnThisRank = [&collectOnIORank](
const int idx) {
161 return collectOnIORank.isCartIdxOnThisRank(idx);
166 if (! Parameters::Get<Parameters::OwnerCellsFirst>()) {
167 const std::string msg =
"The output code does not support --owner-cells-first=false.";
168 if (collectOnIORank.isIORank()) {
171 OPM_THROW_NOLOG(std::runtime_error, msg);
174 if (smryCfg.match(
"[FB]PP[OGW]") || smryCfg.match(
"RPP[OGW]*")) {
175 auto rset = this->
eclState_.fieldProps().fip_regions();
176 rset.push_back(
"PVTNUM");
182 .emplace(this->simulator_.gridView().comm(),
183 FluidSystem::numPhases, rset,
184 [fp = std::cref(this->eclState_.fieldProps())]
185 (
const std::string& rsetName) ->
decltype(
auto)
186 { return fp.get().get_int(rsetName); });
196 const unsigned reportStepNum,
199 const bool isRestart)
205 const auto& problem = this->simulator_.problem();
212 &problem.materialLawManager()->hysteresisConfig(),
213 problem.eclWriter().getOutputNnc().size());
218 const int reportStepNum)
220 this->setupElementExtractors_();
221 this->setupBlockExtractors_(isSubStep, reportStepNum);
227 this->extractors_.clear();
228 this->blockExtractors_.clear();
229 this->extraBlockExtractors_.clear();
243 if (this->extractors_.empty()) {
247 const auto& matLawManager = simulator_.problem().materialLawManager();
250 for (
unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(0); ++dofIdx) {
251 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, 0);
252 const auto& fs = intQuants.fluidState();
255 elemCtx.globalSpaceIndex(dofIdx, 0),
256 elemCtx.primaryVars(dofIdx, 0).pvtRegionIndex(),
263 if (matLawManager->enableHysteresis()) {
264 if (FluidSystem::phaseIsActive(oilPhaseIdx) && FluidSystem::phaseIsActive(waterPhaseIdx)) {
265 matLawManager->oilWaterHysteresisParams(hysterParams.
somax,
270 if (FluidSystem::phaseIsActive(oilPhaseIdx) && FluidSystem::phaseIsActive(gasPhaseIdx)) {
271 matLawManager->gasOilHysteresisParams(hysterParams.
sgmax,
289 if (this->blockExtractors_.empty() && this->extraBlockExtractors_.empty()) {
293 for (
unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(0); ++dofIdx) {
295 const auto globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, 0);
296 const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx);
298 const auto be_it = this->blockExtractors_.find(cartesianIdx);
299 const auto bee_it = this->extraBlockExtractors_.find(cartesianIdx);
300 if (be_it == this->blockExtractors_.end() &&
301 bee_it == this->extraBlockExtractors_.end())
306 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, 0);
307 const auto& fs = intQuants.fluidState();
317 if (be_it != this->blockExtractors_.end()) {
320 if (bee_it != this->extraBlockExtractors_.end()) {
327 const std::size_t reportStepNum,
329 boost::posix_time::ptime currentDate,
334 if (comm.rank() != 0) {
339 std::unique_ptr<FIPConfig> fipSched;
340 if (reportStepNum > 0) {
341 const auto& rpt = this->
schedule_[reportStepNum-1].rpt_config.get();
342 fipSched = std::make_unique<FIPConfig>(rpt);
344 const FIPConfig& fipc = reportStepNum == 0 ? this->
eclState_.getEclipseConfig().fip()
349 this->
logOutput_.timeStamp(
"BALANCE", elapsed, reportStepNum, currentDate);
352 this->
logOutput_.fip(inplace, initial_inplace,
"");
354 if (fipc.output(FIPConfig::OutputField::FIPNUM)) {
355 this->
logOutput_.fip(inplace, initial_inplace,
"FIPNUM");
357 if (fipc.output(FIPConfig::OutputField::RESV))
361 if (fipc.output(FIPConfig::OutputField::FIP)) {
362 for (
const auto& reg : this->regions_) {
363 if (reg.first !=
"FIPNUM") {
364 std::ostringstream ss;
365 ss <<
"BAL" << reg.first.substr(3);
366 this->
logOutput_.timeStamp(ss.str(), elapsed, reportStepNum, currentDate);
367 this->
logOutput_.fip(inplace, initial_inplace, reg.first);
369 if (fipc.output(FIPConfig::OutputField::RESV))
381 if (comm.rank() != 0) {
385 if ((reportStepNum == 0) && (!substep) &&
386 (this->
schedule_.initialReportConfiguration().has_value()) &&
387 (this->schedule_.initialReportConfiguration()->contains(
"CSVFIP"))) {
389 std::ostringstream csv_stream;
395 this->
logOutput_.fip_csv(csv_stream, initial_inplace,
"FIPNUM");
397 for (
const auto& reg : this->regions_) {
398 if (reg.first !=
"FIPNUM") {
399 this->
logOutput_.fip_csv(csv_stream, initial_inplace, reg.first);
403 const IOConfig& io = this->
eclState_.getIOConfig();
404 auto csv_fname = io.getOutputDir() +
"/" + io.getBaseName() +
".CSV";
406 std::ofstream outputFile(csv_fname);
408 outputFile << csv_stream.str();
442 template <
class ActiveIndex,
class CartesianIndex>
444 ActiveIndex&& activeIndex,
445 CartesianIndex&& cartesianIndex)
448 const auto identifyCell = [&activeIndex, &cartesianIndex](
const Element& elem)
451 const auto cellIndex = activeIndex(elem);
454 static_cast<int>(cellIndex),
455 cartesianIndex(cellIndex),
456 elem.partitionType() == Dune::InteriorEntity
460 const auto timeIdx = 0u;
461 const auto& stencil = elemCtx.stencil(timeIdx);
462 const auto numInteriorFaces = elemCtx.numInteriorFaces(timeIdx);
464 for (
auto scvfIdx = 0 * numInteriorFaces; scvfIdx < numInteriorFaces; ++scvfIdx) {
465 const auto& face = stencil.interiorFace(scvfIdx);
466 const auto left = identifyCell(stencil.element(face.interiorIndex()));
467 const auto right = identifyCell(stencil.element(face.exteriorIndex()));
469 const auto rates = this->
470 getComponentSurfaceRates(elemCtx, face.area(), scvfIdx, timeIdx);
503 template <
class Flu
idState>
506 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
510 fs.setSaturation(phaseIdx, this->
saturation_[phaseIdx][elemIdx]);
516 std::array<Scalar, numPhases> pc = {0};
517 const MaterialLawParams& matParams = simulator_.problem().materialLawParams(elemIdx);
518 MaterialLaw::capillaryPressures(pc, matParams, fs);
520 Valgrind::CheckDefined(pc);
522 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
523 if (!FluidSystem::phaseIsActive(phaseIdx))
526 if (Indices::oilEnabled)
527 fs.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx]));
528 else if (Indices::gasEnabled)
529 fs.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx]));
530 else if (Indices::waterEnabled)
532 fs.setPressure(phaseIdx, pressure);
538 if constexpr (enableDissolvedGas) {
539 if (!this->
rs_.empty())
540 fs.setRs(this->rs_[elemIdx]);
541 if (!this->
rv_.empty())
542 fs.setRv(this->rv_[elemIdx]);
544 if constexpr (enableDisgasInWater) {
545 if (!this->
rsw_.empty())
546 fs.setRsw(this->rsw_[elemIdx]);
548 if constexpr (enableVapwat) {
549 if (!this->
rvw_.empty())
550 fs.setRvw(this->rvw_[elemIdx]);
556 if (!this->
soMax_.empty())
557 simulator.problem().setMaxOilSaturation(elemIdx, this->
soMax_[elemIdx]);
559 if (simulator.problem().materialLawManager()->enableHysteresis()) {
560 auto matLawManager = simulator.problem().materialLawManager();
562 if (FluidSystem::phaseIsActive(oilPhaseIdx)
563 && FluidSystem::phaseIsActive(waterPhaseIdx)) {
568 if (matLawManager->enableNonWettingHysteresis()) {
569 if (!this->
soMax_.empty()) {
570 somax = this->
soMax_[elemIdx];
573 if (matLawManager->enableWettingHysteresis()) {
574 if (!this->
swMax_.empty()) {
575 swmax = this->
swMax_[elemIdx];
578 if (matLawManager->enablePCHysteresis()) {
579 if (!this->
swmin_.empty()) {
580 swmin = this->
swmin_[elemIdx];
583 matLawManager->setOilWaterHysteresisParams(
584 somax, swmax, swmin, elemIdx);
586 if (FluidSystem::phaseIsActive(oilPhaseIdx)
587 && FluidSystem::phaseIsActive(gasPhaseIdx)) {
592 if (matLawManager->enableNonWettingHysteresis()) {
593 if (!this->
sgmax_.empty()) {
594 sgmax = this->
sgmax_[elemIdx];
597 if (matLawManager->enableWettingHysteresis()) {
598 if (!this->
shmax_.empty()) {
599 shmax = this->
shmax_[elemIdx];
602 if (matLawManager->enablePCHysteresis()) {
603 if (!this->
somin_.empty()) {
604 somin = this->
somin_[elemIdx];
607 matLawManager->setGasOilHysteresisParams(
608 sgmax, shmax, somin, elemIdx);
613 if (simulator_.vanguard().eclState().fieldProps().has_double(
"SWATINIT")) {
614 simulator.problem().materialLawManager()
615 ->applyRestartSwatInit(elemIdx, this->
ppcw_[elemIdx]);
621 for (
unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(0); ++dofIdx) {
622 updateFluidInPlace_(elemCtx, dofIdx);
627 const IntensiveQuantities& intQuants,
628 const double totVolume)
630 this->updateFluidInPlace_(globalDofIdx, intQuants, totVolume);
634 template <
typename T>
635 using RemoveCVR = std::remove_cv_t<std::remove_reference_t<T>>;
637 template <
typename,
class =
void>
638 struct HasGeoMech :
public std::false_type {};
640 template <
typename Problem>
642 Problem, std::void_t<decltype(std::declval<Problem>().geoMechModel())>
643 > :
public std::true_type {};
645 bool isDefunctParallelWell(
const std::string& wname)
const override
647 if (simulator_.gridView().comm().size() == 1)
649 const auto& parallelWells = simulator_.vanguard().parallelWells();
650 std::pair<std::string, bool> value {wname,
true};
651 auto candidate = std::lower_bound(parallelWells.begin(), parallelWells.end(), value);
652 return candidate == parallelWells.end() || *candidate != value;
655 bool isOwnedByCurrentRank(
const std::string& wname)
const override
657 return this->simulator_.problem().wellModel().isOwner(wname);
660 bool isOnCurrentRank(
const std::string& wname)
const override
662 return this->simulator_.problem().wellModel().hasLocalCells(wname);
665 void updateFluidInPlace_(
const ElementContext& elemCtx,
const unsigned dofIdx)
667 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, 0);
668 const unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, 0);
669 const auto totVolume = elemCtx.simulator().model().dofTotalVolume(globalDofIdx);
671 this->updateFluidInPlace_(globalDofIdx, intQuants, totVolume);
674 void updateFluidInPlace_(
const unsigned globalDofIdx,
675 const IntensiveQuantities& intQuants,
676 const double totVolume)
680 this->updateTotalVolumesAndPressures_(globalDofIdx, intQuants, totVolume);
683 this->updatePhaseInplaceVolumes_(globalDofIdx, intQuants, totVolume);
687 void createLocalRegion_(std::vector<int>& region)
693 region.resize(simulator_.gridView().size(0));
694 std::size_t elemIdx = 0;
695 for (
const auto& elem : elements(simulator_.gridView())) {
696 if (elem.partitionType() != Dune::InteriorEntity) {
704 template <
typename Flu
idState>
705 void aggregateAverageDensityContributions_(
const FluidState& fs,
706 const unsigned int globalDofIdx,
709 auto pvCellValue = RegionPhasePoreVolAverage::CellValue{};
710 pvCellValue.porv = porv;
712 for (
auto phaseIdx = 0*FluidSystem::numPhases;
713 phaseIdx < FluidSystem::numPhases; ++phaseIdx)
715 if (! FluidSystem::phaseIsActive(phaseIdx)) {
719 pvCellValue.value = getValue(fs.density(phaseIdx));
720 pvCellValue.sat = getValue(fs.saturation(phaseIdx));
723 ->addCell(globalDofIdx,
724 RegionPhasePoreVolAverage::Phase { phaseIdx },
744 data::InterRegFlowMap::FlowRates
745 getComponentSurfaceRates(
const ElementContext& elemCtx,
746 const Scalar faceArea,
747 const std::size_t scvfIdx,
748 const std::size_t timeIdx)
const
750 using Component = data::InterRegFlowMap::Component;
752 auto rates = data::InterRegFlowMap::FlowRates {};
754 const auto& extQuant = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
756 const auto alpha = getValue(extQuant.extrusionFactor()) * faceArea;
758 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
759 const auto& up = elemCtx
760 .intensiveQuantities(extQuant.upstreamIndex(oilPhaseIdx), timeIdx);
762 const auto pvtReg = up.pvtRegionIndex();
764 const auto bO = getValue(getInvB_<FluidSystem, FluidState, Scalar>
765 (up.fluidState(), oilPhaseIdx, pvtReg));
767 const auto qO = alpha * bO * getValue(extQuant.volumeFlux(oilPhaseIdx));
769 rates[Component::Oil] += qO;
771 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
772 const auto Rs = getValue(
773 BlackOil::getRs_<FluidSystem, FluidState, Scalar>
774 (up.fluidState(), pvtReg));
776 rates[Component::Gas] += qO * Rs;
777 rates[Component::Disgas] += qO * Rs;
781 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
782 const auto& up = elemCtx
783 .intensiveQuantities(extQuant.upstreamIndex(gasPhaseIdx), timeIdx);
785 const auto pvtReg = up.pvtRegionIndex();
787 const auto bG = getValue(getInvB_<FluidSystem, FluidState, Scalar>
788 (up.fluidState(), gasPhaseIdx, pvtReg));
790 const auto qG = alpha * bG * getValue(extQuant.volumeFlux(gasPhaseIdx));
792 rates[Component::Gas] += qG;
794 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
795 const auto Rv = getValue(
796 BlackOil::getRv_<FluidSystem, FluidState, Scalar>
797 (up.fluidState(), pvtReg));
799 rates[Component::Oil] += qG * Rv;
800 rates[Component::Vapoil] += qG * Rv;
804 if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
805 const auto& up = elemCtx
806 .intensiveQuantities(extQuant.upstreamIndex(waterPhaseIdx), timeIdx);
808 const auto pvtReg = up.pvtRegionIndex();
810 const auto bW = getValue(getInvB_<FluidSystem, FluidState, Scalar>
811 (up.fluidState(), waterPhaseIdx, pvtReg));
813 rates[Component::Water] +=
814 alpha * bW * getValue(extQuant.volumeFlux(waterPhaseIdx));
820 template <
typename Flu
idState>
821 Scalar hydroCarbonFraction(
const FluidState& fs)
const
823 if (this->
eclState_.runspec().co2Storage()) {
830 auto hydrocarbon = Scalar {0};
831 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
832 hydrocarbon += getValue(fs.saturation(oilPhaseIdx));
835 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
836 hydrocarbon += getValue(fs.saturation(gasPhaseIdx));
842 void updateTotalVolumesAndPressures_(
const unsigned globalDofIdx,
843 const IntensiveQuantities& intQuants,
844 const double totVolume)
846 const auto& fs = intQuants.fluidState();
848 const double pv = totVolume * intQuants.porosity().value();
849 const auto hydrocarbon = this->hydroCarbonFraction(fs);
853 totVolume * intQuants.referencePorosity());
860 !this->pressureTimesPoreVolume_.empty())
863 assert(this->
fipC_.
get(Inplace::Phase::PoreVolume).size() == this->pressureTimesPoreVolume_.size());
865 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
867 getValue(fs.pressure(oilPhaseIdx)) * pv;
872 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
874 getValue(fs.pressure(gasPhaseIdx)) * pv;
879 else if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
881 getValue(fs.pressure(waterPhaseIdx)) * pv;
886 void updatePhaseInplaceVolumes_(
const unsigned globalDofIdx,
887 const IntensiveQuantities& intQuants,
888 const double totVolume)
890 std::array<Scalar, FluidSystem::numPhases> fip {};
891 std::array<Scalar, FluidSystem::numPhases> fipr{};
893 const auto& fs = intQuants.fluidState();
894 const auto pv = totVolume * intQuants.porosity().value();
896 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
897 if (!FluidSystem::phaseIsActive(phaseIdx)) {
901 const auto b = getValue(fs.invB(phaseIdx));
902 const auto s = getValue(fs.saturation(phaseIdx));
904 fipr[phaseIdx] = s * pv;
905 fip [phaseIdx] = b * fipr[phaseIdx];
910 fs.saltConcentration().value(),
913 if (FluidSystem::phaseIsActive(oilPhaseIdx) &&
914 FluidSystem::phaseIsActive(gasPhaseIdx))
916 this->updateOilGasDistribution(globalDofIdx, fs, fip);
919 if (FluidSystem::phaseIsActive(waterPhaseIdx) &&
920 FluidSystem::phaseIsActive(gasPhaseIdx))
922 this->updateGasWaterDistribution(globalDofIdx, fs, fip);
925 if (FluidSystem::phaseIsActive(gasPhaseIdx) &&
928 this->updateCO2InGas(globalDofIdx, pv, intQuants);
932 (FluidSystem::phaseIsActive(waterPhaseIdx) ||
933 FluidSystem::phaseIsActive(oilPhaseIdx)))
935 this->updateCO2InWater(globalDofIdx, pv, fs);
938 if constexpr(enableMICP) {
939 const auto surfVolWat = pv * getValue(fs.invB(waterPhaseIdx));
941 this->updateMicrobialMass(globalDofIdx, intQuants, surfVolWat);
944 this->updateOxygenMass(globalDofIdx, intQuants, surfVolWat);
947 this->updateUreaMass(globalDofIdx, intQuants, surfVolWat);
950 this->updateBiofilmMass(globalDofIdx, intQuants, totVolume);
953 this->updateCalciteMass(globalDofIdx, intQuants, totVolume);
959 this->updateWaterMass(globalDofIdx, fs, fip);
963 template <
typename Flu
idState,
typename FIPArray>
964 void updateOilGasDistribution(
const unsigned globalDofIdx,
965 const FluidState& fs,
969 const auto gasInPlaceLiquid = getValue(fs.Rs()) * fip[oilPhaseIdx];
970 const auto oilInPlaceGas = getValue(fs.Rv()) * fip[gasPhaseIdx];
975 template <
typename Flu
idState,
typename FIPArray>
976 void updateGasWaterDistribution(
const unsigned globalDofIdx,
977 const FluidState& fs,
981 const auto gasInPlaceWater = getValue(fs.Rsw()) * fip[waterPhaseIdx];
982 const auto waterInPlaceGas = getValue(fs.Rvw()) * fip[gasPhaseIdx];
987 template <
typename IntensiveQuantities>
988 void updateCO2InGas(
const unsigned globalDofIdx,
990 const IntensiveQuantities& intQuants)
992 const auto& scaledDrainageInfo = this->simulator_.problem().materialLawManager()
993 ->oilWaterScaledEpsInfoDrainage(globalDofIdx);
995 const auto& fs = intQuants.fluidState();
996 Scalar sgcr = scaledDrainageInfo.Sgcr;
997 if (this->simulator_.problem().materialLawManager()->enableHysteresis()) {
998 const auto& matParams = simulator_.problem().materialLawParams(globalDofIdx);
999 sgcr = MaterialLaw::trappedGasSaturation(matParams,
false);
1002 Scalar trappedGasSaturation = scaledDrainageInfo.Sgcr;
1003 if (this->
fipC_.
has(Inplace::Phase::CO2MassInGasPhaseMaximumTrapped) ||
1004 this->fipC_.has(Inplace::Phase::CO2MassInGasPhaseMaximumUnTrapped))
1006 if (this->simulator_.problem().materialLawManager()->enableHysteresis()) {
1007 const auto& matParams = simulator_.problem().materialLawParams(globalDofIdx);
1009 trappedGasSaturation = MaterialLaw::trappedGasSaturation(matParams,
true);
1013 const Scalar sg = getValue(fs.saturation(gasPhaseIdx));
1014 Scalar strandedGasSaturation = scaledDrainageInfo.Sgcr;
1015 if (this->
fipC_.
has(Inplace::Phase::CO2MassInGasPhaseEffectiveTrapped) ||
1016 this->fipC_.has(Inplace::Phase::CO2MassInGasPhaseEffectiveUnTrapped))
1018 if (this->simulator_.problem().materialLawManager()->enableHysteresis()) {
1019 const auto& matParams = simulator_.problem().materialLawParams(globalDofIdx);
1020 const double krg = getValue(intQuants.relativePermeability(gasPhaseIdx));
1021 strandedGasSaturation = MaterialLaw::strandedGasSaturation(matParams, sg, krg);
1025 const typename FIPContainer<FluidSystem>::Co2InGasInput v{
1029 getValue(fs.density(gasPhaseIdx)),
1030 FluidSystem::phaseIsActive(waterPhaseIdx)
1031 ? FluidSystem::convertRvwToXgW(getValue(fs.Rvw()), fs.pvtRegionIndex())
1032 : FluidSystem::convertRvToXgO(getValue(fs.Rv()), fs.pvtRegionIndex()),
1033 FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex()),
1034 trappedGasSaturation,
1035 strandedGasSaturation,
1041 template <
typename Flu
idState>
1042 void updateCO2InWater(
const unsigned globalDofIdx,
1044 const FluidState& fs)
1046 const auto co2InWater = FluidSystem::phaseIsActive(oilPhaseIdx)
1047 ? this->co2InWaterFromOil(fs, pv)
1048 : this->co2InWaterFromWater(fs, pv);
1050 const Scalar mM = FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex());
1055 template <
typename Flu
idState>
1056 Scalar co2InWaterFromWater(
const FluidState& fs,
const double pv)
const
1058 const double rhow = getValue(fs.density(waterPhaseIdx));
1059 const double sw = getValue(fs.saturation(waterPhaseIdx));
1060 const double xwG = FluidSystem::convertRswToXwG(getValue(fs.Rsw()), fs.pvtRegionIndex());
1062 const Scalar mM = FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex());
1064 return xwG * pv * rhow * sw / mM;
1067 template <
typename Flu
idState>
1068 Scalar co2InWaterFromOil(
const FluidState& fs,
const double pv)
const
1070 const double rhoo = getValue(fs.density(oilPhaseIdx));
1071 const double so = getValue(fs.saturation(oilPhaseIdx));
1072 const double xoG = FluidSystem::convertRsToXoG(getValue(fs.Rs()), fs.pvtRegionIndex());
1074 const Scalar mM = FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex());
1076 return xoG * pv * rhoo * so / mM;
1079 template <
typename Flu
idState,
typename FIPArray>
1080 void updateWaterMass(
const unsigned globalDofIdx,
1081 const FluidState& fs,
1085 const Scalar rhoW = FluidSystem::referenceDensity(waterPhaseIdx, fs.pvtRegionIndex());
1090 template <
typename IntensiveQuantities>
1091 void updateMicrobialMass(
const unsigned globalDofIdx,
1092 const IntensiveQuantities& intQuants,
1093 const double surfVolWat)
1095 const Scalar mass = surfVolWat * intQuants.microbialConcentration().value();
1100 template <
typename IntensiveQuantities>
1101 void updateOxygenMass(
const unsigned globalDofIdx,
1102 const IntensiveQuantities& intQuants,
1103 const double surfVolWat)
1105 const Scalar mass = surfVolWat * intQuants.oxygenConcentration().value();
1110 template <
typename IntensiveQuantities>
1111 void updateUreaMass(
const unsigned globalDofIdx,
1112 const IntensiveQuantities& intQuants,
1113 const double surfVolWat)
1115 const Scalar mass = surfVolWat * intQuants.ureaConcentration().value();
1120 template <
typename IntensiveQuantities>
1121 void updateBiofilmMass(
const unsigned globalDofIdx,
1122 const IntensiveQuantities& intQuants,
1123 const double totVolume)
1125 const Scalar mass = totVolume * intQuants.biofilmMass().value();
1130 template <
typename IntensiveQuantities>
1131 void updateCalciteMass(
const unsigned globalDofIdx,
1132 const IntensiveQuantities& intQuants,
1133 const double totVolume)
1135 const Scalar mass = totVolume * intQuants.calciteMass().value();
1141 void setupElementExtractors_()
1143 using Entry =
typename Extractor::Entry;
1144 using Context =
typename Extractor::Context;
1145 using ScalarEntry =
typename Extractor::ScalarEntry;
1146 using PhaseEntry =
typename Extractor::PhaseEntry;
1148 const bool hasResidual = simulator_.model().linearizer().residual().size() > 0;
1149 const auto& hysteresisConfig = simulator_.problem().materialLawManager()->hysteresisConfig();
1151 auto extractors = std::array{
1153 [](
const unsigned phase,
const Context& ectx)
1154 {
return getValue(ectx.fs.saturation(phase)); }
1157 Entry{PhaseEntry{&this->
invB_,
1158 [](
const unsigned phase,
const Context& ectx)
1159 {
return getValue(ectx.fs.invB(phase)); }
1163 [](
const unsigned phase,
const Context& ectx)
1164 {
return getValue(ectx.fs.density(phase)); }
1168 [](
const unsigned phase,
const Context& ectx)
1169 {
return getValue(ectx.intQuants.relativePermeability(phase)); }
1173 [
this](
const unsigned phaseIdx,
const Context& ectx)
1175 if (this->
extboC_.allocated() && phaseIdx == oilPhaseIdx) {
1176 return getValue(ectx.intQuants.oilViscosity());
1178 else if (this->
extboC_.allocated() && phaseIdx == gasPhaseIdx) {
1179 return getValue(ectx.intQuants.gasViscosity());
1182 return getValue(ectx.fs.viscosity(phaseIdx));
1188 [&modelResid = this->simulator_.model().linearizer().residual()]
1189 (
const unsigned phaseIdx,
const Context& ectx)
1191 const unsigned sIdx = FluidSystem::solventComponentIndex(phaseIdx);
1192 const unsigned activeCompIdx = FluidSystem::canonicalToActiveCompIdx(sIdx);
1193 return modelResid[ectx.globalDofIdx][activeCompIdx];
1199 [&problem = this->simulator_.problem()](
const Context& ectx)
1201 return problem.template
1202 rockCompPoroMultiplier<Scalar>(ectx.intQuants,
1208 [&problem = this->simulator_.problem()](
const Context& ectx)
1211 template rockCompTransMultiplier<Scalar>(ectx.intQuants,
1216 [&problem = this->simulator_.problem()](
const Context& ectx)
1218 return std::min(getValue(ectx.fs.pressure(oilPhaseIdx)),
1219 problem.minOilPressure(ectx.globalDofIdx));
1225 &vanguard = this->simulator_.vanguard()](
const Context& ectx)
1229 FluidSystem::bubblePointPressure(ectx.fs,
1230 ectx.intQuants.pvtRegionIndex())
1232 }
catch (
const NumericalProblem&) {
1233 const auto cartesianIdx = vanguard.cartesianIndex(ectx.globalDofIdx);
1234 failedCells.push_back(cartesianIdx);
1242 &vanguard = this->simulator_.vanguard()](
const Context& ectx)
1246 FluidSystem::dewPointPressure(ectx.fs,
1247 ectx.intQuants.pvtRegionIndex())
1249 }
catch (
const NumericalProblem&) {
1250 const auto cartesianIdx = vanguard.cartesianIndex(ectx.globalDofIdx);
1251 failedCells.push_back(cartesianIdx);
1258 [&problem = simulator_.problem()](
const Context& ectx)
1259 {
return problem.overburdenPressure(ectx.globalDofIdx); }
1263 [](
const Context& ectx)
1264 {
return getValue(ectx.fs.temperature(oilPhaseIdx)); }
1267 Entry{ScalarEntry{&this->
sSol_,
1268 [](
const Context& ectx)
1269 {
return getValue(ectx.intQuants.solventSaturation()); }
1272 Entry{ScalarEntry{&this->
rswSol_,
1273 [](
const Context& ectx)
1274 {
return getValue(ectx.intQuants.rsSolw()); }
1278 [](
const Context& ectx)
1279 {
return getValue(ectx.intQuants.polymerConcentration()); }
1282 Entry{ScalarEntry{&this->
cFoam_,
1283 [](
const Context& ectx)
1284 {
return getValue(ectx.intQuants.foamConcentration()); }
1287 Entry{ScalarEntry{&this->
cSalt_,
1288 [](
const Context& ectx)
1289 {
return getValue(ectx.fs.saltConcentration()); }
1292 Entry{ScalarEntry{&this->
pSalt_,
1293 [](
const Context& ectx)
1294 {
return getValue(ectx.fs.saltSaturation()); }
1298 [](
const Context& ectx)
1299 {
return getValue(ectx.intQuants.permFactor()); }
1302 Entry{ScalarEntry{&this->
rPorV_,
1303 [&model = this->simulator_.model()](
const Context& ectx)
1305 const auto totVolume = model.dofTotalVolume(ectx.globalDofIdx);
1306 return totVolume * getValue(ectx.intQuants.porosity());
1310 Entry{ScalarEntry{&this->
rs_,
1311 [](
const Context& ectx)
1312 {
return getValue(ectx.fs.Rs()); }
1315 Entry{ScalarEntry{&this->
rv_,
1316 [](
const Context& ectx)
1317 {
return getValue(ectx.fs.Rv()); }
1320 Entry{ScalarEntry{&this->
rsw_,
1321 [](
const Context& ectx)
1322 {
return getValue(ectx.fs.Rsw()); }
1325 Entry{ScalarEntry{&this->
rvw_,
1326 [](
const Context& ectx)
1327 {
return getValue(ectx.fs.Rvw()); }
1330 Entry{ScalarEntry{&this->
ppcw_,
1331 [&matLawManager = *this->simulator_.problem().materialLawManager()]
1332 (
const Context& ectx)
1334 return matLawManager.
1335 oilWaterScaledEpsInfoDrainage(ectx.globalDofIdx).maxPcow;
1340 [&problem = this->simulator_.problem()](
const Context& ectx)
1342 return problem.drsdtcon(ectx.globalDofIdx,
1347 Entry{ScalarEntry{&this->
pcgw_,
1348 [](
const Context& ectx)
1350 return getValue(ectx.fs.pressure(gasPhaseIdx)) -
1351 getValue(ectx.fs.pressure(waterPhaseIdx));
1355 Entry{ScalarEntry{&this->
pcow_,
1356 [](
const Context& ectx)
1358 return getValue(ectx.fs.pressure(oilPhaseIdx)) -
1359 getValue(ectx.fs.pressure(waterPhaseIdx));
1363 Entry{ScalarEntry{&this->
pcog_,
1364 [](
const Context& ectx)
1366 return getValue(ectx.fs.pressure(gasPhaseIdx)) -
1367 getValue(ectx.fs.pressure(oilPhaseIdx));
1372 [](
const Context& ectx)
1374 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1376 return getValue(ectx.fs.pressure(oilPhaseIdx));
1378 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1380 return getValue(ectx.fs.pressure(gasPhaseIdx));
1384 return getValue(ectx.fs.pressure(waterPhaseIdx));
1390 [&problem = this->simulator_.problem()](
const Context& ectx)
1392 const Scalar SoMax = problem.maxOilSaturation(ectx.globalDofIdx);
1393 return FluidSystem::template
1394 saturatedDissolutionFactor<FluidState, Scalar>(ectx.fs,
1402 [&problem = this->simulator_.problem()](
const Context& ectx)
1404 const Scalar SoMax = problem.maxOilSaturation(ectx.globalDofIdx);
1405 return FluidSystem::template
1406 saturatedDissolutionFactor<FluidState, Scalar>(ectx.fs,
1414 [&problem = this->simulator_.problem()](
const Context& ectx)
1416 const Scalar SwMax = problem.maxWaterSaturation(ectx.globalDofIdx);
1417 return FluidSystem::template
1418 saturatedDissolutionFactor<FluidState, Scalar>(ectx.fs,
1426 [](
const Context& ectx)
1428 return FluidSystem::template
1429 saturatedVaporizationFactor<FluidState, Scalar>(ectx.fs,
1436 [](
const Context& ectx)
1438 return 1.0 / FluidSystem::template
1439 inverseFormationVolumeFactor<FluidState, Scalar>(ectx.fs,
1446 [](
const Context& ectx)
1448 return 1.0 / FluidSystem::template
1449 saturatedInverseFormationVolumeFactor<FluidState, Scalar>(ectx.fs,
1456 [](
const Context& ectx)
1458 return FluidSystem::template
1459 saturationPressure<FluidState, Scalar>(ectx.fs,
1465 Entry{ScalarEntry{&this->
soMax_,
1466 [&problem = this->simulator_.problem()](
const Context& ectx)
1468 return std::max(getValue(ectx.fs.saturation(oilPhaseIdx)),
1469 problem.maxOilSaturation(ectx.globalDofIdx));
1472 !hysteresisConfig.enableHysteresis()
1474 Entry{ScalarEntry{&this->
swMax_,
1475 [&problem = this->simulator_.problem()](
const Context& ectx)
1477 return std::max(getValue(ectx.fs.saturation(waterPhaseIdx)),
1478 problem.maxWaterSaturation(ectx.globalDofIdx));
1481 !hysteresisConfig.enableHysteresis()
1483 Entry{ScalarEntry{&this->
soMax_,
1484 [](
const Context& ectx)
1485 {
return ectx.hParams.somax; }
1487 hysteresisConfig.enableHysteresis() &&
1488 hysteresisConfig.enableNonWettingHysteresis() &&
1489 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1490 FluidSystem::phaseIsActive(waterPhaseIdx)
1492 Entry{ScalarEntry{&this->
swMax_,
1493 [](
const Context& ectx)
1494 {
return ectx.hParams.swmax; }
1496 hysteresisConfig.enableHysteresis() &&
1497 hysteresisConfig.enableWettingHysteresis() &&
1498 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1499 FluidSystem::phaseIsActive(waterPhaseIdx)
1501 Entry{ScalarEntry{&this->
swmin_,
1502 [](
const Context& ectx)
1503 {
return ectx.hParams.swmin; }
1505 hysteresisConfig.enableHysteresis() &&
1506 hysteresisConfig.enablePCHysteresis() &&
1507 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1508 FluidSystem::phaseIsActive(waterPhaseIdx)
1510 Entry{ScalarEntry{&this->
sgmax_,
1511 [](
const Context& ectx)
1512 {
return ectx.hParams.sgmax; }
1514 hysteresisConfig.enableHysteresis() &&
1515 hysteresisConfig.enableNonWettingHysteresis() &&
1516 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1517 FluidSystem::phaseIsActive(gasPhaseIdx)
1519 Entry{ScalarEntry{&this->
shmax_,
1520 [](
const Context& ectx)
1521 {
return ectx.hParams.shmax; }
1523 hysteresisConfig.enableHysteresis() &&
1524 hysteresisConfig.enableWettingHysteresis() &&
1525 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1526 FluidSystem::phaseIsActive(gasPhaseIdx)
1528 Entry{ScalarEntry{&this->
somin_,
1529 [](
const Context& ectx)
1530 {
return ectx.hParams.somin; }
1532 hysteresisConfig.enableHysteresis() &&
1533 hysteresisConfig.enablePCHysteresis() &&
1534 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1535 FluidSystem::phaseIsActive(gasPhaseIdx)
1537 Entry{[&model = this->simulator_.model(),
this](
const Context& ectx)
1541 const auto porv = ectx.intQuants.referencePorosity()
1542 * model.dofTotalVolume(ectx.globalDofIdx);
1544 this->aggregateAverageDensityContributions_(ectx.fs, ectx.globalDofIdx,
1545 static_cast<double>(porv));
1548 Entry{[&extboC = this->
extboC_](
const Context& ectx)
1550 extboC.assignVolumes(ectx.globalDofIdx,
1551 ectx.intQuants.xVolume().value(),
1552 ectx.intQuants.yVolume().value());
1553 extboC.assignZFraction(ectx.globalDofIdx,
1554 ectx.intQuants.zFraction().value());
1556 const Scalar stdVolOil = getValue(ectx.fs.saturation(oilPhaseIdx)) *
1557 getValue(ectx.fs.invB(oilPhaseIdx)) +
1558 getValue(ectx.fs.saturation(gasPhaseIdx)) *
1559 getValue(ectx.fs.invB(gasPhaseIdx)) *
1560 getValue(ectx.fs.Rv());
1561 const Scalar stdVolGas = getValue(ectx.fs.saturation(gasPhaseIdx)) *
1562 getValue(ectx.fs.invB(gasPhaseIdx)) *
1563 (1.0 - ectx.intQuants.yVolume().value()) +
1564 getValue(ectx.fs.saturation(oilPhaseIdx)) *
1565 getValue(ectx.fs.invB(oilPhaseIdx)) *
1566 getValue(ectx.fs.Rs()) *
1567 (1.0 - ectx.intQuants.xVolume().value());
1568 const Scalar stdVolCo2 = getValue(ectx.fs.saturation(gasPhaseIdx)) *
1569 getValue(ectx.fs.invB(gasPhaseIdx)) *
1570 ectx.intQuants.yVolume().value() +
1571 getValue(ectx.fs.saturation(oilPhaseIdx)) *
1572 getValue(ectx.fs.invB(oilPhaseIdx)) *
1573 getValue(ectx.fs.Rs()) *
1574 ectx.intQuants.xVolume().value();
1575 const Scalar rhoO = FluidSystem::referenceDensity(gasPhaseIdx, ectx.pvtRegionIdx);
1576 const Scalar rhoG = FluidSystem::referenceDensity(gasPhaseIdx, ectx.pvtRegionIdx);
1577 const Scalar rhoCO2 = ectx.intQuants.zRefDensity();
1578 const Scalar stdMassTotal = 1.0e-10 + stdVolOil * rhoO + stdVolGas * rhoG + stdVolCo2 * rhoCO2;
1579 extboC.assignMassFractions(ectx.globalDofIdx,
1580 stdVolGas * rhoG / stdMassTotal,
1581 stdVolOil * rhoO / stdMassTotal,
1582 stdVolCo2 * rhoCO2 / stdMassTotal);
1585 Entry{[&micpC = this->
micpC_](
const Context& ectx)
1587 micpC.assign(ectx.globalDofIdx,
1588 ectx.intQuants.microbialConcentration().value(),
1589 ectx.intQuants.oxygenConcentration().value(),
1590 ectx.intQuants.ureaConcentration().value(),
1591 ectx.intQuants.biofilmConcentration().value(),
1592 ectx.intQuants.calciteConcentration().value());
1593 }, this->
micpC_.allocated()
1595 Entry{[&rftC = this->
rftC_,
1596 &vanguard = this->simulator_.vanguard()](
const Context& ectx)
1598 const auto cartesianIdx = vanguard.cartesianIndex(ectx.globalDofIdx);
1599 rftC.assign(cartesianIdx,
1600 [&fs = ectx.fs]() {
return getValue(fs.pressure(oilPhaseIdx)); },
1601 [&fs = ectx.fs]() {
return getValue(fs.saturation(waterPhaseIdx)); },
1602 [&fs = ectx.fs]() {
return getValue(fs.saturation(gasPhaseIdx)); });
1606 &tM = this->simulator_.problem().tracerModel()](
const Context& ectx)
1608 tC.assignFreeConcentrations(ectx.globalDofIdx,
1609 [gIdx = ectx.globalDofIdx, &tM](
const unsigned tracerIdx)
1610 {
return tM.freeTracerConcentration(tracerIdx, gIdx); });
1611 tC.assignSolConcentrations(ectx.globalDofIdx,
1612 [gIdx = ectx.globalDofIdx, &tM](
const unsigned tracerIdx)
1613 {
return tM.solTracerConcentration(tracerIdx, gIdx); });
1616 Entry{[&flowsInf = this->simulator_.problem().model().linearizer().getFlowsInfo(),
1617 &flowsC = this->
flowsC_](
const Context& ectx)
1619 const auto gas_idx = Indices::gasEnabled ?
1620 conti0EqIdx + FluidSystem::canonicalToActiveCompIdx(gasCompIdx) : -1;
1621 const auto oil_idx = Indices::oilEnabled ?
1622 conti0EqIdx + FluidSystem::canonicalToActiveCompIdx(oilCompIdx) : -1;
1623 const auto water_idx = Indices::waterEnabled ?
1624 conti0EqIdx + FluidSystem::canonicalToActiveCompIdx(waterCompIdx) : -1;
1625 const auto& flowsInfos = flowsInf[ectx.globalDofIdx];
1626 for (
const auto& flowsInfo : flowsInfos) {
1627 flowsC.assignFlows(ectx.globalDofIdx,
1630 value_or_zero(gas_idx, flowsInfo.flow),
1631 value_or_zero(oil_idx, flowsInfo.flow),
1632 value_or_zero(water_idx, flowsInfo.flow));
1634 }, !this->simulator_.problem().model().linearizer().getFlowsInfo().empty()
1636 Entry{[&floresInf = this->simulator_.problem().model().linearizer().getFloresInfo(),
1637 &flowsC = this->
flowsC_](
const Context& ectx)
1639 const auto gas_idx = Indices::gasEnabled ?
1640 conti0EqIdx + FluidSystem::canonicalToActiveCompIdx(gasCompIdx) : -1;
1641 const auto oil_idx = Indices::oilEnabled ?
1642 conti0EqIdx + FluidSystem::canonicalToActiveCompIdx(oilCompIdx) : -1;
1643 const auto water_idx = Indices::waterEnabled ?
1644 conti0EqIdx + FluidSystem::canonicalToActiveCompIdx(waterCompIdx) : -1;
1645 const auto& floresInfos = floresInf[ectx.globalDofIdx];
1646 for (
const auto& floresInfo : floresInfos) {
1647 flowsC.assignFlores(ectx.globalDofIdx,
1650 value_or_zero(gas_idx, floresInfo.flow),
1651 value_or_zero(oil_idx, floresInfo.flow),
1652 value_or_zero(water_idx, floresInfo.flow));
1654 }, !this->simulator_.problem().model().linearizer().getFloresInfo().empty()
1661 Entry{ScalarEntry{&this->
rv_,
1662 [&problem = this->simulator_.problem()](
const Context& ectx)
1663 {
return problem.initialFluidState(ectx.globalDofIdx).Rv(); }
1665 simulator_.episodeIndex() < 0 &&
1666 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1667 FluidSystem::phaseIsActive(gasPhaseIdx)
1669 Entry{ScalarEntry{&this->
rs_,
1670 [&problem = this->simulator_.problem()](
const Context& ectx)
1671 {
return problem.initialFluidState(ectx.globalDofIdx).Rs(); }
1673 simulator_.episodeIndex() < 0 &&
1674 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1675 FluidSystem::phaseIsActive(gasPhaseIdx)
1677 Entry{ScalarEntry{&this->
rsw_,
1678 [&problem = this->simulator_.problem()](
const Context& ectx)
1679 {
return problem.initialFluidState(ectx.globalDofIdx).Rsw(); }
1681 simulator_.episodeIndex() < 0 &&
1682 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1683 FluidSystem::phaseIsActive(gasPhaseIdx)
1685 Entry{ScalarEntry{&this->
rvw_,
1686 [&problem = this->simulator_.problem()](
const Context& ectx)
1687 {
return problem.initialFluidState(ectx.globalDofIdx).Rvw(); }
1689 simulator_.episodeIndex() < 0 &&
1690 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1691 FluidSystem::phaseIsActive(gasPhaseIdx)
1695 [&problem = this->simulator_.problem()](
const unsigned phase,
1696 const Context& ectx)
1698 const auto& fsInitial = problem.initialFluidState(ectx.globalDofIdx);
1699 return FluidSystem::density(fsInitial,
1701 ectx.intQuants.pvtRegionIndex());
1704 simulator_.episodeIndex() < 0 &&
1705 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1706 FluidSystem::phaseIsActive(gasPhaseIdx)
1708 Entry{PhaseEntry{&this->
invB_,
1709 [&problem = this->simulator_.problem()](
const unsigned phase,
1710 const Context& ectx)
1712 const auto& fsInitial = problem.initialFluidState(ectx.globalDofIdx);
1713 return FluidSystem::inverseFormationVolumeFactor(fsInitial,
1715 ectx.intQuants.pvtRegionIndex());
1718 simulator_.episodeIndex() < 0 &&
1719 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1720 FluidSystem::phaseIsActive(gasPhaseIdx)
1723 [&problem = this->simulator_.problem()](
const unsigned phase,
1724 const Context& ectx)
1726 const auto& fsInitial = problem.initialFluidState(ectx.globalDofIdx);
1727 return FluidSystem::viscosity(fsInitial,
1729 ectx.intQuants.pvtRegionIndex());
1732 simulator_.episodeIndex() < 0 &&
1733 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1734 FluidSystem::phaseIsActive(gasPhaseIdx)
1741 if constexpr (getPropValue<TypeTag, Properties::EnableMech>()) {
1742 if (this->
mech_.allocated()) {
1743 this->extractors_.push_back(
1744 Entry{[&mech = this->
mech_,
1745 &model = simulator_.problem().geoMechModel()](
const Context& ectx)
1747 mech.assignDelStress(ectx.globalDofIdx,
1748 model.delstress(ectx.globalDofIdx));
1750 mech.assignDisplacement(ectx.globalDofIdx,
1751 model.disp(ectx.globalDofIdx,
true));
1754 mech.assignFracStress(ectx.globalDofIdx,
1755 model.fractureStress(ectx.globalDofIdx));
1757 mech.assignLinStress(ectx.globalDofIdx,
1758 model.linstress(ectx.globalDofIdx));
1760 mech.assignPotentialForces(ectx.globalDofIdx,
1761 model.mechPotentialForce(ectx.globalDofIdx),
1762 model.mechPotentialPressForce(ectx.globalDofIdx),
1763 model.mechPotentialTempForce(ectx.globalDofIdx));
1765 mech.assignStrain(ectx.globalDofIdx,
1766 model.strain(ectx.globalDofIdx,
true));
1769 mech.assignStress(ectx.globalDofIdx,
1770 model.stress(ectx.globalDofIdx,
true));
1779 void setupBlockExtractors_(
const bool isSubStep,
1780 const int reportStepNum)
1783 using Context =
typename BlockExtractor::Context;
1784 using PhaseEntry =
typename BlockExtractor::PhaseEntry;
1785 using ScalarEntry =
typename BlockExtractor::ScalarEntry;
1787 using namespace std::string_view_literals;
1789 const auto pressure_handler =
1790 Entry{ScalarEntry{std::vector{
"BPR"sv,
"BPRESSUR"sv},
1791 [](
const Context& ectx)
1793 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1794 return getValue(ectx.fs.pressure(oilPhaseIdx));
1796 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1797 return getValue(ectx.fs.pressure(gasPhaseIdx));
1800 return getValue(ectx.fs.pressure(waterPhaseIdx));
1806 const auto handlers = std::array{
1808 Entry{PhaseEntry{std::array{
1809 std::array{
"BWSAT"sv,
"BOSAT"sv,
"BGSAT"sv},
1810 std::array{
"BSWAT"sv,
"BSOIL"sv,
"BSGAS"sv}
1812 [](
const unsigned phaseIdx,
const Context& ectx)
1814 return getValue(ectx.fs.saturation(phaseIdx));
1818 Entry{ScalarEntry{
"BNSAT",
1819 [](
const Context& ectx)
1821 return ectx.intQuants.solventSaturation().value();
1825 Entry{ScalarEntry{std::vector{
"BTCNFHEA"sv,
"BTEMP"sv},
1826 [](
const Context& ectx)
1828 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1829 return getValue(ectx.fs.temperature(oilPhaseIdx));
1831 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1832 return getValue(ectx.fs.temperature(gasPhaseIdx));
1835 return getValue(ectx.fs.temperature(waterPhaseIdx));
1840 Entry{PhaseEntry{std::array{
1841 std::array{
"BWKR"sv,
"BOKR"sv,
"BGKR"sv},
1842 std::array{
"BKRW"sv,
"BKRO"sv,
"BKRG"sv}
1844 [](
const unsigned phaseIdx,
const Context& ectx)
1846 return getValue(ectx.intQuants.relativePermeability(phaseIdx));
1850 Entry{ScalarEntry{
"BKROG",
1851 [&problem = this->simulator_.problem()](
const Context& ectx)
1853 const auto& materialParams =
1854 problem.materialLawParams(ectx.elemCtx,
1857 return getValue(MaterialLaw::template
1858 relpermOilInOilGasSystem<Evaluation>(materialParams,
1863 Entry{ScalarEntry{
"BKROW",
1864 [&problem = this->simulator_.problem()](
const Context& ectx)
1866 const auto& materialParams = problem.materialLawParams(ectx.elemCtx,
1869 return getValue(MaterialLaw::template
1870 relpermOilInOilWaterSystem<Evaluation>(materialParams,
1875 Entry{ScalarEntry{
"BWPC",
1876 [](
const Context& ectx)
1878 return getValue(ectx.fs.pressure(oilPhaseIdx)) -
1879 getValue(ectx.fs.pressure(waterPhaseIdx));
1883 Entry{ScalarEntry{
"BGPC",
1884 [](
const Context& ectx)
1886 return getValue(ectx.fs.pressure(gasPhaseIdx)) -
1887 getValue(ectx.fs.pressure(oilPhaseIdx));
1891 Entry{ScalarEntry{
"BWPR",
1892 [](
const Context& ectx)
1894 return getValue(ectx.fs.pressure(waterPhaseIdx));
1898 Entry{ScalarEntry{
"BGPR",
1899 [](
const Context& ectx)
1901 return getValue(ectx.fs.pressure(gasPhaseIdx));
1905 Entry{PhaseEntry{std::array{
1906 std::array{
"BVWAT"sv,
"BVOIL"sv,
"BVGAS"sv},
1907 std::array{
"BWVIS"sv,
"BOVIS"sv,
"BGVIS"sv}
1909 [](
const unsigned phaseIdx,
const Context& ectx)
1911 return getValue(ectx.fs.viscosity(phaseIdx));
1915 Entry{PhaseEntry{std::array{
1916 std::array{
"BWDEN"sv,
"BODEN"sv,
"BGDEN"sv},
1917 std::array{
"BDENW"sv,
"BDENO"sv,
"BDENG"sv}
1919 [](
const unsigned phaseIdx,
const Context& ectx)
1921 return getValue(ectx.fs.density(phaseIdx));
1925 Entry{ScalarEntry{
"BFLOWI",
1926 [&flowsC = this->
flowsC_](
const Context& ectx)
1928 return flowsC.getFlow(ectx.globalDofIdx, Dir::XPlus, waterCompIdx);
1932 Entry{ScalarEntry{
"BFLOWJ",
1933 [&flowsC = this->
flowsC_](
const Context& ectx)
1935 return flowsC.getFlow(ectx.globalDofIdx, Dir::YPlus, waterCompIdx);
1939 Entry{ScalarEntry{
"BFLOWK",
1940 [&flowsC = this->
flowsC_](
const Context& ectx)
1942 return flowsC.getFlow(ectx.globalDofIdx, Dir::ZPlus, waterCompIdx);
1946 Entry{ScalarEntry{
"BRPV",
1947 [&model = this->simulator_.model()](
const Context& ectx)
1949 return getValue(ectx.intQuants.porosity()) *
1950 model.dofTotalVolume(ectx.globalDofIdx);
1954 Entry{PhaseEntry{std::array{
"BWPV"sv,
"BOPV"sv,
"BGPV"sv},
1955 [&model = this->simulator_.model()](
const unsigned phaseIdx,
1956 const Context& ectx)
1958 return getValue(ectx.fs.saturation(phaseIdx)) *
1959 getValue(ectx.intQuants.porosity()) *
1960 model.dofTotalVolume(ectx.globalDofIdx);
1964 Entry{ScalarEntry{
"BRS",
1965 [](
const Context& ectx)
1967 return getValue(ectx.fs.Rs());
1971 Entry{ScalarEntry{
"BRV",
1972 [](
const Context& ectx)
1974 return getValue(ectx.fs.Rv());
1978 Entry{ScalarEntry{
"BOIP",
1979 [&model = this->simulator_.model()](
const Context& ectx)
1981 return (getValue(ectx.fs.invB(oilPhaseIdx)) *
1982 getValue(ectx.fs.saturation(oilPhaseIdx)) +
1983 getValue(ectx.fs.Rv()) *
1984 getValue(ectx.fs.invB(gasPhaseIdx)) *
1985 getValue(ectx.fs.saturation(gasPhaseIdx))) *
1986 model.dofTotalVolume(ectx.globalDofIdx) *
1987 getValue(ectx.intQuants.porosity());
1991 Entry{ScalarEntry{
"BGIP",
1992 [&model = this->simulator_.model()](
const Context& ectx)
1994 Scalar result = getValue(ectx.fs.invB(gasPhaseIdx)) *
1995 getValue(ectx.fs.saturation(gasPhaseIdx));
1997 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1998 result += getValue(ectx.fs.Rs()) *
1999 getValue(ectx.fs.invB(oilPhaseIdx)) *
2000 getValue(ectx.fs.saturation(oilPhaseIdx));
2003 result += getValue(ectx.fs.Rsw()) *
2004 getValue(ectx.fs.invB(waterPhaseIdx)) *
2005 getValue(ectx.fs.saturation(waterPhaseIdx));
2009 model.dofTotalVolume(ectx.globalDofIdx) *
2010 getValue(ectx.intQuants.porosity());
2014 Entry{ScalarEntry{
"BWIP",
2015 [&model = this->simulator_.model()](
const Context& ectx)
2017 return getValue(ectx.fs.invB(waterPhaseIdx)) *
2018 getValue(ectx.fs.saturation(waterPhaseIdx)) *
2019 model.dofTotalVolume(ectx.globalDofIdx) *
2020 getValue(ectx.intQuants.porosity());
2024 Entry{ScalarEntry{
"BOIPL",
2025 [&model = this->simulator_.model()](
const Context& ectx)
2027 return getValue(ectx.fs.invB(oilPhaseIdx)) *
2028 getValue(ectx.fs.saturation(oilPhaseIdx)) *
2029 model.dofTotalVolume(ectx.globalDofIdx) *
2030 getValue(ectx.intQuants.porosity());
2034 Entry{ScalarEntry{
"BGIPL",
2035 [&model = this->simulator_.model()](
const Context& ectx)
2038 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
2039 result = getValue(ectx.fs.Rs()) *
2040 getValue(ectx.fs.invB(oilPhaseIdx)) *
2041 getValue(ectx.fs.saturation(oilPhaseIdx));
2044 result = getValue(ectx.fs.Rsw()) *
2045 getValue(ectx.fs.invB(waterPhaseIdx)) *
2046 getValue(ectx.fs.saturation(waterPhaseIdx));
2049 model.dofTotalVolume(ectx.globalDofIdx) *
2050 getValue(ectx.intQuants.porosity());
2054 Entry{ScalarEntry{
"BGIPG",
2055 [&model = this->simulator_.model()](
const Context& ectx)
2057 return getValue(ectx.fs.invB(gasPhaseIdx)) *
2058 getValue(ectx.fs.saturation(gasPhaseIdx)) *
2059 model.dofTotalVolume(ectx.globalDofIdx) *
2060 getValue(ectx.intQuants.porosity());
2064 Entry{ScalarEntry{
"BOIPG",
2065 [&model = this->simulator_.model()](
const Context& ectx)
2067 return getValue(ectx.fs.Rv()) *
2068 getValue(ectx.fs.invB(gasPhaseIdx)) *
2069 getValue(ectx.fs.saturation(gasPhaseIdx)) *
2070 model.dofTotalVolume(ectx.globalDofIdx) *
2071 getValue(ectx.intQuants.porosity());
2075 Entry{PhaseEntry{std::array{
"BPPW"sv,
"BPPO"sv,
"BPPG"sv},
2076 [&simConfig = this->
eclState_.getSimulationConfig(),
2077 &grav = this->simulator_.problem().gravity(),
2079 &problem = this->simulator_.problem(),
2080 ®ions = this->
regions_](
const unsigned phaseIdx,
const Context& ectx)
2082 auto phase = RegionPhasePoreVolAverage::Phase{};
2083 phase.ix = phaseIdx;
2092 const auto datum = simConfig.datumDepths()(regions[
"FIPNUM"][ectx.dofIdx] - 1);
2095 const auto region = RegionPhasePoreVolAverage::Region {
2096 ectx.elemCtx.primaryVars(ectx.dofIdx, 0).pvtRegionIndex() + 1
2099 const auto density = regionAvgDensity->value(
"PVTNUM", phase, region);
2101 const auto press = getValue(ectx.fs.pressure(phase.ix));
2102 const auto dz = problem.dofCenterDepth(ectx.globalDofIdx) - datum;
2103 return press - density*dz*grav[GridView::dimensionworld - 1];
2107 Entry{ScalarEntry{
"BAMIP",
2108 [&model = this->simulator_.model()](
const Context& ectx)
2110 const Scalar rhoW = FluidSystem::referenceDensity(waterPhaseIdx,
2111 ectx.intQuants.pvtRegionIndex());
2112 return getValue(ectx.fs.invB(waterPhaseIdx)) *
2113 getValue(ectx.fs.saturation(waterPhaseIdx)) *
2115 model.dofTotalVolume(ectx.globalDofIdx) *
2116 getValue(ectx.intQuants.porosity());
2120 Entry{ScalarEntry{
"BMMIP",
2121 [&model = this->simulator_.model()](
const Context& ectx)
2123 return getValue(ectx.intQuants.microbialConcentration()) *
2124 getValue(ectx.intQuants.porosity()) *
2125 model.dofTotalVolume(ectx.globalDofIdx);
2129 Entry{ScalarEntry{
"BMOIP",
2130 [&model = this->simulator_.model()](
const Context& ectx)
2132 return getValue(ectx.intQuants.oxygenConcentration()) *
2133 getValue(ectx.intQuants.porosity()) *
2134 model.dofTotalVolume(ectx.globalDofIdx);
2138 Entry{ScalarEntry{
"BMUIP",
2139 [&model = this->simulator_.model()](
const Context& ectx)
2141 return getValue(ectx.intQuants.ureaConcentration()) *
2142 getValue(ectx.intQuants.porosity()) *
2143 model.dofTotalVolume(ectx.globalDofIdx) * 1;
2147 Entry{ScalarEntry{
"BMBIP",
2148 [&model = this->simulator_.model()](
const Context& ectx)
2150 return model.dofTotalVolume(ectx.globalDofIdx) *
2151 getValue(ectx.intQuants.biofilmMass());
2155 Entry{ScalarEntry{
"BMCIP",
2156 [&model = this->simulator_.model()](
const Context& ectx)
2158 return model.dofTotalVolume(ectx.globalDofIdx) *
2159 getValue(ectx.intQuants.calciteMass());
2163 Entry{ScalarEntry{
"BGMIP",
2164 [&model = this->simulator_.model()](
const Context& ectx)
2166 Scalar result = getValue(ectx.fs.invB(gasPhaseIdx)) *
2167 getValue(ectx.fs.saturation(gasPhaseIdx));
2169 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
2170 result += getValue(ectx.fs.Rs()) *
2171 getValue(ectx.fs.invB(oilPhaseIdx)) *
2172 getValue(ectx.fs.saturation(oilPhaseIdx));
2175 result += getValue(ectx.fs.Rsw()) *
2176 getValue(ectx.fs.invB(waterPhaseIdx)) *
2177 getValue(ectx.fs.saturation(waterPhaseIdx));
2179 const Scalar rhoG = FluidSystem::referenceDensity(gasPhaseIdx,
2180 ectx.intQuants.pvtRegionIndex());
2182 model.dofTotalVolume(ectx.globalDofIdx) *
2183 getValue(ectx.intQuants.porosity()) *
2188 Entry{ScalarEntry{
"BGMGP",
2189 [&model = this->simulator_.model()](
const Context& ectx)
2191 const Scalar rhoG = FluidSystem::referenceDensity(gasPhaseIdx,
2192 ectx.intQuants.pvtRegionIndex());
2193 return getValue(ectx.fs.invB(gasPhaseIdx)) *
2194 getValue(ectx.fs.saturation(gasPhaseIdx)) *
2195 model.dofTotalVolume(ectx.globalDofIdx) *
2196 getValue(ectx.intQuants.porosity()) *
2201 Entry{ScalarEntry{
"BGMDS",
2202 [&model = this->simulator_.model()](
const Context& ectx)
2205 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
2206 result = getValue(ectx.fs.Rs()) *
2207 getValue(ectx.fs.invB(oilPhaseIdx)) *
2208 getValue(ectx.fs.saturation(oilPhaseIdx));
2211 result = getValue(ectx.fs.Rsw()) *
2212 getValue(ectx.fs.invB(waterPhaseIdx)) *
2213 getValue(ectx.fs.saturation(waterPhaseIdx));
2215 const Scalar rhoG = FluidSystem::referenceDensity(gasPhaseIdx,
2216 ectx.intQuants.pvtRegionIndex());
2218 model.dofTotalVolume(ectx.globalDofIdx) *
2219 getValue(ectx.intQuants.porosity()) *
2224 Entry{ScalarEntry{
"BGMST",
2225 [&model = this->simulator_.model(),
2226 &problem = this->simulator_.problem()](
const Context& ectx)
2228 const auto& scaledDrainageInfo = problem.materialLawManager()
2229 ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2230 const Scalar sg = getValue(ectx.fs.saturation(gasPhaseIdx));
2231 Scalar strandedGas = scaledDrainageInfo.Sgcr;
2232 if (problem.materialLawManager()->enableHysteresis()) {
2233 const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2234 const Scalar krg = getValue(ectx.intQuants.relativePermeability(gasPhaseIdx));
2235 strandedGas = MaterialLaw::strandedGasSaturation(matParams, sg, krg);
2237 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2238 FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2239 : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2240 return (1.0 - xgW) *
2241 model.dofTotalVolume(ectx.globalDofIdx) *
2242 getValue(ectx.intQuants.porosity()) *
2243 getValue(ectx.fs.density(gasPhaseIdx)) *
2244 std::min(strandedGas, sg);
2248 Entry{ScalarEntry{
"BGMUS",
2249 [&model = this->simulator_.model(),
2250 &problem = this->simulator_.problem()](
const Context& ectx)
2252 const auto& scaledDrainageInfo = problem.materialLawManager()
2253 ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2254 const Scalar sg = getValue(ectx.fs.saturation(gasPhaseIdx));
2255 Scalar strandedGas = scaledDrainageInfo.Sgcr;
2256 if (problem.materialLawManager()->enableHysteresis()) {
2257 const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2258 const Scalar krg = getValue(ectx.intQuants.relativePermeability(gasPhaseIdx));
2259 strandedGas = MaterialLaw::strandedGasSaturation(matParams, sg, krg);
2261 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2262 FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2263 : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2264 return (1.0 - xgW) *
2265 model.dofTotalVolume(ectx.globalDofIdx) *
2266 getValue(ectx.intQuants.porosity()) *
2267 getValue(ectx.fs.density(gasPhaseIdx)) *
2268 std::max(Scalar{0.0}, sg - strandedGas);
2272 Entry{ScalarEntry{
"BGMTR",
2273 [&model = this->simulator_.model(),
2274 &problem = this->simulator_.problem()](
const Context& ectx)
2276 const auto& scaledDrainageInfo = problem.materialLawManager()
2277 ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2278 Scalar trappedGas = scaledDrainageInfo.Sgcr;
2279 if (problem.materialLawManager()->enableHysteresis()) {
2280 const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2281 trappedGas = MaterialLaw::trappedGasSaturation(matParams,
true);
2283 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2284 FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2285 : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2286 return (1.0 - xgW) *
2287 model.dofTotalVolume(ectx.globalDofIdx) *
2288 getValue(ectx.intQuants.porosity()) *
2289 getValue(ectx.fs.density(gasPhaseIdx)) *
2290 std::min(trappedGas, getValue(ectx.fs.saturation(gasPhaseIdx)));
2294 Entry{ScalarEntry{
"BGMMO",
2295 [&model = this->simulator_.model(),
2296 &problem = this->simulator_.problem()](
const Context& ectx)
2298 const auto& scaledDrainageInfo = problem.materialLawManager()
2299 ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2300 Scalar trappedGas = scaledDrainageInfo.Sgcr;
2301 if (problem.materialLawManager()->enableHysteresis()) {
2302 const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2303 trappedGas = MaterialLaw::trappedGasSaturation(matParams,
true);
2305 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2306 FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2307 : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2308 return (1.0 - xgW) *
2309 model.dofTotalVolume(ectx.globalDofIdx) *
2310 getValue(ectx.intQuants.porosity()) *
2311 getValue(ectx.fs.density(gasPhaseIdx)) *
2312 std::max(Scalar{0.0}, getValue(ectx.fs.saturation(gasPhaseIdx)) - trappedGas);
2316 Entry{ScalarEntry{
"BGKTR",
2317 [&model = this->simulator_.model(),
2318 &problem = this->simulator_.problem()](
const Context& ectx)
2320 const auto& scaledDrainageInfo = problem.materialLawManager()
2321 ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2322 const Scalar sg = getValue(ectx.fs.saturation(gasPhaseIdx));
2323 Scalar sgcr = scaledDrainageInfo.Sgcr;
2324 if (problem.materialLawManager()->enableHysteresis()) {
2325 const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2326 sgcr = MaterialLaw::trappedGasSaturation(matParams,
false);
2332 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2333 FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2334 : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2335 return (1.0 - xgW) *
2336 model.dofTotalVolume(ectx.globalDofIdx) *
2337 getValue(ectx.intQuants.porosity()) *
2338 getValue(ectx.fs.density(gasPhaseIdx)) *
2339 getValue(ectx.fs.saturation(gasPhaseIdx));
2344 Entry{ScalarEntry{
"BGKMO",
2345 [&model = this->simulator_.model(),
2346 &problem = this->simulator_.problem()](
const Context& ectx)
2348 const auto& scaledDrainageInfo = problem.materialLawManager()
2349 ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2350 const Scalar sg = getValue(ectx.fs.saturation(gasPhaseIdx));
2351 Scalar sgcr = scaledDrainageInfo.Sgcr;
2352 if (problem.materialLawManager()->enableHysteresis()) {
2353 const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2354 sgcr = MaterialLaw::trappedGasSaturation(matParams,
false);
2360 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2361 FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2362 : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2363 return (1.0 - xgW) *
2364 model.dofTotalVolume(ectx.globalDofIdx) *
2365 getValue(ectx.intQuants.porosity()) *
2366 getValue(ectx.fs.density(gasPhaseIdx)) *
2367 getValue(ectx.fs.saturation(gasPhaseIdx));
2372 Entry{ScalarEntry{
"BGCDI",
2373 [&model = this->simulator_.model(),
2374 &problem = this->simulator_.problem()](
const Context& ectx)
2376 const auto& scaledDrainageInfo = problem.materialLawManager()
2377 ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2378 Scalar sgcr = scaledDrainageInfo.Sgcr;
2379 if (problem.materialLawManager()->enableHysteresis()) {
2380 const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2381 sgcr = MaterialLaw::trappedGasSaturation(matParams,
false);
2383 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2384 FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2385 : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2386 return (1.0 - xgW) *
2387 model.dofTotalVolume(ectx.globalDofIdx) *
2388 getValue(ectx.intQuants.porosity()) *
2389 getValue(ectx.fs.density(gasPhaseIdx)) *
2390 std::min(sgcr, getValue(ectx.fs.saturation(gasPhaseIdx))) /
2391 FluidSystem::molarMass(gasCompIdx, ectx.intQuants.pvtRegionIndex());
2395 Entry{ScalarEntry{
"BGCDM",
2396 [&model = this->simulator_.model(),
2397 &problem = this->simulator_.problem()](
const Context& ectx)
2399 const auto& scaledDrainageInfo = problem.materialLawManager()
2400 ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2401 Scalar sgcr = scaledDrainageInfo.Sgcr;
2402 if (problem.materialLawManager()->enableHysteresis()) {
2403 const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2404 sgcr = MaterialLaw::trappedGasSaturation(matParams,
false);
2406 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2407 FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2408 : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2409 return (1.0 - xgW) *
2410 model.dofTotalVolume(ectx.globalDofIdx) *
2411 getValue(ectx.intQuants.porosity()) *
2412 getValue(ectx.fs.density(gasPhaseIdx)) *
2413 std::max(Scalar{0.0}, getValue(ectx.fs.saturation(gasPhaseIdx)) - sgcr) /
2414 FluidSystem::molarMass(gasCompIdx, ectx.intQuants.pvtRegionIndex());
2418 Entry{ScalarEntry{
"BGKDI",
2419 [&model = this->simulator_.model(),
2420 &problem = this->simulator_.problem()](
const Context& ectx)
2422 const auto& scaledDrainageInfo = problem.materialLawManager()
2423 ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2424 const Scalar sg = getValue(ectx.fs.saturation(gasPhaseIdx));
2425 Scalar sgcr = scaledDrainageInfo.Sgcr;
2426 if (problem.materialLawManager()->enableHysteresis()) {
2427 const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2428 sgcr = MaterialLaw::trappedGasSaturation(matParams,
false);
2434 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2435 FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2436 : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2437 return (1.0 - xgW) *
2438 model.dofTotalVolume(ectx.globalDofIdx) *
2439 getValue(ectx.intQuants.porosity()) *
2440 getValue(ectx.fs.density(gasPhaseIdx)) *
2441 getValue(ectx.fs.saturation(gasPhaseIdx)) /
2442 FluidSystem::molarMass(gasCompIdx, ectx.intQuants.pvtRegionIndex());
2447 Entry{ScalarEntry{
"BGKDM",
2448 [&model = this->simulator_.model(),
2449 &problem = this->simulator_.problem()](
const Context& ectx)
2451 const auto& scaledDrainageInfo = problem.materialLawManager()
2452 ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2453 const Scalar sg = getValue(ectx.fs.saturation(gasPhaseIdx));
2454 Scalar sgcr = scaledDrainageInfo.Sgcr;
2455 if (problem.materialLawManager()->enableHysteresis()) {
2456 const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2457 sgcr = MaterialLaw::trappedGasSaturation(matParams,
false);
2463 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2464 FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2465 : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2466 return (1.0 - xgW) *
2467 model.dofTotalVolume(ectx.globalDofIdx) *
2468 getValue(ectx.intQuants.porosity()) *
2469 getValue(ectx.fs.density(gasPhaseIdx)) *
2470 getValue(ectx.fs.saturation(gasPhaseIdx)) /
2471 FluidSystem::molarMass(gasCompIdx, ectx.intQuants.pvtRegionIndex());
2476 Entry{ScalarEntry{
"BWCD",
2477 [&model = this->simulator_.model()](
const Context& ectx)
2480 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
2481 result = getValue(ectx.fs.Rs()) *
2482 getValue(ectx.fs.invB(oilPhaseIdx)) *
2483 getValue(ectx.fs.saturation(oilPhaseIdx));
2486 result = getValue(ectx.fs.Rsw()) *
2487 getValue(ectx.fs.invB(waterPhaseIdx)) *
2488 getValue(ectx.fs.saturation(waterPhaseIdx));
2490 const Scalar rhoG = FluidSystem::referenceDensity(gasPhaseIdx,
2491 ectx.intQuants.pvtRegionIndex());
2493 model.dofTotalVolume(ectx.globalDofIdx) *
2494 getValue(ectx.intQuants.porosity()) *
2496 FluidSystem::molarMass(gasCompIdx, ectx.intQuants.pvtRegionIndex());
2500 Entry{ScalarEntry{
"BWIPG",
2501 [&model = this->simulator_.model()](
const Context& ectx)
2503 Scalar result = 0.0;
2504 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
2505 result = getValue(ectx.fs.Rvw()) *
2506 getValue(ectx.fs.invB(gasPhaseIdx)) *
2507 getValue(ectx.fs.saturation(gasPhaseIdx));
2510 model.dofTotalVolume(ectx.globalDofIdx) *
2511 getValue(ectx.intQuants.porosity());
2515 Entry{ScalarEntry{
"BWIPL",
2516 [&model = this->simulator_.model()](
const Context& ectx)
2518 return getValue(ectx.fs.invB(waterPhaseIdx)) *
2519 getValue(ectx.fs.saturation(waterPhaseIdx)) *
2520 model.dofTotalVolume(ectx.globalDofIdx) *
2521 getValue(ectx.intQuants.porosity());
2530 if (reportStepNum > 0 && !isSubStep) {
2532 const auto& rpt = this->
schedule_[reportStepNum - 1].rpt_config.get();
2533 if (rpt.contains(
"WELLS") && rpt.at(
"WELLS") > 1) {
2535 [&c = this->collectOnIORank_](
const int idx)
2536 {
return c.isCartIdxOnThisRank(idx); });
2538 const auto extraHandlers = std::array{
2547 const Simulator& simulator_;
2548 const CollectDataOnIORankType& collectOnIORank_;
2549 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:435
std::array< ScalarBuffer, numPhases > relativePermeability_
Definition: GenericOutputBlackoilModule.hpp:424
ScalarBuffer fluidPressure_
Definition: GenericOutputBlackoilModule.hpp:378
std::array< ScalarBuffer, numPhases > density_
Definition: GenericOutputBlackoilModule.hpp:422
ScalarBuffer saturatedOilFormationVolumeFactor_
Definition: GenericOutputBlackoilModule.hpp:410
ScalarBuffer overburdenPressure_
Definition: GenericOutputBlackoilModule.hpp:384
ScalarBuffer gasDissolutionFactorInWater_
Definition: GenericOutputBlackoilModule.hpp:404
const EclipseState & eclState_
Definition: GenericOutputBlackoilModule.hpp:337
ScalarBuffer swmin_
Definition: GenericOutputBlackoilModule.hpp:400
ScalarBuffer rockCompPorvMultiplier_
Definition: GenericOutputBlackoilModule.hpp:408
MICPContainer< Scalar > micpC_
Definition: GenericOutputBlackoilModule.hpp:412
RFTContainer< GetPropType< TypeTag, Properties::FluidSystem > > rftC_
Definition: GenericOutputBlackoilModule.hpp:432
bool computeFip_
Definition: GenericOutputBlackoilModule.hpp:360
ScalarBuffer dewPointPressure_
Definition: GenericOutputBlackoilModule.hpp:407
LogOutputHelper< Scalar > logOutput_
Definition: GenericOutputBlackoilModule.hpp:344
std::vector< int > failedCellsPb_
Definition: GenericOutputBlackoilModule.hpp:369
ScalarBuffer permFact_
Definition: GenericOutputBlackoilModule.hpp:393
ScalarBuffer rsw_
Definition: GenericOutputBlackoilModule.hpp:381
ScalarBuffer pcog_
Definition: GenericOutputBlackoilModule.hpp:415
std::optional< RegionPhasePoreVolAverage > regionAvgDensity_
Definition: GenericOutputBlackoilModule.hpp:443
std::array< ScalarBuffer, numPhases > invB_
Definition: GenericOutputBlackoilModule.hpp:421
ScalarBuffer pSalt_
Definition: GenericOutputBlackoilModule.hpp:392
ScalarBuffer cFoam_
Definition: GenericOutputBlackoilModule.hpp:390
ScalarBuffer bubblePointPressure_
Definition: GenericOutputBlackoilModule.hpp:406
ScalarBuffer temperature_
Definition: GenericOutputBlackoilModule.hpp:379
ScalarBuffer ppcw_
Definition: GenericOutputBlackoilModule.hpp:401
FIPContainer< GetPropType< TypeTag, Properties::FluidSystem > > fipC_
Definition: GenericOutputBlackoilModule.hpp:362
ScalarBuffer rockCompTransMultiplier_
Definition: GenericOutputBlackoilModule.hpp:411
MechContainer< Scalar > mech_
Definition: GenericOutputBlackoilModule.hpp:418
ScalarBuffer dynamicPoreVolume_
Definition: GenericOutputBlackoilModule.hpp:376
ScalarBuffer minimumOilPressure_
Definition: GenericOutputBlackoilModule.hpp:409
ScalarBuffer gasFormationVolumeFactor_
Definition: GenericOutputBlackoilModule.hpp:372
std::array< ScalarBuffer, numPhases > residual_
Definition: GenericOutputBlackoilModule.hpp:428
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:398
const Schedule & schedule_
Definition: GenericOutputBlackoilModule.hpp:338
FlowsContainer< GetPropType< TypeTag, Properties::FluidSystem > > flowsC_
Definition: GenericOutputBlackoilModule.hpp:430
ExtboContainer< Scalar > extboC_
Definition: GenericOutputBlackoilModule.hpp:394
void setupExtraBlockData(const std::size_t reportStepNum, std::function< bool(int)> isCartIdxOnThisRank)
ScalarBuffer oilSaturationPressure_
Definition: GenericOutputBlackoilModule.hpp:385
InterRegFlowMap interRegionFlows_
Definition: GenericOutputBlackoilModule.hpp:343
ScalarBuffer pcgw_
Definition: GenericOutputBlackoilModule.hpp:413
ScalarBuffer cPolymer_
Definition: GenericOutputBlackoilModule.hpp:389
void setupBlockData(std::function< bool(int)> isCartIdxOnThisRank)
ScalarBuffer rvw_
Definition: GenericOutputBlackoilModule.hpp:383
std::array< ScalarBuffer, numPhases > saturation_
Definition: GenericOutputBlackoilModule.hpp:420
std::unordered_map< std::string, std::vector< int > > regions_
Definition: GenericOutputBlackoilModule.hpp:363
ScalarBuffer rPorV_
Definition: GenericOutputBlackoilModule.hpp:377
ScalarBuffer oilVaporizationFactor_
Definition: GenericOutputBlackoilModule.hpp:403
std::vector< int > failedCellsPd_
Definition: GenericOutputBlackoilModule.hpp:370
ScalarBuffer rs_
Definition: GenericOutputBlackoilModule.hpp:380
ScalarBuffer drsdtcon_
Definition: GenericOutputBlackoilModule.hpp:386
ScalarBuffer sSol_
Definition: GenericOutputBlackoilModule.hpp:387
std::map< std::pair< std::string, int >, double > extraBlockData_
Definition: GenericOutputBlackoilModule.hpp:438
ScalarBuffer pressureTimesPoreVolume_
Definition: GenericOutputBlackoilModule.hpp:374
ScalarBuffer gasDissolutionFactor_
Definition: GenericOutputBlackoilModule.hpp:402
std::array< ScalarBuffer, numPhases > viscosity_
Definition: GenericOutputBlackoilModule.hpp:423
bool forceDisableFipOutput_
Definition: GenericOutputBlackoilModule.hpp:358
ScalarBuffer soMax_
Definition: GenericOutputBlackoilModule.hpp:395
ScalarBuffer sgmax_
Definition: GenericOutputBlackoilModule.hpp:397
ScalarBuffer somin_
Definition: GenericOutputBlackoilModule.hpp:399
ScalarBuffer hydrocarbonPoreVolume_
Definition: GenericOutputBlackoilModule.hpp:373
const std::optional< Inplace > & initialInplace() const
Definition: GenericOutputBlackoilModule.hpp:225
ScalarBuffer waterVaporizationFactor_
Definition: GenericOutputBlackoilModule.hpp:405
ScalarBuffer cSalt_
Definition: GenericOutputBlackoilModule.hpp:391
TracerContainer< GetPropType< TypeTag, Properties::FluidSystem > > tracerC_
Definition: GenericOutputBlackoilModule.hpp:426
ScalarBuffer rv_
Definition: GenericOutputBlackoilModule.hpp:382
ScalarBuffer pcow_
Definition: GenericOutputBlackoilModule.hpp:414
ScalarBuffer swMax_
Definition: GenericOutputBlackoilModule.hpp:396
ScalarBuffer pressureTimesHydrocarbonVolume_
Definition: GenericOutputBlackoilModule.hpp:375
ScalarBuffer rswSol_
Definition: GenericOutputBlackoilModule.hpp:388
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:236
void initializeFluxData()
Prepare for capturing connection fluxes, particularly to account for inter-region flows.
Definition: OutputBlackoilModule.hpp:480
void setupExtractors(const bool isSubStep, const int reportStepNum)
Setup list of active element-level data extractors.
Definition: OutputBlackoilModule.hpp:217
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:195
void processFluxes(const ElementContext &elemCtx, ActiveIndex &&activeIndex, CartesianIndex &&cartesianIndex)
Capture connection fluxes, particularly to account for inter-region flows.
Definition: OutputBlackoilModule.hpp:443
void clearExtractors()
Clear list of active element-level data extractors.
Definition: OutputBlackoilModule.hpp:225
void outputFipAndResvLogToCSV(const std::size_t reportStepNum, const bool substep, const Parallel::Communication &comm)
Definition: OutputBlackoilModule.hpp:377
void assignToFluidState(FluidState &fs, unsigned elemIdx) const
Definition: OutputBlackoilModule.hpp:504
void initHysteresisParams(Simulator &simulator, unsigned elemIdx) const
Definition: OutputBlackoilModule.hpp:554
void updateFluidInPlace(const ElementContext &elemCtx)
Definition: OutputBlackoilModule.hpp:619
OutputBlackOilModule(const Simulator &simulator, const SummaryConfig &smryCfg, const CollectDataOnIORankType &collectOnIORank)
Definition: OutputBlackoilModule.hpp:132
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:326
const InterRegFlowMap & getInterRegFlows() const
Get read-only access to collection of inter-region flows.
Definition: OutputBlackoilModule.hpp:498
void processElementBlockData(const ElementContext &elemCtx)
Definition: OutputBlackoilModule.hpp:282
void finalizeFluxData()
Finalize capturing connection fluxes.
Definition: OutputBlackoilModule.hpp:490
void updateFluidInPlace(const unsigned globalDofIdx, const IntensiveQuantities &intQuants, const double totVolume)
Definition: OutputBlackoilModule.hpp:626
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