My Project
Loading...
Searching...
No Matches
OutputBlackoilModule.hpp
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*
4 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 2 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18
19 Consult the COPYING file in the top-level source directory of this
20 module for the precise wording of the license and the list of
21 copyright holders.
22*/
27#ifndef OPM_OUTPUT_BLACK_OIL_MODULE_HPP
28#define OPM_OUTPUT_BLACK_OIL_MODULE_HPP
29
30#include <dune/common/fvector.hh>
31
32#include <opm/simulators/utils/moduleVersion.hpp>
33
34#include <opm/common/Exceptions.hpp>
35#include <opm/common/TimingMacros.hpp>
36#include <opm/common/OpmLog/OpmLog.hpp>
37#include <opm/common/utility/Visitor.hpp>
38
39#include <opm/input/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
40
41#include <opm/material/common/Valgrind.hpp>
42#include <opm/material/fluidmatrixinteractions/EclEpsScalingPoints.hpp>
43#include <opm/material/fluidstates/BlackOilFluidState.hpp>
44#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
45
50
51#include <opm/output/data/Cells.hpp>
52#include <opm/output/eclipse/EclipseIO.hpp>
53#include <opm/output/eclipse/Inplace.hpp>
54
55#include <opm/simulators/flow/CollectDataOnIORank.hpp>
59
60#include <algorithm>
61#include <array>
62#include <cassert>
63#include <cstddef>
64#include <functional>
65#include <stdexcept>
66#include <string>
67#include <type_traits>
68#include <utility>
69#include <vector>
70
71namespace Opm {
72
73// forward declaration
74template <class TypeTag>
75class EcfvDiscretization;
76
83template <class TypeTag>
84class OutputBlackOilModule : public GenericOutputBlackoilModule<GetPropType<TypeTag, Properties::FluidSystem>>
85{
95 using FluidState = typename IntensiveQuantities::FluidState;
97 using Element = typename GridView::template Codim<0>::Entity;
98 using ElementIterator = typename GridView::template Codim<0>::Iterator;
101 using Dir = FaceDir::DirEnum;
107
108 static constexpr int conti0EqIdx = Indices::conti0EqIdx;
109 static constexpr int numPhases = FluidSystem::numPhases;
110 static constexpr int oilPhaseIdx = FluidSystem::oilPhaseIdx;
111 static constexpr int gasPhaseIdx = FluidSystem::gasPhaseIdx;
112 static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
113 static constexpr int gasCompIdx = FluidSystem::gasCompIdx;
114 static constexpr int oilCompIdx = FluidSystem::oilCompIdx;
115 static constexpr int waterCompIdx = FluidSystem::waterCompIdx;
116 enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
117 enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
118
119 template<int idx, class VectorType>
120 static Scalar value_or_zero(const VectorType& v)
121 {
122 if constexpr (idx == -1) {
123 return 0.0;
124 } else {
125 return v.empty() ? 0.0 : v[idx];
126 }
127 }
128
129public:
130 OutputBlackOilModule(const Simulator& simulator,
131 const SummaryConfig& smryCfg,
132 const CollectDataOnIORankType& collectOnIORank)
133 : BaseType(simulator.vanguard().eclState(),
134 simulator.vanguard().schedule(),
135 smryCfg,
136 simulator.vanguard().summaryState(),
138 [this](const int idx)
139 { return simulator_.problem().eclWriter().collectOnIORank().localIdxToGlobalIdx(idx); },
140 simulator.vanguard().grid().comm(),
151 , simulator_(simulator)
152 , collectOnIORank_(collectOnIORank)
153 {
154 for (auto& region_pair : this->regions_) {
155 this->createLocalRegion_(region_pair.second);
156 }
157
158 auto isCartIdxOnThisRank = [&collectOnIORank](const int idx) {
159 return collectOnIORank.isCartIdxOnThisRank(idx);
160 };
161
162 this->setupBlockData(isCartIdxOnThisRank);
163
164 if (! Parameters::Get<Parameters::OwnerCellsFirst>()) {
165 const std::string msg = "The output code does not support --owner-cells-first=false.";
166 if (collectOnIORank.isIORank()) {
167 OpmLog::error(msg);
168 }
169 OPM_THROW_NOLOG(std::runtime_error, msg);
170 }
171
172 if (smryCfg.match("[FB]PP[OGW]") || smryCfg.match("RPP[OGW]*")) {
173 auto rset = this->eclState_.fieldProps().fip_regions();
174 rset.push_back("PVTNUM");
175
176 // Note: We explicitly use decltype(auto) here because the
177 // default scheme (-> auto) will deduce an undesirable type. We
178 // need the "reference to vector" semantics in this instance.
179 this->regionAvgDensity_
180 .emplace(this->simulator_.gridView().comm(),
181 FluidSystem::numPhases, rset,
182 [fp = std::cref(this->eclState_.fieldProps())]
183 (const std::string& rsetName) -> decltype(auto)
184 { return fp.get().get_int(rsetName); });
185 }
186 }
187
192 void
193 allocBuffers(const unsigned bufferSize,
194 const unsigned reportStepNum,
195 const bool substep,
196 const bool log,
197 const bool isRestart)
198 {
199 if (! std::is_same<Discretization, EcfvDiscretization<TypeTag>>::value) {
200 return;
201 }
202
203 const auto& problem = this->simulator_.problem();
204
205 this->doAllocBuffers(bufferSize,
206 reportStepNum,
207 substep,
208 log,
209 isRestart,
210 &problem.materialLawManager()->hysteresisConfig(),
211 problem.eclWriter().getOutputNnc().size());
212 }
213
215 void setupExtractors(const bool isSubStep,
216 const int reportStepNum)
217 {
218 this->setupElementExtractors_();
219 this->setupBlockExtractors_(isSubStep, reportStepNum);
220 }
221
224 {
225 this->extractors_.clear();
226 this->blockExtractors_.clear();
227 this->extraBlockExtractors_.clear();
228 }
229
234 void processElement(const ElementContext& elemCtx)
235 {
237 if (!std::is_same<Discretization, EcfvDiscretization<TypeTag>>::value) {
238 return;
239 }
240
241 if (this->extractors_.empty()) {
242 assert(0);
243 }
244
245 const auto& matLawManager = simulator_.problem().materialLawManager();
246
248 for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
249 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
250 const auto& fs = intQuants.fluidState();
251
252 const typename Extractor::Context ectx{
253 elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0),
254 elemCtx.primaryVars(dofIdx, /*timeIdx=*/0).pvtRegionIndex(),
255 elemCtx.simulator().episodeIndex(),
256 fs,
257 intQuants,
259 };
260
261 if (matLawManager->enableHysteresis()) {
262 if (FluidSystem::phaseIsActive(oilPhaseIdx) && FluidSystem::phaseIsActive(waterPhaseIdx)) {
263 matLawManager->oilWaterHysteresisParams(hysterParams.somax,
264 hysterParams.swmax,
265 hysterParams.swmin,
266 ectx.globalDofIdx);
267 }
268 if (FluidSystem::phaseIsActive(oilPhaseIdx) && FluidSystem::phaseIsActive(gasPhaseIdx)) {
269 matLawManager->gasOilHysteresisParams(hysterParams.sgmax,
270 hysterParams.shmax,
271 hysterParams.somin,
272 ectx.globalDofIdx);
273 }
274 }
275
276 Extractor::process(ectx, extractors_);
277 }
278 }
279
280 void processElementBlockData(const ElementContext& elemCtx)
281 {
282 OPM_TIMEBLOCK_LOCAL(processElementBlockData);
283 if (!std::is_same<Discretization, EcfvDiscretization<TypeTag>>::value) {
284 return;
285 }
286
287 if (this->blockExtractors_.empty() && this->extraBlockExtractors_.empty()) {
288 return;
289 }
290
291 for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
292 // Adding block data
293 const auto globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
294 const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx);
295
296 const auto be_it = this->blockExtractors_.find(cartesianIdx);
297 const auto bee_it = this->extraBlockExtractors_.find(cartesianIdx);
298 if (be_it == this->blockExtractors_.end() &&
299 bee_it == this->extraBlockExtractors_.end())
300 {
301 continue;
302 }
303
304 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
305 const auto& fs = intQuants.fluidState();
306
307 const typename BlockExtractor::Context ectx{
308 globalDofIdx,
309 dofIdx,
310 fs,
311 intQuants,
312 elemCtx,
313 };
314
315 if (be_it != this->blockExtractors_.end()) {
317 }
318 if (bee_it != this->extraBlockExtractors_.end()) {
320 }
321 }
322 }
323
324 void outputFipAndResvLog(const Inplace& inplace,
325 const std::size_t reportStepNum,
326 double elapsed,
327 boost::posix_time::ptime currentDate,
328 const bool substep,
329 const Parallel::Communication& comm)
330 {
331
332 if (comm.rank() != 0) {
333 return;
334 }
335
336 // For report step 0 we use the RPTSOL config, else derive from RPTSCHED
337 std::unique_ptr<FIPConfig> fipSched;
338 if (reportStepNum > 0) {
339 const auto& rpt = this->schedule_[reportStepNum-1].rpt_config.get();
340 fipSched = std::make_unique<FIPConfig>(rpt);
341 }
342 const FIPConfig& fipc = reportStepNum == 0 ? this->eclState_.getEclipseConfig().fip()
343 : *fipSched;
344
345 if (!substep && !this->forceDisableFipOutput_ && fipc.output(FIPConfig::OutputField::FIELD)) {
346
347 this->logOutput_.timeStamp("BALANCE", elapsed, reportStepNum, currentDate);
348
349 const auto& initial_inplace = this->initialInplace().value();
350 this->logOutput_.fip(inplace, initial_inplace, "");
351
352 if (fipc.output(FIPConfig::OutputField::FIPNUM)) {
353 this->logOutput_.fip(inplace, initial_inplace, "FIPNUM");
354
355 if (fipc.output(FIPConfig::OutputField::RESV))
356 this->logOutput_.fipResv(inplace, "FIPNUM");
357 }
358
359 if (fipc.output(FIPConfig::OutputField::FIP)) {
360 for (const auto& reg : this->regions_) {
361 if (reg.first != "FIPNUM") {
362 std::ostringstream ss;
363 ss << "BAL" << reg.first.substr(3);
364 this->logOutput_.timeStamp(ss.str(), elapsed, reportStepNum, currentDate);
365 this->logOutput_.fip(inplace, initial_inplace, reg.first);
366
367 if (fipc.output(FIPConfig::OutputField::RESV))
368 this->logOutput_.fipResv(inplace, reg.first);
369 }
370 }
371 }
372 }
373 }
374
403 template <class ActiveIndex, class CartesianIndex>
404 void processFluxes(const ElementContext& elemCtx,
405 ActiveIndex&& activeIndex,
406 CartesianIndex&& cartesianIndex)
407 {
409 const auto identifyCell = [&activeIndex, &cartesianIndex](const Element& elem)
411 {
412 const auto cellIndex = activeIndex(elem);
413
414 return {
415 static_cast<int>(cellIndex),
416 cartesianIndex(cellIndex),
417 elem.partitionType() == Dune::InteriorEntity
418 };
419 };
420
421 const auto timeIdx = 0u;
422 const auto& stencil = elemCtx.stencil(timeIdx);
423 const auto numInteriorFaces = elemCtx.numInteriorFaces(timeIdx);
424
425 for (auto scvfIdx = 0 * numInteriorFaces; scvfIdx < numInteriorFaces; ++scvfIdx) {
426 const auto& face = stencil.interiorFace(scvfIdx);
427 const auto left = identifyCell(stencil.element(face.interiorIndex()));
428 const auto right = identifyCell(stencil.element(face.exteriorIndex()));
429
430 const auto rates = this->
431 getComponentSurfaceRates(elemCtx, face.area(), scvfIdx, timeIdx);
432
433 this->interRegionFlows_.addConnection(left, right, rates);
434 }
435 }
436
442 {
443 // Inter-region flow rates. Note: ".clear()" prepares to accumulate
444 // contributions per bulk connection between FIP regions.
445 this->interRegionFlows_.clear();
446 }
447
452 {
453 this->interRegionFlows_.compress();
454 }
455
460 {
461 return this->interRegionFlows_;
462 }
463
464 template <class FluidState>
465 void assignToFluidState(FluidState& fs, unsigned elemIdx) const
466 {
467 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
468 if (this->saturation_[phaseIdx].empty())
469 continue;
470
471 fs.setSaturation(phaseIdx, this->saturation_[phaseIdx][elemIdx]);
472 }
473
474 if (!this->fluidPressure_.empty()) {
475 // this assumes that capillary pressures only depend on the phase saturations
476 // and possibly on temperature. (this is always the case for ECL problems.)
477 std::array<Scalar, numPhases> pc = {0};
478 const MaterialLawParams& matParams = simulator_.problem().materialLawParams(elemIdx);
479 MaterialLaw::capillaryPressures(pc, matParams, fs);
480 Valgrind::CheckDefined(this->fluidPressure_[elemIdx]);
481 Valgrind::CheckDefined(pc);
482 const auto& pressure = this->fluidPressure_[elemIdx];
483 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
484 if (!FluidSystem::phaseIsActive(phaseIdx))
485 continue;
486
487 if (Indices::oilEnabled)
488 fs.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx]));
489 else if (Indices::gasEnabled)
490 fs.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx]));
491 else if (Indices::waterEnabled)
492 //single (water) phase
493 fs.setPressure(phaseIdx, pressure);
494 }
495 }
496
497 if (!this->temperature_.empty())
498 fs.setTemperature(this->temperature_[elemIdx]);
499 if (!this->rs_.empty())
500 fs.setRs(this->rs_[elemIdx]);
501 if (!this->rsw_.empty())
502 fs.setRsw(this->rsw_[elemIdx]);
503 if (!this->rv_.empty())
504 fs.setRv(this->rv_[elemIdx]);
505 if (!this->rvw_.empty())
506 fs.setRvw(this->rvw_[elemIdx]);
507 }
508
509 void initHysteresisParams(Simulator& simulator, unsigned elemIdx) const
510 {
511 if (!this->soMax_.empty())
512 simulator.problem().setMaxOilSaturation(elemIdx, this->soMax_[elemIdx]);
513
514 if (simulator.problem().materialLawManager()->enableHysteresis()) {
515 auto matLawManager = simulator.problem().materialLawManager();
516
517 if (FluidSystem::phaseIsActive(oilPhaseIdx)
518 && FluidSystem::phaseIsActive(waterPhaseIdx)) {
519 Scalar somax = 2.0;
520 Scalar swmax = -2.0;
521 Scalar swmin = 2.0;
522
523 if (matLawManager->enableNonWettingHysteresis()) {
524 if (!this->soMax_.empty()) {
525 somax = this->soMax_[elemIdx];
526 }
527 }
528 if (matLawManager->enableWettingHysteresis()) {
529 if (!this->swMax_.empty()) {
530 swmax = this->swMax_[elemIdx];
531 }
532 }
533 if (matLawManager->enablePCHysteresis()) {
534 if (!this->swmin_.empty()) {
535 swmin = this->swmin_[elemIdx];
536 }
537 }
538 matLawManager->setOilWaterHysteresisParams(
539 somax, swmax, swmin, elemIdx);
540 }
541 if (FluidSystem::phaseIsActive(oilPhaseIdx)
542 && FluidSystem::phaseIsActive(gasPhaseIdx)) {
543 Scalar sgmax = 2.0;
544 Scalar shmax = -2.0;
545 Scalar somin = 2.0;
546
547 if (matLawManager->enableNonWettingHysteresis()) {
548 if (!this->sgmax_.empty()) {
549 sgmax = this->sgmax_[elemIdx];
550 }
551 }
552 if (matLawManager->enableWettingHysteresis()) {
553 if (!this->shmax_.empty()) {
554 shmax = this->shmax_[elemIdx];
555 }
556 }
557 if (matLawManager->enablePCHysteresis()) {
558 if (!this->somin_.empty()) {
559 somin = this->somin_[elemIdx];
560 }
561 }
562 matLawManager->setGasOilHysteresisParams(
563 sgmax, shmax, somin, elemIdx);
564 }
565
566 }
567
568 if (simulator_.vanguard().eclState().fieldProps().has_double("SWATINIT")) {
569 simulator.problem().materialLawManager()
570 ->applyRestartSwatInit(elemIdx, this->ppcw_[elemIdx]);
571 }
572 }
573
574 void updateFluidInPlace(const ElementContext& elemCtx)
575 {
576 for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
577 updateFluidInPlace_(elemCtx, dofIdx);
578 }
579 }
580
581 void updateFluidInPlace(const unsigned globalDofIdx,
582 const IntensiveQuantities& intQuants,
583 const double totVolume)
584 {
585 this->updateFluidInPlace_(globalDofIdx, intQuants, totVolume);
586 }
587
588private:
589 template <typename T>
590 using RemoveCVR = std::remove_cv_t<std::remove_reference_t<T>>;
591
592 template <typename, class = void>
593 struct HasGeoMech : public std::false_type {};
594
595 template <typename Problem>
596 struct HasGeoMech<
597 Problem, std::void_t<decltype(std::declval<Problem>().geoMechModel())>
598 > : public std::true_type {};
599
600 bool isDefunctParallelWell(std::string wname) const override
601 {
602 if (simulator_.gridView().comm().size() == 1)
603 return false;
604 const auto& parallelWells = simulator_.vanguard().parallelWells();
605 std::pair<std::string, bool> value {wname, true};
606 auto candidate = std::lower_bound(parallelWells.begin(), parallelWells.end(), value);
607 return candidate == parallelWells.end() || *candidate != value;
608 }
609
610 void updateFluidInPlace_(const ElementContext& elemCtx, const unsigned dofIdx)
611 {
612 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
613 const unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
614 const auto totVolume = elemCtx.simulator().model().dofTotalVolume(globalDofIdx);
615
616 this->updateFluidInPlace_(globalDofIdx, intQuants, totVolume);
617 }
618
619 void updateFluidInPlace_(const unsigned globalDofIdx,
620 const IntensiveQuantities& intQuants,
621 const double totVolume)
622 {
623 OPM_TIMEBLOCK_LOCAL(updateFluidInPlace);
624
625 this->updateTotalVolumesAndPressures_(globalDofIdx, intQuants, totVolume);
626
627 if (this->computeFip_) {
628 this->updatePhaseInplaceVolumes_(globalDofIdx, intQuants, totVolume);
629 }
630 }
631
632 void createLocalRegion_(std::vector<int>& region)
633 {
634 std::size_t elemIdx = 0;
635 for (const auto& elem : elements(simulator_.gridView())) {
636 if (elem.partitionType() != Dune::InteriorEntity) {
637 region[elemIdx] = 0;
638 }
639
640 ++elemIdx;
641 }
642 }
643
644 template <typename FluidState>
645 void aggregateAverageDensityContributions_(const FluidState& fs,
646 const unsigned int globalDofIdx,
647 const double porv)
648 {
649 auto pvCellValue = RegionPhasePoreVolAverage::CellValue{};
650 pvCellValue.porv = porv;
651
652 for (auto phaseIdx = 0*FluidSystem::numPhases;
653 phaseIdx < FluidSystem::numPhases; ++phaseIdx)
654 {
655 if (! FluidSystem::phaseIsActive(phaseIdx)) {
656 continue;
657 }
658
659 pvCellValue.value = getValue(fs.density(phaseIdx));
660 pvCellValue.sat = getValue(fs.saturation(phaseIdx));
661
662 this->regionAvgDensity_
663 ->addCell(globalDofIdx,
664 RegionPhasePoreVolAverage::Phase { phaseIdx },
666 }
667 }
668
684 data::InterRegFlowMap::FlowRates
685 getComponentSurfaceRates(const ElementContext& elemCtx,
686 const Scalar faceArea,
687 const std::size_t scvfIdx,
688 const std::size_t timeIdx) const
689 {
690 using Component = data::InterRegFlowMap::Component;
691
692 auto rates = data::InterRegFlowMap::FlowRates {};
693
694 const auto& extQuant = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
695
696 const auto alpha = getValue(extQuant.extrusionFactor()) * faceArea;
697
698 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
699 const auto& up = elemCtx
700 .intensiveQuantities(extQuant.upstreamIndex(oilPhaseIdx), timeIdx);
701
702 const auto pvtReg = up.pvtRegionIndex();
703
705 (up.fluidState(), oilPhaseIdx, pvtReg));
706
707 const auto qO = alpha * bO * getValue(extQuant.volumeFlux(oilPhaseIdx));
708
709 rates[Component::Oil] += qO;
710
711 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
712 const auto Rs = getValue(
713 BlackOil::getRs_<FluidSystem, FluidState, Scalar>
714 (up.fluidState(), pvtReg));
715
716 rates[Component::Gas] += qO * Rs;
717 rates[Component::Disgas] += qO * Rs;
718 }
719 }
720
721 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
722 const auto& up = elemCtx
723 .intensiveQuantities(extQuant.upstreamIndex(gasPhaseIdx), timeIdx);
724
725 const auto pvtReg = up.pvtRegionIndex();
726
728 (up.fluidState(), gasPhaseIdx, pvtReg));
729
730 const auto qG = alpha * bG * getValue(extQuant.volumeFlux(gasPhaseIdx));
731
732 rates[Component::Gas] += qG;
733
734 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
735 const auto Rv = getValue(
736 BlackOil::getRv_<FluidSystem, FluidState, Scalar>
737 (up.fluidState(), pvtReg));
738
739 rates[Component::Oil] += qG * Rv;
740 rates[Component::Vapoil] += qG * Rv;
741 }
742 }
743
744 if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
745 const auto& up = elemCtx
746 .intensiveQuantities(extQuant.upstreamIndex(waterPhaseIdx), timeIdx);
747
748 const auto pvtReg = up.pvtRegionIndex();
749
751 (up.fluidState(), waterPhaseIdx, pvtReg));
752
753 rates[Component::Water] +=
754 alpha * bW * getValue(extQuant.volumeFlux(waterPhaseIdx));
755 }
756
757 return rates;
758 }
759
760 template <typename FluidState>
761 Scalar hydroCarbonFraction(const FluidState& fs) const
762 {
763 if (this->eclState_.runspec().co2Storage()) {
764 // CO2 storage: Hydrocarbon volume is full pore-volume.
765 return 1.0;
766 }
767
768 // Common case. Hydrocarbon volume is fraction occupied by actual
769 // hydrocarbons.
770 auto hydrocarbon = Scalar {0};
771 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
772 hydrocarbon += getValue(fs.saturation(oilPhaseIdx));
773 }
774
775 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
776 hydrocarbon += getValue(fs.saturation(gasPhaseIdx));
777 }
778
779 return hydrocarbon;
780 }
781
782 void updateTotalVolumesAndPressures_(const unsigned globalDofIdx,
783 const IntensiveQuantities& intQuants,
784 const double totVolume)
785 {
786 const auto& fs = intQuants.fluidState();
787
788 const double pv = totVolume * intQuants.porosity().value();
789 const auto hydrocarbon = this->hydroCarbonFraction(fs);
790
791 if (! this->hydrocarbonPoreVolume_.empty()) {
792 this->fipC_.assignPoreVolume(globalDofIdx,
793 totVolume * intQuants.referencePorosity());
794
795 this->dynamicPoreVolume_[globalDofIdx] = pv;
796 this->hydrocarbonPoreVolume_[globalDofIdx] = pv * hydrocarbon;
797 }
798
799 if (!this->pressureTimesHydrocarbonVolume_.empty() &&
800 !this->pressureTimesPoreVolume_.empty())
801 {
802 assert(this->hydrocarbonPoreVolume_.size() == this->pressureTimesHydrocarbonVolume_.size());
803 assert(this->fipC_.get(Inplace::Phase::PoreVolume).size() == this->pressureTimesPoreVolume_.size());
804
805 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
806 this->pressureTimesPoreVolume_[globalDofIdx] =
807 getValue(fs.pressure(oilPhaseIdx)) * pv;
808
809 this->pressureTimesHydrocarbonVolume_[globalDofIdx] =
810 this->pressureTimesPoreVolume_[globalDofIdx] * hydrocarbon;
811 }
812 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
813 this->pressureTimesPoreVolume_[globalDofIdx] =
814 getValue(fs.pressure(gasPhaseIdx)) * pv;
815
816 this->pressureTimesHydrocarbonVolume_[globalDofIdx] =
817 this->pressureTimesPoreVolume_[globalDofIdx] * hydrocarbon;
818 }
819 else if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
820 this->pressureTimesPoreVolume_[globalDofIdx] =
821 getValue(fs.pressure(waterPhaseIdx)) * pv;
822 }
823 }
824 }
825
826 void updatePhaseInplaceVolumes_(const unsigned globalDofIdx,
827 const IntensiveQuantities& intQuants,
828 const double totVolume)
829 {
830 std::array<Scalar, FluidSystem::numPhases> fip {};
831 std::array<Scalar, FluidSystem::numPhases> fipr{}; // at reservoir condition
832
833 const auto& fs = intQuants.fluidState();
834 const auto pv = totVolume * intQuants.porosity().value();
835
836 for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
837 if (!FluidSystem::phaseIsActive(phaseIdx)) {
838 continue;
839 }
840
841 const auto b = getValue(fs.invB(phaseIdx));
842 const auto s = getValue(fs.saturation(phaseIdx));
843
844 fipr[phaseIdx] = s * pv;
845 fip [phaseIdx] = b * fipr[phaseIdx];
846 }
847
848 this->fipC_.assignVolumesSurface(globalDofIdx, fip);
849 this->fipC_.assignVolumesReservoir(globalDofIdx,
850 fs.saltConcentration().value(),
851 fipr);
852
853 if (FluidSystem::phaseIsActive(oilPhaseIdx) &&
854 FluidSystem::phaseIsActive(gasPhaseIdx))
855 {
856 this->updateOilGasDistribution(globalDofIdx, fs, fip);
857 }
858
859 if (FluidSystem::phaseIsActive(waterPhaseIdx) &&
860 FluidSystem::phaseIsActive(gasPhaseIdx))
861 {
862 this->updateGasWaterDistribution(globalDofIdx, fs, fip);
863 }
864
865 if (FluidSystem::phaseIsActive(gasPhaseIdx) &&
866 this->fipC_.hasCo2InGas())
867 {
868 this->updateCO2InGas(globalDofIdx, pv, intQuants);
869 }
870
871 if (this->fipC_.hasCo2InWater() &&
872 (FluidSystem::phaseIsActive(waterPhaseIdx) ||
873 FluidSystem::phaseIsActive(oilPhaseIdx)))
874 {
875 this->updateCO2InWater(globalDofIdx, pv, fs);
876 }
877
878 if constexpr(enableMICP) {
879 const auto surfVolWat = pv * getValue(fs.invB(waterPhaseIdx));
880 if (this->fipC_.hasMicrobialMass()) {
881 this->updateMicrobialMass(globalDofIdx, intQuants, surfVolWat);
882 }
883 if (this->fipC_.hasOxygenMass()) {
884 this->updateOxygenMass(globalDofIdx, intQuants, surfVolWat);
885 }
886 if (this->fipC_.hasUreaMass()) {
887 this->updateUreaMass(globalDofIdx, intQuants, surfVolWat);
888 }
889 if (this->fipC_.hasBiofilmMass()) {
890 this->updateBiofilmMass(globalDofIdx, intQuants, totVolume);
891 }
892 if (this->fipC_.hasCalciteMass()) {
893 this->updateCalciteMass(globalDofIdx, intQuants, totVolume);
894 }
895 }
896 }
897
898 template <typename FluidState, typename FIPArray>
899 void updateOilGasDistribution(const unsigned globalDofIdx,
900 const FluidState& fs,
901 const FIPArray& fip)
902 {
903 // Gas dissolved in oil and vaporized oil
904 const auto gasInPlaceLiquid = getValue(fs.Rs()) * fip[oilPhaseIdx];
905 const auto oilInPlaceGas = getValue(fs.Rv()) * fip[gasPhaseIdx];
906
907 this->fipC_.assignOilGasDistribution(globalDofIdx, gasInPlaceLiquid, oilInPlaceGas);
908 }
909
910 template <typename FluidState, typename FIPArray>
911 void updateGasWaterDistribution(const unsigned globalDofIdx,
912 const FluidState& fs,
913 const FIPArray& fip)
914 {
915 // Gas dissolved in water and vaporized water
916 const auto gasInPlaceWater = getValue(fs.Rsw()) * fip[waterPhaseIdx];
917 const auto waterInPlaceGas = getValue(fs.Rvw()) * fip[gasPhaseIdx];
918
919 this->fipC_.assignGasWater(globalDofIdx, fip, gasInPlaceWater, waterInPlaceGas);
920 }
921
922 template <typename IntensiveQuantities>
923 void updateCO2InGas(const unsigned globalDofIdx,
924 const double pv,
925 const IntensiveQuantities& intQuants)
926 {
927 const auto& scaledDrainageInfo = this->simulator_.problem().materialLawManager()
928 ->oilWaterScaledEpsInfoDrainage(globalDofIdx);
929
930 const auto& fs = intQuants.fluidState();
931 Scalar sgcr = scaledDrainageInfo.Sgcr;
932 if (this->simulator_.problem().materialLawManager()->enableHysteresis()) {
933 const auto& matParams = simulator_.problem().materialLawParams(globalDofIdx);
934 sgcr = MaterialLaw::trappedGasSaturation(matParams, /*maximumTrapping*/false);
935 }
936
938 if (this->fipC_.has(Inplace::Phase::CO2MassInGasPhaseMaximumTrapped) ||
939 this->fipC_.has(Inplace::Phase::CO2MassInGasPhaseMaximumUnTrapped))
940 {
941 if (this->simulator_.problem().materialLawManager()->enableHysteresis()) {
942 const auto& matParams = simulator_.problem().materialLawParams(globalDofIdx);
943 // Get the maximum trapped gas saturation
944 trappedGasSaturation = MaterialLaw::trappedGasSaturation(matParams, /*maximumTrapping*/true);
945 }
946 }
947
948 const Scalar sg = getValue(fs.saturation(gasPhaseIdx));
950 if (this->fipC_.has(Inplace::Phase::CO2MassInGasPhaseEffectiveTrapped) ||
951 this->fipC_.has(Inplace::Phase::CO2MassInGasPhaseEffectiveUnTrapped))
952 {
953 if (this->simulator_.problem().materialLawManager()->enableHysteresis()) {
954 const auto& matParams = simulator_.problem().materialLawParams(globalDofIdx);
955 const double krg = getValue(intQuants.relativePermeability(gasPhaseIdx));
956 strandedGasSaturation = MaterialLaw::strandedGasSaturation(matParams, sg, krg);
957 }
958 }
959
960 const typename FIPContainer<FluidSystem>::Co2InGasInput v{
961 pv,
962 sg,
963 sgcr,
964 getValue(fs.density(gasPhaseIdx)),
965 FluidSystem::phaseIsActive(waterPhaseIdx)
966 ? FluidSystem::convertRvwToXgW(getValue(fs.Rvw()), fs.pvtRegionIndex())
967 : FluidSystem::convertRvToXgO(getValue(fs.Rv()), fs.pvtRegionIndex()),
968 FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex()),
971 };
972
973 this->fipC_.assignCo2InGas(globalDofIdx, v);
974 }
975
976 template <typename FluidState>
977 void updateCO2InWater(const unsigned globalDofIdx,
978 const double pv,
979 const FluidState& fs)
980 {
981 const auto co2InWater = FluidSystem::phaseIsActive(oilPhaseIdx)
982 ? this->co2InWaterFromOil(fs, pv)
983 : this->co2InWaterFromWater(fs, pv);
984
985 const Scalar mM = FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex());
986
987 this->fipC_.assignCo2InWater(globalDofIdx, co2InWater, mM);
988 }
989
990 template <typename FluidState>
991 Scalar co2InWaterFromWater(const FluidState& fs, const double pv) const
992 {
993 const double rhow = getValue(fs.density(waterPhaseIdx));
994 const double sw = getValue(fs.saturation(waterPhaseIdx));
995 const double xwG = FluidSystem::convertRswToXwG(getValue(fs.Rsw()), fs.pvtRegionIndex());
996
997 const Scalar mM = FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex());
998
999 return xwG * pv * rhow * sw / mM;
1000 }
1001
1002 template <typename FluidState>
1003 Scalar co2InWaterFromOil(const FluidState& fs, const double pv) const
1004 {
1005 const double rhoo = getValue(fs.density(oilPhaseIdx));
1006 const double so = getValue(fs.saturation(oilPhaseIdx));
1007 const double xoG = FluidSystem::convertRsToXoG(getValue(fs.Rs()), fs.pvtRegionIndex());
1008
1009 const Scalar mM = FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex());
1010
1011 return xoG * pv * rhoo * so / mM;
1012 }
1013
1014 template <typename IntensiveQuantities>
1015 void updateMicrobialMass(const unsigned globalDofIdx,
1016 const IntensiveQuantities& intQuants,
1017 const double surfVolWat)
1018 {
1019 const Scalar mass = surfVolWat * intQuants.microbialConcentration().value();
1020
1021 this->fipC_.assignMicrobialMass(globalDofIdx, mass);
1022 }
1023
1024 template <typename IntensiveQuantities>
1025 void updateOxygenMass(const unsigned globalDofIdx,
1026 const IntensiveQuantities& intQuants,
1027 const double surfVolWat)
1028 {
1029 const Scalar mass = surfVolWat * intQuants.oxygenConcentration().value();
1030
1031 this->fipC_.assignOxygenMass(globalDofIdx, mass);
1032 }
1033
1034 template <typename IntensiveQuantities>
1035 void updateUreaMass(const unsigned globalDofIdx,
1036 const IntensiveQuantities& intQuants,
1037 const double surfVolWat)
1038 {
1039 const Scalar mass = surfVolWat * intQuants.ureaConcentration().value();
1040
1041 this->fipC_.assignUreaMass(globalDofIdx, mass);
1042 }
1043
1044 template <typename IntensiveQuantities>
1045 void updateBiofilmMass(const unsigned globalDofIdx,
1046 const IntensiveQuantities& intQuants,
1047 const double totVolume)
1048 {
1049 const Scalar mass = totVolume * intQuants.biofilmMass().value();
1050
1051 this->fipC_.assignBiofilmMass(globalDofIdx, mass);
1052 }
1053
1054 template <typename IntensiveQuantities>
1055 void updateCalciteMass(const unsigned globalDofIdx,
1056 const IntensiveQuantities& intQuants,
1057 const double totVolume)
1058 {
1059 const Scalar mass = totVolume * intQuants.calciteMass().value();
1060
1061 this->fipC_.assignCalciteMass(globalDofIdx, mass);
1062 }
1063
1065 void setupElementExtractors_()
1066 {
1067 using Entry = typename Extractor::Entry;
1068 using Context = typename Extractor::Context;
1069 using ScalarEntry = typename Extractor::ScalarEntry;
1070 using PhaseEntry = typename Extractor::PhaseEntry;
1071
1072 const bool hasResidual = simulator_.model().linearizer().residual().size() > 0;
1073 const auto& hysteresisConfig = simulator_.problem().materialLawManager()->hysteresisConfig();
1074
1075 auto extractors = std::array{
1076 Entry{PhaseEntry{&this->saturation_,
1077 [](const unsigned phase, const Context& ectx)
1078 { return getValue(ectx.fs.saturation(phase)); }
1079 }
1080 },
1081 Entry{PhaseEntry{&this->invB_,
1082 [](const unsigned phase, const Context& ectx)
1083 { return getValue(ectx.fs.invB(phase)); }
1084 }
1085 },
1086 Entry{PhaseEntry{&this->density_,
1087 [](const unsigned phase, const Context& ectx)
1088 { return getValue(ectx.fs.density(phase)); }
1089 }
1090 },
1091 Entry{PhaseEntry{&this->relativePermeability_,
1092 [](const unsigned phase, const Context& ectx)
1093 { return getValue(ectx.intQuants.relativePermeability(phase)); }
1094 }
1095 },
1096 Entry{PhaseEntry{&this->viscosity_,
1097 [this](const unsigned phaseIdx, const Context& ectx)
1098 {
1099 if (this->extboC_.allocated() && phaseIdx == oilPhaseIdx) {
1100 return getValue(ectx.intQuants.oilViscosity());
1101 }
1102 else if (this->extboC_.allocated() && phaseIdx == gasPhaseIdx) {
1103 return getValue(ectx.intQuants.gasViscosity());
1104 }
1105 else {
1106 return getValue(ectx.fs.viscosity(phaseIdx));
1107 }
1108 }
1109 }
1110 },
1111 Entry{PhaseEntry{&this->residual_,
1112 [&modelResid = this->simulator_.model().linearizer().residual()]
1113 (const unsigned phaseIdx, const Context& ectx)
1114 {
1115 const unsigned sIdx = FluidSystem::solventComponentIndex(phaseIdx);
1116 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(sIdx);
1117 return modelResid[ectx.globalDofIdx][activeCompIdx];
1118 }
1119 },
1121 },
1122 Entry{ScalarEntry{&this->rockCompPorvMultiplier_,
1123 [&problem = this->simulator_.problem()](const Context& ectx)
1124 {
1125 return problem.template
1127 ectx.globalDofIdx);
1128 }
1129 }
1130 },
1131 Entry{ScalarEntry{&this->rockCompTransMultiplier_,
1132 [&problem = this->simulator_.problem()](const Context& ectx)
1133 {
1134 return problem.
1135 template rockCompTransMultiplier<Scalar>(ectx.intQuants,
1136 ectx.globalDofIdx);
1137 }}
1138 },
1139 Entry{ScalarEntry{&this->minimumOilPressure_,
1140 [&problem = this->simulator_.problem()](const Context& ectx)
1141 {
1142 return std::min(getValue(ectx.fs.pressure(oilPhaseIdx)),
1143 problem.minOilPressure(ectx.globalDofIdx));
1144 }
1145 }
1146 },
1147 Entry{ScalarEntry{&this->bubblePointPressure_,
1148 [&failedCells = this->failedCellsPb_,
1149 &vanguard = this->simulator_.vanguard()](const Context& ectx)
1150 {
1151 try {
1152 return getValue(
1153 FluidSystem::bubblePointPressure(ectx.fs,
1154 ectx.intQuants.pvtRegionIndex())
1155 );
1156 } catch (const NumericalProblem&) {
1157 const auto cartesianIdx = vanguard.cartesianIndex(ectx.globalDofIdx);
1158 failedCells.push_back(cartesianIdx);
1159 return Scalar{0};
1160 }
1161 }
1162 }
1163 },
1164 Entry{ScalarEntry{&this->dewPointPressure_,
1165 [&failedCells = this->failedCellsPd_,
1166 &vanguard = this->simulator_.vanguard()](const Context& ectx)
1167 {
1168 try {
1169 return getValue(
1170 FluidSystem::dewPointPressure(ectx.fs,
1171 ectx.intQuants.pvtRegionIndex())
1172 );
1173 } catch (const NumericalProblem&) {
1174 const auto cartesianIdx = vanguard.cartesianIndex(ectx.globalDofIdx);
1175 failedCells.push_back(cartesianIdx);
1176 return Scalar{0};
1177 }
1178 }
1179 }
1180 },
1181 Entry{ScalarEntry{&this->overburdenPressure_,
1182 [&problem = simulator_.problem()](const Context& ectx)
1183 { return problem.overburdenPressure(ectx.globalDofIdx); }
1184 }
1185 },
1186 Entry{ScalarEntry{&this->temperature_,
1187 [](const Context& ectx)
1188 { return getValue(ectx.fs.temperature(oilPhaseIdx)); }
1189 }
1190 },
1191 Entry{ScalarEntry{&this->sSol_,
1192 [](const Context& ectx)
1193 { return getValue(ectx.intQuants.solventSaturation()); }
1194 }
1195 },
1196 Entry{ScalarEntry{&this->rswSol_,
1197 [](const Context& ectx)
1198 { return getValue(ectx.intQuants.rsSolw()); }
1199 }
1200 },
1201 Entry{ScalarEntry{&this->cPolymer_,
1202 [](const Context& ectx)
1203 { return getValue(ectx.intQuants.polymerConcentration()); }
1204 }
1205 },
1206 Entry{ScalarEntry{&this->cFoam_,
1207 [](const Context& ectx)
1208 { return getValue(ectx.intQuants.foamConcentration()); }
1209 }
1210 },
1211 Entry{ScalarEntry{&this->cSalt_,
1212 [](const Context& ectx)
1213 { return getValue(ectx.fs.saltConcentration()); }
1214 }
1215 },
1216 Entry{ScalarEntry{&this->pSalt_,
1217 [](const Context& ectx)
1218 { return getValue(ectx.fs.saltSaturation()); }
1219 }
1220 },
1221 Entry{ScalarEntry{&this->permFact_,
1222 [](const Context& ectx)
1223 { return getValue(ectx.intQuants.permFactor()); }
1224 }
1225 },
1226 Entry{ScalarEntry{&this->rPorV_,
1227 [&model = this->simulator_.model()](const Context& ectx)
1228 {
1229 const auto totVolume = model.dofTotalVolume(ectx.globalDofIdx);
1230 return totVolume * getValue(ectx.intQuants.porosity());
1231 }
1232 }
1233 },
1234 Entry{ScalarEntry{&this->rs_,
1235 [](const Context& ectx)
1236 { return getValue(ectx.fs.Rs()); }
1237 }
1238 },
1239 Entry{ScalarEntry{&this->rv_,
1240 [](const Context& ectx)
1241 { return getValue(ectx.fs.Rv()); }
1242 }
1243 },
1244 Entry{ScalarEntry{&this->rsw_,
1245 [](const Context& ectx)
1246 { return getValue(ectx.fs.Rsw()); }
1247 }
1248 },
1249 Entry{ScalarEntry{&this->rvw_,
1250 [](const Context& ectx)
1251 { return getValue(ectx.fs.Rvw()); }
1252 }
1253 },
1254 Entry{ScalarEntry{&this->ppcw_,
1255 [&matLawManager = *this->simulator_.problem().materialLawManager()]
1256 (const Context& ectx)
1257 {
1258 return matLawManager.
1259 oilWaterScaledEpsInfoDrainage(ectx.globalDofIdx).maxPcow;
1260 }
1261 }
1262 },
1263 Entry{ScalarEntry{&this->drsdtcon_,
1264 [&problem = this->simulator_.problem()](const Context& ectx)
1265 {
1266 return problem.drsdtcon(ectx.globalDofIdx,
1267 ectx.episodeIndex);
1268 }
1269 }
1270 },
1271 Entry{ScalarEntry{&this->pcgw_,
1272 [](const Context& ectx)
1273 {
1274 return getValue(ectx.fs.pressure(gasPhaseIdx)) -
1275 getValue(ectx.fs.pressure(waterPhaseIdx));
1276 }
1277 }
1278 },
1279 Entry{ScalarEntry{&this->pcow_,
1280 [](const Context& ectx)
1281 {
1282 return getValue(ectx.fs.pressure(oilPhaseIdx)) -
1283 getValue(ectx.fs.pressure(waterPhaseIdx));
1284 }
1285 }
1286 },
1287 Entry{ScalarEntry{&this->pcog_,
1288 [](const Context& ectx)
1289 {
1290 return getValue(ectx.fs.pressure(gasPhaseIdx)) -
1291 getValue(ectx.fs.pressure(oilPhaseIdx));
1292 }
1293 }
1294 },
1295 Entry{ScalarEntry{&this->fluidPressure_,
1296 [](const Context& ectx)
1297 {
1298 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1299 // Output oil pressure as default
1300 return getValue(ectx.fs.pressure(oilPhaseIdx));
1301 }
1302 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1303 // Output gas if oil is not present
1304 return getValue(ectx.fs.pressure(gasPhaseIdx));
1305 }
1306 else {
1307 // Output water if neither oil nor gas is present
1308 return getValue(ectx.fs.pressure(waterPhaseIdx));
1309 }
1310 }
1311 }
1312 },
1313 Entry{ScalarEntry{&this->gasDissolutionFactor_,
1314 [&problem = this->simulator_.problem()](const Context& ectx)
1315 {
1316 const Scalar SoMax = problem.maxOilSaturation(ectx.globalDofIdx);
1317 return FluidSystem::template
1319 oilPhaseIdx,
1320 ectx.pvtRegionIdx,
1321 SoMax);
1322 }
1323 }
1324 },
1325 Entry{ScalarEntry{&this->oilVaporizationFactor_,
1326 [&problem = this->simulator_.problem()](const Context& ectx)
1327 {
1328 const Scalar SoMax = problem.maxOilSaturation(ectx.globalDofIdx);
1329 return FluidSystem::template
1331 gasPhaseIdx,
1332 ectx.pvtRegionIdx,
1333 SoMax);
1334 }
1335 }
1336 },
1337 Entry{ScalarEntry{&this->gasDissolutionFactorInWater_,
1338 [&problem = this->simulator_.problem()](const Context& ectx)
1339 {
1340 const Scalar SwMax = problem.maxWaterSaturation(ectx.globalDofIdx);
1341 return FluidSystem::template
1343 waterPhaseIdx,
1344 ectx.pvtRegionIdx,
1345 SwMax);
1346 }
1347 }
1348 },
1349 Entry{ScalarEntry{&this->waterVaporizationFactor_,
1350 [](const Context& ectx)
1351 {
1352 return FluidSystem::template
1354 gasPhaseIdx,
1355 ectx.pvtRegionIdx);
1356 }
1357 }
1358 },
1359 Entry{ScalarEntry{&this->gasFormationVolumeFactor_,
1360 [](const Context& ectx)
1361 {
1362 return 1.0 / FluidSystem::template
1364 gasPhaseIdx,
1365 ectx.pvtRegionIdx);
1366 }
1367 }
1368 },
1369 Entry{ScalarEntry{&this->saturatedOilFormationVolumeFactor_,
1370 [](const Context& ectx)
1371 {
1372 return 1.0 / FluidSystem::template
1374 oilPhaseIdx,
1375 ectx.pvtRegionIdx);
1376 }
1377 }
1378 },
1379 Entry{ScalarEntry{&this->oilSaturationPressure_,
1380 [](const Context& ectx)
1381 {
1382 return FluidSystem::template
1384 oilPhaseIdx,
1385 ectx.pvtRegionIdx);
1386 }
1387 }
1388 },
1389 Entry{ScalarEntry{&this->soMax_,
1390 [&problem = this->simulator_.problem()](const Context& ectx)
1391 {
1392 return std::max(getValue(ectx.fs.saturation(oilPhaseIdx)),
1393 problem.maxOilSaturation(ectx.globalDofIdx));
1394 }
1395 },
1396 !hysteresisConfig.enableHysteresis()
1397 },
1398 Entry{ScalarEntry{&this->swMax_,
1399 [&problem = this->simulator_.problem()](const Context& ectx)
1400 {
1401 return std::max(getValue(ectx.fs.saturation(waterPhaseIdx)),
1402 problem.maxWaterSaturation(ectx.globalDofIdx));
1403 }
1404 },
1405 !hysteresisConfig.enableHysteresis()
1406 },
1407 Entry{ScalarEntry{&this->soMax_,
1408 [](const Context& ectx)
1409 { return ectx.hParams.somax; }
1410 },
1411 hysteresisConfig.enableHysteresis() &&
1412 hysteresisConfig.enableNonWettingHysteresis() &&
1413 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1414 FluidSystem::phaseIsActive(waterPhaseIdx)
1415 },
1416 Entry{ScalarEntry{&this->swMax_,
1417 [](const Context& ectx)
1418 { return ectx.hParams.swmax; }
1419 },
1420 hysteresisConfig.enableHysteresis() &&
1421 hysteresisConfig.enableWettingHysteresis() &&
1422 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1423 FluidSystem::phaseIsActive(waterPhaseIdx)
1424 },
1425 Entry{ScalarEntry{&this->swmin_,
1426 [](const Context& ectx)
1427 { return ectx.hParams.swmin; }
1428 },
1429 hysteresisConfig.enableHysteresis() &&
1430 hysteresisConfig.enablePCHysteresis() &&
1431 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1432 FluidSystem::phaseIsActive(waterPhaseIdx)
1433 },
1434 Entry{ScalarEntry{&this->sgmax_,
1435 [](const Context& ectx)
1436 { return ectx.hParams.sgmax; }
1437 },
1438 hysteresisConfig.enableHysteresis() &&
1439 hysteresisConfig.enableNonWettingHysteresis() &&
1440 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1441 FluidSystem::phaseIsActive(gasPhaseIdx)
1442 },
1443 Entry{ScalarEntry{&this->shmax_,
1444 [](const Context& ectx)
1445 { return ectx.hParams.shmax; }
1446 },
1447 hysteresisConfig.enableHysteresis() &&
1448 hysteresisConfig.enableWettingHysteresis() &&
1449 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1450 FluidSystem::phaseIsActive(gasPhaseIdx)
1451 },
1452 Entry{ScalarEntry{&this->somin_,
1453 [](const Context& ectx)
1454 { return ectx.hParams.somin; }
1455 },
1456 hysteresisConfig.enableHysteresis() &&
1457 hysteresisConfig.enablePCHysteresis() &&
1458 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1459 FluidSystem::phaseIsActive(gasPhaseIdx)
1460 },
1461 Entry{[&model = this->simulator_.model(), this](const Context& ectx)
1462 {
1463 // Note: We intentionally exclude effects of rock
1464 // compressibility by using referencePorosity() here.
1465 const auto porv = ectx.intQuants.referencePorosity()
1466 * model.dofTotalVolume(ectx.globalDofIdx);
1467
1468 this->aggregateAverageDensityContributions_(ectx.fs, ectx.globalDofIdx,
1469 static_cast<double>(porv));
1470 }, this->regionAvgDensity_.has_value()
1471 },
1472 Entry{[&extboC = this->extboC_](const Context& ectx)
1473 {
1474 extboC.assignVolumes(ectx.globalDofIdx,
1475 ectx.intQuants.xVolume().value(),
1476 ectx.intQuants.yVolume().value());
1477 extboC.assignZFraction(ectx.globalDofIdx,
1478 ectx.intQuants.zFraction().value());
1479
1480 const Scalar stdVolOil = getValue(ectx.fs.saturation(oilPhaseIdx)) *
1481 getValue(ectx.fs.invB(oilPhaseIdx)) +
1482 getValue(ectx.fs.saturation(gasPhaseIdx)) *
1483 getValue(ectx.fs.invB(gasPhaseIdx)) *
1484 getValue(ectx.fs.Rv());
1485 const Scalar stdVolGas = getValue(ectx.fs.saturation(gasPhaseIdx)) *
1486 getValue(ectx.fs.invB(gasPhaseIdx)) *
1487 (1.0 - ectx.intQuants.yVolume().value()) +
1488 getValue(ectx.fs.saturation(oilPhaseIdx)) *
1489 getValue(ectx.fs.invB(oilPhaseIdx)) *
1490 getValue(ectx.fs.Rs()) *
1491 (1.0 - ectx.intQuants.xVolume().value());
1492 const Scalar stdVolCo2 = getValue(ectx.fs.saturation(gasPhaseIdx)) *
1493 getValue(ectx.fs.invB(gasPhaseIdx)) *
1494 ectx.intQuants.yVolume().value() +
1495 getValue(ectx.fs.saturation(oilPhaseIdx)) *
1496 getValue(ectx.fs.invB(oilPhaseIdx)) *
1497 getValue(ectx.fs.Rs()) *
1498 ectx.intQuants.xVolume().value();
1499 const Scalar rhoO = FluidSystem::referenceDensity(gasPhaseIdx, ectx.pvtRegionIdx);
1500 const Scalar rhoG = FluidSystem::referenceDensity(gasPhaseIdx, ectx.pvtRegionIdx);
1501 const Scalar rhoCO2 = ectx.intQuants.zRefDensity();
1502 const Scalar stdMassTotal = 1.0e-10 + stdVolOil * rhoO + stdVolGas * rhoG + stdVolCo2 * rhoCO2;
1503 extboC.assignMassFractions(ectx.globalDofIdx,
1507 }, this->extboC_.allocated()
1508 },
1509 Entry{[&micpC = this->micpC_](const Context& ectx)
1510 {
1511 micpC.assign(ectx.globalDofIdx,
1512 ectx.intQuants.microbialConcentration().value(),
1513 ectx.intQuants.oxygenConcentration().value(),
1514 ectx.intQuants.ureaConcentration().value(),
1515 ectx.intQuants.biofilmConcentration().value(),
1516 ectx.intQuants.calciteConcentration().value());
1517 }, this->micpC_.allocated()
1518 },
1519 Entry{[&rftC = this->rftC_,
1520 &vanguard = this->simulator_.vanguard()](const Context& ectx)
1521 {
1522 const auto cartesianIdx = vanguard.cartesianIndex(ectx.globalDofIdx);
1523 rftC.assign(cartesianIdx,
1524 [&fs = ectx.fs]() { return getValue(fs.pressure(oilPhaseIdx)); },
1525 [&fs = ectx.fs]() { return getValue(fs.saturation(waterPhaseIdx)); },
1526 [&fs = ectx.fs]() { return getValue(fs.saturation(gasPhaseIdx)); });
1527 }
1528 },
1529 Entry{[&tC = this->tracerC_,
1530 &tM = this->simulator_.problem().tracerModel()](const Context& ectx)
1531 {
1532 tC.assignFreeConcentrations(ectx.globalDofIdx,
1533 [gIdx = ectx.globalDofIdx, &tM](const unsigned tracerIdx)
1534 { return tM.freeTracerConcentration(tracerIdx, gIdx); });
1535 tC.assignSolConcentrations(ectx.globalDofIdx,
1536 [gIdx = ectx.globalDofIdx, &tM](const unsigned tracerIdx)
1537 { return tM.solTracerConcentration(tracerIdx, gIdx); });
1538 }
1539 },
1540 Entry{[&flowsInf = this->simulator_.problem().model().linearizer().getFlowsInfo(),
1541 &flowsC = this->flowsC_](const Context& ectx)
1542 {
1543 constexpr auto gas_idx = Indices::gasEnabled ?
1544 conti0EqIdx + Indices::canonicalToActiveComponentIndex(gasCompIdx) : -1;
1545 constexpr auto oil_idx = Indices::oilEnabled ?
1546 conti0EqIdx + Indices::canonicalToActiveComponentIndex(oilCompIdx) : -1;
1547 constexpr auto water_idx = Indices::waterEnabled ?
1548 conti0EqIdx + Indices::canonicalToActiveComponentIndex(waterCompIdx) : -1;
1549 const auto& flowsInfos = flowsInf[ectx.globalDofIdx];
1550 for (const auto& flowsInfo : flowsInfos) {
1551 flowsC.assignFlows(ectx.globalDofIdx,
1552 flowsInfo.faceId,
1553 flowsInfo.nncId,
1557 }
1558 }, !this->simulator_.problem().model().linearizer().getFlowsInfo().empty()
1559 },
1560 Entry{[&floresInf = this->simulator_.problem().model().linearizer().getFloresInfo(),
1561 &flowsC = this->flowsC_](const Context& ectx)
1562 {
1563 constexpr auto gas_idx = Indices::gasEnabled ?
1564 conti0EqIdx + Indices::canonicalToActiveComponentIndex(gasCompIdx) : -1;
1565 constexpr auto oil_idx = Indices::oilEnabled ?
1566 conti0EqIdx + Indices::canonicalToActiveComponentIndex(oilCompIdx) : -1;
1567 constexpr auto water_idx = Indices::waterEnabled ?
1568 conti0EqIdx + Indices::canonicalToActiveComponentIndex(waterCompIdx) : -1;
1569 const auto& floresInfos = floresInf[ectx.globalDofIdx];
1570 for (const auto& floresInfo : floresInfos) {
1571 flowsC.assignFlores(ectx.globalDofIdx,
1572 floresInfo.faceId,
1573 floresInfo.nncId,
1577 }
1578 }, !this->simulator_.problem().model().linearizer().getFloresInfo().empty()
1579 },
1580 // hack to make the intial output of rs and rv Ecl compatible.
1581 // For cells with swat == 1 Ecl outputs; rs = rsSat and rv=rvSat, in all but the initial step
1582 // where it outputs rs and rv values calculated by the initialization. To be compatible we overwrite
1583 // rs and rv with the values computed in the initially.
1584 // Volume factors, densities and viscosities need to be recalculated with the updated rs and rv values.
1585 Entry{ScalarEntry{&this->rv_,
1586 [&problem = this->simulator_.problem()](const Context& ectx)
1587 { return problem.initialFluidState(ectx.globalDofIdx).Rv(); }
1588 },
1589 simulator_.episodeIndex() < 0 &&
1590 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1591 FluidSystem::phaseIsActive(gasPhaseIdx)
1592 },
1593 Entry{ScalarEntry{&this->rs_,
1594 [&problem = this->simulator_.problem()](const Context& ectx)
1595 { return problem.initialFluidState(ectx.globalDofIdx).Rs(); }
1596 },
1597 simulator_.episodeIndex() < 0 &&
1598 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1599 FluidSystem::phaseIsActive(gasPhaseIdx)
1600 },
1601 Entry{ScalarEntry{&this->rsw_,
1602 [&problem = this->simulator_.problem()](const Context& ectx)
1603 { return problem.initialFluidState(ectx.globalDofIdx).Rsw(); }
1604 },
1605 simulator_.episodeIndex() < 0 &&
1606 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1607 FluidSystem::phaseIsActive(gasPhaseIdx)
1608 },
1609 Entry{ScalarEntry{&this->rvw_,
1610 [&problem = this->simulator_.problem()](const Context& ectx)
1611 { return problem.initialFluidState(ectx.globalDofIdx).Rvw(); }
1612 },
1613 simulator_.episodeIndex() < 0 &&
1614 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1615 FluidSystem::phaseIsActive(gasPhaseIdx)
1616 },
1617 // re-compute the volume factors, viscosities and densities if asked for
1618 Entry{PhaseEntry{&this->density_,
1619 [&problem = this->simulator_.problem()](const unsigned phase,
1620 const Context& ectx)
1621 {
1622 const auto& fsInitial = problem.initialFluidState(ectx.globalDofIdx);
1623 return FluidSystem::density(fsInitial,
1624 phase,
1625 ectx.intQuants.pvtRegionIndex());
1626 }
1627 },
1628 simulator_.episodeIndex() < 0 &&
1629 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1630 FluidSystem::phaseIsActive(gasPhaseIdx)
1631 },
1632 Entry{PhaseEntry{&this->invB_,
1633 [&problem = this->simulator_.problem()](const unsigned phase,
1634 const Context& ectx)
1635 {
1636 const auto& fsInitial = problem.initialFluidState(ectx.globalDofIdx);
1637 return FluidSystem::inverseFormationVolumeFactor(fsInitial,
1638 phase,
1639 ectx.intQuants.pvtRegionIndex());
1640 }
1641 },
1642 simulator_.episodeIndex() < 0 &&
1643 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1644 FluidSystem::phaseIsActive(gasPhaseIdx)
1645 },
1646 Entry{PhaseEntry{&this->viscosity_,
1647 [&problem = this->simulator_.problem()](const unsigned phase,
1648 const Context& ectx)
1649 {
1650 const auto& fsInitial = problem.initialFluidState(ectx.globalDofIdx);
1651 return FluidSystem::viscosity(fsInitial,
1652 phase,
1653 ectx.intQuants.pvtRegionIndex());
1654 }
1655 },
1656 simulator_.episodeIndex() < 0 &&
1657 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1658 FluidSystem::phaseIsActive(gasPhaseIdx)
1659 },
1660 };
1661
1662 // Setup active extractors
1663 this->extractors_ = Extractor::removeInactive(extractors);
1664
1666 if (this->mech_.allocated()) {
1667 this->extractors_.emplace(
1668 [&mech = this->mech_,
1669 &model = simulator_.problem().geoMechModel()](const Context& ectx)
1670 {
1671 mech.assignDelStress(ectx.globalDofIdx,
1672 model.delstress(ectx.globalDofIdx));
1673
1674 mech.assignDisplacement(ectx.globalDofIdx,
1675 model.disp(ectx.globalDofIdx, /*include_fracture*/true));
1676
1677 // is the tresagii stress which make rock fracture
1678 mech.assignFracStress(ectx.globalDofIdx,
1679 model.fractureStress(ectx.globalDofIdx));
1680
1681 mech.assignLinStress(ectx.globalDofIdx,
1682 model.linstress(ectx.globalDofIdx));
1683
1684 mech.assignPotentialForces(ectx.globalDofIdx,
1685 model.mechPotentialForce(ectx.globalDofIdx),
1686 model.mechPotentialPressForce(ectx.globalDofIdx),
1687 model.mechPotentialTempForce(ectx.globalDofIdx));
1688
1689 mech.assignStrain(ectx.globalDofIdx,
1690 model.strain(ectx.globalDofIdx, /*include_fracture*/true));
1691
1692 // Total stress is not stored but calculated result is Voigt notation
1693 mech.assignStress(ectx.globalDofIdx,
1694 model.stress(ectx.globalDofIdx, /*include_fracture*/true));
1695 }, true
1696 );
1697 }
1698 }
1699 }
1700
1702 void setupBlockExtractors_(const bool isSubStep,
1703 const int reportStepNum)
1704 {
1705 using Entry = typename BlockExtractor::Entry;
1706 using Context = typename BlockExtractor::Context;
1707 using PhaseEntry = typename BlockExtractor::PhaseEntry;
1708 using ScalarEntry = typename BlockExtractor::ScalarEntry;
1709
1710 using namespace std::string_view_literals;
1711
1712 const auto pressure_handler =
1713 Entry{ScalarEntry{std::vector{"BPR"sv, "BPRESSUR"sv},
1714 [](const Context& ectx)
1715 {
1716 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1717 return getValue(ectx.fs.pressure(oilPhaseIdx));
1718 }
1719 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1720 return getValue(ectx.fs.pressure(gasPhaseIdx));
1721 }
1722 else { //if (FluidSystem::phaseIsActive(waterPhaseIdx))
1723 return getValue(ectx.fs.pressure(waterPhaseIdx));
1724 }
1725 }
1726 }
1727 };
1728
1729 const auto handlers = std::array{
1731 Entry{PhaseEntry{std::array{
1732 std::array{"BWSAT"sv, "BOSAT"sv, "BGSAT"sv},
1733 std::array{"BSWAT"sv, "BSOIL"sv, "BSGAS"sv}
1734 },
1735 [](const unsigned phaseIdx, const Context& ectx)
1736 {
1737 return getValue(ectx.fs.saturation(phaseIdx));
1738 }
1739 }
1740 },
1741 Entry{ScalarEntry{"BNSAT",
1742 [](const Context& ectx)
1743 {
1744 return ectx.intQuants.solventSaturation().value();
1745 }
1746 }
1747 },
1748 Entry{ScalarEntry{std::vector{"BTCNFHEA"sv, "BTEMP"sv},
1749 [](const Context& ectx)
1750 {
1751 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1752 return getValue(ectx.fs.temperature(oilPhaseIdx));
1753 }
1754 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1755 return getValue(ectx.fs.temperature(gasPhaseIdx));
1756 }
1757 else { //if (FluidSystem::phaseIsActive(waterPhaseIdx))
1758 return getValue(ectx.fs.temperature(waterPhaseIdx));
1759 }
1760 }
1761 }
1762 },
1763 Entry{PhaseEntry{std::array{
1764 std::array{"BWKR"sv, "BOKR"sv, "BGKR"sv},
1765 std::array{"BKRW"sv, "BKRO"sv, "BKRG"sv}
1766 },
1767 [](const unsigned phaseIdx, const Context& ectx)
1768 {
1769 return getValue(ectx.intQuants.relativePermeability(phaseIdx));
1770 }
1771 }
1772 },
1773 Entry{ScalarEntry{"BKROG",
1774 [&problem = this->simulator_.problem()](const Context& ectx)
1775 {
1776 const auto& materialParams =
1777 problem.materialLawParams(ectx.elemCtx,
1778 ectx.dofIdx,
1779 /* timeIdx = */ 0);
1780 return getValue(MaterialLaw::template
1782 ectx.fs));
1783 }
1784 }
1785 },
1786 Entry{ScalarEntry{"BKROW",
1787 [&problem = this->simulator_.problem()](const Context& ectx)
1788 {
1789 const auto& materialParams = problem.materialLawParams(ectx.elemCtx,
1790 ectx.dofIdx,
1791 /* timeIdx = */ 0);
1792 return getValue(MaterialLaw::template
1794 ectx.fs));
1795 }
1796 }
1797 },
1798 Entry{ScalarEntry{"BWPC",
1799 [](const Context& ectx)
1800 {
1801 return getValue(ectx.fs.pressure(oilPhaseIdx)) -
1802 getValue(ectx.fs.pressure(waterPhaseIdx));
1803 }
1804 }
1805 },
1806 Entry{ScalarEntry{"BGPC",
1807 [](const Context& ectx)
1808 {
1809 return getValue(ectx.fs.pressure(gasPhaseIdx)) -
1810 getValue(ectx.fs.pressure(oilPhaseIdx));
1811 }
1812 }
1813 },
1814 Entry{ScalarEntry{"BWPR",
1815 [](const Context& ectx)
1816 {
1817 return getValue(ectx.fs.pressure(waterPhaseIdx));
1818 }
1819 }
1820 },
1821 Entry{ScalarEntry{"BGPR",
1822 [](const Context& ectx)
1823 {
1824 return getValue(ectx.fs.pressure(gasPhaseIdx));
1825 }
1826 }
1827 },
1828 Entry{PhaseEntry{std::array{
1829 std::array{"BVWAT"sv, "BVOIL"sv, "BVGAS"sv},
1830 std::array{"BWVIS"sv, "BOVIS"sv, "BGVIS"sv}
1831 },
1832 [](const unsigned phaseIdx, const Context& ectx)
1833 {
1834 return getValue(ectx.fs.viscosity(phaseIdx));
1835 }
1836 }
1837 },
1838 Entry{PhaseEntry{std::array{
1839 std::array{"BWDEN"sv, "BODEN"sv, "BGDEN"sv},
1840 std::array{"BDENW"sv, "BDENO"sv, "BDENG"sv}
1841 },
1842 [](const unsigned phaseIdx, const Context& ectx)
1843 {
1844 return getValue(ectx.fs.density(phaseIdx));
1845 }
1846 }
1847 },
1848 Entry{ScalarEntry{"BFLOWI",
1849 [&flowsC = this->flowsC_](const Context& ectx)
1850 {
1851 return flowsC.getFlow(ectx.globalDofIdx, Dir::XPlus, waterCompIdx);
1852 }
1853 }
1854 },
1855 Entry{ScalarEntry{"BFLOWJ",
1856 [&flowsC = this->flowsC_](const Context& ectx)
1857 {
1858 return flowsC.getFlow(ectx.globalDofIdx, Dir::YPlus, waterCompIdx);
1859 }
1860 }
1861 },
1862 Entry{ScalarEntry{"BFLOWK",
1863 [&flowsC = this->flowsC_](const Context& ectx)
1864 {
1865 return flowsC.getFlow(ectx.globalDofIdx, Dir::ZPlus, waterCompIdx);
1866 }
1867 }
1868 },
1869 Entry{ScalarEntry{"BRPV",
1870 [&model = this->simulator_.model()](const Context& ectx)
1871 {
1872 return getValue(ectx.intQuants.porosity()) *
1873 model.dofTotalVolume(ectx.globalDofIdx);
1874 }
1875 }
1876 },
1877 Entry{PhaseEntry{std::array{"BWPV"sv, "BOPV"sv, "BGPV"sv},
1878 [&model = this->simulator_.model()](const unsigned phaseIdx,
1879 const Context& ectx)
1880 {
1881 return getValue(ectx.fs.saturation(phaseIdx)) *
1882 getValue(ectx.intQuants.porosity()) *
1883 model.dofTotalVolume(ectx.globalDofIdx);
1884 }
1885 }
1886 },
1887 Entry{ScalarEntry{"BRS",
1888 [](const Context& ectx)
1889 {
1890 return getValue(ectx.fs.Rs());
1891 }
1892 }
1893 },
1894 Entry{ScalarEntry{"BRV",
1895 [](const Context& ectx)
1896 {
1897 return getValue(ectx.fs.Rv());
1898 }
1899 }
1900 },
1901 Entry{ScalarEntry{"BOIP",
1902 [&model = this->simulator_.model()](const Context& ectx)
1903 {
1904 return (getValue(ectx.fs.invB(oilPhaseIdx)) *
1905 getValue(ectx.fs.saturation(oilPhaseIdx)) +
1906 getValue(ectx.fs.Rv()) *
1907 getValue(ectx.fs.invB(gasPhaseIdx)) *
1908 getValue(ectx.fs.saturation(gasPhaseIdx))) *
1909 model.dofTotalVolume(ectx.globalDofIdx) *
1910 getValue(ectx.intQuants.porosity());
1911 }
1912 }
1913 },
1914 Entry{ScalarEntry{"BGIP",
1915 [&model = this->simulator_.model()](const Context& ectx)
1916 {
1917 Scalar result = getValue(ectx.fs.invB(gasPhaseIdx)) *
1918 getValue(ectx.fs.saturation(gasPhaseIdx));
1919
1920 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1921 result += getValue(ectx.fs.Rs()) *
1922 getValue(ectx.fs.invB(oilPhaseIdx)) *
1923 getValue(ectx.fs.saturation(oilPhaseIdx));
1924 }
1925 else {
1926 result += getValue(ectx.fs.Rsw()) *
1927 getValue(ectx.fs.invB(waterPhaseIdx)) *
1928 getValue(ectx.fs.saturation(waterPhaseIdx));
1929 }
1930
1931 return result *
1932 model.dofTotalVolume(ectx.globalDofIdx) *
1933 getValue(ectx.intQuants.porosity());
1934 }
1935 }
1936 },
1937 Entry{ScalarEntry{"BWIP",
1938 [&model = this->simulator_.model()](const Context& ectx)
1939 {
1940 return getValue(ectx.fs.invB(waterPhaseIdx)) *
1941 getValue(ectx.fs.saturation(waterPhaseIdx)) *
1942 model.dofTotalVolume(ectx.globalDofIdx) *
1943 getValue(ectx.intQuants.porosity());
1944 }
1945 }
1946 },
1947 Entry{ScalarEntry{"BOIPL",
1948 [&model = this->simulator_.model()](const Context& ectx)
1949 {
1950 return getValue(ectx.fs.invB(oilPhaseIdx)) *
1951 getValue(ectx.fs.saturation(oilPhaseIdx)) *
1952 model.dofTotalVolume(ectx.globalDofIdx) *
1953 getValue(ectx.intQuants.porosity());
1954 }
1955 }
1956 },
1957 Entry{ScalarEntry{"BGIPL",
1958 [&model = this->simulator_.model()](const Context& ectx)
1959 {
1960 Scalar result;
1961 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1962 result = getValue(ectx.fs.Rs()) *
1963 getValue(ectx.fs.invB(oilPhaseIdx)) *
1964 getValue(ectx.fs.saturation(oilPhaseIdx));
1965 }
1966 else {
1967 result = getValue(ectx.fs.Rsw()) *
1968 getValue(ectx.fs.invB(waterPhaseIdx)) *
1969 getValue(ectx.fs.saturation(waterPhaseIdx));
1970 }
1971 return result *
1972 model.dofTotalVolume(ectx.globalDofIdx) *
1973 getValue(ectx.intQuants.porosity());
1974 }
1975 }
1976 },
1977 Entry{ScalarEntry{"BGIPG",
1978 [&model = this->simulator_.model()](const Context& ectx)
1979 {
1980 return getValue(ectx.fs.invB(gasPhaseIdx)) *
1981 getValue(ectx.fs.saturation(gasPhaseIdx)) *
1982 model.dofTotalVolume(ectx.globalDofIdx) *
1983 getValue(ectx.intQuants.porosity());
1984 }
1985 }
1986 },
1987 Entry{ScalarEntry{"BOIPG",
1988 [&model = this->simulator_.model()](const Context& ectx)
1989 {
1990 return getValue(ectx.fs.Rv()) *
1991 getValue(ectx.fs.invB(gasPhaseIdx)) *
1992 getValue(ectx.fs.saturation(gasPhaseIdx)) *
1993 model.dofTotalVolume(ectx.globalDofIdx) *
1994 getValue(ectx.intQuants.porosity());
1995 }
1996 }
1997 },
1998 Entry{PhaseEntry{std::array{"BPPW"sv, "BPPO"sv, "BPPG"sv},
1999 [&simConfig = this->eclState_.getSimulationConfig(),
2000 &grav = this->simulator_.problem().gravity(),
2001 &regionAvgDensity = this->regionAvgDensity_,
2002 &problem = this->simulator_.problem(),
2003 &regions = this->regions_](const unsigned phaseIdx, const Context& ectx)
2004 {
2005 auto phase = RegionPhasePoreVolAverage::Phase{};
2006 phase.ix = phaseIdx;
2007
2008 // Note different region handling here. FIPNUM is
2009 // one-based, but we need zero-based lookup in
2010 // DatumDepth. On the other hand, pvtRegionIndex is
2011 // zero-based but we need one-based lookup in
2012 // RegionPhasePoreVolAverage.
2013
2014 // Subtract one to convert FIPNUM to region index.
2015 const auto datum = simConfig.datumDepths()(regions["FIPNUM"][ectx.dofIdx] - 1);
2016
2017 // Add one to convert region index to region ID.
2018 const auto region = RegionPhasePoreVolAverage::Region {
2019 ectx.elemCtx.primaryVars(ectx.dofIdx, /*timeIdx=*/0).pvtRegionIndex() + 1
2020 };
2021
2022 const auto density = regionAvgDensity->value("PVTNUM", phase, region);
2023
2024 const auto press = getValue(ectx.fs.pressure(phase.ix));
2025 const auto dz = problem.dofCenterDepth(ectx.globalDofIdx) - datum;
2026 return press - density*dz*grav[GridView::dimensionworld - 1];
2027 }
2028 }
2029 },
2030 Entry{ScalarEntry{"BMMIP",
2031 [&model = this->simulator_.model()](const Context& ectx)
2032 {
2033 return getValue(ectx.intQuants.microbialConcentration()) *
2034 getValue(ectx.intQuants.porosity()) *
2035 model.dofTotalVolume(ectx.globalDofIdx);
2036 }
2037 }
2038 },
2039 Entry{ScalarEntry{"BMOIP",
2040 [&model = this->simulator_.model()](const Context& ectx)
2041 {
2042 return getValue(ectx.intQuants.oxygenConcentration()) *
2043 getValue(ectx.intQuants.porosity()) *
2044 model.dofTotalVolume(ectx.globalDofIdx);
2045 }
2046 }
2047 },
2048 Entry{ScalarEntry{"BMUIP",
2049 [&model = this->simulator_.model()](const Context& ectx)
2050 {
2051 return getValue(ectx.intQuants.ureaConcentration()) *
2052 getValue(ectx.intQuants.porosity()) *
2053 model.dofTotalVolume(ectx.globalDofIdx) * 1;
2054 }
2055 }
2056 },
2057 Entry{ScalarEntry{"BMBIP",
2058 [&model = this->simulator_.model()](const Context& ectx)
2059 {
2060 return model.dofTotalVolume(ectx.globalDofIdx) *
2061 getValue(ectx.intQuants.biofilmMass());
2062 }
2063 }
2064 },
2065 Entry{ScalarEntry{"BMCIP",
2066 [&model = this->simulator_.model()](const Context& ectx)
2067 {
2068 return model.dofTotalVolume(ectx.globalDofIdx) *
2069 getValue(ectx.intQuants.calciteMass());
2070 }
2071 }
2072 },
2073 Entry{ScalarEntry{"BGMIP",
2074 [&model = this->simulator_.model()](const Context& ectx)
2075 {
2076 Scalar result = getValue(ectx.fs.invB(gasPhaseIdx)) *
2077 getValue(ectx.fs.saturation(gasPhaseIdx));
2078
2079 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
2080 result += getValue(ectx.fs.Rs()) *
2081 getValue(ectx.fs.invB(oilPhaseIdx)) *
2082 getValue(ectx.fs.saturation(oilPhaseIdx));
2083 }
2084 else {
2085 result += getValue(ectx.fs.Rsw()) *
2086 getValue(ectx.fs.invB(waterPhaseIdx)) *
2087 getValue(ectx.fs.saturation(waterPhaseIdx));
2088 }
2089 const Scalar rhoG = FluidSystem::referenceDensity(gasPhaseIdx,
2090 ectx.intQuants.pvtRegionIndex());
2091 return result *
2092 model.dofTotalVolume(ectx.globalDofIdx) *
2093 getValue(ectx.intQuants.porosity()) *
2094 rhoG;
2095 }
2096 }
2097 },
2098 Entry{ScalarEntry{"BGMGP",
2099 [&model = this->simulator_.model()](const Context& ectx)
2100 {
2101 const Scalar rhoG = FluidSystem::referenceDensity(gasPhaseIdx,
2102 ectx.intQuants.pvtRegionIndex());
2103 return getValue(ectx.fs.invB(gasPhaseIdx)) *
2104 getValue(ectx.fs.saturation(gasPhaseIdx)) *
2105 model.dofTotalVolume(ectx.globalDofIdx) *
2106 getValue(ectx.intQuants.porosity()) *
2107 rhoG;
2108 }
2109 }
2110 },
2111 Entry{ScalarEntry{"BGMDS",
2112 [&model = this->simulator_.model()](const Context& ectx)
2113 {
2114 Scalar result;
2115 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
2116 result = getValue(ectx.fs.Rs()) *
2117 getValue(ectx.fs.invB(oilPhaseIdx)) *
2118 getValue(ectx.fs.saturation(oilPhaseIdx));
2119 }
2120 else {
2121 result = getValue(ectx.fs.Rsw()) *
2122 getValue(ectx.fs.invB(waterPhaseIdx)) *
2123 getValue(ectx.fs.saturation(waterPhaseIdx));
2124 }
2125 const Scalar rhoG = FluidSystem::referenceDensity(gasPhaseIdx,
2126 ectx.intQuants.pvtRegionIndex());
2127 return result *
2128 model.dofTotalVolume(ectx.globalDofIdx) *
2129 getValue(ectx.intQuants.porosity()) *
2130 rhoG;
2131 }
2132 }
2133 },
2134 Entry{ScalarEntry{"BGMST",
2135 [&model = this->simulator_.model(),
2136 &problem = this->simulator_.problem()](const Context& ectx)
2137 {
2138 const auto& scaledDrainageInfo = problem.materialLawManager()
2139 ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2140 const Scalar sg = getValue(ectx.fs.saturation(gasPhaseIdx));
2141 Scalar strandedGas = scaledDrainageInfo.Sgcr;
2142 if (problem.materialLawManager()->enableHysteresis()) {
2143 const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2144 const Scalar krg = getValue(ectx.intQuants.relativePermeability(gasPhaseIdx));
2145 strandedGas = MaterialLaw::strandedGasSaturation(matParams, sg, krg);
2146 }
2147 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2148 FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2149 : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2150 return (1.0 - xgW) *
2151 model.dofTotalVolume(ectx.globalDofIdx) *
2152 getValue(ectx.intQuants.porosity()) *
2153 getValue(ectx.fs.density(gasPhaseIdx)) *
2154 std::min(strandedGas, sg);
2155 }
2156 }
2157 },
2158 Entry{ScalarEntry{"BGMUS",
2159 [&model = this->simulator_.model(),
2160 &problem = this->simulator_.problem()](const Context& ectx)
2161 {
2162 const auto& scaledDrainageInfo = problem.materialLawManager()
2163 ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2164 const Scalar sg = getValue(ectx.fs.saturation(gasPhaseIdx));
2165 Scalar strandedGas = scaledDrainageInfo.Sgcr;
2166 if (problem.materialLawManager()->enableHysteresis()) {
2167 const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2168 const Scalar krg = getValue(ectx.intQuants.relativePermeability(gasPhaseIdx));
2169 strandedGas = MaterialLaw::strandedGasSaturation(matParams, sg, krg);
2170 }
2171 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2172 FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2173 : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2174 return (1.0 - xgW) *
2175 model.dofTotalVolume(ectx.globalDofIdx) *
2176 getValue(ectx.intQuants.porosity()) *
2177 getValue(ectx.fs.density(gasPhaseIdx)) *
2178 std::max(Scalar{0.0}, sg - strandedGas);
2179 }
2180 }
2181 },
2182 Entry{ScalarEntry{"BGMTR",
2183 [&model = this->simulator_.model(),
2184 &problem = this->simulator_.problem()](const Context& ectx)
2185 {
2186 const auto& scaledDrainageInfo = problem.materialLawManager()
2187 ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2188 Scalar trappedGas = scaledDrainageInfo.Sgcr;
2189 if (problem.materialLawManager()->enableHysteresis()) {
2190 const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2191 trappedGas = MaterialLaw::trappedGasSaturation(matParams, /*maxTrapping*/true);
2192 }
2193 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2194 FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2195 : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2196 return (1.0 - xgW) *
2197 model.dofTotalVolume(ectx.globalDofIdx) *
2198 getValue(ectx.intQuants.porosity()) *
2199 getValue(ectx.fs.density(gasPhaseIdx)) *
2200 std::min(trappedGas, getValue(ectx.fs.saturation(gasPhaseIdx)));
2201 }
2202 }
2203 },
2204 Entry{ScalarEntry{"BGMMO",
2205 [&model = this->simulator_.model(),
2206 &problem = this->simulator_.problem()](const Context& ectx)
2207 {
2208 const auto& scaledDrainageInfo = problem.materialLawManager()
2209 ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2210 Scalar trappedGas = scaledDrainageInfo.Sgcr;
2211 if (problem.materialLawManager()->enableHysteresis()) {
2212 const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2213 trappedGas = MaterialLaw::trappedGasSaturation(matParams, /*maxTrapping*/true);
2214 }
2215 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2216 FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2217 : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2218 return (1.0 - xgW) *
2219 model.dofTotalVolume(ectx.globalDofIdx) *
2220 getValue(ectx.intQuants.porosity()) *
2221 getValue(ectx.fs.density(gasPhaseIdx)) *
2222 std::max(Scalar{0.0}, getValue(ectx.fs.saturation(gasPhaseIdx)) - trappedGas);
2223 }
2224 }
2225 },
2226 Entry{ScalarEntry{"BGKTR",
2227 [&model = this->simulator_.model(),
2228 &problem = this->simulator_.problem()](const Context& ectx)
2229 {
2230 const auto& scaledDrainageInfo = problem.materialLawManager()
2231 ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2232 const Scalar sg = getValue(ectx.fs.saturation(gasPhaseIdx));
2233 Scalar sgcr = scaledDrainageInfo.Sgcr;
2234 if (problem.materialLawManager()->enableHysteresis()) {
2235 const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2236 sgcr = MaterialLaw::trappedGasSaturation(matParams, /*maxTrapping*/false);
2237 }
2238 if (sg > sgcr) {
2239 return 0.0;
2240 }
2241 else {
2242 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2243 FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2244 : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2245 return (1.0 - xgW) *
2246 model.dofTotalVolume(ectx.globalDofIdx) *
2247 getValue(ectx.intQuants.porosity()) *
2248 getValue(ectx.fs.density(gasPhaseIdx)) *
2249 getValue(ectx.fs.saturation(gasPhaseIdx));
2250 }
2251 }
2252 }
2253 },
2254 Entry{ScalarEntry{"BGKMO",
2255 [&model = this->simulator_.model(),
2256 &problem = this->simulator_.problem()](const Context& ectx)
2257 {
2258 const auto& scaledDrainageInfo = problem.materialLawManager()
2259 ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2260 const Scalar sg = getValue(ectx.fs.saturation(gasPhaseIdx));
2261 Scalar sgcr = scaledDrainageInfo.Sgcr;
2262 if (problem.materialLawManager()->enableHysteresis()) {
2263 const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2264 sgcr = MaterialLaw::trappedGasSaturation(matParams, /*maxTrapping*/false);
2265 }
2266 if (sgcr >= sg) {
2267 return 0.0;
2268 }
2269 else {
2270 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2271 FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2272 : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2273 return (1.0 - xgW) *
2274 model.dofTotalVolume(ectx.globalDofIdx) *
2275 getValue(ectx.intQuants.porosity()) *
2276 getValue(ectx.fs.density(gasPhaseIdx)) *
2277 getValue(ectx.fs.saturation(gasPhaseIdx));
2278 }
2279 }
2280 }
2281 },
2282 Entry{ScalarEntry{"BGCDI",
2283 [&model = this->simulator_.model(),
2284 &problem = this->simulator_.problem()](const Context& ectx)
2285 {
2286 const auto& scaledDrainageInfo = problem.materialLawManager()
2287 ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2288 Scalar sgcr = scaledDrainageInfo.Sgcr;
2289 if (problem.materialLawManager()->enableHysteresis()) {
2290 const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2291 sgcr = MaterialLaw::trappedGasSaturation(matParams, /*maxTrapping*/false);
2292 }
2293 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2294 FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2295 : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2296 return (1.0 - xgW) *
2297 model.dofTotalVolume(ectx.globalDofIdx) *
2298 getValue(ectx.intQuants.porosity()) *
2299 getValue(ectx.fs.density(gasPhaseIdx)) *
2300 std::min(sgcr, getValue(ectx.fs.saturation(gasPhaseIdx))) /
2301 FluidSystem::molarMass(gasCompIdx, ectx.intQuants.pvtRegionIndex());
2302 }
2303 }
2304 },
2305 Entry{ScalarEntry{"BGCDM",
2306 [&model = this->simulator_.model(),
2307 &problem = this->simulator_.problem()](const Context& ectx)
2308 {
2309 const auto& scaledDrainageInfo = problem.materialLawManager()
2310 ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2311 Scalar sgcr = scaledDrainageInfo.Sgcr;
2312 if (problem.materialLawManager()->enableHysteresis()) {
2313 const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2314 sgcr = MaterialLaw::trappedGasSaturation(matParams, /*maxTrapping*/false);
2315 }
2316 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2317 FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2318 : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2319 return (1.0 - xgW) *
2320 model.dofTotalVolume(ectx.globalDofIdx) *
2321 getValue(ectx.intQuants.porosity()) *
2322 getValue(ectx.fs.density(gasPhaseIdx)) *
2323 std::max(Scalar{0.0}, getValue(ectx.fs.saturation(gasPhaseIdx)) - sgcr) /
2324 FluidSystem::molarMass(gasCompIdx, ectx.intQuants.pvtRegionIndex());
2325 }
2326 }
2327 },
2328 Entry{ScalarEntry{"BGKDI",
2329 [&model = this->simulator_.model(),
2330 &problem = this->simulator_.problem()](const Context& ectx)
2331 {
2332 const auto& scaledDrainageInfo = problem.materialLawManager()
2333 ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2334 const Scalar sg = getValue(ectx.fs.saturation(gasPhaseIdx));
2335 Scalar sgcr = scaledDrainageInfo.Sgcr;
2336 if (problem.materialLawManager()->enableHysteresis()) {
2337 const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2338 sgcr = MaterialLaw::trappedGasSaturation(matParams, /*maxTrapping*/false);
2339 }
2340 if (sg > sgcr) {
2341 return 0.0;
2342 }
2343 else {
2344 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2345 FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2346 : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2347 return (1.0 - xgW) *
2348 model.dofTotalVolume(ectx.globalDofIdx) *
2349 getValue(ectx.intQuants.porosity()) *
2350 getValue(ectx.fs.density(gasPhaseIdx)) *
2351 getValue(ectx.fs.saturation(gasPhaseIdx)) /
2352 FluidSystem::molarMass(gasCompIdx, ectx.intQuants.pvtRegionIndex());
2353 }
2354 }
2355 }
2356 },
2357 Entry{ScalarEntry{"BGKDM",
2358 [&model = this->simulator_.model(),
2359 &problem = this->simulator_.problem()](const Context& ectx)
2360 {
2361 const auto& scaledDrainageInfo = problem.materialLawManager()
2362 ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2363 const Scalar sg = getValue(ectx.fs.saturation(gasPhaseIdx));
2364 Scalar sgcr = scaledDrainageInfo.Sgcr;
2365 if (problem.materialLawManager()->enableHysteresis()) {
2366 const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2367 sgcr = MaterialLaw::trappedGasSaturation(matParams, /*maxTrapping*/false);
2368 }
2369 if (sgcr >= sg) {
2370 return 0.0;
2371 }
2372 else {
2373 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2374 FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2375 : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2376 return (1.0 - xgW) *
2377 model.dofTotalVolume(ectx.globalDofIdx) *
2378 getValue(ectx.intQuants.porosity()) *
2379 getValue(ectx.fs.density(gasPhaseIdx)) *
2380 getValue(ectx.fs.saturation(gasPhaseIdx)) /
2381 FluidSystem::molarMass(gasCompIdx, ectx.intQuants.pvtRegionIndex());
2382 }
2383 }
2384 }
2385 },
2386 Entry{ScalarEntry{"BWCD",
2387 [&model = this->simulator_.model()](const Context& ectx)
2388 {
2389 Scalar result;
2390 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
2391 result = getValue(ectx.fs.Rs()) *
2392 getValue(ectx.fs.invB(oilPhaseIdx)) *
2393 getValue(ectx.fs.saturation(oilPhaseIdx));
2394 }
2395 else {
2396 result = getValue(ectx.fs.Rsw()) *
2397 getValue(ectx.fs.invB(waterPhaseIdx)) *
2398 getValue(ectx.fs.saturation(waterPhaseIdx));
2399 }
2400 const Scalar rhoG = FluidSystem::referenceDensity(gasPhaseIdx,
2401 ectx.intQuants.pvtRegionIndex());
2402 return result *
2403 model.dofTotalVolume(ectx.globalDofIdx) *
2404 getValue(ectx.intQuants.porosity()) *
2405 rhoG /
2406 FluidSystem::molarMass(gasCompIdx, ectx.intQuants.pvtRegionIndex());
2407 }
2408 }
2409 },
2410 Entry{ScalarEntry{"BWIPG",
2411 [&model = this->simulator_.model()](const Context& ectx)
2412 {
2413 Scalar result = 0.0;
2414 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
2415 result = getValue(ectx.fs.Rvw()) *
2416 getValue(ectx.fs.invB(gasPhaseIdx)) *
2417 getValue(ectx.fs.saturation(gasPhaseIdx));
2418 }
2419 return result *
2420 model.dofTotalVolume(ectx.globalDofIdx) *
2421 getValue(ectx.intQuants.porosity());
2422 }
2423 }
2424 },
2425 Entry{ScalarEntry{"BWIPL",
2426 [&model = this->simulator_.model()](const Context& ectx)
2427 {
2428 return getValue(ectx.fs.invB(waterPhaseIdx)) *
2429 getValue(ectx.fs.saturation(waterPhaseIdx)) *
2430 model.dofTotalVolume(ectx.globalDofIdx) *
2431 getValue(ectx.intQuants.porosity());
2432 }
2433 }
2434 },
2435 };
2436
2437 this->blockExtractors_ = BlockExtractor::setupExecMap(this->blockData_, handlers);
2438
2439 this->extraBlockData_.clear();
2440 if (reportStepNum > 0 && !isSubStep) {
2441 // check we need extra block pressures for RPTSCHED
2442 const auto& rpt = this->schedule_[reportStepNum - 1].rpt_config.get();
2443 if (rpt.contains("WELLS") && rpt.at("WELLS") > 1) {
2444 this->setupExtraBlockData(reportStepNum,
2445 [&c = this->collectOnIORank_](const int idx)
2446 { return c.isCartIdxOnThisRank(idx); });
2447
2448 const auto extraHandlers = std::array{
2450 };
2451
2452 this->extraBlockExtractors_ = BlockExtractor::setupExecMap(this->extraBlockData_, extraHandlers);
2453 }
2454 }
2455 }
2456
2457 const Simulator& simulator_;
2458 const CollectDataOnIORankType& collectOnIORank_;
2459 std::vector<typename Extractor::Entry> extractors_;
2460 typename BlockExtractor::ExecMap blockExtractors_;
2461 typename BlockExtractor::ExecMap extraBlockExtractors_;
2462};
2463
2464} // namespace Opm
2465
2466#endif // OPM_OUTPUT_BLACK_OIL_MODULE_HPP
Helper class for grid instantiation of ECL file-format using problems.
Output module for the results black oil model writing in ECL binary format.
Output module for the results black oil model writing in ECL binary format.
Declares the properties required by the black oil model.
Definition CollectDataOnIORank.hpp:56
The base class for the element-centered finite-volume discretization scheme.
Definition ecfvdiscretization.hh:147
Definition GenericOutputBlackoilModule.hpp:76
Inter-region flow accumulation maps for all region definition arrays.
Definition InterRegFlows.hpp:179
void compress()
Form CSR adjacency matrix representation of input graph from connections established in previous call...
Definition InterRegFlows.cpp:165
void addConnection(const Cell &source, const Cell &destination, const data::InterRegFlowMap::FlowRates &rates)
Add flow rate connection between regions for all region definitions.
Definition InterRegFlows.cpp:156
void clear()
Clear all internal buffers, but preserve allocated capacity.
Definition InterRegFlows.cpp:172
Output module for the results black oil model writing in ECL binary format.
Definition OutputBlackoilModule.hpp:85
void processElement(const ElementContext &elemCtx)
Modify the internal buffers according to the intensive quanties relevant for an element.
Definition OutputBlackoilModule.hpp:234
void initializeFluxData()
Prepare for capturing connection fluxes, particularly to account for inter-region flows.
Definition OutputBlackoilModule.hpp:441
void setupExtractors(const bool isSubStep, const int reportStepNum)
Setup list of active element-level data extractors.
Definition OutputBlackoilModule.hpp:215
void allocBuffers(const unsigned bufferSize, const unsigned reportStepNum, const bool substep, const bool log, const bool isRestart)
Allocate memory for the scalar fields we would like to write to ECL output files.
Definition OutputBlackoilModule.hpp:193
void processFluxes(const ElementContext &elemCtx, ActiveIndex &&activeIndex, CartesianIndex &&cartesianIndex)
Capture connection fluxes, particularly to account for inter-region flows.
Definition OutputBlackoilModule.hpp:404
void clearExtractors()
Clear list of active element-level data extractors.
Definition OutputBlackoilModule.hpp:223
const InterRegFlowMap & getInterRegFlows() const
Get read-only access to collection of inter-region flows.
Definition OutputBlackoilModule.hpp:459
void finalizeFluxData()
Finalize capturing connection fluxes.
Definition OutputBlackoilModule.hpp:451
Declare the properties used by the infrastructure code of the finite volume discretizations.
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
std::string moduleVersionName()
Return the version name of the module, for example "2015.10" (for a release branch) or "2016....
Definition moduleVersion.cpp:34
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition propertysystem.hh:235
This file provides the infrastructure to retrieve run-time parameters.
The Opm property system, traits with inheritance.
Minimal characteristics of a cell from a simulation grid.
Definition InterRegFlows.hpp:50
Wrapping struct holding types used for block-level data extraction.
Definition OutputExtractor.hpp:193
std::unordered_map< int, std::vector< Exec > > ExecMap
A map of extraction executors, keyed by cartesian cell index.
Definition OutputExtractor.hpp:260
std::variant< ScalarEntry, PhaseEntry > Entry
Descriptor for extractors.
Definition OutputExtractor.hpp:245
static ExecMap setupExecMap(std::map< std::pair< std::string, int >, double > &blockData, const std::array< Entry, size > &handlers)
Setup an extractor executor map from a map of evaluations to perform.
Definition OutputExtractor.hpp:264
static void process(const std::vector< Exec > &blockExtractors, const Context &ectx)
Process a list of block extractors.
Definition OutputExtractor.hpp:367
Context passed to extractor functions.
Definition OutputExtractor.hpp:74
int episodeIndex
Current report step.
Definition OutputExtractor.hpp:77
Struct holding hysteresis parameters.
Definition OutputExtractor.hpp:63
Wrapping struct holding types used for element-level data extraction.
Definition OutputExtractor.hpp:54
static void process(const Context &ectx, const std::vector< Entry > &extractors)
Process the given extractor entries.
Definition OutputExtractor.hpp:158
static std::vector< Entry > removeInactive(std::array< Entry, size > &input)
Obtain vector of active extractors from an array of extractors.
Definition OutputExtractor.hpp:120