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