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;
118 template<
int idx,
class VectorType>
119 static Scalar value_or_zero(
const VectorType&
v)
121 if constexpr (
idx == -1) {
124 return v.empty() ? 0.0 :
v[
idx];
132 :
BaseType(simulator.vanguard().eclState(),
133 simulator.vanguard().schedule(),
135 simulator.vanguard().summaryState(),
137 [
this](
const int idx)
138 {
return simulator_.problem().eclWriter().collectOnIORank().localIdxToGlobalIdx(
idx); },
139 simulator.vanguard().grid().comm(),
150 , simulator_(simulator)
151 , collectOnIORank_(collectOnIORank)
157 auto isCartIdxOnThisRank = [&collectOnIORank](
const int idx) {
158 return collectOnIORank.isCartIdxOnThisRank(
idx);
161 this->setupBlockData(isCartIdxOnThisRank);
163 if (! Parameters::Get<Parameters::OwnerCellsFirst>()) {
164 const std::string
msg =
"The output code does not support --owner-cells-first=false.";
165 if (collectOnIORank.isIORank()) {
172 auto rset = this->eclState_.fieldProps().fip_regions();
173 rset.push_back(
"PVTNUM");
178 this->regionAvgDensity_
179 .emplace(this->simulator_.gridView().comm(),
180 FluidSystem::numPhases,
rset,
181 [
fp = std::cref(
this->eclState_.fieldProps())]
182 (
const std::string&
rsetName) ->
decltype(
auto)
183 { return fp.get().get_int(rsetName); });
193 const unsigned reportStepNum,
196 const bool isRestart)
202 const auto& problem = this->simulator_.problem();
209 &problem.materialLawManager()->hysteresisConfig(),
210 problem.eclWriter().getOutputNnc().size());
215 const int reportStepNum)
217 this->setupElementExtractors_();
218 this->setupBlockExtractors_(
isSubStep, reportStepNum);
224 this->extractors_.clear();
225 this->blockExtractors_.clear();
226 this->extraBlockExtractors_.clear();
240 if (this->extractors_.empty()) {
244 const auto&
matLawManager = simulator_.problem().materialLawManager();
247 for (
unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(0); ++dofIdx) {
248 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, 0);
249 const auto& fs = intQuants.fluidState();
252 elemCtx.globalSpaceIndex(dofIdx, 0),
253 elemCtx.primaryVars(dofIdx, 0).pvtRegionIndex(),
261 if (FluidSystem::phaseIsActive(oilPhaseIdx) && FluidSystem::phaseIsActive(waterPhaseIdx)) {
267 if (FluidSystem::phaseIsActive(oilPhaseIdx) && FluidSystem::phaseIsActive(gasPhaseIdx)) {
279 void processElementBlockData(
const ElementContext& elemCtx)
286 if (this->blockExtractors_.empty() &&
this->extraBlockExtractors_.empty()) {
290 for (
unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(0); ++dofIdx) {
292 const auto globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, 0);
293 const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx);
297 if (
be_it == this->blockExtractors_.end() &&
303 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, 0);
304 const auto& fs = intQuants.fluidState();
306 const typename BlockExtractor::Context
ectx{
314 if (
be_it != this->blockExtractors_.end()) {
317 if (
bee_it != this->extraBlockExtractors_.end()) {
351 template <
class ActiveIndex,
class CartesianIndex>
357 const auto identifyCell = [&activeIndex, &cartesianIndex](
const Element&
elem)
365 elem.partitionType() == Dune::InteriorEntity
370 const auto& stencil = elemCtx.stencil(
timeIdx);
371 const auto numInteriorFaces = elemCtx.numInteriorFaces(
timeIdx);
378 const auto rates = this->
393 this->interRegionFlows_.
clear();
409 return this->interRegionFlows_;
412 template <
class Flu
idState>
413 void assignToFluidState(FluidState& fs,
unsigned elemIdx)
const
416 if (this->saturation_[
phaseIdx].empty())
422 if (!this->fluidPressure_.empty()) {
425 std::array<Scalar, numPhases>
pc = {0};
426 const MaterialLawParams&
matParams = simulator_.problem().materialLawParams(elemIdx);
427 MaterialLaw::capillaryPressures(
pc,
matParams, fs);
428 Valgrind::CheckDefined(this->fluidPressure_[elemIdx]);
429 Valgrind::CheckDefined(
pc);
430 const auto& pressure = this->fluidPressure_[elemIdx];
432 if (!FluidSystem::phaseIsActive(
phaseIdx))
435 if (Indices::oilEnabled)
437 else if (Indices::gasEnabled)
439 else if (Indices::waterEnabled)
445 if (!this->temperature_.empty())
446 fs.setTemperature(this->temperature_[elemIdx]);
447 if (!this->rs_.empty())
448 fs.setRs(this->rs_[elemIdx]);
449 if (!this->rsw_.empty())
450 fs.setRsw(this->rsw_[elemIdx]);
451 if (!this->rv_.empty())
452 fs.setRv(this->rv_[elemIdx]);
453 if (!this->rvw_.empty())
454 fs.setRvw(this->rvw_[elemIdx]);
457 void initHysteresisParams(Simulator& simulator,
unsigned elemIdx)
const
459 if (!this->soMax_.empty())
460 simulator.problem().setMaxOilSaturation(elemIdx, this->soMax_[elemIdx]);
462 if (simulator.problem().materialLawManager()->enableHysteresis()) {
463 auto matLawManager = simulator.problem().materialLawManager();
465 if (FluidSystem::phaseIsActive(oilPhaseIdx)
466 && FluidSystem::phaseIsActive(waterPhaseIdx)) {
472 if (!this->soMax_.empty()) {
473 somax = this->soMax_[elemIdx];
477 if (!this->swMax_.empty()) {
478 swmax = this->swMax_[elemIdx];
482 if (!this->swmin_.empty()) {
483 swmin = this->swmin_[elemIdx];
487 somax, swmax, swmin, elemIdx);
489 if (FluidSystem::phaseIsActive(oilPhaseIdx)
490 && FluidSystem::phaseIsActive(gasPhaseIdx)) {
496 if (!this->sgmax_.empty()) {
497 sgmax = this->sgmax_[elemIdx];
501 if (!this->shmax_.empty()) {
502 shmax = this->shmax_[elemIdx];
506 if (!this->somin_.empty()) {
507 somin = this->somin_[elemIdx];
511 sgmax, shmax, somin, elemIdx);
516 if (simulator_.vanguard().eclState().fieldProps().has_double(
"SWATINIT")) {
517 simulator.problem().materialLawManager()
518 ->applyRestartSwatInit(elemIdx, this->ppcw_[elemIdx]);
522 void updateFluidInPlace(
const ElementContext& elemCtx)
524 for (
unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(0); ++dofIdx) {
525 updateFluidInPlace_(elemCtx, dofIdx);
529 void updateFluidInPlace(
const unsigned globalDofIdx,
530 const IntensiveQuantities& intQuants,
533 this->updateFluidInPlace_(globalDofIdx, intQuants,
totVolume);
537 template <
typename T>
538 using RemoveCVR = std::remove_cv_t<std::remove_reference_t<T>>;
540 template <
typename,
class =
void>
541 struct HasGeoMech :
public std::false_type {};
543 template <
typename Problem>
545 Problem, std::
void_t<decltype(std::declval<Problem>().geoMechModel())>
546 > :
public std::true_type {};
548 bool isDefunctParallelWell(std::string
wname)
const override
550 if (simulator_.gridView().comm().size() == 1)
552 const auto& parallelWells = simulator_.vanguard().parallelWells();
553 std::pair<std::string, bool> value {
wname,
true};
554 auto candidate = std::lower_bound(parallelWells.begin(), parallelWells.end(), value);
558 void updateFluidInPlace_(
const ElementContext& elemCtx,
const unsigned dofIdx)
560 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, 0);
561 const unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, 0);
562 const auto totVolume = elemCtx.simulator().model().dofTotalVolume(globalDofIdx);
564 this->updateFluidInPlace_(globalDofIdx, intQuants,
totVolume);
567 void updateFluidInPlace_(
const unsigned globalDofIdx,
568 const IntensiveQuantities& intQuants,
573 this->updateTotalVolumesAndPressures_(globalDofIdx, intQuants,
totVolume);
575 if (this->computeFip_) {
576 this->updatePhaseInplaceVolumes_(globalDofIdx, intQuants,
totVolume);
580 void createLocalRegion_(std::vector<int>& region)
582 std::size_t elemIdx = 0;
583 for (
const auto&
elem :
elements(simulator_.gridView())) {
584 if (
elem.partitionType() != Dune::InteriorEntity) {
592 template <
typename Flu
idState>
593 void aggregateAverageDensityContributions_(
const FluidState& fs,
594 const unsigned int globalDofIdx,
597 auto pvCellValue = RegionPhasePoreVolAverage::CellValue{};
600 for (
auto phaseIdx = 0*FluidSystem::numPhases;
603 if (! FluidSystem::phaseIsActive(
phaseIdx)) {
610 this->regionAvgDensity_
611 ->addCell(globalDofIdx,
612 RegionPhasePoreVolAverage::Phase {
phaseIdx },
632 data::InterRegFlowMap::FlowRates
633 getComponentSurfaceRates(
const ElementContext& elemCtx,
634 const Scalar faceArea,
636 const std::size_t
timeIdx)
const
638 using Component = data::InterRegFlowMap::Component;
640 auto rates = data::InterRegFlowMap::FlowRates {};
646 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
647 const auto&
up = elemCtx
650 const auto pvtReg =
up.pvtRegionIndex();
653 (
up.fluidState(), oilPhaseIdx,
pvtReg));
657 rates[Component::Oil] +=
qO;
659 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
661 BlackOil::getRs_<FluidSystem, FluidState, Scalar>
664 rates[Component::Gas] +=
qO * Rs;
665 rates[Component::Disgas] +=
qO * Rs;
669 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
670 const auto&
up = elemCtx
673 const auto pvtReg =
up.pvtRegionIndex();
676 (
up.fluidState(), gasPhaseIdx,
pvtReg));
680 rates[Component::Gas] +=
qG;
682 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
684 BlackOil::getRv_<FluidSystem, FluidState, Scalar>
687 rates[Component::Oil] +=
qG * Rv;
688 rates[Component::Vapoil] +=
qG * Rv;
692 if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
693 const auto&
up = elemCtx
696 const auto pvtReg =
up.pvtRegionIndex();
699 (
up.fluidState(), waterPhaseIdx,
pvtReg));
701 rates[Component::Water] +=
708 template <
typename Flu
idState>
709 Scalar hydroCarbonFraction(
const FluidState& fs)
const
711 if (this->eclState_.runspec().co2Storage()) {
719 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
723 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
730 void updateTotalVolumesAndPressures_(
const unsigned globalDofIdx,
731 const IntensiveQuantities& intQuants,
734 const auto& fs = intQuants.fluidState();
736 const double pv =
totVolume * intQuants.porosity().value();
737 const auto hydrocarbon = this->hydroCarbonFraction(fs);
739 if (! this->hydrocarbonPoreVolume_.empty()) {
740 this->fipC_.assignPoreVolume(globalDofIdx,
741 totVolume * intQuants.referencePorosity());
743 this->dynamicPoreVolume_[globalDofIdx] = pv;
744 this->hydrocarbonPoreVolume_[globalDofIdx] = pv *
hydrocarbon;
747 if (!this->pressureTimesHydrocarbonVolume_.empty() &&
748 !
this->pressureTimesPoreVolume_.empty())
750 assert(this->hydrocarbonPoreVolume_.size() ==
this->pressureTimesHydrocarbonVolume_.size());
751 assert(this->fipC_.get(Inplace::Phase::PoreVolume).size() ==
this->pressureTimesPoreVolume_.size());
753 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
754 this->pressureTimesPoreVolume_[globalDofIdx] =
755 getValue(fs.pressure(oilPhaseIdx)) * pv;
757 this->pressureTimesHydrocarbonVolume_[globalDofIdx] =
758 this->pressureTimesPoreVolume_[globalDofIdx] *
hydrocarbon;
760 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
761 this->pressureTimesPoreVolume_[globalDofIdx] =
762 getValue(fs.pressure(gasPhaseIdx)) * pv;
764 this->pressureTimesHydrocarbonVolume_[globalDofIdx] =
765 this->pressureTimesPoreVolume_[globalDofIdx] *
hydrocarbon;
767 else if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
768 this->pressureTimesPoreVolume_[globalDofIdx] =
769 getValue(fs.pressure(waterPhaseIdx)) * pv;
774 void updatePhaseInplaceVolumes_(
const unsigned globalDofIdx,
775 const IntensiveQuantities& intQuants,
778 std::array<Scalar, FluidSystem::numPhases> fip {};
779 std::array<Scalar, FluidSystem::numPhases>
fipr{};
781 const auto& fs = intQuants.fluidState();
782 const auto pv =
totVolume * intQuants.porosity().value();
785 if (!FluidSystem::phaseIsActive(
phaseIdx)) {
796 this->fipC_.assignVolumesSurface(globalDofIdx, fip);
797 this->fipC_.assignVolumesReservoir(globalDofIdx,
798 fs.saltConcentration().value(),
801 if (FluidSystem::phaseIsActive(oilPhaseIdx) &&
802 FluidSystem::phaseIsActive(gasPhaseIdx))
804 this->updateOilGasDistribution(globalDofIdx, fs, fip);
807 if (FluidSystem::phaseIsActive(waterPhaseIdx) &&
808 FluidSystem::phaseIsActive(gasPhaseIdx))
810 this->updateGasWaterDistribution(globalDofIdx, fs, fip);
813 if (FluidSystem::phaseIsActive(gasPhaseIdx) &&
814 this->fipC_.hasCo2InGas())
816 this->updateCO2InGas(globalDofIdx, pv, intQuants);
819 if (this->fipC_.hasCo2InWater() &&
820 (FluidSystem::phaseIsActive(waterPhaseIdx) ||
821 FluidSystem::phaseIsActive(oilPhaseIdx)))
823 this->updateCO2InWater(globalDofIdx, pv, fs);
827 template <
typename Flu
idState,
typename FIPArray>
828 void updateOilGasDistribution(
const unsigned globalDofIdx,
829 const FluidState& fs,
839 template <
typename Flu
idState,
typename FIPArray>
840 void updateGasWaterDistribution(
const unsigned globalDofIdx,
841 const FluidState& fs,
851 template <
typename IntensiveQuantities>
852 void updateCO2InGas(
const unsigned globalDofIdx,
854 const IntensiveQuantities& intQuants)
857 ->oilWaterScaledEpsInfoDrainage(globalDofIdx);
859 const auto& fs = intQuants.fluidState();
861 if (this->simulator_.problem().materialLawManager()->enableHysteresis()) {
862 const auto&
matParams = simulator_.problem().materialLawParams(globalDofIdx);
863 sgcr = MaterialLaw::trappedGasSaturation(
matParams,
false);
867 if (this->fipC_.has(Inplace::Phase::CO2MassInGasPhaseMaximumTrapped) ||
868 this->fipC_.has(Inplace::Phase::CO2MassInGasPhaseMaximumUnTrapped))
870 if (this->simulator_.problem().materialLawManager()->enableHysteresis()) {
871 const auto&
matParams = simulator_.problem().materialLawParams(globalDofIdx);
877 const Scalar sg =
getValue(fs.saturation(gasPhaseIdx));
879 if (this->fipC_.has(Inplace::Phase::CO2MassInGasPhaseEffectiveTrapped) ||
880 this->fipC_.has(Inplace::Phase::CO2MassInGasPhaseEffectiveUnTrapped))
882 if (this->simulator_.problem().materialLawManager()->enableHysteresis()) {
883 const auto&
matParams = simulator_.problem().materialLawParams(globalDofIdx);
884 const double krg =
getValue(intQuants.relativePermeability(gasPhaseIdx));
889 const typename FIPContainer<FluidSystem>::Co2InGasInput
v{
894 FluidSystem::phaseIsActive(waterPhaseIdx)
895 ? FluidSystem::convertRvwToXgW(
getValue(fs.Rvw()), fs.pvtRegionIndex())
897 FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex()),
902 this->fipC_.assignCo2InGas(globalDofIdx,
v);
905 template <
typename Flu
idState>
906 void updateCO2InWater(
const unsigned globalDofIdx,
908 const FluidState& fs)
910 const auto co2InWater = FluidSystem::phaseIsActive(oilPhaseIdx)
911 ? this->co2InWaterFromOil(fs, pv)
912 :
this->co2InWaterFromWater(fs, pv);
914 const Scalar mM = FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex());
916 this->fipC_.assignCo2InWater(globalDofIdx,
co2InWater, mM);
919 template <
typename Flu
idState>
920 Scalar co2InWaterFromWater(
const FluidState& fs,
const double pv)
const
922 const double rhow =
getValue(fs.density(waterPhaseIdx));
923 const double sw =
getValue(fs.saturation(waterPhaseIdx));
924 const double xwG = FluidSystem::convertRswToXwG(
getValue(fs.Rsw()), fs.pvtRegionIndex());
926 const Scalar mM = FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex());
931 template <
typename Flu
idState>
932 Scalar co2InWaterFromOil(
const FluidState& fs,
const double pv)
const
935 const double so =
getValue(fs.saturation(oilPhaseIdx));
936 const double xoG = FluidSystem::convertRsToXoG(
getValue(fs.Rs()), fs.pvtRegionIndex());
938 const Scalar mM = FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex());
944 void setupElementExtractors_()
946 using Entry =
typename Extractor::Entry;
947 using Context =
typename Extractor::Context;
948 using ScalarEntry =
typename Extractor::ScalarEntry;
949 using PhaseEntry =
typename Extractor::PhaseEntry;
951 const bool hasResidual = simulator_.model().linearizer().residual().size() > 0;
952 const auto&
hysteresisConfig = simulator_.problem().materialLawManager()->hysteresisConfig();
955 Entry{PhaseEntry{&this->saturation_,
956 [](
const unsigned phase,
const Context&
ectx)
960 Entry{PhaseEntry{&this->invB_,
961 [](
const unsigned phase,
const Context&
ectx)
965 Entry{PhaseEntry{&this->density_,
966 [](
const unsigned phase,
const Context&
ectx)
970 Entry{PhaseEntry{&this->relativePermeability_,
971 [](
const unsigned phase,
const Context&
ectx)
972 {
return getValue(
ectx.intQuants.relativePermeability(phase)); }
975 Entry{PhaseEntry{&this->viscosity_,
978 if (this->extboC_.allocated() &&
phaseIdx == oilPhaseIdx) {
981 else if (this->extboC_.allocated() &&
phaseIdx == gasPhaseIdx) {
990 Entry{PhaseEntry{&this->residual_,
991 [&
modelResid = this->simulator_.model().linearizer().residual()]
994 const unsigned sIdx = FluidSystem::solventComponentIndex(
phaseIdx);
1001 Entry{ScalarEntry{&this->rockCompPorvMultiplier_,
1002 [&problem = this->simulator_.problem()](
const Context&
ectx)
1004 return problem.template
1010 Entry{ScalarEntry{&this->rockCompTransMultiplier_,
1011 [&problem = this->simulator_.problem()](
const Context&
ectx)
1018 Entry{ScalarEntry{&this->minimumOilPressure_,
1019 [&problem = this->simulator_.problem()](
const Context&
ectx)
1021 return std::min(
getValue(
ectx.fs.pressure(oilPhaseIdx)),
1022 problem.minOilPressure(
ectx.globalDofIdx));
1026 Entry{ScalarEntry{&this->bubblePointPressure_,
1028 &vanguard = this->simulator_.vanguard()](
const Context&
ectx)
1032 FluidSystem::bubblePointPressure(
ectx.fs,
1033 ectx.intQuants.pvtRegionIndex())
1043 Entry{ScalarEntry{&this->dewPointPressure_,
1045 &vanguard = this->simulator_.vanguard()](
const Context&
ectx)
1049 FluidSystem::dewPointPressure(
ectx.fs,
1050 ectx.intQuants.pvtRegionIndex())
1060 Entry{ScalarEntry{&this->overburdenPressure_,
1061 [&problem = simulator_.problem()](
const Context&
ectx)
1062 {
return problem.overburdenPressure(
ectx.globalDofIdx); }
1065 Entry{ScalarEntry{&this->temperature_,
1066 [](
const Context&
ectx)
1070 Entry{ScalarEntry{&this->sSol_,
1071 [](
const Context&
ectx)
1072 {
return getValue(
ectx.intQuants.solventSaturation()); }
1075 Entry{ScalarEntry{&this->rswSol_,
1076 [](
const Context&
ectx)
1080 Entry{ScalarEntry{&this->cPolymer_,
1081 [](
const Context&
ectx)
1082 {
return getValue(
ectx.intQuants.polymerConcentration()); }
1085 Entry{ScalarEntry{&this->cFoam_,
1086 [](
const Context&
ectx)
1087 {
return getValue(
ectx.intQuants.foamConcentration()); }
1090 Entry{ScalarEntry{&this->cSalt_,
1091 [](
const Context&
ectx)
1095 Entry{ScalarEntry{&this->pSalt_,
1096 [](
const Context&
ectx)
1100 Entry{ScalarEntry{&this->permFact_,
1101 [](
const Context&
ectx)
1105 Entry{ScalarEntry{&this->rPorV_,
1106 [&model = this->simulator_.model()](
const Context&
ectx)
1108 const auto totVolume = model.dofTotalVolume(
ectx.globalDofIdx);
1113 Entry{ScalarEntry{&this->rs_,
1114 [](
const Context&
ectx)
1118 Entry{ScalarEntry{&this->rv_,
1119 [](
const Context&
ectx)
1123 Entry{ScalarEntry{&this->rsw_,
1124 [](
const Context&
ectx)
1128 Entry{ScalarEntry{&this->rvw_,
1129 [](
const Context&
ectx)
1133 Entry{ScalarEntry{&this->ppcw_,
1134 [&
matLawManager = *this->simulator_.problem().materialLawManager()]
1135 (
const Context&
ectx)
1142 Entry{ScalarEntry{&this->drsdtcon_,
1143 [&problem = this->simulator_.problem()](
const Context&
ectx)
1145 return problem.drsdtcon(
ectx.globalDofIdx,
1150 Entry{ScalarEntry{&this->pcgw_,
1151 [](
const Context&
ectx)
1158 Entry{ScalarEntry{&this->pcow_,
1159 [](
const Context&
ectx)
1166 Entry{ScalarEntry{&this->pcog_,
1167 [](
const Context&
ectx)
1174 Entry{ScalarEntry{&this->fluidPressure_,
1175 [](
const Context&
ectx)
1177 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1181 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1192 Entry{ScalarEntry{&this->gasDissolutionFactor_,
1193 [&problem = this->simulator_.problem()](
const Context&
ectx)
1195 const Scalar
SoMax = problem.maxOilSaturation(
ectx.globalDofIdx);
1196 return FluidSystem::template
1204 Entry{ScalarEntry{&this->oilVaporizationFactor_,
1205 [&problem = this->simulator_.problem()](
const Context&
ectx)
1207 const Scalar
SoMax = problem.maxOilSaturation(
ectx.globalDofIdx);
1208 return FluidSystem::template
1216 Entry{ScalarEntry{&this->gasDissolutionFactorInWater_,
1217 [&problem = this->simulator_.problem()](
const Context&
ectx)
1219 const Scalar
SwMax = problem.maxWaterSaturation(
ectx.globalDofIdx);
1220 return FluidSystem::template
1228 Entry{ScalarEntry{&this->waterVaporizationFactor_,
1229 [](
const Context&
ectx)
1231 return FluidSystem::template
1238 Entry{ScalarEntry{&this->gasFormationVolumeFactor_,
1239 [](
const Context&
ectx)
1241 return 1.0 / FluidSystem::template
1248 Entry{ScalarEntry{&this->saturatedOilFormationVolumeFactor_,
1249 [](
const Context&
ectx)
1251 return 1.0 / FluidSystem::template
1258 Entry{ScalarEntry{&this->oilSaturationPressure_,
1259 [](
const Context&
ectx)
1261 return FluidSystem::template
1268 Entry{ScalarEntry{&this->soMax_,
1269 [&problem = this->simulator_.problem()](
const Context&
ectx)
1271 return std::max(
getValue(
ectx.fs.saturation(oilPhaseIdx)),
1272 problem.maxOilSaturation(
ectx.globalDofIdx));
1277 Entry{ScalarEntry{&this->swMax_,
1278 [&problem = this->simulator_.problem()](
const Context&
ectx)
1280 return std::max(
getValue(
ectx.fs.saturation(waterPhaseIdx)),
1281 problem.maxWaterSaturation(
ectx.globalDofIdx));
1286 Entry{ScalarEntry{&this->soMax_,
1287 [](
const Context&
ectx)
1288 {
return ectx.hParams.somax; }
1292 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1293 FluidSystem::phaseIsActive(waterPhaseIdx)
1295 Entry{ScalarEntry{&this->swMax_,
1296 [](
const Context&
ectx)
1297 {
return ectx.hParams.swmax; }
1301 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1302 FluidSystem::phaseIsActive(waterPhaseIdx)
1304 Entry{ScalarEntry{&this->swmin_,
1305 [](
const Context&
ectx)
1306 {
return ectx.hParams.swmin; }
1310 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1311 FluidSystem::phaseIsActive(waterPhaseIdx)
1313 Entry{ScalarEntry{&this->sgmax_,
1314 [](
const Context&
ectx)
1315 {
return ectx.hParams.sgmax; }
1319 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1320 FluidSystem::phaseIsActive(gasPhaseIdx)
1322 Entry{ScalarEntry{&this->shmax_,
1323 [](
const Context&
ectx)
1324 {
return ectx.hParams.shmax; }
1328 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1329 FluidSystem::phaseIsActive(gasPhaseIdx)
1331 Entry{ScalarEntry{&this->somin_,
1332 [](
const Context&
ectx)
1333 {
return ectx.hParams.somin; }
1337 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1338 FluidSystem::phaseIsActive(gasPhaseIdx)
1340 Entry{[&model = this->simulator_.model(),
this](
const Context&
ectx)
1344 const auto porv =
ectx.intQuants.referencePorosity()
1345 * model.dofTotalVolume(
ectx.globalDofIdx);
1347 this->aggregateAverageDensityContributions_(
ectx.fs,
ectx.globalDofIdx,
1348 static_cast<double>(porv));
1349 }, this->regionAvgDensity_.has_value()
1351 Entry{[&
extboC = this->extboC_](
const Context&
ectx)
1354 ectx.intQuants.xVolume().value(),
1355 ectx.intQuants.yVolume().value());
1357 ectx.intQuants.zFraction().value());
1366 (1.0 -
ectx.intQuants.yVolume().value()) +
1370 (1.0 -
ectx.intQuants.xVolume().value());
1373 ectx.intQuants.yVolume().value() +
1377 ectx.intQuants.xVolume().value();
1378 const Scalar
rhoO = FluidSystem::referenceDensity(gasPhaseIdx,
ectx.pvtRegionIdx);
1379 const Scalar
rhoG = FluidSystem::referenceDensity(gasPhaseIdx,
ectx.pvtRegionIdx);
1380 const Scalar
rhoCO2 =
ectx.intQuants.zRefDensity();
1382 extboC.assignMassFractions(
ectx.globalDofIdx,
1386 }, this->extboC_.allocated()
1388 Entry{[&
micpC = this->micpC_](
const Context&
ectx)
1391 ectx.intQuants.microbialConcentration().value(),
1392 ectx.intQuants.oxygenConcentration().value(),
1394 10 *
ectx.intQuants.ureaConcentration().value(),
1395 ectx.intQuants.biofilmConcentration().value(),
1396 ectx.intQuants.calciteConcentration().value());
1397 }, this->micpC_.allocated()
1399 Entry{[&
rftC = this->rftC_,
1400 &vanguard = this->simulator_.vanguard()](
const Context&
ectx)
1404 [&fs =
ectx.fs]() {
return getValue(fs.pressure(oilPhaseIdx)); },
1405 [&fs =
ectx.fs]() {
return getValue(fs.saturation(waterPhaseIdx)); },
1406 [&fs =
ectx.fs]() {
return getValue(fs.saturation(gasPhaseIdx)); });
1409 Entry{[&
tC = this->tracerC_,
1410 &
tM = this->simulator_.problem().tracerModel()](
const Context&
ectx)
1412 tC.assignFreeConcentrations(
ectx.globalDofIdx,
1415 tC.assignSolConcentrations(
ectx.globalDofIdx,
1420 Entry{[&
flowsInf = this->simulator_.problem().model().linearizer().getFlowsInfo(),
1421 &
flowsC = this->flowsC_](
const Context&
ectx)
1423 constexpr auto gas_idx = Indices::gasEnabled ?
1424 conti0EqIdx + Indices::canonicalToActiveComponentIndex(gasCompIdx) : -1;
1425 constexpr auto oil_idx = Indices::oilEnabled ?
1426 conti0EqIdx + Indices::canonicalToActiveComponentIndex(oilCompIdx) : -1;
1427 constexpr auto water_idx = Indices::waterEnabled ?
1428 conti0EqIdx + Indices::canonicalToActiveComponentIndex(waterCompIdx) : -1;
1438 }, !this->simulator_.problem().model().linearizer().getFlowsInfo().empty()
1440 Entry{[&
floresInf = this->simulator_.problem().model().linearizer().getFloresInfo(),
1441 &
flowsC = this->flowsC_](
const Context&
ectx)
1443 constexpr auto gas_idx = Indices::gasEnabled ?
1444 conti0EqIdx + Indices::canonicalToActiveComponentIndex(gasCompIdx) : -1;
1445 constexpr auto oil_idx = Indices::oilEnabled ?
1446 conti0EqIdx + Indices::canonicalToActiveComponentIndex(oilCompIdx) : -1;
1447 constexpr auto water_idx = Indices::waterEnabled ?
1448 conti0EqIdx + Indices::canonicalToActiveComponentIndex(waterCompIdx) : -1;
1458 }, !this->simulator_.problem().model().linearizer().getFloresInfo().empty()
1465 Entry{ScalarEntry{&this->rv_,
1466 [&problem = this->simulator_.problem()](
const Context&
ectx)
1467 {
return problem.initialFluidState(
ectx.globalDofIdx).Rv(); }
1469 simulator_.episodeIndex() < 0 &&
1470 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1471 FluidSystem::phaseIsActive(gasPhaseIdx)
1473 Entry{ScalarEntry{&this->rs_,
1474 [&problem = this->simulator_.problem()](
const Context&
ectx)
1475 {
return problem.initialFluidState(
ectx.globalDofIdx).Rs(); }
1477 simulator_.episodeIndex() < 0 &&
1478 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1479 FluidSystem::phaseIsActive(gasPhaseIdx)
1481 Entry{ScalarEntry{&this->rsw_,
1482 [&problem = this->simulator_.problem()](
const Context&
ectx)
1483 {
return problem.initialFluidState(
ectx.globalDofIdx).Rsw(); }
1485 simulator_.episodeIndex() < 0 &&
1486 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1487 FluidSystem::phaseIsActive(gasPhaseIdx)
1489 Entry{ScalarEntry{&this->rvw_,
1490 [&problem = this->simulator_.problem()](
const Context&
ectx)
1491 {
return problem.initialFluidState(
ectx.globalDofIdx).Rvw(); }
1493 simulator_.episodeIndex() < 0 &&
1494 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1495 FluidSystem::phaseIsActive(gasPhaseIdx)
1498 Entry{PhaseEntry{&this->density_,
1499 [&problem = this->simulator_.problem()](
const unsigned phase,
1500 const Context&
ectx)
1502 const auto&
fsInitial = problem.initialFluidState(
ectx.globalDofIdx);
1505 ectx.intQuants.pvtRegionIndex());
1508 simulator_.episodeIndex() < 0 &&
1509 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1510 FluidSystem::phaseIsActive(gasPhaseIdx)
1512 Entry{PhaseEntry{&this->invB_,
1513 [&problem = this->simulator_.problem()](
const unsigned phase,
1514 const Context&
ectx)
1516 const auto&
fsInitial = problem.initialFluidState(
ectx.globalDofIdx);
1517 return FluidSystem::inverseFormationVolumeFactor(
fsInitial,
1519 ectx.intQuants.pvtRegionIndex());
1522 simulator_.episodeIndex() < 0 &&
1523 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1524 FluidSystem::phaseIsActive(gasPhaseIdx)
1526 Entry{PhaseEntry{&this->viscosity_,
1527 [&problem = this->simulator_.problem()](
const unsigned phase,
1528 const Context&
ectx)
1530 const auto&
fsInitial = problem.initialFluidState(
ectx.globalDofIdx);
1531 return FluidSystem::viscosity(
fsInitial,
1533 ectx.intQuants.pvtRegionIndex());
1536 simulator_.episodeIndex() < 0 &&
1537 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1538 FluidSystem::phaseIsActive(gasPhaseIdx)
1546 if (this->mech_.allocated()) {
1547 this->extractors_.emplace(
1548 [&
mech = this->mech_,
1549 &model = simulator_.problem().geoMechModel()](
const Context&
ectx)
1551 mech.assignDelStress(
ectx.globalDofIdx,
1552 model.delstress(
ectx.globalDofIdx));
1554 mech.assignDisplacement(
ectx.globalDofIdx,
1555 model.disp(
ectx.globalDofIdx,
true));
1558 mech.assignFracStress(
ectx.globalDofIdx,
1559 model.fractureStress(
ectx.globalDofIdx));
1561 mech.assignLinStress(
ectx.globalDofIdx,
1562 model.linstress(
ectx.globalDofIdx));
1564 mech.assignPotentialForces(
ectx.globalDofIdx,
1565 model.mechPotentialForce(
ectx.globalDofIdx),
1566 model.mechPotentialPressForce(
ectx.globalDofIdx),
1567 model.mechPotentialTempForce(
ectx.globalDofIdx));
1569 mech.assignStrain(
ectx.globalDofIdx,
1570 model.strain(
ectx.globalDofIdx,
true));
1573 mech.assignStress(
ectx.globalDofIdx,
1574 model.stress(
ectx.globalDofIdx,
true));
1582 void setupBlockExtractors_(
const bool isSubStep,
1583 const int reportStepNum)
1586 using Context =
typename BlockExtractor::Context;
1587 using PhaseEntry =
typename BlockExtractor::PhaseEntry;
1588 using ScalarEntry =
typename BlockExtractor::ScalarEntry;
1590 using namespace std::string_view_literals;
1593 Entry{ScalarEntry{std::vector{
"BPR"sv,
"BPRESSUR"sv},
1594 [](
const Context&
ectx)
1596 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1599 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1611 Entry{PhaseEntry{std::array{
1612 std::array{
"BWSAT"sv,
"BOSAT"sv,
"BGSAT"sv},
1613 std::array{
"BSWAT"sv,
"BSOIL"sv,
"BSGAS"sv}
1621 Entry{ScalarEntry{
"BNSAT",
1622 [](
const Context&
ectx)
1624 return ectx.intQuants.solventSaturation().value();
1628 Entry{ScalarEntry{std::vector{
"BTCNFHEA"sv,
"BTEMP"sv},
1629 [](
const Context&
ectx)
1631 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1634 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1643 Entry{PhaseEntry{std::array{
1644 std::array{
"BWKR"sv,
"BOKR"sv,
"BGKR"sv},
1645 std::array{
"BKRW"sv,
"BKRO"sv,
"BKRG"sv}
1653 Entry{ScalarEntry{
"BKROG",
1654 [&problem = this->simulator_.problem()](
const Context&
ectx)
1657 problem.materialLawParams(
ectx.elemCtx,
1660 return getValue(MaterialLaw::template
1666 Entry{ScalarEntry{
"BKROW",
1667 [&problem = this->simulator_.problem()](
const Context&
ectx)
1672 return getValue(MaterialLaw::template
1678 Entry{ScalarEntry{
"BWPC",
1679 [](
const Context&
ectx)
1686 Entry{ScalarEntry{
"BGPC",
1687 [](
const Context&
ectx)
1694 Entry{ScalarEntry{
"BWPR",
1695 [](
const Context&
ectx)
1701 Entry{ScalarEntry{
"BGPR",
1702 [](
const Context&
ectx)
1708 Entry{PhaseEntry{std::array{
1709 std::array{
"BVWAT"sv,
"BVOIL"sv,
"BVGAS"sv},
1710 std::array{
"BWVIS"sv,
"BOVIS"sv,
"BGVIS"sv}
1718 Entry{PhaseEntry{std::array{
1719 std::array{
"BWDEN"sv,
"BODEN"sv,
"BGDEN"sv},
1720 std::array{
"BDENW"sv,
"BDENO"sv,
"BDENG"sv}
1728 Entry{ScalarEntry{
"BFLOWI",
1729 [&
flowsC = this->flowsC_](
const Context&
ectx)
1731 return flowsC.getFlow(
ectx.globalDofIdx, Dir::XPlus, waterCompIdx);
1735 Entry{ScalarEntry{
"BFLOWJ",
1736 [&
flowsC = this->flowsC_](
const Context&
ectx)
1738 return flowsC.getFlow(
ectx.globalDofIdx, Dir::YPlus, waterCompIdx);
1742 Entry{ScalarEntry{
"BFLOWK",
1743 [&
flowsC = this->flowsC_](
const Context&
ectx)
1745 return flowsC.getFlow(
ectx.globalDofIdx, Dir::ZPlus, waterCompIdx);
1749 Entry{ScalarEntry{
"BRPV",
1750 [&model = this->simulator_.model()](
const Context&
ectx)
1753 model.dofTotalVolume(
ectx.globalDofIdx);
1757 Entry{PhaseEntry{std::array{
"BWPV"sv,
"BOPV"sv,
"BGPV"sv},
1758 [&model = this->simulator_.model()](
const unsigned phaseIdx,
1759 const Context&
ectx)
1763 model.dofTotalVolume(
ectx.globalDofIdx);
1767 Entry{ScalarEntry{
"BRS",
1768 [](
const Context&
ectx)
1774 Entry{ScalarEntry{
"BRV",
1775 [](
const Context&
ectx)
1781 Entry{ScalarEntry{
"BOIP",
1782 [&model = this->simulator_.model()](
const Context&
ectx)
1789 model.dofTotalVolume(
ectx.globalDofIdx) *
1794 Entry{ScalarEntry{
"BGIP",
1795 [&model = this->simulator_.model()](
const Context&
ectx)
1800 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1812 model.dofTotalVolume(
ectx.globalDofIdx) *
1817 Entry{ScalarEntry{
"BWIP",
1818 [&model = this->simulator_.model()](
const Context&
ectx)
1822 model.dofTotalVolume(
ectx.globalDofIdx) *
1827 Entry{ScalarEntry{
"BOIPL",
1828 [&model = this->simulator_.model()](
const Context&
ectx)
1832 model.dofTotalVolume(
ectx.globalDofIdx) *
1837 Entry{ScalarEntry{
"BGIPL",
1838 [&model = this->simulator_.model()](
const Context&
ectx)
1841 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1852 model.dofTotalVolume(
ectx.globalDofIdx) *
1857 Entry{ScalarEntry{
"BGIPG",
1858 [&model = this->simulator_.model()](
const Context&
ectx)
1862 model.dofTotalVolume(
ectx.globalDofIdx) *
1867 Entry{ScalarEntry{
"BOIPG",
1868 [&model = this->simulator_.model()](
const Context&
ectx)
1873 model.dofTotalVolume(
ectx.globalDofIdx) *
1878 Entry{PhaseEntry{std::array{
"BPPW"sv,
"BPPO"sv,
"BPPG"sv},
1879 [&
simConfig = this->eclState_.getSimulationConfig(),
1880 &
grav = this->simulator_.problem().gravity(),
1882 &problem = this->simulator_.problem(),
1885 auto phase = RegionPhasePoreVolAverage::Phase{};
1898 const auto region = RegionPhasePoreVolAverage::Region {
1899 ectx.elemCtx.primaryVars(
ectx.dofIdx, 0).pvtRegionIndex() + 1
1904 const auto press =
getValue(
ectx.fs.pressure(phase.ix));
1905 const auto dz = problem.dofCenterDepth(
ectx.globalDofIdx) - datum;
1906 return press - density*
dz*
grav[GridView::dimensionworld - 1];
1914 this->extraBlockData_.clear();
1917 const auto&
rpt = this->schedule_[reportStepNum - 1].rpt_config.get();
1918 if (
rpt.contains(
"WELLS") &&
rpt.at(
"WELLS") > 1) {
1919 this->setupExtraBlockData(reportStepNum,
1920 [&
c = this->collectOnIORank_](
const int idx)
1921 {
return c.isCartIdxOnThisRank(
idx); });
1932 const Simulator& simulator_;
1933 const CollectDataOnIORankType& collectOnIORank_;
1934 std::vector<typename Extractor::Entry> extractors_;