117 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
118 using Discretization = GetPropType<TypeTag, Properties::Discretization>;
119 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
120 using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
121 using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
122 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
123 using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
124 using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
125 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
126 using GridView = GetPropType<TypeTag, Properties::GridView>;
127 using Element =
typename GridView::template Codim<0>::Entity;
128 using ElementIterator =
typename GridView::template Codim<0>::Iterator;
130 using Indices = GetPropType<TypeTag, Properties::Indices>;
131 using Dir = FaceDir::DirEnum;
133 enum { conti0EqIdx = Indices::conti0EqIdx };
134 enum { numPhases = FluidSystem::numPhases };
135 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
136 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
137 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
138 enum { gasCompIdx = FluidSystem::gasCompIdx };
139 enum { oilCompIdx = FluidSystem::oilCompIdx };
140 enum { waterCompIdx = FluidSystem::waterCompIdx };
141 enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
144 template <
class CollectDataToIORankType>
146 const SummaryConfig& smryCfg,
147 const CollectDataToIORankType& collectToIORank)
148 : BaseType(simulator.vanguard().eclState(),
149 simulator.vanguard().schedule(),
151 simulator.vanguard().summaryState(),
153 getPropValue<TypeTag, Properties::EnableEnergy>(),
154 getPropValue<TypeTag, Properties::EnableTemperature>(),
155 getPropValue<TypeTag, Properties::EnableMech>(),
156 getPropValue<TypeTag, Properties::EnableSolvent>(),
157 getPropValue<TypeTag, Properties::EnablePolymer>(),
158 getPropValue<TypeTag, Properties::EnableFoam>(),
159 getPropValue<TypeTag, Properties::EnableBrine>(),
160 getPropValue<TypeTag, Properties::EnableSaltPrecipitation>(),
161 getPropValue<TypeTag, Properties::EnableExtbo>(),
162 getPropValue<TypeTag, Properties::EnableMICP>())
163 , simulator_(simulator)
165 for (
auto& region_pair : this->regions_) {
166 this->createLocalRegion_(region_pair.second);
169 auto isCartIdxOnThisRank = [&collectToIORank](
const int idx) {
170 return collectToIORank.isCartIdxOnThisRank(idx);
173 this->setupBlockData(isCartIdxOnThisRank);
175 this->forceDisableFipOutput_ =
176 Parameters::get<TypeTag, Properties::ForceDisableFluidInPlaceOutput>();
178 this->forceDisableFipresvOutput_ =
179 Parameters::get<TypeTag, Properties::ForceDisableResvFluidInPlaceOutput>();
181 if (! Parameters::get<TypeTag, Properties::OwnerCellsFirst>()) {
182 const std::string msg =
"The output code does not support --owner-cells-first=false.";
183 if (collectToIORank.isIORank()) {
186 OPM_THROW_NOLOG(std::runtime_error, msg);
189 if (smryCfg.match(
"[FB]PP[OGW]") || smryCfg.match(
"RPP[OGW]*")) {
190 auto rset = this->eclState_.fieldProps().fip_regions();
191 rset.push_back(
"PVTNUM");
196 this->regionAvgDensity_
197 .emplace(this->simulator_.gridView().comm(),
198 FluidSystem::numPhases, rset,
199 [fp = std::cref(this->eclState_.fieldProps())]
200 (
const std::string& rsetName) ->
decltype(
auto)
201 { return fp.get().get_int(rsetName); });
210 Parameters::registerParam<TypeTag, Properties::ForceDisableFluidInPlaceOutput>
211 (
"Do not print fluid-in-place values after each report step "
212 "even if requested by the deck.");
213 Parameters::registerParam<TypeTag, Properties::ForceDisableResvFluidInPlaceOutput>
214 (
"Do not print reservoir volumes values after each report step "
215 "even if requested by the deck.");
224 const unsigned reportStepNum,
227 const bool isRestart)
233 const auto& problem = this->simulator_.problem();
235 this->doAllocBuffers(bufferSize,
240 problem.vapparsActive(std::max(simulator_.episodeIndex(), 0)),
241 problem.materialLawManager()->enableHysteresis(),
242 problem.tracerModel().numTracers(),
243 problem.eclWriter()->getOutputNnc().size());
246 void processElementMech(
const ElementContext& elemCtx)
248 if constexpr (getPropValue<TypeTag, Properties::EnableMech>()) {
249 const auto& problem = elemCtx.simulator().problem();
250 const auto& model = problem.geoMechModel();
251 for (
unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(0); ++dofIdx) {
252 unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, 0);
253 if (!this->mechPotentialForce_.empty()) {
255 this->mechPotentialForce_[globalDofIdx] = model.mechPotentialForce(globalDofIdx);
256 this->mechPotentialPressForce_[globalDofIdx] = model.mechPotentialPressForce(globalDofIdx);
257 this->mechPotentialTempForce_[globalDofIdx] = model.mechPotentialTempForce(globalDofIdx);
259 this->dispX_[globalDofIdx] = model.disp(globalDofIdx, 0);
260 this->dispY_[globalDofIdx] = model.disp(globalDofIdx, 1);
261 this->dispZ_[globalDofIdx] = model.disp(globalDofIdx, 2);
262 this->stressXX_[globalDofIdx] = model.stress(globalDofIdx, 0);
263 this->stressYY_[globalDofIdx] = model.stress(globalDofIdx, 1);
264 this->stressZZ_[globalDofIdx] = model.stress(globalDofIdx, 2);
266 this->stressXY_[globalDofIdx] = model.stress(globalDofIdx, 5);
267 this->stressXZ_[globalDofIdx] = model.stress(globalDofIdx, 4);
268 this->stressYZ_[globalDofIdx] = model.stress(globalDofIdx, 3);
270 this->strainXX_[globalDofIdx] = model.strain(globalDofIdx, 0);
271 this->strainYY_[globalDofIdx] = model.strain(globalDofIdx, 1);
272 this->strainZZ_[globalDofIdx] = model.strain(globalDofIdx, 2);
274 this->strainXY_[globalDofIdx] = model.strain(globalDofIdx, 5);
275 this->strainXZ_[globalDofIdx] = model.strain(globalDofIdx, 4);
276 this->strainYZ_[globalDofIdx] = model.strain(globalDofIdx, 3);
279 this->delstressXX_[globalDofIdx] = model.delstress(globalDofIdx, 0);
280 this->delstressYY_[globalDofIdx] = model.delstress(globalDofIdx, 1);
281 this->delstressZZ_[globalDofIdx] = model.delstress(globalDofIdx, 2);
283 this->delstressXY_[globalDofIdx] = model.delstress(globalDofIdx, 5);
284 this->delstressXZ_[globalDofIdx] = model.delstress(globalDofIdx, 4);
285 this->delstressYZ_[globalDofIdx] = model.delstress(globalDofIdx, 3);
301 const auto& problem = elemCtx.simulator().problem();
302 const auto& modelResid = elemCtx.simulator().model().linearizer().residual();
303 for (
unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(0); ++dofIdx) {
304 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, 0);
305 const auto& fs = intQuants.fluidState();
307 using FluidState = std::remove_cv_t<std::remove_reference_t<
decltype(fs)>>;
309 const unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, 0);
310 const unsigned pvtRegionIdx = elemCtx.primaryVars(dofIdx, 0).pvtRegionIndex();
312 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
313 if (this->saturation_[phaseIdx].empty())
316 this->saturation_[phaseIdx][globalDofIdx] = getValue(fs.saturation(phaseIdx));
317 Valgrind::CheckDefined(this->saturation_[phaseIdx][globalDofIdx]);
320 if (this->regionAvgDensity_.has_value()) {
323 const auto porv = intQuants.referencePorosity()
324 * elemCtx.simulator().model().dofTotalVolume(globalDofIdx);
326 this->aggregateAverageDensityContributions_(fs, globalDofIdx,
327 static_cast<double>(porv));
330 if (!this->fluidPressure_.empty()) {
331 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
333 this->fluidPressure_[globalDofIdx] = getValue(fs.pressure(oilPhaseIdx));
334 }
else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
336 this->fluidPressure_[globalDofIdx] = getValue(fs.pressure(gasPhaseIdx));
339 this->fluidPressure_[globalDofIdx] = getValue(fs.pressure(waterPhaseIdx));
341 Valgrind::CheckDefined(this->fluidPressure_[globalDofIdx]);
344 if (!this->temperature_.empty()) {
345 this->temperature_[globalDofIdx] = getValue(fs.temperature(oilPhaseIdx));
346 Valgrind::CheckDefined(this->temperature_[globalDofIdx]);
348 if (!this->gasDissolutionFactor_.empty()) {
349 Scalar SoMax = elemCtx.problem().maxOilSaturation(globalDofIdx);
350 this->gasDissolutionFactor_[globalDofIdx]
351 = FluidSystem::template saturatedDissolutionFactor<FluidState, Scalar>(
352 fs, oilPhaseIdx, pvtRegionIdx, SoMax);
353 Valgrind::CheckDefined(this->gasDissolutionFactor_[globalDofIdx]);
355 if (!this->oilVaporizationFactor_.empty()) {
356 Scalar SoMax = elemCtx.problem().maxOilSaturation(globalDofIdx);
357 this->oilVaporizationFactor_[globalDofIdx]
358 = FluidSystem::template saturatedDissolutionFactor<FluidState, Scalar>(
359 fs, gasPhaseIdx, pvtRegionIdx, SoMax);
360 Valgrind::CheckDefined(this->oilVaporizationFactor_[globalDofIdx]);
362 if (!this->gasFormationVolumeFactor_.empty()) {
363 this->gasFormationVolumeFactor_[globalDofIdx] = 1.0
364 / FluidSystem::template inverseFormationVolumeFactor<FluidState, Scalar>(
365 fs, gasPhaseIdx, pvtRegionIdx);
366 Valgrind::CheckDefined(this->gasFormationVolumeFactor_[globalDofIdx]);
368 if (!this->saturatedOilFormationVolumeFactor_.empty()) {
369 this->saturatedOilFormationVolumeFactor_[globalDofIdx] = 1.0
370 / FluidSystem::template saturatedInverseFormationVolumeFactor<FluidState, Scalar>(
371 fs, oilPhaseIdx, pvtRegionIdx);
372 Valgrind::CheckDefined(this->saturatedOilFormationVolumeFactor_[globalDofIdx]);
374 if (!this->oilSaturationPressure_.empty()) {
375 this->oilSaturationPressure_[globalDofIdx]
376 = FluidSystem::template saturationPressure<FluidState, Scalar>(fs, oilPhaseIdx, pvtRegionIdx);
377 Valgrind::CheckDefined(this->oilSaturationPressure_[globalDofIdx]);
380 if (!this->rs_.empty()) {
381 this->rs_[globalDofIdx] = getValue(fs.Rs());
382 Valgrind::CheckDefined(this->rs_[globalDofIdx]);
384 if (!this->rsw_.empty()) {
385 this->rsw_[globalDofIdx] = getValue(fs.Rsw());
386 Valgrind::CheckDefined(this->rsw_[globalDofIdx]);
389 if (!this->rv_.empty()) {
390 this->rv_[globalDofIdx] = getValue(fs.Rv());
391 Valgrind::CheckDefined(this->rv_[globalDofIdx]);
393 if (!this->pcgw_.empty()) {
394 this->pcgw_[globalDofIdx] = getValue(fs.pressure(gasPhaseIdx)) - getValue(fs.pressure(waterPhaseIdx));
395 Valgrind::CheckDefined(this->pcgw_[globalDofIdx]);
397 if (!this->pcow_.empty()) {
398 this->pcow_[globalDofIdx] = getValue(fs.pressure(oilPhaseIdx)) - getValue(fs.pressure(waterPhaseIdx));
399 Valgrind::CheckDefined(this->pcow_[globalDofIdx]);
401 if (!this->pcog_.empty()) {
402 this->pcog_[globalDofIdx] = getValue(fs.pressure(gasPhaseIdx)) - getValue(fs.pressure(oilPhaseIdx));
403 Valgrind::CheckDefined(this->pcog_[globalDofIdx]);
406 if (!this->rvw_.empty()) {
407 this->rvw_[globalDofIdx] = getValue(fs.Rvw());
408 Valgrind::CheckDefined(this->rvw_[globalDofIdx]);
411 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
412 if (this->invB_[phaseIdx].empty())
415 this->invB_[phaseIdx][globalDofIdx] = getValue(fs.invB(phaseIdx));
416 Valgrind::CheckDefined(this->invB_[phaseIdx][globalDofIdx]);
419 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
420 if (this->density_[phaseIdx].empty())
423 this->density_[phaseIdx][globalDofIdx] = getValue(fs.density(phaseIdx));
424 Valgrind::CheckDefined(this->density_[phaseIdx][globalDofIdx]);
427 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
428 if (this->viscosity_[phaseIdx].empty())
431 if (!this->extboX_.empty() && phaseIdx == oilPhaseIdx)
432 this->viscosity_[phaseIdx][globalDofIdx] = getValue(intQuants.oilViscosity());
433 else if (!this->extboX_.empty() && phaseIdx == gasPhaseIdx)
434 this->viscosity_[phaseIdx][globalDofIdx] = getValue(intQuants.gasViscosity());
436 this->viscosity_[phaseIdx][globalDofIdx] = getValue(fs.viscosity(phaseIdx));
437 Valgrind::CheckDefined(this->viscosity_[phaseIdx][globalDofIdx]);
440 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
441 if (this->relativePermeability_[phaseIdx].empty())
444 this->relativePermeability_[phaseIdx][globalDofIdx]
445 = getValue(intQuants.relativePermeability(phaseIdx));
446 Valgrind::CheckDefined(this->relativePermeability_[phaseIdx][globalDofIdx]);
449 if (!this->drsdtcon_.empty()) {
450 this->drsdtcon_[globalDofIdx] = problem.drsdtcon(globalDofIdx, elemCtx.simulator().episodeIndex());
453 if (!this->sSol_.empty()) {
454 this->sSol_[globalDofIdx] = intQuants.solventSaturation().value();
457 if (!this->rswSol_.empty()) {
458 this->rswSol_[globalDofIdx] = intQuants.rsSolw().value();
461 if (!this->cPolymer_.empty()) {
462 this->cPolymer_[globalDofIdx] = intQuants.polymerConcentration().value();
465 if (!this->cFoam_.empty()) {
466 this->cFoam_[globalDofIdx] = intQuants.foamConcentration().value();
469 if (!this->cSalt_.empty()) {
470 this->cSalt_[globalDofIdx] = fs.saltConcentration().value();
473 if (!this->pSalt_.empty()) {
474 this->pSalt_[globalDofIdx] = intQuants.saltSaturation().value();
477 if (!this->permFact_.empty()) {
478 this->permFact_[globalDofIdx] = intQuants.permFactor().value();
481 if (!this->extboX_.empty()) {
482 this->extboX_[globalDofIdx] = intQuants.xVolume().value();
485 if (!this->extboY_.empty()) {
486 this->extboY_[globalDofIdx] = intQuants.yVolume().value();
489 if (!this->extboZ_.empty()) {
490 this->extboZ_[globalDofIdx] = intQuants.zFraction().value();
493 if (!this->rPorV_.empty()) {
494 const auto totVolume = elemCtx.simulator().model().dofTotalVolume(globalDofIdx);
495 this->rPorV_[globalDofIdx] = totVolume * intQuants.porosity().value();
498 if (!this->mFracCo2_.empty()) {
499 const Scalar stdVolOil = getValue(fs.saturation(oilPhaseIdx)) * getValue(fs.invB(oilPhaseIdx))
500 + getValue(fs.saturation(gasPhaseIdx)) * getValue(fs.invB(gasPhaseIdx)) * getValue(fs.Rv());
501 const Scalar stdVolGas = getValue(fs.saturation(gasPhaseIdx)) * getValue(fs.invB(gasPhaseIdx))
502 * (1.0 - intQuants.yVolume().value())
503 + getValue(fs.saturation(oilPhaseIdx)) * getValue(fs.invB(oilPhaseIdx)) * getValue(fs.Rs())
504 * (1.0 - intQuants.xVolume().value());
505 const Scalar stdVolCo2 = getValue(fs.saturation(gasPhaseIdx)) * getValue(fs.invB(gasPhaseIdx))
506 * intQuants.yVolume().value()
507 + getValue(fs.saturation(oilPhaseIdx)) * getValue(fs.invB(oilPhaseIdx)) * getValue(fs.Rs())
508 * intQuants.xVolume().value();
509 const Scalar rhoO = FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
510 const Scalar rhoG = FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
511 const Scalar rhoCO2 = intQuants.zRefDensity();
512 const Scalar stdMassTotal = 1.0e-10 + stdVolOil * rhoO + stdVolGas * rhoG + stdVolCo2 * rhoCO2;
513 this->mFracOil_[globalDofIdx] = stdVolOil * rhoO / stdMassTotal;
514 this->mFracGas_[globalDofIdx] = stdVolGas * rhoG / stdMassTotal;
515 this->mFracCo2_[globalDofIdx] = stdVolCo2 * rhoCO2 / stdMassTotal;
518 if (!this->cMicrobes_.empty()) {
519 this->cMicrobes_[globalDofIdx] = intQuants.microbialConcentration().value();
522 if (!this->cOxygen_.empty()) {
523 this->cOxygen_[globalDofIdx] = intQuants.oxygenConcentration().value();
526 if (!this->cUrea_.empty()) {
527 this->cUrea_[globalDofIdx] = 10
528 * intQuants.ureaConcentration()
532 if (!this->cBiofilm_.empty()) {
533 this->cBiofilm_[globalDofIdx] = intQuants.biofilmConcentration().value();
536 if (!this->cCalcite_.empty()) {
537 this->cCalcite_[globalDofIdx] = intQuants.calciteConcentration().value();
540 if (!this->bubblePointPressure_.empty()) {
542 this->bubblePointPressure_[globalDofIdx]
543 = getValue(FluidSystem::bubblePointPressure(fs, intQuants.pvtRegionIndex()));
544 }
catch (
const NumericalProblem&) {
545 const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx);
546 this->failedCellsPb_.push_back(cartesianIdx);
550 if (!this->dewPointPressure_.empty()) {
552 this->dewPointPressure_[globalDofIdx]
553 = getValue(FluidSystem::dewPointPressure(fs, intQuants.pvtRegionIndex()));
554 }
catch (
const NumericalProblem&) {
555 const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx);
556 this->failedCellsPd_.push_back(cartesianIdx);
560 if (!this->soMax_.empty())
561 this->soMax_[globalDofIdx]
562 = std::max(getValue(fs.saturation(oilPhaseIdx)), problem.maxOilSaturation(globalDofIdx));
564 if (!this->swMax_.empty())
565 this->swMax_[globalDofIdx]
566 = std::max(getValue(fs.saturation(waterPhaseIdx)), problem.maxWaterSaturation(globalDofIdx));
568 if (!this->minimumOilPressure_.empty())
569 this->minimumOilPressure_[globalDofIdx]
570 = std::min(getValue(fs.pressure(oilPhaseIdx)), problem.minOilPressure(globalDofIdx));
572 if (!this->overburdenPressure_.empty())
573 this->overburdenPressure_[globalDofIdx] = problem.overburdenPressure(globalDofIdx);
575 if (!this->rockCompPorvMultiplier_.empty())
576 this->rockCompPorvMultiplier_[globalDofIdx]
577 = problem.template rockCompPoroMultiplier<Scalar>(intQuants, globalDofIdx);
579 if (!this->rockCompTransMultiplier_.empty())
580 this->rockCompTransMultiplier_[globalDofIdx]
581 = problem.template rockCompTransMultiplier<Scalar>(intQuants, globalDofIdx);
583 const auto& matLawManager = problem.materialLawManager();
584 if (matLawManager->enableHysteresis()) {
585 if (!this->pcSwMdcOw_.empty() && !this->krnSwMdcOw_.empty()) {
586 matLawManager->oilWaterHysteresisParams(
587 this->pcSwMdcOw_[globalDofIdx], this->krnSwMdcOw_[globalDofIdx], globalDofIdx);
589 if (!this->pcSwMdcGo_.empty() && !this->krnSwMdcGo_.empty()) {
590 matLawManager->gasOilHysteresisParams(
591 this->pcSwMdcGo_[globalDofIdx], this->krnSwMdcGo_[globalDofIdx], globalDofIdx);
595 if (!this->ppcw_.empty()) {
596 this->ppcw_[globalDofIdx] = matLawManager->oilWaterScaledEpsInfoDrainage(globalDofIdx).maxPcow;
605 if ((elemCtx.simulator().episodeIndex() < 0) &&
606 FluidSystem::phaseIsActive(oilPhaseIdx) &&
607 FluidSystem::phaseIsActive(gasPhaseIdx))
609 const auto& fsInitial = problem.initialFluidState(globalDofIdx);
612 if (!this->rv_.empty())
613 this->rv_[globalDofIdx] = fsInitial.Rv();
615 if (!this->rs_.empty())
616 this->rs_[globalDofIdx] = fsInitial.Rs();
618 if (!this->rsw_.empty())
619 this->rsw_[globalDofIdx] = fsInitial.Rsw();
621 if (!this->rvw_.empty())
622 this->rvw_[globalDofIdx] = fsInitial.Rvw();
625 if (!this->density_[oilPhaseIdx].empty())
626 this->density_[oilPhaseIdx][globalDofIdx]
627 = FluidSystem::density(fsInitial, oilPhaseIdx, intQuants.pvtRegionIndex());
629 if (!this->density_[gasPhaseIdx].empty())
630 this->density_[gasPhaseIdx][globalDofIdx]
631 = FluidSystem::density(fsInitial, gasPhaseIdx, intQuants.pvtRegionIndex());
633 if (!this->invB_[oilPhaseIdx].empty())
634 this->invB_[oilPhaseIdx][globalDofIdx]
635 = FluidSystem::inverseFormationVolumeFactor(fsInitial, oilPhaseIdx, intQuants.pvtRegionIndex());
637 if (!this->invB_[gasPhaseIdx].empty())
638 this->invB_[gasPhaseIdx][globalDofIdx]
639 = FluidSystem::inverseFormationVolumeFactor(fsInitial, gasPhaseIdx, intQuants.pvtRegionIndex());
641 if (!this->viscosity_[oilPhaseIdx].empty())
642 this->viscosity_[oilPhaseIdx][globalDofIdx]
643 = FluidSystem::viscosity(fsInitial, oilPhaseIdx, intQuants.pvtRegionIndex());
645 if (!this->viscosity_[gasPhaseIdx].empty())
646 this->viscosity_[gasPhaseIdx][globalDofIdx]
647 = FluidSystem::viscosity(fsInitial, gasPhaseIdx, intQuants.pvtRegionIndex());
651 const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx);
652 if (this->oilConnectionPressures_.count(cartesianIdx) > 0) {
653 this->oilConnectionPressures_[cartesianIdx] = getValue(fs.pressure(oilPhaseIdx));
655 if (this->waterConnectionSaturations_.count(cartesianIdx) > 0) {
656 this->waterConnectionSaturations_[cartesianIdx] = getValue(fs.saturation(waterPhaseIdx));
658 if (this->gasConnectionSaturations_.count(cartesianIdx) > 0) {
659 this->gasConnectionSaturations_[cartesianIdx] = getValue(fs.saturation(gasPhaseIdx));
663 const auto& tracerModel = simulator_.problem().tracerModel();
664 if (! this->tracerConcentrations_.empty()) {
665 for (
int tracerIdx = 0; tracerIdx < tracerModel.numTracers(); ++tracerIdx) {
666 if (this->tracerConcentrations_[tracerIdx].empty()) {
670 this->tracerConcentrations_[tracerIdx][globalDofIdx] =
671 tracerModel.tracerConcentration(tracerIdx, globalDofIdx);
676 for (
int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx )
678 if (!this->residual_[phaseIdx].empty()) {
679 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
680 this->residual_[phaseIdx][globalDofIdx] = modelResid[globalDofIdx][activeCompIdx];
686 void processElementFlows(
const ElementContext& elemCtx)
688 OPM_TIMEBLOCK_LOCAL(processElementBlockData);
692 const auto& problem = elemCtx.simulator().problem();
693 for (
unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(0); ++dofIdx) {
695 unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, 0);
696 if (!problem.model().linearizer().getFlowsInfo().empty()) {
697 const auto& flowsInf = problem.model().linearizer().getFlowsInfo();
698 auto flowsInfos = flowsInf[globalDofIdx];
699 for (
auto& flowsInfo : flowsInfos) {
700 if (flowsInfo.faceId >= 0) {
701 if (!this->flows_[flowsInfo.faceId][gasCompIdx].empty()) {
702 this->flows_[flowsInfo.faceId][gasCompIdx][globalDofIdx]
703 = flowsInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(gasCompIdx)];
705 if (!this->flows_[flowsInfo.faceId][oilCompIdx].empty()) {
706 this->flows_[flowsInfo.faceId][oilCompIdx][globalDofIdx]
707 = flowsInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(oilCompIdx)];
709 if (!this->flows_[flowsInfo.faceId][waterCompIdx].empty()) {
710 this->flows_[flowsInfo.faceId][waterCompIdx][globalDofIdx]
711 = flowsInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(waterCompIdx)];
714 if (flowsInfo.faceId == -2) {
715 if (!this->flowsn_[gasCompIdx].indices.empty()) {
716 this->flowsn_[gasCompIdx].indices[flowsInfo.nncId] = flowsInfo.nncId;
717 this->flowsn_[gasCompIdx].values[flowsInfo.nncId]
718 = flowsInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(gasCompIdx)];
720 if (!this->flowsn_[oilCompIdx].indices.empty()) {
721 this->flowsn_[oilCompIdx].indices[flowsInfo.nncId] = flowsInfo.nncId;
722 this->flowsn_[oilCompIdx].values[flowsInfo.nncId]
723 = flowsInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(oilCompIdx)];
725 if (!this->flowsn_[waterCompIdx].indices.empty()) {
726 this->flowsn_[waterCompIdx].indices[flowsInfo.nncId] = flowsInfo.nncId;
727 this->flowsn_[waterCompIdx].values[flowsInfo.nncId]
728 = flowsInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(waterCompIdx)];
735 if (!problem.model().linearizer().getFloresInfo().empty()) {
736 const auto& floresInf = problem.model().linearizer().getFloresInfo();
737 auto floresInfos =floresInf[globalDofIdx];
738 for (
auto& floresInfo : floresInfos) {
739 if (floresInfo.faceId >= 0) {
740 if (!this->flores_[floresInfo.faceId][gasCompIdx].empty()) {
741 this->flores_[floresInfo.faceId][gasCompIdx][globalDofIdx]
742 = floresInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(gasCompIdx)];
744 if (!this->flores_[floresInfo.faceId][oilCompIdx].empty()) {
745 this->flores_[floresInfo.faceId][oilCompIdx][globalDofIdx]
746 = floresInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(oilCompIdx)];
748 if (!this->flores_[floresInfo.faceId][waterCompIdx].empty()) {
749 this->flores_[floresInfo.faceId][waterCompIdx][globalDofIdx]
750 = floresInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(waterCompIdx)];
754 if (floresInfo.faceId == -2) {
755 if (!this->floresn_[gasCompIdx].indices.empty()) {
756 this->floresn_[gasCompIdx].indices[floresInfo.nncId] = floresInfo.nncId;
757 this->floresn_[gasCompIdx].values[floresInfo.nncId]
758 = floresInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(gasCompIdx)];
760 if (!this->floresn_[oilCompIdx].indices.empty()) {
761 this->floresn_[oilCompIdx].indices[floresInfo.nncId] = floresInfo.nncId;
762 this->floresn_[oilCompIdx].values[floresInfo.nncId]
763 = floresInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(oilCompIdx)];
765 if (!this->floresn_[waterCompIdx].indices.empty()) {
766 this->floresn_[waterCompIdx].indices[floresInfo.nncId] = floresInfo.nncId;
767 this->floresn_[waterCompIdx].values[floresInfo.nncId]
768 = floresInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(waterCompIdx)];
776 void processElementBlockData(
const ElementContext& elemCtx)
778 OPM_TIMEBLOCK_LOCAL(processElementBlockData);
779 if (!std::is_same<Discretization, EcfvDiscretization<TypeTag>>::value)
782 const auto& problem = elemCtx.simulator().problem();
783 for (
unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(0); ++dofIdx) {
785 const auto globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, 0);
786 const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx);
787 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, 0);
788 const auto& fs = intQuants.fluidState();
789 for (
auto& val : this->blockData_) {
790 const auto& key = val.first;
791 assert(key.second > 0);
793 const auto cartesianIdxBlock =
static_cast<std::remove_cv_t<
794 std::remove_reference_t<decltype(cartesianIdx)
>>>(key.second - 1);
796 if (cartesianIdx == cartesianIdxBlock) {
797 if ((key.first ==
"BWSAT") || (key.first ==
"BSWAT"))
798 val.second = getValue(fs.saturation(waterPhaseIdx));
799 else if ((key.first ==
"BGSAT") || (key.first ==
"BSGAS"))
800 val.second = getValue(fs.saturation(gasPhaseIdx));
801 else if ((key.first ==
"BOSAT") || (key.first ==
"BSOIL"))
802 val.second = getValue(fs.saturation(oilPhaseIdx));
803 else if (key.first ==
"BNSAT")
804 val.second = intQuants.solventSaturation().value();
805 else if ((key.first ==
"BPR") || (key.first ==
"BPRESSUR")) {
806 if (FluidSystem::phaseIsActive(oilPhaseIdx))
807 val.second = getValue(fs.pressure(oilPhaseIdx));
808 else if (FluidSystem::phaseIsActive(gasPhaseIdx))
809 val.second = getValue(fs.pressure(gasPhaseIdx));
810 else if (FluidSystem::phaseIsActive(waterPhaseIdx))
811 val.second = getValue(fs.pressure(waterPhaseIdx));
813 else if ((key.first ==
"BTCNFHEA") || (key.first ==
"BTEMP")) {
814 if (FluidSystem::phaseIsActive(oilPhaseIdx))
815 val.second = getValue(fs.temperature(oilPhaseIdx));
816 else if (FluidSystem::phaseIsActive(gasPhaseIdx))
817 val.second = getValue(fs.temperature(gasPhaseIdx));
818 else if (FluidSystem::phaseIsActive(waterPhaseIdx))
819 val.second = getValue(fs.temperature(waterPhaseIdx));
821 else if (key.first ==
"BWKR" || key.first ==
"BKRW")
822 val.second = getValue(intQuants.relativePermeability(waterPhaseIdx));
823 else if (key.first ==
"BGKR" || key.first ==
"BKRG")
824 val.second = getValue(intQuants.relativePermeability(gasPhaseIdx));
825 else if (key.first ==
"BOKR" || key.first ==
"BKRO")
826 val.second = getValue(intQuants.relativePermeability(oilPhaseIdx));
827 else if (key.first ==
"BKROG") {
828 const auto& materialParams = problem.materialLawParams(elemCtx, dofIdx, 0);
830 = MaterialLaw::template relpermOilInOilGasSystem<Evaluation>(materialParams, fs);
831 val.second = getValue(krog);
833 else if (key.first ==
"BKROW") {
834 const auto& materialParams = problem.materialLawParams(elemCtx, dofIdx, 0);
836 = MaterialLaw::template relpermOilInOilWaterSystem<Evaluation>(materialParams, fs);
837 val.second = getValue(krow);
839 else if (key.first ==
"BWPC")
840 val.second = getValue(fs.pressure(oilPhaseIdx)) - getValue(fs.pressure(waterPhaseIdx));
841 else if (key.first ==
"BGPC")
842 val.second = getValue(fs.pressure(gasPhaseIdx)) - getValue(fs.pressure(oilPhaseIdx));
843 else if (key.first ==
"BWPR")
844 val.second = getValue(fs.pressure(waterPhaseIdx));
845 else if (key.first ==
"BGPR")
846 val.second = getValue(fs.pressure(gasPhaseIdx));
847 else if (key.first ==
"BVWAT" || key.first ==
"BWVIS")
848 val.second = getValue(fs.viscosity(waterPhaseIdx));
849 else if (key.first ==
"BVGAS" || key.first ==
"BGVIS")
850 val.second = getValue(fs.viscosity(gasPhaseIdx));
851 else if (key.first ==
"BVOIL" || key.first ==
"BOVIS")
852 val.second = getValue(fs.viscosity(oilPhaseIdx));
853 else if ((key.first ==
"BODEN") || (key.first ==
"BDENO"))
854 val.second = getValue(fs.density(oilPhaseIdx));
855 else if ((key.first ==
"BGDEN") || (key.first ==
"BDENG"))
856 val.second = getValue(fs.density(gasPhaseIdx));
857 else if ((key.first ==
"BWDEN") || (key.first ==
"BDENW"))
858 val.second = getValue(fs.density(waterPhaseIdx));
859 else if ((key.first ==
"BRPV") ||
860 (key.first ==
"BOPV") ||
861 (key.first ==
"BWPV") ||
862 (key.first ==
"BGPV"))
864 if (key.first ==
"BRPV") {
867 else if (key.first ==
"BOPV") {
868 val.second = getValue(fs.saturation(oilPhaseIdx));
870 else if (key.first ==
"BWPV") {
871 val.second = getValue(fs.saturation(waterPhaseIdx));
874 val.second = getValue(fs.saturation(gasPhaseIdx));
878 val.second *= getValue(intQuants.porosity())
879 * elemCtx.simulator().model().dofTotalVolume(globalDofIdx);
881 else if (key.first ==
"BRS")
882 val.second = getValue(fs.Rs());
883 else if (key.first ==
"BRV")
884 val.second = getValue(fs.Rv());
885 else if ((key.first ==
"BOIP") || (key.first ==
"BOIPL") || (key.first ==
"BOIPG") ||
886 (key.first ==
"BGIP") || (key.first ==
"BGIPL") || (key.first ==
"BGIPG") ||
887 (key.first ==
"BWIP"))
889 if ((key.first ==
"BOIP") || (key.first ==
"BOIPL")) {
890 val.second = getValue(fs.invB(oilPhaseIdx)) * getValue(fs.saturation(oilPhaseIdx));
892 if (key.first ==
"BOIP") {
893 val.second += getValue(fs.Rv()) * getValue(fs.invB(gasPhaseIdx))
894 * getValue(fs.saturation(gasPhaseIdx));
897 else if (key.first ==
"BOIPG") {
898 val.second = getValue(fs.Rv()) * getValue(fs.invB(gasPhaseIdx))
899 * getValue(fs.saturation(gasPhaseIdx));
901 else if ((key.first ==
"BGIP") || (key.first ==
"BGIPG")) {
902 val.second = getValue(fs.invB(gasPhaseIdx)) * getValue(fs.saturation(gasPhaseIdx));
904 if (key.first ==
"BGIP") {
905 if (!FluidSystem::phaseIsActive(oilPhaseIdx)) {
906 val.second += getValue(fs.Rsw()) * getValue(fs.invB(waterPhaseIdx))
907 * getValue(fs.saturation(waterPhaseIdx));
910 val.second += getValue(fs.Rs()) * getValue(fs.invB(oilPhaseIdx))
911 * getValue(fs.saturation(oilPhaseIdx));
915 else if (key.first ==
"BGIPL") {
916 if (!FluidSystem::phaseIsActive(oilPhaseIdx)) {
917 val.second = getValue(fs.Rsw()) * getValue(fs.invB(waterPhaseIdx))
918 * getValue(fs.saturation(waterPhaseIdx));
921 val.second = getValue(fs.Rs()) * getValue(fs.invB(oilPhaseIdx))
922 * getValue(fs.saturation(oilPhaseIdx));
926 val.second = getValue(fs.invB(waterPhaseIdx)) * getValue(fs.saturation(waterPhaseIdx));
930 val.second *= elemCtx.simulator().model().dofTotalVolume(globalDofIdx)
931 * getValue(intQuants.porosity());
933 else if ((key.first ==
"BPPO") ||
934 (key.first ==
"BPPG") ||
935 (key.first ==
"BPPW"))
937 auto phase = RegionPhasePoreVolAverage::Phase{};
939 if (key.first ==
"BPPO") {
940 phase.ix = oilPhaseIdx;
942 else if (key.first ==
"BPPG") {
943 phase.ix = gasPhaseIdx;
946 phase.ix = waterPhaseIdx;
956 const auto datum = this->eclState_.getSimulationConfig()
957 .datumDepths()(this->regions_[
"FIPNUM"][dofIdx] - 1);
960 const auto region = RegionPhasePoreVolAverage::Region {
961 elemCtx.primaryVars(dofIdx, 0).pvtRegionIndex() + 1
964 const auto density = this->regionAvgDensity_
965 ->value(
"PVTNUM", phase, region);
967 const auto press = getValue(fs.pressure(phase.ix));
969 elemCtx.problem().gravity()[GridView::dimensionworld - 1];
970 const auto dz = problem.dofCenterDepth(globalDofIdx) - datum;
972 val.second = press - density*dz*grav;
974 else if ((key.first ==
"BFLOWI") ||
975 (key.first ==
"BLFOWJ") ||
976 (key.first ==
"BFLOWK"))
978 auto dir = FaceDir::ToIntersectionIndex(Dir::XPlus);
980 if (key.first ==
"BFLOWJ") {
981 dir = FaceDir::ToIntersectionIndex(Dir::YPlus);
983 else if (key.first ==
"BFLOWK") {
984 dir = FaceDir::ToIntersectionIndex(Dir::ZPlus);
987 val.second = this->flows_[dir][waterCompIdx][globalDofIdx];
990 std::string logstring =
"Keyword '";
991 logstring.append(key.first);
992 logstring.append(
"' is unhandled for output to summary file.");
993 OpmLog::warning(
"Unhandled output keyword", logstring);
1028 template <
class ActiveIndex,
class CartesianIndex>
1030 ActiveIndex&& activeIndex,
1031 CartesianIndex&& cartesianIndex)
1034 const auto identifyCell = [&activeIndex, &cartesianIndex](
const Element& elem)
1037 const auto cellIndex = activeIndex(elem);
1040 static_cast<int>(cellIndex),
1041 cartesianIndex(cellIndex),
1042 elem.partitionType() == Dune::InteriorEntity
1046 const auto timeIdx = 0u;
1047 const auto& stencil = elemCtx.stencil(timeIdx);
1048 const auto numInteriorFaces = elemCtx.numInteriorFaces(timeIdx);
1050 for (
auto scvfIdx = 0 * numInteriorFaces; scvfIdx < numInteriorFaces; ++scvfIdx) {
1051 const auto& face = stencil.interiorFace(scvfIdx);
1052 const auto left = identifyCell(stencil.element(face.interiorIndex()));
1053 const auto right = identifyCell(stencil.element(face.exteriorIndex()));
1055 const auto rates = this->
1056 getComponentSurfaceRates(elemCtx, face.area(), scvfIdx, timeIdx);
1070 this->interRegionFlows_.
clear();
1078 this->interRegionFlows_.
compress();
1086 return this->interRegionFlows_;
1089 template <
class Flu
idState>
1090 void assignToFluidState(FluidState& fs,
unsigned elemIdx)
const
1092 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1093 if (this->saturation_[phaseIdx].empty())
1096 fs.setSaturation(phaseIdx, this->saturation_[phaseIdx][elemIdx]);
1099 if (!this->fluidPressure_.empty()) {
1102 std::array<Scalar, numPhases> pc = {0};
1103 const MaterialLawParams& matParams = simulator_.problem().materialLawParams(elemIdx);
1104 MaterialLaw::capillaryPressures(pc, matParams, fs);
1105 Valgrind::CheckDefined(this->fluidPressure_[elemIdx]);
1106 Valgrind::CheckDefined(pc);
1107 assert(FluidSystem::phaseIsActive(oilPhaseIdx));
1108 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1109 if (!FluidSystem::phaseIsActive(phaseIdx))
1112 fs.setPressure(phaseIdx, this->fluidPressure_[elemIdx] + (pc[phaseIdx] - pc[oilPhaseIdx]));
1116 if (!this->temperature_.empty())
1117 fs.setTemperature(this->temperature_[elemIdx]);
1118 if (!this->rs_.empty())
1119 fs.setRs(this->rs_[elemIdx]);
1120 if (!this->rsw_.empty())
1121 fs.setRsw(this->rsw_[elemIdx]);
1122 if (!this->rv_.empty())
1123 fs.setRv(this->rv_[elemIdx]);
1124 if (!this->rvw_.empty())
1125 fs.setRvw(this->rvw_[elemIdx]);
1128 void initHysteresisParams(Simulator& simulator,
unsigned elemIdx)
const
1130 if (!this->soMax_.empty())
1131 simulator.problem().setMaxOilSaturation(elemIdx, this->soMax_[elemIdx]);
1133 if (simulator.problem().materialLawManager()->enableHysteresis()) {
1134 auto matLawManager = simulator.problem().materialLawManager();
1136 if (!this->pcSwMdcOw_.empty() && !this->krnSwMdcOw_.empty()) {
1137 matLawManager->setOilWaterHysteresisParams(
1138 this->pcSwMdcOw_[elemIdx], this->krnSwMdcOw_[elemIdx], elemIdx);
1140 if (!this->pcSwMdcGo_.empty() && !this->krnSwMdcGo_.empty()) {
1141 matLawManager->setGasOilHysteresisParams(
1142 this->pcSwMdcGo_[elemIdx], this->krnSwMdcGo_[elemIdx], elemIdx);
1146 if (simulator_.vanguard().eclState().fieldProps().has_double(
"SWATINIT")) {
1147 simulator.problem().materialLawManager()
1148 ->applyRestartSwatInit(elemIdx, this->ppcw_[elemIdx]);
1152 void updateFluidInPlace(
const ElementContext& elemCtx)
1154 for (
unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(0); ++dofIdx) {
1155 updateFluidInPlace_(elemCtx, dofIdx);
1159 void updateFluidInPlace(
const unsigned globalDofIdx,
1160 const IntensiveQuantities& intQuants,
1161 const double totVolume)
1163 this->updateFluidInPlace_(globalDofIdx, intQuants, totVolume);
1167 bool isDefunctParallelWell(std::string wname)
const override
1169 if (simulator_.gridView().comm().size() == 1)
1171 const auto& parallelWells = simulator_.vanguard().parallelWells();
1172 std::pair<std::string, bool> value {wname,
true};
1173 auto candidate = std::lower_bound(parallelWells.begin(), parallelWells.end(), value);
1174 return candidate == parallelWells.end() || *candidate != value;
1177 void updateFluidInPlace_(
const ElementContext& elemCtx,
const unsigned dofIdx)
1179 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, 0);
1180 const unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, 0);
1181 const auto totVolume = elemCtx.simulator().model().dofTotalVolume(globalDofIdx);
1183 this->updateFluidInPlace_(globalDofIdx, intQuants, totVolume);
1186 void updateFluidInPlace_(
const unsigned globalDofIdx,
1187 const IntensiveQuantities& intQuants,
1188 const double totVolume)
1190 OPM_TIMEBLOCK_LOCAL(updateFluidInPlace);
1192 this->updateTotalVolumesAndPressures_(globalDofIdx, intQuants, totVolume);
1194 if (this->computeFip_) {
1195 this->updatePhaseInplaceVolumes_(globalDofIdx, intQuants, totVolume);
1199 void createLocalRegion_(std::vector<int>& region)
1201 std::size_t elemIdx = 0;
1202 for (
const auto& elem : elements(simulator_.gridView())) {
1203 if (elem.partitionType() != Dune::InteriorEntity) {
1204 region[elemIdx] = 0;
1211 template <
typename Flu
idState>
1212 void aggregateAverageDensityContributions_(
const FluidState& fs,
1213 const unsigned int globalDofIdx,
1216 auto pvCellValue = RegionPhasePoreVolAverage::CellValue{};
1217 pvCellValue.porv = porv;
1219 for (
auto phaseIdx = 0*FluidSystem::numPhases;
1220 phaseIdx < FluidSystem::numPhases; ++phaseIdx)
1222 if (! FluidSystem::phaseIsActive(phaseIdx)) {
1226 pvCellValue.value = getValue(fs.density(phaseIdx));
1227 pvCellValue.sat = getValue(fs.saturation(phaseIdx));
1229 this->regionAvgDensity_
1230 ->addCell(globalDofIdx,
1231 RegionPhasePoreVolAverage::Phase { phaseIdx },
1251 data::InterRegFlowMap::FlowRates
1252 getComponentSurfaceRates(
const ElementContext& elemCtx,
1253 const Scalar faceArea,
1254 const std::size_t scvfIdx,
1255 const std::size_t timeIdx)
const
1257 using Component = data::InterRegFlowMap::Component;
1259 auto rates = data::InterRegFlowMap::FlowRates {};
1261 const auto& extQuant = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
1263 const auto alpha = getValue(extQuant.extrusionFactor()) * faceArea;
1265 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1266 const auto& up = elemCtx
1267 .intensiveQuantities(extQuant.upstreamIndex(oilPhaseIdx), timeIdx);
1269 using FluidState = std::remove_cv_t<std::remove_reference_t<
1270 decltype(up.fluidState())>>;
1272 const auto pvtReg = up.pvtRegionIndex();
1274 const auto bO = getValue(getInvB_<FluidSystem, FluidState, Scalar>
1275 (up.fluidState(), oilPhaseIdx, pvtReg));
1277 const auto qO = alpha * bO * getValue(extQuant.volumeFlux(oilPhaseIdx));
1279 rates[Component::Oil] += qO;
1281 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1282 const auto Rs = getValue(
1283 BlackOil::getRs_<FluidSystem, FluidState, Scalar>
1284 (up.fluidState(), pvtReg));
1286 rates[Component::Gas] += qO * Rs;
1287 rates[Component::Disgas] += qO * Rs;
1291 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1292 const auto& up = elemCtx
1293 .intensiveQuantities(extQuant.upstreamIndex(gasPhaseIdx), timeIdx);
1295 using FluidState = std::remove_cv_t<std::remove_reference_t<
1296 decltype(up.fluidState())>>;
1298 const auto pvtReg = up.pvtRegionIndex();
1300 const auto bG = getValue(getInvB_<FluidSystem, FluidState, Scalar>
1301 (up.fluidState(), gasPhaseIdx, pvtReg));
1303 const auto qG = alpha * bG * getValue(extQuant.volumeFlux(gasPhaseIdx));
1305 rates[Component::Gas] += qG;
1307 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1308 const auto Rv = getValue(
1309 BlackOil::getRv_<FluidSystem, FluidState, Scalar>
1310 (up.fluidState(), pvtReg));
1312 rates[Component::Oil] += qG * Rv;
1313 rates[Component::Vapoil] += qG * Rv;
1317 if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
1318 const auto& up = elemCtx
1319 .intensiveQuantities(extQuant.upstreamIndex(waterPhaseIdx), timeIdx);
1321 using FluidState = std::remove_cv_t<std::remove_reference_t<
1322 decltype(up.fluidState())>>;
1324 const auto pvtReg = up.pvtRegionIndex();
1326 const auto bW = getValue(getInvB_<FluidSystem, FluidState, Scalar>
1327 (up.fluidState(), waterPhaseIdx, pvtReg));
1329 rates[Component::Water] +=
1330 alpha * bW * getValue(extQuant.volumeFlux(waterPhaseIdx));
1336 template <
typename Flu
idState>
1337 Scalar hydroCarbonFraction(
const FluidState& fs)
const
1339 if (this->eclState_.runspec().co2Storage()) {
1346 auto hydrocarbon = Scalar {0};
1347 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1348 hydrocarbon += getValue(fs.saturation(oilPhaseIdx));
1351 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1352 hydrocarbon += getValue(fs.saturation(gasPhaseIdx));
1358 void updateTotalVolumesAndPressures_(
const unsigned globalDofIdx,
1359 const IntensiveQuantities& intQuants,
1360 const double totVolume)
1362 const auto& fs = intQuants.fluidState();
1364 const double pv = totVolume * intQuants.porosity().value();
1365 const auto hydrocarbon = this->hydroCarbonFraction(fs);
1367 if (! this->hydrocarbonPoreVolume_.empty()) {
1368 this->fip_[Inplace::Phase::PoreVolume][globalDofIdx] =
1369 totVolume * intQuants.referencePorosity();
1371 this->dynamicPoreVolume_[globalDofIdx] = pv;
1372 this->hydrocarbonPoreVolume_[globalDofIdx] = pv * hydrocarbon;
1375 if (!this->pressureTimesHydrocarbonVolume_.empty() &&
1376 !this->pressureTimesPoreVolume_.empty())
1378 assert(this->hydrocarbonPoreVolume_.size() == this->pressureTimesHydrocarbonVolume_.size());
1379 assert(this->fip_[Inplace::Phase::PoreVolume].size() == this->pressureTimesPoreVolume_.size());
1381 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1382 this->pressureTimesPoreVolume_[globalDofIdx] =
1383 getValue(fs.pressure(oilPhaseIdx)) * pv;
1385 this->pressureTimesHydrocarbonVolume_[globalDofIdx] =
1386 this->pressureTimesPoreVolume_[globalDofIdx] * hydrocarbon;
1388 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1389 this->pressureTimesPoreVolume_[globalDofIdx] =
1390 getValue(fs.pressure(gasPhaseIdx)) * pv;
1392 this->pressureTimesHydrocarbonVolume_[globalDofIdx] =
1393 this->pressureTimesPoreVolume_[globalDofIdx] * hydrocarbon;
1395 else if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
1396 this->pressureTimesPoreVolume_[globalDofIdx] =
1397 getValue(fs.pressure(waterPhaseIdx)) * pv;
1402 void updatePhaseInplaceVolumes_(
const unsigned globalDofIdx,
1403 const IntensiveQuantities& intQuants,
1404 const double totVolume)
1406 std::array<Scalar, FluidSystem::numPhases> fip {};
1407 std::array<Scalar, FluidSystem::numPhases> fipr{};
1409 const auto& fs = intQuants.fluidState();
1410 const auto pv = totVolume * intQuants.porosity().value();
1412 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1413 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1417 const auto b = getValue(fs.invB(phaseIdx));
1418 const auto s = getValue(fs.saturation(phaseIdx));
1420 fipr[phaseIdx] = s * pv;
1421 fip [phaseIdx] = b * fipr[phaseIdx];
1424 this->updateInplaceVolumesSurface(globalDofIdx, fip);
1425 this->updateInplaceVolumesReservoir(globalDofIdx, fs, fipr);
1427 if (FluidSystem::phaseIsActive(oilPhaseIdx) &&
1428 FluidSystem::phaseIsActive(gasPhaseIdx))
1430 this->updateOilGasDistribution(globalDofIdx, fs, fip);
1433 if (FluidSystem::phaseIsActive(waterPhaseIdx) &&
1434 FluidSystem::phaseIsActive(gasPhaseIdx))
1436 this->updateGasWaterDistribution(globalDofIdx, fs, fip);
1439 if (FluidSystem::phaseIsActive(gasPhaseIdx) &&
1440 (!this->fip_[Inplace::Phase::CO2InGasPhaseInMob].empty() ||
1441 !this->fip_[Inplace::Phase::CO2InGasPhaseMob].empty() ||
1442 !this->fip_[Inplace::Phase::CO2MassInGasPhaseInMob].empty() ||
1443 !this->fip_[Inplace::Phase::CO2MassInGasPhaseMob].empty() ||
1444 !this->fip_[Inplace::Phase::CO2Mass].empty() ||
1445 !this->fip_[Inplace::Phase::CO2InGasPhaseInMobKrg].empty() ||
1446 !this->fip_[Inplace::Phase::CO2InGasPhaseMobKrg].empty() ||
1447 !this->fip_[Inplace::Phase::CO2MassInGasPhaseInMobKrg].empty() ||
1448 !this->fip_[Inplace::Phase::CO2MassInGasPhaseMobKrg].empty()))
1450 this->updateCO2InGas(globalDofIdx, pv, fs);
1453 if ((!this->fip_[Inplace::Phase::CO2InWaterPhase].empty() ||
1454 !this->fip_[Inplace::Phase::CO2MassInWaterPhase].empty() ||
1455 !this->fip_[Inplace::Phase::CO2Mass].empty()) &&
1456 (FluidSystem::phaseIsActive(waterPhaseIdx) ||
1457 FluidSystem::phaseIsActive(oilPhaseIdx)))
1459 this->updateCO2InWater(globalDofIdx, pv, fs);
1463 template <
typename FIPArray>
1464 void updateInplaceVolumesSurface(
const unsigned globalDofIdx,
1465 const FIPArray& fip)
1467 if (FluidSystem::phaseIsActive(oilPhaseIdx) &&
1468 !this->fip_[Inplace::Phase::OIL].empty())
1470 this->fip_[Inplace::Phase::OIL][globalDofIdx] = fip[oilPhaseIdx];
1473 if (FluidSystem::phaseIsActive(oilPhaseIdx) &&
1474 !this->fip_[Inplace::Phase::OilInLiquidPhase].empty())
1476 this->fip_[Inplace::Phase::OilInLiquidPhase][globalDofIdx] = fip[oilPhaseIdx];
1479 if (FluidSystem::phaseIsActive(gasPhaseIdx) &&
1480 !this->fip_[Inplace::Phase::GAS].empty())
1482 this->fip_[Inplace::Phase::GAS][globalDofIdx] = fip[gasPhaseIdx];
1485 if (FluidSystem::phaseIsActive(gasPhaseIdx) &&
1486 !this->fip_[Inplace::Phase::GasInGasPhase].empty())
1488 this->fip_[Inplace::Phase::GasInGasPhase][globalDofIdx] = fip[gasPhaseIdx];
1491 if (FluidSystem::phaseIsActive(waterPhaseIdx) &&
1492 !this->fip_[Inplace::Phase::WATER].empty())
1494 this->fip_[Inplace::Phase::WATER][globalDofIdx] = fip[waterPhaseIdx];
1498 template <
typename Flu
idState,
typename FIPArray>
1499 void updateInplaceVolumesReservoir(
const unsigned globalDofIdx,
1500 const FluidState& fs,
1501 const FIPArray& fipr)
1503 if (FluidSystem::phaseIsActive(oilPhaseIdx) &&
1504 !this->fip_[Inplace::Phase::OilResVolume].empty())
1506 this->fip_[Inplace::Phase::OilResVolume][globalDofIdx] = fipr[oilPhaseIdx];
1509 if (FluidSystem::phaseIsActive(gasPhaseIdx) &&
1510 !this->fip_[Inplace::Phase::GasResVolume].empty())
1512 this->fip_[Inplace::Phase::GasResVolume][globalDofIdx] = fipr[gasPhaseIdx];
1515 if (FluidSystem::phaseIsActive(waterPhaseIdx) &&
1516 !this->fip_[Inplace::Phase::WaterResVolume].empty())
1518 this->fip_[Inplace::Phase::WaterResVolume][globalDofIdx] = fipr[waterPhaseIdx];
1521 if (FluidSystem::phaseIsActive(waterPhaseIdx) &&
1522 !this->fip_[Inplace::Phase::SALT].empty())
1524 this->fip_[Inplace::Phase::SALT][globalDofIdx] =
1525 fipr[waterPhaseIdx] * fs.saltConcentration().value();
1529 template <
typename Flu
idState,
typename FIPArray>
1530 void updateOilGasDistribution(
const unsigned globalDofIdx,
1531 const FluidState& fs,
1532 const FIPArray& fip)
1535 const auto gasInPlaceLiquid = getValue(fs.Rs()) * fip[oilPhaseIdx];
1536 const auto oilInPlaceGas = getValue(fs.Rv()) * fip[gasPhaseIdx];
1538 if (!this->fip_[Inplace::Phase::GasInLiquidPhase].empty()) {
1539 this->fip_[Inplace::Phase::GasInLiquidPhase][globalDofIdx] = gasInPlaceLiquid;
1542 if (!this->fip_[Inplace::Phase::OilInGasPhase].empty()) {
1543 this->fip_[Inplace::Phase::OilInGasPhase][globalDofIdx] = oilInPlaceGas;
1547 if (!this->fip_[Inplace::Phase::OIL].empty()) {
1548 this->fip_[Inplace::Phase::OIL][globalDofIdx] += oilInPlaceGas;
1551 if (!this->fip_[Inplace::Phase::GAS].empty()) {
1552 this->fip_[Inplace::Phase::GAS][globalDofIdx] += gasInPlaceLiquid;
1556 template <
typename Flu
idState,
typename FIPArray>
1557 void updateGasWaterDistribution(
const unsigned globalDofIdx,
1558 const FluidState& fs,
1559 const FIPArray& fip)
1562 const auto gasInPlaceWater = getValue(fs.Rsw()) * fip[waterPhaseIdx];
1563 const auto waterInPlaceGas = getValue(fs.Rvw()) * fip[gasPhaseIdx];
1565 if (!this->fip_[Inplace::Phase::WaterInGasPhase].empty()) {
1566 this->fip_[Inplace::Phase::WaterInGasPhase][globalDofIdx] = waterInPlaceGas;
1569 if (!this->fip_[Inplace::Phase::WaterInWaterPhase].empty()) {
1570 this->fip_[Inplace::Phase::WaterInWaterPhase][globalDofIdx] = fip[waterPhaseIdx];
1574 if (!this->fip_[Inplace::Phase::GasInLiquidPhase].empty() &&
1575 !FluidSystem::phaseIsActive(oilPhaseIdx))
1577 this->fip_[Inplace::Phase::GasInLiquidPhase][globalDofIdx] = gasInPlaceWater;
1581 if (!this->fip_[Inplace::Phase::WATER].empty()) {
1582 this->fip_[Inplace::Phase::WATER][globalDofIdx] += waterInPlaceGas;
1585 if (!this->fip_[Inplace::Phase::GAS].empty()) {
1586 this->fip_[Inplace::Phase::GAS][globalDofIdx] += gasInPlaceWater;
1590 template <
typename Flu
idState>
1591 void updateCO2InGas(
const unsigned globalDofIdx,
1593 const FluidState& fs)
1595 const auto& scaledDrainageInfo = this->simulator_.problem().materialLawManager()
1596 ->oilWaterScaledEpsInfoDrainage(globalDofIdx);
1598 Scalar sgcr = scaledDrainageInfo.Sgcr;
1599 if (this->simulator_.problem().materialLawManager()->enableHysteresis()) {
1600 const auto& matParams = simulator_.problem().materialLawParams(globalDofIdx);
1601 sgcr = MaterialLaw::trappedGasSaturation(matParams);
1604 const double sg = getValue(fs.saturation(gasPhaseIdx));
1605 const double rhog = getValue(fs.density(gasPhaseIdx));
1606 const double xgW = FluidSystem::phaseIsActive(waterPhaseIdx)
1607 ? FluidSystem::convertRvwToXgW(getValue(fs.Rvw()), fs.pvtRegionIndex())
1608 : FluidSystem::convertRvToXgO(getValue(fs.Rv()), fs.pvtRegionIndex());
1610 const Scalar mM = FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex());
1611 const Scalar massGas = (1 - xgW) * pv * rhog;
1612 if (!this->fip_[Inplace::Phase::CO2Mass].empty()) {
1613 this->fip_[Inplace::Phase::CO2Mass][globalDofIdx] = massGas;
1616 if (!this->fip_[Inplace::Phase::CO2InGasPhaseInMob].empty()) {
1617 const Scalar imMobileGas = massGas / mM * std::min(sgcr , sg);
1618 this->fip_[Inplace::Phase::CO2InGasPhaseInMob][globalDofIdx] = imMobileGas;
1621 if (!this->fip_[Inplace::Phase::CO2InGasPhaseMob].empty()) {
1622 const Scalar mobileGas = massGas / mM * std::max(0.0, sg - sgcr);
1623 this->fip_[Inplace::Phase::CO2InGasPhaseMob][globalDofIdx] = mobileGas;
1626 if (!this->fip_[Inplace::Phase::CO2InGasPhaseInMobKrg].empty()) {
1628 const Scalar imMobileGasKrg = massGas / mM * sg;
1629 this->fip_[Inplace::Phase::CO2InGasPhaseInMobKrg][globalDofIdx] = imMobileGasKrg;
1631 this->fip_[Inplace::Phase::CO2InGasPhaseInMobKrg][globalDofIdx] = 0;
1635 if (!this->fip_[Inplace::Phase::CO2InGasPhaseMobKrg].empty()) {
1637 const Scalar mobileGasKrg = massGas / mM * sg;
1638 this->fip_[Inplace::Phase::CO2InGasPhaseMobKrg][globalDofIdx] = mobileGasKrg;
1640 this->fip_[Inplace::Phase::CO2InGasPhaseMobKrg][globalDofIdx] = 0;
1644 if (!this->fip_[Inplace::Phase::CO2MassInGasPhaseInMob].empty()) {
1645 const Scalar imMobileMassGas = massGas * std::min(sgcr , sg);
1646 this->fip_[Inplace::Phase::CO2MassInGasPhaseInMob][globalDofIdx] = imMobileMassGas;
1649 if (!this->fip_[Inplace::Phase::CO2MassInGasPhaseMob].empty()) {
1650 const Scalar mobileMassGas = massGas * std::max(0.0, sg - sgcr);
1651 this->fip_[Inplace::Phase::CO2MassInGasPhaseMob][globalDofIdx] = mobileMassGas;
1654 if (!this->fip_[Inplace::Phase::CO2MassInGasPhaseInMobKrg].empty()) {
1656 const Scalar imMobileMassGasKrg = massGas * sg;
1657 this->fip_[Inplace::Phase::CO2MassInGasPhaseInMobKrg][globalDofIdx] = imMobileMassGasKrg;
1659 this->fip_[Inplace::Phase::CO2MassInGasPhaseInMobKrg][globalDofIdx] = 0;
1663 if (!this->fip_[Inplace::Phase::CO2MassInGasPhaseMobKrg].empty()) {
1665 const Scalar mobileMassGasKrg = massGas * sg;
1666 this->fip_[Inplace::Phase::CO2MassInGasPhaseMobKrg][globalDofIdx] = mobileMassGasKrg;
1668 this->fip_[Inplace::Phase::CO2MassInGasPhaseMobKrg][globalDofIdx] = 0;
1673 template <
typename Flu
idState>
1674 void updateCO2InWater(
const unsigned globalDofIdx,
1676 const FluidState& fs)
1678 const auto co2InWater = FluidSystem::phaseIsActive(oilPhaseIdx)
1679 ? this->co2InWaterFromOil(fs, pv)
1680 : this->co2InWaterFromWater(fs, pv);
1682 const Scalar mM = FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex());
1683 if (!this->fip_[Inplace::Phase::CO2Mass].empty()) {
1684 this->fip_[Inplace::Phase::CO2Mass][globalDofIdx] += co2InWater * mM;
1686 if (!this->fip_[Inplace::Phase::CO2MassInWaterPhase].empty()) {
1687 this->fip_[Inplace::Phase::CO2MassInWaterPhase][globalDofIdx] = co2InWater * mM;
1689 if (!this->fip_[Inplace::Phase::CO2InWaterPhase].empty()) {
1690 this->fip_[Inplace::Phase::CO2InWaterPhase][globalDofIdx] = co2InWater;
1694 template <
typename Flu
idState>
1695 Scalar co2InWaterFromWater(
const FluidState& fs,
const double pv)
const
1697 const double rhow = getValue(fs.density(waterPhaseIdx));
1698 const double sw = getValue(fs.saturation(waterPhaseIdx));
1699 const double xwG = FluidSystem::convertRswToXwG(getValue(fs.Rsw()), fs.pvtRegionIndex());
1701 const Scalar mM = FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex());
1703 return xwG * pv * rhow * sw / mM;
1706 template <
typename Flu
idState>
1707 Scalar co2InWaterFromOil(
const FluidState& fs,
const double pv)
const
1709 const double rhoo = getValue(fs.density(oilPhaseIdx));
1710 const double so = getValue(fs.saturation(oilPhaseIdx));
1711 const double xoG = FluidSystem::convertRsToXoG(getValue(fs.Rs()), fs.pvtRegionIndex());
1713 const Scalar mM = FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex());
1715 return xoG * pv * rhoo * so / mM;
1718 const Simulator& simulator_;