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