101 using Element =
typename GridView::template Codim<0>::Entity;
102 using ElementIterator =
typename GridView::template Codim<0>::Iterator;
105 using Dir = FaceDir::DirEnum;
107 enum { conti0EqIdx = Indices::conti0EqIdx };
108 enum { numPhases = FluidSystem::numPhases };
109 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
110 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
111 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
112 enum { gasCompIdx = FluidSystem::gasCompIdx };
113 enum { oilCompIdx = FluidSystem::oilCompIdx };
114 enum { waterCompIdx = FluidSystem::waterCompIdx };
118 template <
class CollectDataToIORankType>
120 const SummaryConfig& smryCfg,
121 const CollectDataToIORankType& collectToIORank)
122 :
BaseType(simulator.vanguard().eclState(),
123 simulator.vanguard().schedule(),
125 simulator.vanguard().summaryState(),
137 , simulator_(simulator)
139 for (
auto& region_pair : this->regions_) {
140 this->createLocalRegion_(region_pair.second);
143 auto isCartIdxOnThisRank = [&collectToIORank](
const int idx) {
144 return collectToIORank.isCartIdxOnThisRank(idx);
147 this->setupBlockData(isCartIdxOnThisRank);
149 this->forceDisableFipOutput_ =
150 Parameters::Get<Parameters::ForceDisableFluidInPlaceOutput>();
152 this->forceDisableFipresvOutput_ =
153 Parameters::Get<Parameters::ForceDisableResvFluidInPlaceOutput>();
155 if (! Parameters::Get<Parameters::OwnerCellsFirst>()) {
156 const std::string msg =
"The output code does not support --owner-cells-first=false.";
157 if (collectToIORank.isIORank()) {
160 OPM_THROW_NOLOG(std::runtime_error, msg);
163 if (smryCfg.match(
"[FB]PP[OGW]") || smryCfg.match(
"RPP[OGW]*")) {
164 auto rset = this->eclState_.fieldProps().fip_regions();
165 rset.push_back(
"PVTNUM");
170 this->regionAvgDensity_
171 .emplace(this->simulator_.gridView().comm(),
172 FluidSystem::numPhases, rset,
173 [fp = std::cref(this->eclState_.fieldProps())]
174 (
const std::string& rsetName) ->
decltype(
auto)
175 { return fp.get().get_int(rsetName); });
184 Parameters::Register<Parameters::ForceDisableFluidInPlaceOutput>
185 (
"Do not print fluid-in-place values after each report step "
186 "even if requested by the deck.");
187 Parameters::Register<Parameters::ForceDisableResvFluidInPlaceOutput>
188 (
"Do not print reservoir volumes values after each report step "
189 "even if requested by the deck.");
198 const unsigned reportStepNum,
201 const bool isRestart)
207 const auto& problem = this->simulator_.problem();
209 this->doAllocBuffers(bufferSize,
214 problem.vapparsActive(std::max(simulator_.episodeIndex(), 0)),
215 problem.materialLawManager()->enablePCHysteresis(),
216 problem.materialLawManager()->enableNonWettingHysteresis(),
217 problem.materialLawManager()->enableWettingHysteresis(),
218 problem.tracerModel().numTracers(),
219 problem.tracerModel().enableSolTracers(),
220 problem.eclWriter()->getOutputNnc().size());
223 void processElementMech(
const ElementContext& elemCtx)
226 const auto& problem = elemCtx.simulator().problem();
227 const auto& model = problem.geoMechModel();
228 for (
unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(0); ++dofIdx) {
229 unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, 0);
230 if (!this->mechPotentialForce_.empty()) {
232 this->mechPotentialForce_[globalDofIdx] = model.mechPotentialForce(globalDofIdx);
233 this->mechPotentialPressForce_[globalDofIdx] = model.mechPotentialPressForce(globalDofIdx);
234 this->mechPotentialTempForce_[globalDofIdx] = model.mechPotentialTempForce(globalDofIdx);
236 this->dispX_[globalDofIdx] = model.disp(globalDofIdx, 0);
237 this->dispY_[globalDofIdx] = model.disp(globalDofIdx, 1);
238 this->dispZ_[globalDofIdx] = model.disp(globalDofIdx, 2);
239 this->stressXX_[globalDofIdx] = model.stress(globalDofIdx, 0);
240 this->stressYY_[globalDofIdx] = model.stress(globalDofIdx, 1);
241 this->stressZZ_[globalDofIdx] = model.stress(globalDofIdx, 2);
243 this->stressXY_[globalDofIdx] = model.stress(globalDofIdx, 5);
244 this->stressXZ_[globalDofIdx] = model.stress(globalDofIdx, 4);
245 this->stressYZ_[globalDofIdx] = model.stress(globalDofIdx, 3);
247 this->strainXX_[globalDofIdx] = model.strain(globalDofIdx, 0);
248 this->strainYY_[globalDofIdx] = model.strain(globalDofIdx, 1);
249 this->strainZZ_[globalDofIdx] = model.strain(globalDofIdx, 2);
251 this->strainXY_[globalDofIdx] = model.strain(globalDofIdx, 5);
252 this->strainXZ_[globalDofIdx] = model.strain(globalDofIdx, 4);
253 this->strainYZ_[globalDofIdx] = model.strain(globalDofIdx, 3);
256 this->delstressXX_[globalDofIdx] = model.delstress(globalDofIdx, 0);
257 this->delstressYY_[globalDofIdx] = model.delstress(globalDofIdx, 1);
258 this->delstressZZ_[globalDofIdx] = model.delstress(globalDofIdx, 2);
260 this->delstressXY_[globalDofIdx] = model.delstress(globalDofIdx, 5);
261 this->delstressXZ_[globalDofIdx] = model.delstress(globalDofIdx, 4);
262 this->delstressYZ_[globalDofIdx] = model.delstress(globalDofIdx, 3);
278 const auto& problem = elemCtx.simulator().problem();
279 const auto& modelResid = elemCtx.simulator().model().linearizer().residual();
280 for (
unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(0); ++dofIdx) {
281 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, 0);
282 const auto& fs = intQuants.fluidState();
284 using FluidState = std::remove_cv_t<std::remove_reference_t<
decltype(fs)>>;
286 const unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, 0);
287 const unsigned pvtRegionIdx = elemCtx.primaryVars(dofIdx, 0).pvtRegionIndex();
289 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
290 if (this->saturation_[phaseIdx].empty())
293 this->saturation_[phaseIdx][globalDofIdx] = getValue(fs.saturation(phaseIdx));
294 Valgrind::CheckDefined(this->saturation_[phaseIdx][globalDofIdx]);
297 if (this->regionAvgDensity_.has_value()) {
300 const auto porv = intQuants.referencePorosity()
301 * elemCtx.simulator().model().dofTotalVolume(globalDofIdx);
303 this->aggregateAverageDensityContributions_(fs, globalDofIdx,
304 static_cast<double>(porv));
307 if (!this->fluidPressure_.empty()) {
308 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
310 this->fluidPressure_[globalDofIdx] = getValue(fs.pressure(oilPhaseIdx));
311 }
else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
313 this->fluidPressure_[globalDofIdx] = getValue(fs.pressure(gasPhaseIdx));
316 this->fluidPressure_[globalDofIdx] = getValue(fs.pressure(waterPhaseIdx));
318 Valgrind::CheckDefined(this->fluidPressure_[globalDofIdx]);
321 if (!this->temperature_.empty()) {
322 this->temperature_[globalDofIdx] = getValue(fs.temperature(oilPhaseIdx));
323 Valgrind::CheckDefined(this->temperature_[globalDofIdx]);
325 if (!this->gasDissolutionFactor_.empty()) {
326 Scalar SoMax = elemCtx.problem().maxOilSaturation(globalDofIdx);
327 this->gasDissolutionFactor_[globalDofIdx]
328 = FluidSystem::template saturatedDissolutionFactor<FluidState, Scalar>(
329 fs, oilPhaseIdx, pvtRegionIdx, SoMax);
330 Valgrind::CheckDefined(this->gasDissolutionFactor_[globalDofIdx]);
332 if (!this->oilVaporizationFactor_.empty()) {
333 Scalar SoMax = elemCtx.problem().maxOilSaturation(globalDofIdx);
334 this->oilVaporizationFactor_[globalDofIdx]
335 = FluidSystem::template saturatedDissolutionFactor<FluidState, Scalar>(
336 fs, gasPhaseIdx, pvtRegionIdx, SoMax);
337 Valgrind::CheckDefined(this->oilVaporizationFactor_[globalDofIdx]);
339 if (!this->gasDissolutionFactorInWater_.empty()) {
340 Scalar SwMax = elemCtx.problem().maxWaterSaturation(globalDofIdx);
341 this->gasDissolutionFactorInWater_[globalDofIdx]
342 = FluidSystem::template saturatedDissolutionFactor<FluidState, Scalar>(
343 fs, waterPhaseIdx, pvtRegionIdx, SwMax);
344 Valgrind::CheckDefined(this->gasDissolutionFactorInWater_[globalDofIdx]);
346 if (!this->waterVaporizationFactor_.empty()) {
347 this->waterVaporizationFactor_[globalDofIdx]
348 = FluidSystem::template saturatedVaporizationFactor<FluidState, Scalar>(
349 fs, gasPhaseIdx, pvtRegionIdx);
350 Valgrind::CheckDefined(this->waterVaporizationFactor_[globalDofIdx]);
352 if (!this->gasFormationVolumeFactor_.empty()) {
353 this->gasFormationVolumeFactor_[globalDofIdx] = 1.0
354 / FluidSystem::template inverseFormationVolumeFactor<FluidState, Scalar>(
355 fs, gasPhaseIdx, pvtRegionIdx);
356 Valgrind::CheckDefined(this->gasFormationVolumeFactor_[globalDofIdx]);
358 if (!this->saturatedOilFormationVolumeFactor_.empty()) {
359 this->saturatedOilFormationVolumeFactor_[globalDofIdx] = 1.0
360 / FluidSystem::template saturatedInverseFormationVolumeFactor<FluidState, Scalar>(
361 fs, oilPhaseIdx, pvtRegionIdx);
362 Valgrind::CheckDefined(this->saturatedOilFormationVolumeFactor_[globalDofIdx]);
364 if (!this->oilSaturationPressure_.empty()) {
365 this->oilSaturationPressure_[globalDofIdx]
366 = FluidSystem::template saturationPressure<FluidState, Scalar>(fs, oilPhaseIdx, pvtRegionIdx);
367 Valgrind::CheckDefined(this->oilSaturationPressure_[globalDofIdx]);
370 if (!this->rs_.empty()) {
371 this->rs_[globalDofIdx] = getValue(fs.Rs());
372 Valgrind::CheckDefined(this->rs_[globalDofIdx]);
374 if (!this->rsw_.empty()) {
375 this->rsw_[globalDofIdx] = getValue(fs.Rsw());
376 Valgrind::CheckDefined(this->rsw_[globalDofIdx]);
379 if (!this->rv_.empty()) {
380 this->rv_[globalDofIdx] = getValue(fs.Rv());
381 Valgrind::CheckDefined(this->rv_[globalDofIdx]);
383 if (!this->pcgw_.empty()) {
384 this->pcgw_[globalDofIdx] = getValue(fs.pressure(gasPhaseIdx)) - getValue(fs.pressure(waterPhaseIdx));
385 Valgrind::CheckDefined(this->pcgw_[globalDofIdx]);
387 if (!this->pcow_.empty()) {
388 this->pcow_[globalDofIdx] = getValue(fs.pressure(oilPhaseIdx)) - getValue(fs.pressure(waterPhaseIdx));
389 Valgrind::CheckDefined(this->pcow_[globalDofIdx]);
391 if (!this->pcog_.empty()) {
392 this->pcog_[globalDofIdx] = getValue(fs.pressure(gasPhaseIdx)) - getValue(fs.pressure(oilPhaseIdx));
393 Valgrind::CheckDefined(this->pcog_[globalDofIdx]);
396 if (!this->rvw_.empty()) {
397 this->rvw_[globalDofIdx] = getValue(fs.Rvw());
398 Valgrind::CheckDefined(this->rvw_[globalDofIdx]);
401 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
402 if (this->invB_[phaseIdx].empty())
405 this->invB_[phaseIdx][globalDofIdx] = getValue(fs.invB(phaseIdx));
406 Valgrind::CheckDefined(this->invB_[phaseIdx][globalDofIdx]);
409 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
410 if (this->density_[phaseIdx].empty())
413 this->density_[phaseIdx][globalDofIdx] = getValue(fs.density(phaseIdx));
414 Valgrind::CheckDefined(this->density_[phaseIdx][globalDofIdx]);
417 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
418 if (this->viscosity_[phaseIdx].empty())
421 if (!this->extboX_.empty() && phaseIdx == oilPhaseIdx)
422 this->viscosity_[phaseIdx][globalDofIdx] = getValue(intQuants.oilViscosity());
423 else if (!this->extboX_.empty() && phaseIdx == gasPhaseIdx)
424 this->viscosity_[phaseIdx][globalDofIdx] = getValue(intQuants.gasViscosity());
426 this->viscosity_[phaseIdx][globalDofIdx] = getValue(fs.viscosity(phaseIdx));
427 Valgrind::CheckDefined(this->viscosity_[phaseIdx][globalDofIdx]);
430 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
431 if (this->relativePermeability_[phaseIdx].empty())
434 this->relativePermeability_[phaseIdx][globalDofIdx]
435 = getValue(intQuants.relativePermeability(phaseIdx));
436 Valgrind::CheckDefined(this->relativePermeability_[phaseIdx][globalDofIdx]);
439 if (!this->drsdtcon_.empty()) {
440 this->drsdtcon_[globalDofIdx] = problem.drsdtcon(globalDofIdx, elemCtx.simulator().episodeIndex());
443 if (!this->sSol_.empty()) {
444 this->sSol_[globalDofIdx] = intQuants.solventSaturation().value();
447 if (!this->rswSol_.empty()) {
448 this->rswSol_[globalDofIdx] = intQuants.rsSolw().value();
451 if (!this->cPolymer_.empty()) {
452 this->cPolymer_[globalDofIdx] = intQuants.polymerConcentration().value();
455 if (!this->cFoam_.empty()) {
456 this->cFoam_[globalDofIdx] = intQuants.foamConcentration().value();
459 if (!this->cSalt_.empty()) {
460 this->cSalt_[globalDofIdx] = fs.saltConcentration().value();
463 if (!this->pSalt_.empty()) {
464 this->pSalt_[globalDofIdx] = intQuants.saltSaturation().value();
467 if (!this->permFact_.empty()) {
468 this->permFact_[globalDofIdx] = intQuants.permFactor().value();
471 if (!this->extboX_.empty()) {
472 this->extboX_[globalDofIdx] = intQuants.xVolume().value();
475 if (!this->extboY_.empty()) {
476 this->extboY_[globalDofIdx] = intQuants.yVolume().value();
479 if (!this->extboZ_.empty()) {
480 this->extboZ_[globalDofIdx] = intQuants.zFraction().value();
483 if (!this->rPorV_.empty()) {
484 const auto totVolume = elemCtx.simulator().model().dofTotalVolume(globalDofIdx);
485 this->rPorV_[globalDofIdx] = totVolume * intQuants.porosity().value();
488 if (!this->mFracCo2_.empty()) {
489 const Scalar stdVolOil = getValue(fs.saturation(oilPhaseIdx)) * getValue(fs.invB(oilPhaseIdx))
490 + getValue(fs.saturation(gasPhaseIdx)) * getValue(fs.invB(gasPhaseIdx)) * getValue(fs.Rv());
491 const Scalar stdVolGas = getValue(fs.saturation(gasPhaseIdx)) * getValue(fs.invB(gasPhaseIdx))
492 * (1.0 - intQuants.yVolume().value())
493 + getValue(fs.saturation(oilPhaseIdx)) * getValue(fs.invB(oilPhaseIdx)) * getValue(fs.Rs())
494 * (1.0 - intQuants.xVolume().value());
495 const Scalar stdVolCo2 = getValue(fs.saturation(gasPhaseIdx)) * getValue(fs.invB(gasPhaseIdx))
496 * intQuants.yVolume().value()
497 + getValue(fs.saturation(oilPhaseIdx)) * getValue(fs.invB(oilPhaseIdx)) * getValue(fs.Rs())
498 * intQuants.xVolume().value();
499 const Scalar rhoO = FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
500 const Scalar rhoG = FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
501 const Scalar rhoCO2 = intQuants.zRefDensity();
502 const Scalar stdMassTotal = 1.0e-10 + stdVolOil * rhoO + stdVolGas * rhoG + stdVolCo2 * rhoCO2;
503 this->mFracOil_[globalDofIdx] = stdVolOil * rhoO / stdMassTotal;
504 this->mFracGas_[globalDofIdx] = stdVolGas * rhoG / stdMassTotal;
505 this->mFracCo2_[globalDofIdx] = stdVolCo2 * rhoCO2 / stdMassTotal;
508 if (!this->cMicrobes_.empty()) {
509 this->cMicrobes_[globalDofIdx] = intQuants.microbialConcentration().value();
512 if (!this->cOxygen_.empty()) {
513 this->cOxygen_[globalDofIdx] = intQuants.oxygenConcentration().value();
516 if (!this->cUrea_.empty()) {
517 this->cUrea_[globalDofIdx] = 10
518 * intQuants.ureaConcentration()
522 if (!this->cBiofilm_.empty()) {
523 this->cBiofilm_[globalDofIdx] = intQuants.biofilmConcentration().value();
526 if (!this->cCalcite_.empty()) {
527 this->cCalcite_[globalDofIdx] = intQuants.calciteConcentration().value();
530 if (!this->bubblePointPressure_.empty()) {
532 this->bubblePointPressure_[globalDofIdx]
533 = getValue(FluidSystem::bubblePointPressure(fs, intQuants.pvtRegionIndex()));
534 }
catch (
const NumericalProblem&) {
535 const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx);
536 this->failedCellsPb_.push_back(cartesianIdx);
540 if (!this->dewPointPressure_.empty()) {
542 this->dewPointPressure_[globalDofIdx]
543 = getValue(FluidSystem::dewPointPressure(fs, intQuants.pvtRegionIndex()));
544 }
catch (
const NumericalProblem&) {
545 const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx);
546 this->failedCellsPd_.push_back(cartesianIdx);
550 if (!this->minimumOilPressure_.empty())
551 this->minimumOilPressure_[globalDofIdx]
552 = std::min(getValue(fs.pressure(oilPhaseIdx)), problem.minOilPressure(globalDofIdx));
554 if (!this->overburdenPressure_.empty())
555 this->overburdenPressure_[globalDofIdx] = problem.overburdenPressure(globalDofIdx);
557 if (!this->rockCompPorvMultiplier_.empty())
558 this->rockCompPorvMultiplier_[globalDofIdx]
559 = problem.template rockCompPoroMultiplier<Scalar>(intQuants, globalDofIdx);
561 if (!this->rockCompTransMultiplier_.empty())
562 this->rockCompTransMultiplier_[globalDofIdx]
563 = problem.template rockCompTransMultiplier<Scalar>(intQuants, globalDofIdx);
565 const auto& matLawManager = problem.materialLawManager();
566 if (matLawManager->enableHysteresis()) {
567 if (FluidSystem::phaseIsActive(oilPhaseIdx)
568 && FluidSystem::phaseIsActive(waterPhaseIdx)) {
573 matLawManager->oilWaterHysteresisParams(
574 somax, swmax, swmin, globalDofIdx);
576 if (matLawManager->enableNonWettingHysteresis()) {
577 if (!this->soMax_.empty()) {
578 this->soMax_[globalDofIdx] = somax;
581 if (matLawManager->enableWettingHysteresis()) {
582 if (!this->swMax_.empty()) {
583 this->swMax_[globalDofIdx] = swmax;
586 if (matLawManager->enablePCHysteresis()) {
587 if (!this->swmin_.empty()) {
588 this->swmin_[globalDofIdx] = swmin;
593 if (FluidSystem::phaseIsActive(oilPhaseIdx)
594 && FluidSystem::phaseIsActive(gasPhaseIdx)) {
598 matLawManager->gasOilHysteresisParams(
599 sgmax, shmax, somin, globalDofIdx);
601 if (matLawManager->enableNonWettingHysteresis()) {
602 if (!this->sgmax_.empty()) {
603 this->sgmax_[globalDofIdx] = sgmax;
606 if (matLawManager->enableWettingHysteresis()) {
607 if (!this->shmax_.empty()) {
608 this->shmax_[globalDofIdx] = shmax;
611 if (matLawManager->enablePCHysteresis()) {
612 if (!this->somin_.empty()) {
613 this->somin_[globalDofIdx] = somin;
619 if (!this->soMax_.empty())
620 this->soMax_[globalDofIdx]
621 = std::max(getValue(fs.saturation(oilPhaseIdx)), problem.maxOilSaturation(globalDofIdx));
623 if (!this->swMax_.empty())
624 this->swMax_[globalDofIdx]
625 = std::max(getValue(fs.saturation(waterPhaseIdx)), problem.maxWaterSaturation(globalDofIdx));
628 if (!this->ppcw_.empty()) {
629 this->ppcw_[globalDofIdx] = matLawManager->oilWaterScaledEpsInfoDrainage(globalDofIdx).maxPcow;
638 if ((elemCtx.simulator().episodeIndex() < 0) &&
639 FluidSystem::phaseIsActive(oilPhaseIdx) &&
640 FluidSystem::phaseIsActive(gasPhaseIdx))
642 const auto& fsInitial = problem.initialFluidState(globalDofIdx);
645 if (!this->rv_.empty())
646 this->rv_[globalDofIdx] = fsInitial.Rv();
648 if (!this->rs_.empty())
649 this->rs_[globalDofIdx] = fsInitial.Rs();
651 if (!this->rsw_.empty())
652 this->rsw_[globalDofIdx] = fsInitial.Rsw();
654 if (!this->rvw_.empty())
655 this->rvw_[globalDofIdx] = fsInitial.Rvw();
658 if (!this->density_[oilPhaseIdx].empty())
659 this->density_[oilPhaseIdx][globalDofIdx]
660 = FluidSystem::density(fsInitial, oilPhaseIdx, intQuants.pvtRegionIndex());
662 if (!this->density_[gasPhaseIdx].empty())
663 this->density_[gasPhaseIdx][globalDofIdx]
664 = FluidSystem::density(fsInitial, gasPhaseIdx, intQuants.pvtRegionIndex());
666 if (!this->invB_[oilPhaseIdx].empty())
667 this->invB_[oilPhaseIdx][globalDofIdx]
668 = FluidSystem::inverseFormationVolumeFactor(fsInitial, oilPhaseIdx, intQuants.pvtRegionIndex());
670 if (!this->invB_[gasPhaseIdx].empty())
671 this->invB_[gasPhaseIdx][globalDofIdx]
672 = FluidSystem::inverseFormationVolumeFactor(fsInitial, gasPhaseIdx, intQuants.pvtRegionIndex());
674 if (!this->viscosity_[oilPhaseIdx].empty())
675 this->viscosity_[oilPhaseIdx][globalDofIdx]
676 = FluidSystem::viscosity(fsInitial, oilPhaseIdx, intQuants.pvtRegionIndex());
678 if (!this->viscosity_[gasPhaseIdx].empty())
679 this->viscosity_[gasPhaseIdx][globalDofIdx]
680 = FluidSystem::viscosity(fsInitial, gasPhaseIdx, intQuants.pvtRegionIndex());
684 const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx);
685 if (this->oilConnectionPressures_.count(cartesianIdx) > 0) {
686 this->oilConnectionPressures_[cartesianIdx] = getValue(fs.pressure(oilPhaseIdx));
688 if (this->waterConnectionSaturations_.count(cartesianIdx) > 0) {
689 this->waterConnectionSaturations_[cartesianIdx] = getValue(fs.saturation(waterPhaseIdx));
691 if (this->gasConnectionSaturations_.count(cartesianIdx) > 0) {
692 this->gasConnectionSaturations_[cartesianIdx] = getValue(fs.saturation(gasPhaseIdx));
696 const auto& tracerModel = simulator_.problem().tracerModel();
697 if (! this->freeTracerConcentrations_.empty()) {
698 for (
int tracerIdx = 0; tracerIdx < tracerModel.numTracers(); ++tracerIdx) {
699 if (this->freeTracerConcentrations_[tracerIdx].empty()) {
702 this->freeTracerConcentrations_[tracerIdx][globalDofIdx] =
703 tracerModel.freeTracerConcentration(tracerIdx, globalDofIdx);
706 if (! this->solTracerConcentrations_.empty()) {
707 for (
int tracerIdx = 0; tracerIdx < tracerModel.numTracers(); ++tracerIdx) {
708 if (this->solTracerConcentrations_[tracerIdx].empty()) {
711 this->solTracerConcentrations_[tracerIdx][globalDofIdx] =
712 tracerModel.solTracerConcentration(tracerIdx, globalDofIdx);
718 for (
int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx )
720 if (!this->residual_[phaseIdx].empty()) {
721 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
722 this->residual_[phaseIdx][globalDofIdx] = modelResid[globalDofIdx][activeCompIdx];
728 void processElementFlows(
const ElementContext& elemCtx)
730 OPM_TIMEBLOCK_LOCAL(processElementBlockData);
734 const auto& problem = elemCtx.simulator().problem();
735 for (
unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(0); ++dofIdx) {
737 unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, 0);
738 if (!problem.model().linearizer().getFlowsInfo().empty()) {
739 const auto& flowsInf = problem.model().linearizer().getFlowsInfo();
740 auto flowsInfos = flowsInf[globalDofIdx];
741 for (
auto& flowsInfo : flowsInfos) {
742 if (flowsInfo.faceId >= 0) {
743 if (!this->flows_[flowsInfo.faceId][gasCompIdx].empty()) {
744 this->flows_[flowsInfo.faceId][gasCompIdx][globalDofIdx]
745 = flowsInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(gasCompIdx)];
747 if (!this->flows_[flowsInfo.faceId][oilCompIdx].empty()) {
748 this->flows_[flowsInfo.faceId][oilCompIdx][globalDofIdx]
749 = flowsInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(oilCompIdx)];
751 if (!this->flows_[flowsInfo.faceId][waterCompIdx].empty()) {
752 this->flows_[flowsInfo.faceId][waterCompIdx][globalDofIdx]
753 = flowsInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(waterCompIdx)];
756 if (flowsInfo.faceId == -2) {
757 if (!this->flowsn_[gasCompIdx].indices.empty()) {
758 this->flowsn_[gasCompIdx].indices[flowsInfo.nncId] = flowsInfo.nncId;
759 this->flowsn_[gasCompIdx].values[flowsInfo.nncId]
760 = flowsInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(gasCompIdx)];
762 if (!this->flowsn_[oilCompIdx].indices.empty()) {
763 this->flowsn_[oilCompIdx].indices[flowsInfo.nncId] = flowsInfo.nncId;
764 this->flowsn_[oilCompIdx].values[flowsInfo.nncId]
765 = flowsInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(oilCompIdx)];
767 if (!this->flowsn_[waterCompIdx].indices.empty()) {
768 this->flowsn_[waterCompIdx].indices[flowsInfo.nncId] = flowsInfo.nncId;
769 this->flowsn_[waterCompIdx].values[flowsInfo.nncId]
770 = flowsInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(waterCompIdx)];
777 if (!problem.model().linearizer().getFloresInfo().empty()) {
778 const auto& floresInf = problem.model().linearizer().getFloresInfo();
779 auto floresInfos =floresInf[globalDofIdx];
780 for (
auto& floresInfo : floresInfos) {
781 if (floresInfo.faceId >= 0) {
782 if (!this->flores_[floresInfo.faceId][gasCompIdx].empty()) {
783 this->flores_[floresInfo.faceId][gasCompIdx][globalDofIdx]
784 = floresInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(gasCompIdx)];
786 if (!this->flores_[floresInfo.faceId][oilCompIdx].empty()) {
787 this->flores_[floresInfo.faceId][oilCompIdx][globalDofIdx]
788 = floresInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(oilCompIdx)];
790 if (!this->flores_[floresInfo.faceId][waterCompIdx].empty()) {
791 this->flores_[floresInfo.faceId][waterCompIdx][globalDofIdx]
792 = floresInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(waterCompIdx)];
796 if (floresInfo.faceId == -2) {
797 if (!this->floresn_[gasCompIdx].indices.empty()) {
798 this->floresn_[gasCompIdx].indices[floresInfo.nncId] = floresInfo.nncId;
799 this->floresn_[gasCompIdx].values[floresInfo.nncId]
800 = floresInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(gasCompIdx)];
802 if (!this->floresn_[oilCompIdx].indices.empty()) {
803 this->floresn_[oilCompIdx].indices[floresInfo.nncId] = floresInfo.nncId;
804 this->floresn_[oilCompIdx].values[floresInfo.nncId]
805 = floresInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(oilCompIdx)];
807 if (!this->floresn_[waterCompIdx].indices.empty()) {
808 this->floresn_[waterCompIdx].indices[floresInfo.nncId] = floresInfo.nncId;
809 this->floresn_[waterCompIdx].values[floresInfo.nncId]
810 = floresInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(waterCompIdx)];
818 void processElementBlockData(
const ElementContext& elemCtx)
820 OPM_TIMEBLOCK_LOCAL(processElementBlockData);
821 if (!std::is_same<Discretization, EcfvDiscretization<TypeTag>>::value)
824 const auto& problem = elemCtx.simulator().problem();
825 for (
unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(0); ++dofIdx) {
827 const auto globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, 0);
828 const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx);
829 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, 0);
830 const auto& fs = intQuants.fluidState();
831 for (
auto& val : this->blockData_) {
832 const auto& key = val.first;
833 assert(key.second > 0);
835 const auto cartesianIdxBlock =
static_cast<std::remove_cv_t<
836 std::remove_reference_t<decltype(cartesianIdx)
>>>(key.second - 1);
838 if (cartesianIdx == cartesianIdxBlock) {
839 if ((key.first ==
"BWSAT") || (key.first ==
"BSWAT"))
840 val.second = getValue(fs.saturation(waterPhaseIdx));
841 else if ((key.first ==
"BGSAT") || (key.first ==
"BSGAS"))
842 val.second = getValue(fs.saturation(gasPhaseIdx));
843 else if ((key.first ==
"BOSAT") || (key.first ==
"BSOIL"))
844 val.second = getValue(fs.saturation(oilPhaseIdx));
845 else if (key.first ==
"BNSAT")
846 val.second = intQuants.solventSaturation().value();
847 else if ((key.first ==
"BPR") || (key.first ==
"BPRESSUR")) {
848 if (FluidSystem::phaseIsActive(oilPhaseIdx))
849 val.second = getValue(fs.pressure(oilPhaseIdx));
850 else if (FluidSystem::phaseIsActive(gasPhaseIdx))
851 val.second = getValue(fs.pressure(gasPhaseIdx));
852 else if (FluidSystem::phaseIsActive(waterPhaseIdx))
853 val.second = getValue(fs.pressure(waterPhaseIdx));
855 else if ((key.first ==
"BTCNFHEA") || (key.first ==
"BTEMP")) {
856 if (FluidSystem::phaseIsActive(oilPhaseIdx))
857 val.second = getValue(fs.temperature(oilPhaseIdx));
858 else if (FluidSystem::phaseIsActive(gasPhaseIdx))
859 val.second = getValue(fs.temperature(gasPhaseIdx));
860 else if (FluidSystem::phaseIsActive(waterPhaseIdx))
861 val.second = getValue(fs.temperature(waterPhaseIdx));
863 else if (key.first ==
"BWKR" || key.first ==
"BKRW")
864 val.second = getValue(intQuants.relativePermeability(waterPhaseIdx));
865 else if (key.first ==
"BGKR" || key.first ==
"BKRG")
866 val.second = getValue(intQuants.relativePermeability(gasPhaseIdx));
867 else if (key.first ==
"BOKR" || key.first ==
"BKRO")
868 val.second = getValue(intQuants.relativePermeability(oilPhaseIdx));
869 else if (key.first ==
"BKROG") {
870 const auto& materialParams = problem.materialLawParams(elemCtx, dofIdx, 0);
872 = MaterialLaw::template relpermOilInOilGasSystem<Evaluation>(materialParams, fs);
873 val.second = getValue(krog);
875 else if (key.first ==
"BKROW") {
876 const auto& materialParams = problem.materialLawParams(elemCtx, dofIdx, 0);
878 = MaterialLaw::template relpermOilInOilWaterSystem<Evaluation>(materialParams, fs);
879 val.second = getValue(krow);
881 else if (key.first ==
"BWPC")
882 val.second = getValue(fs.pressure(oilPhaseIdx)) - getValue(fs.pressure(waterPhaseIdx));
883 else if (key.first ==
"BGPC")
884 val.second = getValue(fs.pressure(gasPhaseIdx)) - getValue(fs.pressure(oilPhaseIdx));
885 else if (key.first ==
"BWPR")
886 val.second = getValue(fs.pressure(waterPhaseIdx));
887 else if (key.first ==
"BGPR")
888 val.second = getValue(fs.pressure(gasPhaseIdx));
889 else if (key.first ==
"BVWAT" || key.first ==
"BWVIS")
890 val.second = getValue(fs.viscosity(waterPhaseIdx));
891 else if (key.first ==
"BVGAS" || key.first ==
"BGVIS")
892 val.second = getValue(fs.viscosity(gasPhaseIdx));
893 else if (key.first ==
"BVOIL" || key.first ==
"BOVIS")
894 val.second = getValue(fs.viscosity(oilPhaseIdx));
895 else if ((key.first ==
"BODEN") || (key.first ==
"BDENO"))
896 val.second = getValue(fs.density(oilPhaseIdx));
897 else if ((key.first ==
"BGDEN") || (key.first ==
"BDENG"))
898 val.second = getValue(fs.density(gasPhaseIdx));
899 else if ((key.first ==
"BWDEN") || (key.first ==
"BDENW"))
900 val.second = getValue(fs.density(waterPhaseIdx));
901 else if ((key.first ==
"BRPV") ||
902 (key.first ==
"BOPV") ||
903 (key.first ==
"BWPV") ||
904 (key.first ==
"BGPV"))
906 if (key.first ==
"BRPV") {
909 else if (key.first ==
"BOPV") {
910 val.second = getValue(fs.saturation(oilPhaseIdx));
912 else if (key.first ==
"BWPV") {
913 val.second = getValue(fs.saturation(waterPhaseIdx));
916 val.second = getValue(fs.saturation(gasPhaseIdx));
920 val.second *= getValue(intQuants.porosity())
921 * elemCtx.simulator().model().dofTotalVolume(globalDofIdx);
923 else if (key.first ==
"BRS")
924 val.second = getValue(fs.Rs());
925 else if (key.first ==
"BRV")
926 val.second = getValue(fs.Rv());
927 else if ((key.first ==
"BOIP") || (key.first ==
"BOIPL") || (key.first ==
"BOIPG") ||
928 (key.first ==
"BGIP") || (key.first ==
"BGIPL") || (key.first ==
"BGIPG") ||
929 (key.first ==
"BWIP"))
931 if ((key.first ==
"BOIP") || (key.first ==
"BOIPL")) {
932 val.second = getValue(fs.invB(oilPhaseIdx)) * getValue(fs.saturation(oilPhaseIdx));
934 if (key.first ==
"BOIP") {
935 val.second += getValue(fs.Rv()) * getValue(fs.invB(gasPhaseIdx))
936 * getValue(fs.saturation(gasPhaseIdx));
939 else if (key.first ==
"BOIPG") {
940 val.second = getValue(fs.Rv()) * getValue(fs.invB(gasPhaseIdx))
941 * getValue(fs.saturation(gasPhaseIdx));
943 else if ((key.first ==
"BGIP") || (key.first ==
"BGIPG")) {
944 val.second = getValue(fs.invB(gasPhaseIdx)) * getValue(fs.saturation(gasPhaseIdx));
946 if (key.first ==
"BGIP") {
947 if (!FluidSystem::phaseIsActive(oilPhaseIdx)) {
948 val.second += getValue(fs.Rsw()) * getValue(fs.invB(waterPhaseIdx))
949 * getValue(fs.saturation(waterPhaseIdx));
952 val.second += getValue(fs.Rs()) * getValue(fs.invB(oilPhaseIdx))
953 * getValue(fs.saturation(oilPhaseIdx));
957 else if (key.first ==
"BGIPL") {
958 if (!FluidSystem::phaseIsActive(oilPhaseIdx)) {
959 val.second = getValue(fs.Rsw()) * getValue(fs.invB(waterPhaseIdx))
960 * getValue(fs.saturation(waterPhaseIdx));
963 val.second = getValue(fs.Rs()) * getValue(fs.invB(oilPhaseIdx))
964 * getValue(fs.saturation(oilPhaseIdx));
968 val.second = getValue(fs.invB(waterPhaseIdx)) * getValue(fs.saturation(waterPhaseIdx));
972 val.second *= elemCtx.simulator().model().dofTotalVolume(globalDofIdx)
973 * getValue(intQuants.porosity());
975 else if ((key.first ==
"BPPO") ||
976 (key.first ==
"BPPG") ||
977 (key.first ==
"BPPW"))
979 auto phase = RegionPhasePoreVolAverage::Phase{};
981 if (key.first ==
"BPPO") {
982 phase.ix = oilPhaseIdx;
984 else if (key.first ==
"BPPG") {
985 phase.ix = gasPhaseIdx;
988 phase.ix = waterPhaseIdx;
998 const auto datum = this->eclState_.getSimulationConfig()
999 .datumDepths()(this->regions_[
"FIPNUM"][dofIdx] - 1);
1002 const auto region = RegionPhasePoreVolAverage::Region {
1003 elemCtx.primaryVars(dofIdx, 0).pvtRegionIndex() + 1
1006 const auto density = this->regionAvgDensity_
1007 ->value(
"PVTNUM", phase, region);
1009 const auto press = getValue(fs.pressure(phase.ix));
1011 elemCtx.problem().gravity()[GridView::dimensionworld - 1];
1012 const auto dz = problem.dofCenterDepth(globalDofIdx) - datum;
1014 val.second = press - density*dz*grav;
1016 else if ((key.first ==
"BFLOWI") ||
1017 (key.first ==
"BFLOWJ") ||
1018 (key.first ==
"BFLOWK"))
1020 auto dir = FaceDir::ToIntersectionIndex(Dir::XPlus);
1022 if (key.first ==
"BFLOWJ") {
1023 dir = FaceDir::ToIntersectionIndex(Dir::YPlus);
1025 else if (key.first ==
"BFLOWK") {
1026 dir = FaceDir::ToIntersectionIndex(Dir::ZPlus);
1029 val.second = this->flows_[dir][waterCompIdx][globalDofIdx];
1032 std::string logstring =
"Keyword '";
1033 logstring.append(key.first);
1034 logstring.append(
"' is unhandled for output to summary file.");
1035 OpmLog::warning(
"Unhandled output keyword", logstring);
1070 template <
class ActiveIndex,
class CartesianIndex>
1072 ActiveIndex&& activeIndex,
1073 CartesianIndex&& cartesianIndex)
1076 const auto identifyCell = [&activeIndex, &cartesianIndex](
const Element& elem)
1079 const auto cellIndex = activeIndex(elem);
1082 static_cast<int>(cellIndex),
1083 cartesianIndex(cellIndex),
1084 elem.partitionType() == Dune::InteriorEntity
1088 const auto timeIdx = 0u;
1089 const auto& stencil = elemCtx.stencil(timeIdx);
1090 const auto numInteriorFaces = elemCtx.numInteriorFaces(timeIdx);
1092 for (
auto scvfIdx = 0 * numInteriorFaces; scvfIdx < numInteriorFaces; ++scvfIdx) {
1093 const auto& face = stencil.interiorFace(scvfIdx);
1094 const auto left = identifyCell(stencil.element(face.interiorIndex()));
1095 const auto right = identifyCell(stencil.element(face.exteriorIndex()));
1097 const auto rates = this->
1098 getComponentSurfaceRates(elemCtx, face.area(), scvfIdx, timeIdx);
1112 this->interRegionFlows_.
clear();
1120 this->interRegionFlows_.
compress();
1128 return this->interRegionFlows_;
1131 template <
class Flu
idState>
1132 void assignToFluidState(FluidState& fs,
unsigned elemIdx)
const
1134 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1135 if (this->saturation_[phaseIdx].empty())
1138 fs.setSaturation(phaseIdx, this->saturation_[phaseIdx][elemIdx]);
1141 if (!this->fluidPressure_.empty()) {
1144 std::array<Scalar, numPhases> pc = {0};
1145 const MaterialLawParams& matParams = simulator_.problem().materialLawParams(elemIdx);
1146 MaterialLaw::capillaryPressures(pc, matParams, fs);
1147 Valgrind::CheckDefined(this->fluidPressure_[elemIdx]);
1148 Valgrind::CheckDefined(pc);
1149 const auto& pressure = this->fluidPressure_[elemIdx];
1150 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1151 if (!FluidSystem::phaseIsActive(phaseIdx))
1154 if (Indices::oilEnabled)
1155 fs.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx]));
1156 else if (Indices::gasEnabled)
1157 fs.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx]));
1158 else if (Indices::waterEnabled)
1160 fs.setPressure(phaseIdx, pressure);
1164 if (!this->temperature_.empty())
1165 fs.setTemperature(this->temperature_[elemIdx]);
1166 if (!this->rs_.empty())
1167 fs.setRs(this->rs_[elemIdx]);
1168 if (!this->rsw_.empty())
1169 fs.setRsw(this->rsw_[elemIdx]);
1170 if (!this->rv_.empty())
1171 fs.setRv(this->rv_[elemIdx]);
1172 if (!this->rvw_.empty())
1173 fs.setRvw(this->rvw_[elemIdx]);
1176 void initHysteresisParams(Simulator& simulator,
unsigned elemIdx)
const
1178 if (!this->soMax_.empty())
1179 simulator.problem().setMaxOilSaturation(elemIdx, this->soMax_[elemIdx]);
1181 if (simulator.problem().materialLawManager()->enableHysteresis()) {
1182 auto matLawManager = simulator.problem().materialLawManager();
1184 if (FluidSystem::phaseIsActive(oilPhaseIdx)
1185 && FluidSystem::phaseIsActive(waterPhaseIdx)) {
1187 Scalar swmax = -2.0;
1190 if (matLawManager->enableNonWettingHysteresis()) {
1191 if (!this->soMax_.empty()) {
1192 somax = this->soMax_[elemIdx];
1195 if (matLawManager->enableWettingHysteresis()) {
1196 if (!this->swMax_.empty()) {
1197 swmax = this->swMax_[elemIdx];
1200 if (matLawManager->enablePCHysteresis()) {
1201 if (!this->swmin_.empty()) {
1202 swmin = this->swmin_[elemIdx];
1205 matLawManager->setOilWaterHysteresisParams(
1206 somax, swmax, swmin, elemIdx);
1208 if (FluidSystem::phaseIsActive(oilPhaseIdx)
1209 && FluidSystem::phaseIsActive(gasPhaseIdx)) {
1211 Scalar shmax = -2.0;
1214 if (matLawManager->enableNonWettingHysteresis()) {
1215 if (!this->sgmax_.empty()) {
1216 sgmax = this->sgmax_[elemIdx];
1219 if (matLawManager->enableWettingHysteresis()) {
1220 if (!this->shmax_.empty()) {
1221 shmax = this->shmax_[elemIdx];
1224 if (matLawManager->enablePCHysteresis()) {
1225 if (!this->somin_.empty()) {
1226 somin = this->somin_[elemIdx];
1229 matLawManager->setGasOilHysteresisParams(
1230 sgmax, shmax, somin, elemIdx);
1235 if (simulator_.vanguard().eclState().fieldProps().has_double(
"SWATINIT")) {
1236 simulator.problem().materialLawManager()
1237 ->applyRestartSwatInit(elemIdx, this->ppcw_[elemIdx]);
1241 void updateFluidInPlace(
const ElementContext& elemCtx)
1243 for (
unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(0); ++dofIdx) {
1244 updateFluidInPlace_(elemCtx, dofIdx);
1248 void updateFluidInPlace(
const unsigned globalDofIdx,
1249 const IntensiveQuantities& intQuants,
1250 const double totVolume)
1252 this->updateFluidInPlace_(globalDofIdx, intQuants, totVolume);
1256 bool isDefunctParallelWell(std::string wname)
const override
1258 if (simulator_.gridView().comm().size() == 1)
1260 const auto& parallelWells = simulator_.vanguard().parallelWells();
1261 std::pair<std::string, bool> value {wname,
true};
1262 auto candidate = std::lower_bound(parallelWells.begin(), parallelWells.end(), value);
1263 return candidate == parallelWells.end() || *candidate != value;
1266 void updateFluidInPlace_(
const ElementContext& elemCtx,
const unsigned dofIdx)
1268 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, 0);
1269 const unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, 0);
1270 const auto totVolume = elemCtx.simulator().model().dofTotalVolume(globalDofIdx);
1272 this->updateFluidInPlace_(globalDofIdx, intQuants, totVolume);
1275 void updateFluidInPlace_(
const unsigned globalDofIdx,
1276 const IntensiveQuantities& intQuants,
1277 const double totVolume)
1279 OPM_TIMEBLOCK_LOCAL(updateFluidInPlace);
1281 this->updateTotalVolumesAndPressures_(globalDofIdx, intQuants, totVolume);
1283 if (this->computeFip_) {
1284 this->updatePhaseInplaceVolumes_(globalDofIdx, intQuants, totVolume);
1288 void createLocalRegion_(std::vector<int>& region)
1290 std::size_t elemIdx = 0;
1291 for (
const auto& elem : elements(simulator_.gridView())) {
1292 if (elem.partitionType() != Dune::InteriorEntity) {
1293 region[elemIdx] = 0;
1300 template <
typename Flu
idState>
1301 void aggregateAverageDensityContributions_(
const FluidState& fs,
1302 const unsigned int globalDofIdx,
1305 auto pvCellValue = RegionPhasePoreVolAverage::CellValue{};
1306 pvCellValue.porv = porv;
1308 for (
auto phaseIdx = 0*FluidSystem::numPhases;
1309 phaseIdx < FluidSystem::numPhases; ++phaseIdx)
1311 if (! FluidSystem::phaseIsActive(phaseIdx)) {
1315 pvCellValue.value = getValue(fs.density(phaseIdx));
1316 pvCellValue.sat = getValue(fs.saturation(phaseIdx));
1318 this->regionAvgDensity_
1319 ->addCell(globalDofIdx,
1320 RegionPhasePoreVolAverage::Phase { phaseIdx },
1340 data::InterRegFlowMap::FlowRates
1341 getComponentSurfaceRates(
const ElementContext& elemCtx,
1342 const Scalar faceArea,
1343 const std::size_t scvfIdx,
1344 const std::size_t timeIdx)
const
1346 using Component = data::InterRegFlowMap::Component;
1348 auto rates = data::InterRegFlowMap::FlowRates {};
1350 const auto& extQuant = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
1352 const auto alpha = getValue(extQuant.extrusionFactor()) * faceArea;
1354 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1355 const auto& up = elemCtx
1356 .intensiveQuantities(extQuant.upstreamIndex(oilPhaseIdx), timeIdx);
1358 using FluidState = std::remove_cv_t<std::remove_reference_t<
1359 decltype(up.fluidState())>>;
1361 const auto pvtReg = up.pvtRegionIndex();
1363 const auto bO = getValue(getInvB_<FluidSystem, FluidState, Scalar>
1364 (up.fluidState(), oilPhaseIdx, pvtReg));
1366 const auto qO = alpha * bO * getValue(extQuant.volumeFlux(oilPhaseIdx));
1368 rates[Component::Oil] += qO;
1370 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1371 const auto Rs = getValue(
1372 BlackOil::getRs_<FluidSystem, FluidState, Scalar>
1373 (up.fluidState(), pvtReg));
1375 rates[Component::Gas] += qO * Rs;
1376 rates[Component::Disgas] += qO * Rs;
1380 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1381 const auto& up = elemCtx
1382 .intensiveQuantities(extQuant.upstreamIndex(gasPhaseIdx), timeIdx);
1384 using FluidState = std::remove_cv_t<std::remove_reference_t<
1385 decltype(up.fluidState())>>;
1387 const auto pvtReg = up.pvtRegionIndex();
1389 const auto bG = getValue(getInvB_<FluidSystem, FluidState, Scalar>
1390 (up.fluidState(), gasPhaseIdx, pvtReg));
1392 const auto qG = alpha * bG * getValue(extQuant.volumeFlux(gasPhaseIdx));
1394 rates[Component::Gas] += qG;
1396 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1397 const auto Rv = getValue(
1398 BlackOil::getRv_<FluidSystem, FluidState, Scalar>
1399 (up.fluidState(), pvtReg));
1401 rates[Component::Oil] += qG * Rv;
1402 rates[Component::Vapoil] += qG * Rv;
1406 if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
1407 const auto& up = elemCtx
1408 .intensiveQuantities(extQuant.upstreamIndex(waterPhaseIdx), timeIdx);
1410 using FluidState = std::remove_cv_t<std::remove_reference_t<
1411 decltype(up.fluidState())>>;
1413 const auto pvtReg = up.pvtRegionIndex();
1415 const auto bW = getValue(getInvB_<FluidSystem, FluidState, Scalar>
1416 (up.fluidState(), waterPhaseIdx, pvtReg));
1418 rates[Component::Water] +=
1419 alpha * bW * getValue(extQuant.volumeFlux(waterPhaseIdx));
1425 template <
typename Flu
idState>
1426 Scalar hydroCarbonFraction(
const FluidState& fs)
const
1428 if (this->eclState_.runspec().co2Storage()) {
1435 auto hydrocarbon = Scalar {0};
1436 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1437 hydrocarbon += getValue(fs.saturation(oilPhaseIdx));
1440 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1441 hydrocarbon += getValue(fs.saturation(gasPhaseIdx));
1447 void updateTotalVolumesAndPressures_(
const unsigned globalDofIdx,
1448 const IntensiveQuantities& intQuants,
1449 const double totVolume)
1451 const auto& fs = intQuants.fluidState();
1453 const double pv = totVolume * intQuants.porosity().value();
1454 const auto hydrocarbon = this->hydroCarbonFraction(fs);
1456 if (! this->hydrocarbonPoreVolume_.empty()) {
1457 this->fip_[Inplace::Phase::PoreVolume][globalDofIdx] =
1458 totVolume * intQuants.referencePorosity();
1460 this->dynamicPoreVolume_[globalDofIdx] = pv;
1461 this->hydrocarbonPoreVolume_[globalDofIdx] = pv * hydrocarbon;
1464 if (!this->pressureTimesHydrocarbonVolume_.empty() &&
1465 !this->pressureTimesPoreVolume_.empty())
1467 assert(this->hydrocarbonPoreVolume_.size() == this->pressureTimesHydrocarbonVolume_.size());
1468 assert(this->fip_[Inplace::Phase::PoreVolume].size() == this->pressureTimesPoreVolume_.size());
1470 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1471 this->pressureTimesPoreVolume_[globalDofIdx] =
1472 getValue(fs.pressure(oilPhaseIdx)) * pv;
1474 this->pressureTimesHydrocarbonVolume_[globalDofIdx] =
1475 this->pressureTimesPoreVolume_[globalDofIdx] * hydrocarbon;
1477 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1478 this->pressureTimesPoreVolume_[globalDofIdx] =
1479 getValue(fs.pressure(gasPhaseIdx)) * pv;
1481 this->pressureTimesHydrocarbonVolume_[globalDofIdx] =
1482 this->pressureTimesPoreVolume_[globalDofIdx] * hydrocarbon;
1484 else if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
1485 this->pressureTimesPoreVolume_[globalDofIdx] =
1486 getValue(fs.pressure(waterPhaseIdx)) * pv;
1491 void updatePhaseInplaceVolumes_(
const unsigned globalDofIdx,
1492 const IntensiveQuantities& intQuants,
1493 const double totVolume)
1495 std::array<Scalar, FluidSystem::numPhases> fip {};
1496 std::array<Scalar, FluidSystem::numPhases> fipr{};
1498 const auto& fs = intQuants.fluidState();
1499 const auto pv = totVolume * intQuants.porosity().value();
1501 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1502 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1506 const auto b = getValue(fs.invB(phaseIdx));
1507 const auto s = getValue(fs.saturation(phaseIdx));
1509 fipr[phaseIdx] = s * pv;
1510 fip [phaseIdx] = b * fipr[phaseIdx];
1513 this->updateInplaceVolumesSurface(globalDofIdx, fip);
1514 this->updateInplaceVolumesReservoir(globalDofIdx, fs, fipr);
1516 if (FluidSystem::phaseIsActive(oilPhaseIdx) &&
1517 FluidSystem::phaseIsActive(gasPhaseIdx))
1519 this->updateOilGasDistribution(globalDofIdx, fs, fip);
1522 if (FluidSystem::phaseIsActive(waterPhaseIdx) &&
1523 FluidSystem::phaseIsActive(gasPhaseIdx))
1525 this->updateGasWaterDistribution(globalDofIdx, fs, fip);
1528 if (FluidSystem::phaseIsActive(gasPhaseIdx) &&
1529 (!this->fip_[Inplace::Phase::CO2InGasPhaseInMob].empty() ||
1530 !this->fip_[Inplace::Phase::CO2InGasPhaseMob].empty() ||
1531 !this->fip_[Inplace::Phase::CO2MassInGasPhaseInMob].empty() ||
1532 !this->fip_[Inplace::Phase::CO2MassInGasPhaseMob].empty() ||
1533 !this->fip_[Inplace::Phase::CO2Mass].empty() ||
1534 !this->fip_[Inplace::Phase::CO2MassInGasPhase].empty() ||
1535 !this->fip_[Inplace::Phase::CO2InGasPhaseInMobKrg].empty() ||
1536 !this->fip_[Inplace::Phase::CO2InGasPhaseMobKrg].empty() ||
1537 !this->fip_[Inplace::Phase::CO2MassInGasPhaseInMobKrg].empty() ||
1538 !this->fip_[Inplace::Phase::CO2MassInGasPhaseMobKrg].empty() ||
1539 !this->fip_[Inplace::Phase::CO2MassInGasPhaseEffectiveTrapped].empty() ||
1540 !this->fip_[Inplace::Phase::CO2MassInGasPhaseEffectiveUnTrapped].empty() ||
1541 !this->fip_[Inplace::Phase::CO2MassInGasPhaseMaximumTrapped].empty() ||
1542 !this->fip_[Inplace::Phase::CO2MassInGasPhaseMaximumUnTrapped].empty()))
1544 this->updateCO2InGas(globalDofIdx, pv, intQuants);
1547 if ((!this->fip_[Inplace::Phase::CO2InWaterPhase].empty() ||
1548 !this->fip_[Inplace::Phase::CO2MassInWaterPhase].empty() ||
1549 !this->fip_[Inplace::Phase::CO2Mass].empty()) &&
1550 (FluidSystem::phaseIsActive(waterPhaseIdx) ||
1551 FluidSystem::phaseIsActive(oilPhaseIdx)))
1553 this->updateCO2InWater(globalDofIdx, pv, fs);
1557 template <
typename FIPArray>
1558 void updateInplaceVolumesSurface(
const unsigned globalDofIdx,
1559 const FIPArray& fip)
1561 if (FluidSystem::phaseIsActive(oilPhaseIdx) &&
1562 !this->fip_[Inplace::Phase::OIL].empty())
1564 this->fip_[Inplace::Phase::OIL][globalDofIdx] = fip[oilPhaseIdx];
1567 if (FluidSystem::phaseIsActive(oilPhaseIdx) &&
1568 !this->fip_[Inplace::Phase::OilInLiquidPhase].empty())
1570 this->fip_[Inplace::Phase::OilInLiquidPhase][globalDofIdx] = fip[oilPhaseIdx];
1573 if (FluidSystem::phaseIsActive(gasPhaseIdx) &&
1574 !this->fip_[Inplace::Phase::GAS].empty())
1576 this->fip_[Inplace::Phase::GAS][globalDofIdx] = fip[gasPhaseIdx];
1579 if (FluidSystem::phaseIsActive(gasPhaseIdx) &&
1580 !this->fip_[Inplace::Phase::GasInGasPhase].empty())
1582 this->fip_[Inplace::Phase::GasInGasPhase][globalDofIdx] = fip[gasPhaseIdx];
1585 if (FluidSystem::phaseIsActive(waterPhaseIdx) &&
1586 !this->fip_[Inplace::Phase::WATER].empty())
1588 this->fip_[Inplace::Phase::WATER][globalDofIdx] = fip[waterPhaseIdx];
1592 template <
typename Flu
idState,
typename FIPArray>
1593 void updateInplaceVolumesReservoir(
const unsigned globalDofIdx,
1594 const FluidState& fs,
1595 const FIPArray& fipr)
1597 if (FluidSystem::phaseIsActive(oilPhaseIdx) &&
1598 !this->fip_[Inplace::Phase::OilResVolume].empty())
1600 this->fip_[Inplace::Phase::OilResVolume][globalDofIdx] = fipr[oilPhaseIdx];
1603 if (FluidSystem::phaseIsActive(gasPhaseIdx) &&
1604 !this->fip_[Inplace::Phase::GasResVolume].empty())
1606 this->fip_[Inplace::Phase::GasResVolume][globalDofIdx] = fipr[gasPhaseIdx];
1609 if (FluidSystem::phaseIsActive(waterPhaseIdx) &&
1610 !this->fip_[Inplace::Phase::WaterResVolume].empty())
1612 this->fip_[Inplace::Phase::WaterResVolume][globalDofIdx] = fipr[waterPhaseIdx];
1615 if (FluidSystem::phaseIsActive(waterPhaseIdx) &&
1616 !this->fip_[Inplace::Phase::SALT].empty())
1618 this->fip_[Inplace::Phase::SALT][globalDofIdx] =
1619 fipr[waterPhaseIdx] * fs.saltConcentration().value();
1623 template <
typename Flu
idState,
typename FIPArray>
1624 void updateOilGasDistribution(
const unsigned globalDofIdx,
1625 const FluidState& fs,
1626 const FIPArray& fip)
1629 const auto gasInPlaceLiquid = getValue(fs.Rs()) * fip[oilPhaseIdx];
1630 const auto oilInPlaceGas = getValue(fs.Rv()) * fip[gasPhaseIdx];
1632 if (!this->fip_[Inplace::Phase::GasInLiquidPhase].empty()) {
1633 this->fip_[Inplace::Phase::GasInLiquidPhase][globalDofIdx] = gasInPlaceLiquid;
1636 if (!this->fip_[Inplace::Phase::OilInGasPhase].empty()) {
1637 this->fip_[Inplace::Phase::OilInGasPhase][globalDofIdx] = oilInPlaceGas;
1641 if (!this->fip_[Inplace::Phase::OIL].empty()) {
1642 this->fip_[Inplace::Phase::OIL][globalDofIdx] += oilInPlaceGas;
1645 if (!this->fip_[Inplace::Phase::GAS].empty()) {
1646 this->fip_[Inplace::Phase::GAS][globalDofIdx] += gasInPlaceLiquid;
1650 template <
typename Flu
idState,
typename FIPArray>
1651 void updateGasWaterDistribution(
const unsigned globalDofIdx,
1652 const FluidState& fs,
1653 const FIPArray& fip)
1656 const auto gasInPlaceWater = getValue(fs.Rsw()) * fip[waterPhaseIdx];
1657 const auto waterInPlaceGas = getValue(fs.Rvw()) * fip[gasPhaseIdx];
1659 if (!this->fip_[Inplace::Phase::WaterInGasPhase].empty()) {
1660 this->fip_[Inplace::Phase::WaterInGasPhase][globalDofIdx] = waterInPlaceGas;
1663 if (!this->fip_[Inplace::Phase::WaterInWaterPhase].empty()) {
1664 this->fip_[Inplace::Phase::WaterInWaterPhase][globalDofIdx] = fip[waterPhaseIdx];
1668 if (!this->fip_[Inplace::Phase::GasInLiquidPhase].empty() &&
1669 !FluidSystem::phaseIsActive(oilPhaseIdx))
1671 this->fip_[Inplace::Phase::GasInLiquidPhase][globalDofIdx] = gasInPlaceWater;
1675 if (!this->fip_[Inplace::Phase::WATER].empty()) {
1676 this->fip_[Inplace::Phase::WATER][globalDofIdx] += waterInPlaceGas;
1679 if (!this->fip_[Inplace::Phase::GAS].empty()) {
1680 this->fip_[Inplace::Phase::GAS][globalDofIdx] += gasInPlaceWater;
1684 template <
typename IntensiveQuantities>
1685 void updateCO2InGas(
const unsigned globalDofIdx,
1687 const IntensiveQuantities& intQuants)
1689 const auto& scaledDrainageInfo = this->simulator_.problem().materialLawManager()
1690 ->oilWaterScaledEpsInfoDrainage(globalDofIdx);
1692 const auto& fs = intQuants.fluidState();
1693 Scalar sgcr = scaledDrainageInfo.Sgcr;
1694 if (this->simulator_.problem().materialLawManager()->enableHysteresis()) {
1695 const auto& matParams = simulator_.problem().materialLawParams(globalDofIdx);
1696 sgcr = MaterialLaw::trappedGasSaturation(matParams,
false);
1699 const Scalar sg = getValue(fs.saturation(gasPhaseIdx));
1700 const Scalar rhog = getValue(fs.density(gasPhaseIdx));
1701 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx)
1702 ? FluidSystem::convertRvwToXgW(getValue(fs.Rvw()), fs.pvtRegionIndex())
1703 : FluidSystem::convertRvToXgO(getValue(fs.Rv()), fs.pvtRegionIndex());
1705 const Scalar mM = FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex());
1706 const Scalar massGas = (1 - xgW) * pv * rhog;
1707 if (!this->fip_[Inplace::Phase::CO2Mass].empty()) {
1708 this->fip_[Inplace::Phase::CO2Mass][globalDofIdx] = massGas * sg;
1711 if (!this->fip_[Inplace::Phase::CO2MassInGasPhase].empty()) {
1712 this->fip_[Inplace::Phase::CO2MassInGasPhase][globalDofIdx] = massGas * sg;
1715 if (!this->fip_[Inplace::Phase::CO2InGasPhaseInMob].empty()) {
1716 const Scalar imMobileGas = massGas / mM * std::min(sgcr , sg);
1717 this->fip_[Inplace::Phase::CO2InGasPhaseInMob][globalDofIdx] = imMobileGas;
1720 if (!this->fip_[Inplace::Phase::CO2InGasPhaseMob].empty()) {
1721 const Scalar mobileGas = massGas / mM * std::max(Scalar{0.0}, sg - sgcr);
1722 this->fip_[Inplace::Phase::CO2InGasPhaseMob][globalDofIdx] = mobileGas;
1725 if (!this->fip_[Inplace::Phase::CO2InGasPhaseInMobKrg].empty()) {
1727 const Scalar imMobileGasKrg = massGas / mM * sg;
1728 this->fip_[Inplace::Phase::CO2InGasPhaseInMobKrg][globalDofIdx] = imMobileGasKrg;
1730 this->fip_[Inplace::Phase::CO2InGasPhaseInMobKrg][globalDofIdx] = 0;
1734 if (!this->fip_[Inplace::Phase::CO2InGasPhaseMobKrg].empty()) {
1736 const Scalar mobileGasKrg = massGas / mM * sg;
1737 this->fip_[Inplace::Phase::CO2InGasPhaseMobKrg][globalDofIdx] = mobileGasKrg;
1739 this->fip_[Inplace::Phase::CO2InGasPhaseMobKrg][globalDofIdx] = 0;
1743 if (!this->fip_[Inplace::Phase::CO2MassInGasPhaseInMob].empty()) {
1744 const Scalar imMobileMassGas = massGas * std::min(sgcr , sg);
1745 this->fip_[Inplace::Phase::CO2MassInGasPhaseInMob][globalDofIdx] = imMobileMassGas;
1748 if (!this->fip_[Inplace::Phase::CO2MassInGasPhaseMob].empty()) {
1749 const Scalar mobileMassGas = massGas * std::max(Scalar{0.0}, sg - sgcr);
1750 this->fip_[Inplace::Phase::CO2MassInGasPhaseMob][globalDofIdx] = mobileMassGas;
1753 if (!this->fip_[Inplace::Phase::CO2MassInGasPhaseInMobKrg].empty()) {
1755 const Scalar imMobileMassGasKrg = massGas * sg;
1756 this->fip_[Inplace::Phase::CO2MassInGasPhaseInMobKrg][globalDofIdx] = imMobileMassGasKrg;
1758 this->fip_[Inplace::Phase::CO2MassInGasPhaseInMobKrg][globalDofIdx] = 0;
1762 if (!this->fip_[Inplace::Phase::CO2MassInGasPhaseMobKrg].empty()) {
1764 const Scalar mobileMassGasKrg = massGas * sg;
1765 this->fip_[Inplace::Phase::CO2MassInGasPhaseMobKrg][globalDofIdx] = mobileMassGasKrg;
1767 this->fip_[Inplace::Phase::CO2MassInGasPhaseMobKrg][globalDofIdx] = 0;
1771 if (!this->fip_[Inplace::Phase::CO2MassInGasPhaseMaximumTrapped].empty() ||
1772 !this->fip_[Inplace::Phase::CO2MassInGasPhaseMaximumUnTrapped].empty() ) {
1773 Scalar trappedGasSaturation = scaledDrainageInfo.Sgcr;
1774 if (this->simulator_.problem().materialLawManager()->enableHysteresis()) {
1775 const auto& matParams = simulator_.problem().materialLawParams(globalDofIdx);
1777 trappedGasSaturation = MaterialLaw::trappedGasSaturation(matParams,
true);
1779 if (!this->fip_[Inplace::Phase::CO2MassInGasPhaseMaximumTrapped].empty()) {
1780 const Scalar imMobileMassGas = massGas * std::min(trappedGasSaturation , sg);
1781 this->fip_[Inplace::Phase::CO2MassInGasPhaseMaximumTrapped][globalDofIdx] = imMobileMassGas;
1783 if (!this->fip_[Inplace::Phase::CO2MassInGasPhaseMaximumUnTrapped].empty()) {
1784 const Scalar mobileMassGas = massGas * std::max(Scalar{0.0}, sg - trappedGasSaturation);
1785 this->fip_[Inplace::Phase::CO2MassInGasPhaseMaximumUnTrapped][globalDofIdx] = mobileMassGas;
1789 if (!this->fip_[Inplace::Phase::CO2MassInGasPhaseEffectiveTrapped].empty() ||
1790 !this->fip_[Inplace::Phase::CO2MassInGasPhaseEffectiveUnTrapped].empty()) {
1791 Scalar trappedGasSaturation = scaledDrainageInfo.Sgcr;
1792 if (this->simulator_.problem().materialLawManager()->enableHysteresis()) {
1793 const auto& matParams = simulator_.problem().materialLawParams(globalDofIdx);
1794 const double krg = getValue(intQuants.relativePermeability(gasPhaseIdx));
1795 trappedGasSaturation = MaterialLaw::strandedGasSaturation(matParams, sg, krg);
1797 if (!this->fip_[Inplace::Phase::CO2MassInGasPhaseEffectiveTrapped].empty()) {
1798 const Scalar imMobileMassGas = massGas * std::min(trappedGasSaturation , sg);
1799 this->fip_[Inplace::Phase::CO2MassInGasPhaseEffectiveTrapped][globalDofIdx] = imMobileMassGas;
1801 if (!this->fip_[Inplace::Phase::CO2MassInGasPhaseEffectiveUnTrapped].empty()) {
1802 const Scalar mobileMassGas = massGas * std::max(Scalar{0.0}, sg - trappedGasSaturation);
1803 this->fip_[Inplace::Phase::CO2MassInGasPhaseEffectiveUnTrapped][globalDofIdx] = mobileMassGas;
1808 template <
typename Flu
idState>
1809 void updateCO2InWater(
const unsigned globalDofIdx,
1811 const FluidState& fs)
1813 const auto co2InWater = FluidSystem::phaseIsActive(oilPhaseIdx)
1814 ? this->co2InWaterFromOil(fs, pv)
1815 : this->co2InWaterFromWater(fs, pv);
1817 const Scalar mM = FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex());
1818 if (!this->fip_[Inplace::Phase::CO2Mass].empty()) {
1819 this->fip_[Inplace::Phase::CO2Mass][globalDofIdx] += co2InWater * mM;
1821 if (!this->fip_[Inplace::Phase::CO2MassInWaterPhase].empty()) {
1822 this->fip_[Inplace::Phase::CO2MassInWaterPhase][globalDofIdx] = co2InWater * mM;
1824 if (!this->fip_[Inplace::Phase::CO2InWaterPhase].empty()) {
1825 this->fip_[Inplace::Phase::CO2InWaterPhase][globalDofIdx] = co2InWater;
1829 template <
typename Flu
idState>
1830 Scalar co2InWaterFromWater(
const FluidState& fs,
const double pv)
const
1832 const double rhow = getValue(fs.density(waterPhaseIdx));
1833 const double sw = getValue(fs.saturation(waterPhaseIdx));
1834 const double xwG = FluidSystem::convertRswToXwG(getValue(fs.Rsw()), fs.pvtRegionIndex());
1836 const Scalar mM = FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex());
1838 return xwG * pv * rhow * sw / mM;
1841 template <
typename Flu
idState>
1842 Scalar co2InWaterFromOil(
const FluidState& fs,
const double pv)
const
1844 const double rhoo = getValue(fs.density(oilPhaseIdx));
1845 const double so = getValue(fs.saturation(oilPhaseIdx));
1846 const double xoG = FluidSystem::convertRsToXoG(getValue(fs.Rs()), fs.pvtRegionIndex());
1848 const Scalar mM = FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex());
1850 return xoG * pv * rhoo * so / mM;
1853 const Simulator& simulator_;