My Project
Loading...
Searching...
No Matches
OutputBlackoilModule.hpp
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*
4 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 2 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18
19 Consult the COPYING file in the top-level source directory of this
20 module for the precise wording of the license and the list of
21 copyright holders.
22*/
27#ifndef OPM_OUTPUT_BLACK_OIL_MODULE_HPP
28#define OPM_OUTPUT_BLACK_OIL_MODULE_HPP
29
30#include <dune/common/fvector.hh>
31
32#include <opm/simulators/utils/moduleVersion.hpp>
33
34#include <opm/common/Exceptions.hpp>
35#include <opm/common/TimingMacros.hpp>
36#include <opm/common/OpmLog/OpmLog.hpp>
37#include <opm/common/utility/Visitor.hpp>
38
39#include <opm/input/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
40
41#include <opm/material/common/Valgrind.hpp>
42#include <opm/material/fluidmatrixinteractions/EclEpsScalingPoints.hpp>
43#include <opm/material/fluidstates/BlackOilFluidState.hpp>
44#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
45
50
51#include <opm/output/data/Cells.hpp>
52#include <opm/output/eclipse/EclipseIO.hpp>
53#include <opm/output/eclipse/Inplace.hpp>
54
55#include <opm/simulators/flow/CollectDataOnIORank.hpp>
59
60#include <algorithm>
61#include <array>
62#include <cassert>
63#include <cstddef>
64#include <functional>
65#include <stdexcept>
66#include <string>
67#include <type_traits>
68#include <utility>
69#include <vector>
70
71namespace Opm {
72
73// forward declaration
74template <class TypeTag>
75class EcfvDiscretization;
76
83template <class TypeTag>
84class OutputBlackOilModule : public GenericOutputBlackoilModule<GetPropType<TypeTag, Properties::FluidSystem>>
85{
95 using FluidState = typename IntensiveQuantities::FluidState;
97 using Element = typename GridView::template Codim<0>::Entity;
98 using ElementIterator = typename GridView::template Codim<0>::Iterator;
101 using Dir = FaceDir::DirEnum;
107
108 static constexpr int conti0EqIdx = Indices::conti0EqIdx;
109 static constexpr int numPhases = FluidSystem::numPhases;
110 static constexpr int oilPhaseIdx = FluidSystem::oilPhaseIdx;
111 static constexpr int gasPhaseIdx = FluidSystem::gasPhaseIdx;
112 static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
113 static constexpr int gasCompIdx = FluidSystem::gasCompIdx;
114 static constexpr int oilCompIdx = FluidSystem::oilCompIdx;
115 static constexpr int waterCompIdx = FluidSystem::waterCompIdx;
116 enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
117
118 template<int idx, class VectorType>
119 static Scalar value_or_zero(const VectorType& v)
120 {
121 if constexpr (idx == -1) {
122 return 0.0;
123 } else {
124 return v.empty() ? 0.0 : v[idx];
125 }
126 }
127
128public:
129 OutputBlackOilModule(const Simulator& simulator,
130 const SummaryConfig& smryCfg,
131 const CollectDataOnIORankType& collectOnIORank)
132 : BaseType(simulator.vanguard().eclState(),
133 simulator.vanguard().schedule(),
134 smryCfg,
135 simulator.vanguard().summaryState(),
137 [this](const int idx)
138 { return simulator_.problem().eclWriter().collectOnIORank().localIdxToGlobalIdx(idx); },
139 simulator.vanguard().grid().comm(),
150 , simulator_(simulator)
151 , collectOnIORank_(collectOnIORank)
152 {
153 for (auto& region_pair : this->regions_) {
154 this->createLocalRegion_(region_pair.second);
155 }
156
157 auto isCartIdxOnThisRank = [&collectOnIORank](const int idx) {
158 return collectOnIORank.isCartIdxOnThisRank(idx);
159 };
160
161 this->setupBlockData(isCartIdxOnThisRank);
162
163 if (! Parameters::Get<Parameters::OwnerCellsFirst>()) {
164 const std::string msg = "The output code does not support --owner-cells-first=false.";
165 if (collectOnIORank.isIORank()) {
166 OpmLog::error(msg);
167 }
168 OPM_THROW_NOLOG(std::runtime_error, msg);
169 }
170
171 if (smryCfg.match("[FB]PP[OGW]") || smryCfg.match("RPP[OGW]*")) {
172 auto rset = this->eclState_.fieldProps().fip_regions();
173 rset.push_back("PVTNUM");
174
175 // Note: We explicitly use decltype(auto) here because the
176 // default scheme (-> auto) will deduce an undesirable type. We
177 // need the "reference to vector" semantics in this instance.
178 this->regionAvgDensity_
179 .emplace(this->simulator_.gridView().comm(),
180 FluidSystem::numPhases, rset,
181 [fp = std::cref(this->eclState_.fieldProps())]
182 (const std::string& rsetName) -> decltype(auto)
183 { return fp.get().get_int(rsetName); });
184 }
185 }
186
191 void
192 allocBuffers(const unsigned bufferSize,
193 const unsigned reportStepNum,
194 const bool substep,
195 const bool log,
196 const bool isRestart)
197 {
198 if (! std::is_same<Discretization, EcfvDiscretization<TypeTag>>::value) {
199 return;
200 }
201
202 const auto& problem = this->simulator_.problem();
203
204 this->doAllocBuffers(bufferSize,
205 reportStepNum,
206 substep,
207 log,
208 isRestart,
209 &problem.materialLawManager()->hysteresisConfig(),
210 problem.eclWriter().getOutputNnc().size());
211 }
212
214 void setupExtractors(const bool isSubStep,
215 const int reportStepNum)
216 {
217 this->setupElementExtractors_();
218 this->setupBlockExtractors_(isSubStep, reportStepNum);
219 }
220
223 {
224 this->extractors_.clear();
225 this->blockExtractors_.clear();
226 this->extraBlockExtractors_.clear();
227 }
228
233 void processElement(const ElementContext& elemCtx)
234 {
236 if (!std::is_same<Discretization, EcfvDiscretization<TypeTag>>::value) {
237 return;
238 }
239
240 if (this->extractors_.empty()) {
241 assert(0);
242 }
243
244 const auto& matLawManager = simulator_.problem().materialLawManager();
245
247 for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
248 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
249 const auto& fs = intQuants.fluidState();
250
251 const typename Extractor::Context ectx{
252 elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0),
253 elemCtx.primaryVars(dofIdx, /*timeIdx=*/0).pvtRegionIndex(),
254 elemCtx.simulator().episodeIndex(),
255 fs,
256 intQuants,
258 };
259
260 if (matLawManager->enableHysteresis()) {
261 if (FluidSystem::phaseIsActive(oilPhaseIdx) && FluidSystem::phaseIsActive(waterPhaseIdx)) {
262 matLawManager->oilWaterHysteresisParams(hysterParams.somax,
263 hysterParams.swmax,
264 hysterParams.swmin,
265 ectx.globalDofIdx);
266 }
267 if (FluidSystem::phaseIsActive(oilPhaseIdx) && FluidSystem::phaseIsActive(gasPhaseIdx)) {
268 matLawManager->gasOilHysteresisParams(hysterParams.sgmax,
269 hysterParams.shmax,
270 hysterParams.somin,
271 ectx.globalDofIdx);
272 }
273 }
274
275 Extractor::process(ectx, extractors_);
276 }
277 }
278
279 void processElementBlockData(const ElementContext& elemCtx)
280 {
281 OPM_TIMEBLOCK_LOCAL(processElementBlockData);
282 if (!std::is_same<Discretization, EcfvDiscretization<TypeTag>>::value) {
283 return;
284 }
285
286 if (this->blockExtractors_.empty() && this->extraBlockExtractors_.empty()) {
287 return;
288 }
289
290 for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
291 // Adding block data
292 const auto globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
293 const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx);
294
295 const auto be_it = this->blockExtractors_.find(cartesianIdx);
296 const auto bee_it = this->extraBlockExtractors_.find(cartesianIdx);
297 if (be_it == this->blockExtractors_.end() &&
298 bee_it == this->extraBlockExtractors_.end())
299 {
300 continue;
301 }
302
303 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
304 const auto& fs = intQuants.fluidState();
305
306 const typename BlockExtractor::Context ectx{
307 globalDofIdx,
308 dofIdx,
309 fs,
310 intQuants,
311 elemCtx,
312 };
313
314 if (be_it != this->blockExtractors_.end()) {
316 }
317 if (bee_it != this->extraBlockExtractors_.end()) {
319 }
320 }
321 }
322
351 template <class ActiveIndex, class CartesianIndex>
352 void processFluxes(const ElementContext& elemCtx,
353 ActiveIndex&& activeIndex,
354 CartesianIndex&& cartesianIndex)
355 {
357 const auto identifyCell = [&activeIndex, &cartesianIndex](const Element& elem)
359 {
360 const auto cellIndex = activeIndex(elem);
361
362 return {
363 static_cast<int>(cellIndex),
364 cartesianIndex(cellIndex),
365 elem.partitionType() == Dune::InteriorEntity
366 };
367 };
368
369 const auto timeIdx = 0u;
370 const auto& stencil = elemCtx.stencil(timeIdx);
371 const auto numInteriorFaces = elemCtx.numInteriorFaces(timeIdx);
372
373 for (auto scvfIdx = 0 * numInteriorFaces; scvfIdx < numInteriorFaces; ++scvfIdx) {
374 const auto& face = stencil.interiorFace(scvfIdx);
375 const auto left = identifyCell(stencil.element(face.interiorIndex()));
376 const auto right = identifyCell(stencil.element(face.exteriorIndex()));
377
378 const auto rates = this->
379 getComponentSurfaceRates(elemCtx, face.area(), scvfIdx, timeIdx);
380
381 this->interRegionFlows_.addConnection(left, right, rates);
382 }
383 }
384
390 {
391 // Inter-region flow rates. Note: ".clear()" prepares to accumulate
392 // contributions per bulk connection between FIP regions.
393 this->interRegionFlows_.clear();
394 }
395
400 {
401 this->interRegionFlows_.compress();
402 }
403
408 {
409 return this->interRegionFlows_;
410 }
411
412 template <class FluidState>
413 void assignToFluidState(FluidState& fs, unsigned elemIdx) const
414 {
415 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
416 if (this->saturation_[phaseIdx].empty())
417 continue;
418
419 fs.setSaturation(phaseIdx, this->saturation_[phaseIdx][elemIdx]);
420 }
421
422 if (!this->fluidPressure_.empty()) {
423 // this assumes that capillary pressures only depend on the phase saturations
424 // and possibly on temperature. (this is always the case for ECL problems.)
425 std::array<Scalar, numPhases> pc = {0};
426 const MaterialLawParams& matParams = simulator_.problem().materialLawParams(elemIdx);
427 MaterialLaw::capillaryPressures(pc, matParams, fs);
428 Valgrind::CheckDefined(this->fluidPressure_[elemIdx]);
429 Valgrind::CheckDefined(pc);
430 const auto& pressure = this->fluidPressure_[elemIdx];
431 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
432 if (!FluidSystem::phaseIsActive(phaseIdx))
433 continue;
434
435 if (Indices::oilEnabled)
436 fs.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx]));
437 else if (Indices::gasEnabled)
438 fs.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx]));
439 else if (Indices::waterEnabled)
440 //single (water) phase
441 fs.setPressure(phaseIdx, pressure);
442 }
443 }
444
445 if (!this->temperature_.empty())
446 fs.setTemperature(this->temperature_[elemIdx]);
447 if (!this->rs_.empty())
448 fs.setRs(this->rs_[elemIdx]);
449 if (!this->rsw_.empty())
450 fs.setRsw(this->rsw_[elemIdx]);
451 if (!this->rv_.empty())
452 fs.setRv(this->rv_[elemIdx]);
453 if (!this->rvw_.empty())
454 fs.setRvw(this->rvw_[elemIdx]);
455 }
456
457 void initHysteresisParams(Simulator& simulator, unsigned elemIdx) const
458 {
459 if (!this->soMax_.empty())
460 simulator.problem().setMaxOilSaturation(elemIdx, this->soMax_[elemIdx]);
461
462 if (simulator.problem().materialLawManager()->enableHysteresis()) {
463 auto matLawManager = simulator.problem().materialLawManager();
464
465 if (FluidSystem::phaseIsActive(oilPhaseIdx)
466 && FluidSystem::phaseIsActive(waterPhaseIdx)) {
467 Scalar somax = 2.0;
468 Scalar swmax = -2.0;
469 Scalar swmin = 2.0;
470
471 if (matLawManager->enableNonWettingHysteresis()) {
472 if (!this->soMax_.empty()) {
473 somax = this->soMax_[elemIdx];
474 }
475 }
476 if (matLawManager->enableWettingHysteresis()) {
477 if (!this->swMax_.empty()) {
478 swmax = this->swMax_[elemIdx];
479 }
480 }
481 if (matLawManager->enablePCHysteresis()) {
482 if (!this->swmin_.empty()) {
483 swmin = this->swmin_[elemIdx];
484 }
485 }
486 matLawManager->setOilWaterHysteresisParams(
487 somax, swmax, swmin, elemIdx);
488 }
489 if (FluidSystem::phaseIsActive(oilPhaseIdx)
490 && FluidSystem::phaseIsActive(gasPhaseIdx)) {
491 Scalar sgmax = 2.0;
492 Scalar shmax = -2.0;
493 Scalar somin = 2.0;
494
495 if (matLawManager->enableNonWettingHysteresis()) {
496 if (!this->sgmax_.empty()) {
497 sgmax = this->sgmax_[elemIdx];
498 }
499 }
500 if (matLawManager->enableWettingHysteresis()) {
501 if (!this->shmax_.empty()) {
502 shmax = this->shmax_[elemIdx];
503 }
504 }
505 if (matLawManager->enablePCHysteresis()) {
506 if (!this->somin_.empty()) {
507 somin = this->somin_[elemIdx];
508 }
509 }
510 matLawManager->setGasOilHysteresisParams(
511 sgmax, shmax, somin, elemIdx);
512 }
513
514 }
515
516 if (simulator_.vanguard().eclState().fieldProps().has_double("SWATINIT")) {
517 simulator.problem().materialLawManager()
518 ->applyRestartSwatInit(elemIdx, this->ppcw_[elemIdx]);
519 }
520 }
521
522 void updateFluidInPlace(const ElementContext& elemCtx)
523 {
524 for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
525 updateFluidInPlace_(elemCtx, dofIdx);
526 }
527 }
528
529 void updateFluidInPlace(const unsigned globalDofIdx,
530 const IntensiveQuantities& intQuants,
531 const double totVolume)
532 {
533 this->updateFluidInPlace_(globalDofIdx, intQuants, totVolume);
534 }
535
536private:
537 template <typename T>
538 using RemoveCVR = std::remove_cv_t<std::remove_reference_t<T>>;
539
540 template <typename, class = void>
541 struct HasGeoMech : public std::false_type {};
542
543 template <typename Problem>
544 struct HasGeoMech<
545 Problem, std::void_t<decltype(std::declval<Problem>().geoMechModel())>
546 > : public std::true_type {};
547
548 bool isDefunctParallelWell(std::string wname) const override
549 {
550 if (simulator_.gridView().comm().size() == 1)
551 return false;
552 const auto& parallelWells = simulator_.vanguard().parallelWells();
553 std::pair<std::string, bool> value {wname, true};
554 auto candidate = std::lower_bound(parallelWells.begin(), parallelWells.end(), value);
555 return candidate == parallelWells.end() || *candidate != value;
556 }
557
558 void updateFluidInPlace_(const ElementContext& elemCtx, const unsigned dofIdx)
559 {
560 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
561 const unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
562 const auto totVolume = elemCtx.simulator().model().dofTotalVolume(globalDofIdx);
563
564 this->updateFluidInPlace_(globalDofIdx, intQuants, totVolume);
565 }
566
567 void updateFluidInPlace_(const unsigned globalDofIdx,
568 const IntensiveQuantities& intQuants,
569 const double totVolume)
570 {
571 OPM_TIMEBLOCK_LOCAL(updateFluidInPlace);
572
573 this->updateTotalVolumesAndPressures_(globalDofIdx, intQuants, totVolume);
574
575 if (this->computeFip_) {
576 this->updatePhaseInplaceVolumes_(globalDofIdx, intQuants, totVolume);
577 }
578 }
579
580 void createLocalRegion_(std::vector<int>& region)
581 {
582 std::size_t elemIdx = 0;
583 for (const auto& elem : elements(simulator_.gridView())) {
584 if (elem.partitionType() != Dune::InteriorEntity) {
585 region[elemIdx] = 0;
586 }
587
588 ++elemIdx;
589 }
590 }
591
592 template <typename FluidState>
593 void aggregateAverageDensityContributions_(const FluidState& fs,
594 const unsigned int globalDofIdx,
595 const double porv)
596 {
597 auto pvCellValue = RegionPhasePoreVolAverage::CellValue{};
598 pvCellValue.porv = porv;
599
600 for (auto phaseIdx = 0*FluidSystem::numPhases;
601 phaseIdx < FluidSystem::numPhases; ++phaseIdx)
602 {
603 if (! FluidSystem::phaseIsActive(phaseIdx)) {
604 continue;
605 }
606
607 pvCellValue.value = getValue(fs.density(phaseIdx));
608 pvCellValue.sat = getValue(fs.saturation(phaseIdx));
609
610 this->regionAvgDensity_
611 ->addCell(globalDofIdx,
612 RegionPhasePoreVolAverage::Phase { phaseIdx },
614 }
615 }
616
632 data::InterRegFlowMap::FlowRates
633 getComponentSurfaceRates(const ElementContext& elemCtx,
634 const Scalar faceArea,
635 const std::size_t scvfIdx,
636 const std::size_t timeIdx) const
637 {
638 using Component = data::InterRegFlowMap::Component;
639
640 auto rates = data::InterRegFlowMap::FlowRates {};
641
642 const auto& extQuant = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
643
644 const auto alpha = getValue(extQuant.extrusionFactor()) * faceArea;
645
646 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
647 const auto& up = elemCtx
648 .intensiveQuantities(extQuant.upstreamIndex(oilPhaseIdx), timeIdx);
649
650 const auto pvtReg = up.pvtRegionIndex();
651
653 (up.fluidState(), oilPhaseIdx, pvtReg));
654
655 const auto qO = alpha * bO * getValue(extQuant.volumeFlux(oilPhaseIdx));
656
657 rates[Component::Oil] += qO;
658
659 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
660 const auto Rs = getValue(
661 BlackOil::getRs_<FluidSystem, FluidState, Scalar>
662 (up.fluidState(), pvtReg));
663
664 rates[Component::Gas] += qO * Rs;
665 rates[Component::Disgas] += qO * Rs;
666 }
667 }
668
669 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
670 const auto& up = elemCtx
671 .intensiveQuantities(extQuant.upstreamIndex(gasPhaseIdx), timeIdx);
672
673 const auto pvtReg = up.pvtRegionIndex();
674
676 (up.fluidState(), gasPhaseIdx, pvtReg));
677
678 const auto qG = alpha * bG * getValue(extQuant.volumeFlux(gasPhaseIdx));
679
680 rates[Component::Gas] += qG;
681
682 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
683 const auto Rv = getValue(
684 BlackOil::getRv_<FluidSystem, FluidState, Scalar>
685 (up.fluidState(), pvtReg));
686
687 rates[Component::Oil] += qG * Rv;
688 rates[Component::Vapoil] += qG * Rv;
689 }
690 }
691
692 if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
693 const auto& up = elemCtx
694 .intensiveQuantities(extQuant.upstreamIndex(waterPhaseIdx), timeIdx);
695
696 const auto pvtReg = up.pvtRegionIndex();
697
699 (up.fluidState(), waterPhaseIdx, pvtReg));
700
701 rates[Component::Water] +=
702 alpha * bW * getValue(extQuant.volumeFlux(waterPhaseIdx));
703 }
704
705 return rates;
706 }
707
708 template <typename FluidState>
709 Scalar hydroCarbonFraction(const FluidState& fs) const
710 {
711 if (this->eclState_.runspec().co2Storage()) {
712 // CO2 storage: Hydrocarbon volume is full pore-volume.
713 return 1.0;
714 }
715
716 // Common case. Hydrocarbon volume is fraction occupied by actual
717 // hydrocarbons.
718 auto hydrocarbon = Scalar {0};
719 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
720 hydrocarbon += getValue(fs.saturation(oilPhaseIdx));
721 }
722
723 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
724 hydrocarbon += getValue(fs.saturation(gasPhaseIdx));
725 }
726
727 return hydrocarbon;
728 }
729
730 void updateTotalVolumesAndPressures_(const unsigned globalDofIdx,
731 const IntensiveQuantities& intQuants,
732 const double totVolume)
733 {
734 const auto& fs = intQuants.fluidState();
735
736 const double pv = totVolume * intQuants.porosity().value();
737 const auto hydrocarbon = this->hydroCarbonFraction(fs);
738
739 if (! this->hydrocarbonPoreVolume_.empty()) {
740 this->fipC_.assignPoreVolume(globalDofIdx,
741 totVolume * intQuants.referencePorosity());
742
743 this->dynamicPoreVolume_[globalDofIdx] = pv;
744 this->hydrocarbonPoreVolume_[globalDofIdx] = pv * hydrocarbon;
745 }
746
747 if (!this->pressureTimesHydrocarbonVolume_.empty() &&
748 !this->pressureTimesPoreVolume_.empty())
749 {
750 assert(this->hydrocarbonPoreVolume_.size() == this->pressureTimesHydrocarbonVolume_.size());
751 assert(this->fipC_.get(Inplace::Phase::PoreVolume).size() == this->pressureTimesPoreVolume_.size());
752
753 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
754 this->pressureTimesPoreVolume_[globalDofIdx] =
755 getValue(fs.pressure(oilPhaseIdx)) * pv;
756
757 this->pressureTimesHydrocarbonVolume_[globalDofIdx] =
758 this->pressureTimesPoreVolume_[globalDofIdx] * hydrocarbon;
759 }
760 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
761 this->pressureTimesPoreVolume_[globalDofIdx] =
762 getValue(fs.pressure(gasPhaseIdx)) * pv;
763
764 this->pressureTimesHydrocarbonVolume_[globalDofIdx] =
765 this->pressureTimesPoreVolume_[globalDofIdx] * hydrocarbon;
766 }
767 else if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
768 this->pressureTimesPoreVolume_[globalDofIdx] =
769 getValue(fs.pressure(waterPhaseIdx)) * pv;
770 }
771 }
772 }
773
774 void updatePhaseInplaceVolumes_(const unsigned globalDofIdx,
775 const IntensiveQuantities& intQuants,
776 const double totVolume)
777 {
778 std::array<Scalar, FluidSystem::numPhases> fip {};
779 std::array<Scalar, FluidSystem::numPhases> fipr{}; // at reservoir condition
780
781 const auto& fs = intQuants.fluidState();
782 const auto pv = totVolume * intQuants.porosity().value();
783
784 for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
785 if (!FluidSystem::phaseIsActive(phaseIdx)) {
786 continue;
787 }
788
789 const auto b = getValue(fs.invB(phaseIdx));
790 const auto s = getValue(fs.saturation(phaseIdx));
791
792 fipr[phaseIdx] = s * pv;
793 fip [phaseIdx] = b * fipr[phaseIdx];
794 }
795
796 this->fipC_.assignVolumesSurface(globalDofIdx, fip);
797 this->fipC_.assignVolumesReservoir(globalDofIdx,
798 fs.saltConcentration().value(),
799 fipr);
800
801 if (FluidSystem::phaseIsActive(oilPhaseIdx) &&
802 FluidSystem::phaseIsActive(gasPhaseIdx))
803 {
804 this->updateOilGasDistribution(globalDofIdx, fs, fip);
805 }
806
807 if (FluidSystem::phaseIsActive(waterPhaseIdx) &&
808 FluidSystem::phaseIsActive(gasPhaseIdx))
809 {
810 this->updateGasWaterDistribution(globalDofIdx, fs, fip);
811 }
812
813 if (FluidSystem::phaseIsActive(gasPhaseIdx) &&
814 this->fipC_.hasCo2InGas())
815 {
816 this->updateCO2InGas(globalDofIdx, pv, intQuants);
817 }
818
819 if (this->fipC_.hasCo2InWater() &&
820 (FluidSystem::phaseIsActive(waterPhaseIdx) ||
821 FluidSystem::phaseIsActive(oilPhaseIdx)))
822 {
823 this->updateCO2InWater(globalDofIdx, pv, fs);
824 }
825 }
826
827 template <typename FluidState, typename FIPArray>
828 void updateOilGasDistribution(const unsigned globalDofIdx,
829 const FluidState& fs,
830 const FIPArray& fip)
831 {
832 // Gas dissolved in oil and vaporized oil
833 const auto gasInPlaceLiquid = getValue(fs.Rs()) * fip[oilPhaseIdx];
834 const auto oilInPlaceGas = getValue(fs.Rv()) * fip[gasPhaseIdx];
835
836 this->fipC_.assignOilGasDistribution(globalDofIdx, gasInPlaceLiquid, oilInPlaceGas);
837 }
838
839 template <typename FluidState, typename FIPArray>
840 void updateGasWaterDistribution(const unsigned globalDofIdx,
841 const FluidState& fs,
842 const FIPArray& fip)
843 {
844 // Gas dissolved in water and vaporized water
845 const auto gasInPlaceWater = getValue(fs.Rsw()) * fip[waterPhaseIdx];
846 const auto waterInPlaceGas = getValue(fs.Rvw()) * fip[gasPhaseIdx];
847
848 this->fipC_.assignGasWater(globalDofIdx, fip, gasInPlaceWater, waterInPlaceGas);
849 }
850
851 template <typename IntensiveQuantities>
852 void updateCO2InGas(const unsigned globalDofIdx,
853 const double pv,
854 const IntensiveQuantities& intQuants)
855 {
856 const auto& scaledDrainageInfo = this->simulator_.problem().materialLawManager()
857 ->oilWaterScaledEpsInfoDrainage(globalDofIdx);
858
859 const auto& fs = intQuants.fluidState();
860 Scalar sgcr = scaledDrainageInfo.Sgcr;
861 if (this->simulator_.problem().materialLawManager()->enableHysteresis()) {
862 const auto& matParams = simulator_.problem().materialLawParams(globalDofIdx);
863 sgcr = MaterialLaw::trappedGasSaturation(matParams, /*maximumTrapping*/false);
864 }
865
867 if (this->fipC_.has(Inplace::Phase::CO2MassInGasPhaseMaximumTrapped) ||
868 this->fipC_.has(Inplace::Phase::CO2MassInGasPhaseMaximumUnTrapped))
869 {
870 if (this->simulator_.problem().materialLawManager()->enableHysteresis()) {
871 const auto& matParams = simulator_.problem().materialLawParams(globalDofIdx);
872 // Get the maximum trapped gas saturation
873 trappedGasSaturation = MaterialLaw::trappedGasSaturation(matParams, /*maximumTrapping*/true);
874 }
875 }
876
877 const Scalar sg = getValue(fs.saturation(gasPhaseIdx));
879 if (this->fipC_.has(Inplace::Phase::CO2MassInGasPhaseEffectiveTrapped) ||
880 this->fipC_.has(Inplace::Phase::CO2MassInGasPhaseEffectiveUnTrapped))
881 {
882 if (this->simulator_.problem().materialLawManager()->enableHysteresis()) {
883 const auto& matParams = simulator_.problem().materialLawParams(globalDofIdx);
884 const double krg = getValue(intQuants.relativePermeability(gasPhaseIdx));
885 strandedGasSaturation = MaterialLaw::strandedGasSaturation(matParams, sg, krg);
886 }
887 }
888
889 const typename FIPContainer<FluidSystem>::Co2InGasInput v{
890 pv,
891 sg,
892 sgcr,
893 getValue(fs.density(gasPhaseIdx)),
894 FluidSystem::phaseIsActive(waterPhaseIdx)
895 ? FluidSystem::convertRvwToXgW(getValue(fs.Rvw()), fs.pvtRegionIndex())
896 : FluidSystem::convertRvToXgO(getValue(fs.Rv()), fs.pvtRegionIndex()),
897 FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex()),
900 };
901
902 this->fipC_.assignCo2InGas(globalDofIdx, v);
903 }
904
905 template <typename FluidState>
906 void updateCO2InWater(const unsigned globalDofIdx,
907 const double pv,
908 const FluidState& fs)
909 {
910 const auto co2InWater = FluidSystem::phaseIsActive(oilPhaseIdx)
911 ? this->co2InWaterFromOil(fs, pv)
912 : this->co2InWaterFromWater(fs, pv);
913
914 const Scalar mM = FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex());
915
916 this->fipC_.assignCo2InWater(globalDofIdx, co2InWater, mM);
917 }
918
919 template <typename FluidState>
920 Scalar co2InWaterFromWater(const FluidState& fs, const double pv) const
921 {
922 const double rhow = getValue(fs.density(waterPhaseIdx));
923 const double sw = getValue(fs.saturation(waterPhaseIdx));
924 const double xwG = FluidSystem::convertRswToXwG(getValue(fs.Rsw()), fs.pvtRegionIndex());
925
926 const Scalar mM = FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex());
927
928 return xwG * pv * rhow * sw / mM;
929 }
930
931 template <typename FluidState>
932 Scalar co2InWaterFromOil(const FluidState& fs, const double pv) const
933 {
934 const double rhoo = getValue(fs.density(oilPhaseIdx));
935 const double so = getValue(fs.saturation(oilPhaseIdx));
936 const double xoG = FluidSystem::convertRsToXoG(getValue(fs.Rs()), fs.pvtRegionIndex());
937
938 const Scalar mM = FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex());
939
940 return xoG * pv * rhoo * so / mM;
941 }
942
944 void setupElementExtractors_()
945 {
946 using Entry = typename Extractor::Entry;
947 using Context = typename Extractor::Context;
948 using ScalarEntry = typename Extractor::ScalarEntry;
949 using PhaseEntry = typename Extractor::PhaseEntry;
950
951 const bool hasResidual = simulator_.model().linearizer().residual().size() > 0;
952 const auto& hysteresisConfig = simulator_.problem().materialLawManager()->hysteresisConfig();
953
954 auto extractors = std::array{
955 Entry{PhaseEntry{&this->saturation_,
956 [](const unsigned phase, const Context& ectx)
957 { return getValue(ectx.fs.saturation(phase)); }
958 }
959 },
960 Entry{PhaseEntry{&this->invB_,
961 [](const unsigned phase, const Context& ectx)
962 { return getValue(ectx.fs.invB(phase)); }
963 }
964 },
965 Entry{PhaseEntry{&this->density_,
966 [](const unsigned phase, const Context& ectx)
967 { return getValue(ectx.fs.density(phase)); }
968 }
969 },
970 Entry{PhaseEntry{&this->relativePermeability_,
971 [](const unsigned phase, const Context& ectx)
972 { return getValue(ectx.intQuants.relativePermeability(phase)); }
973 }
974 },
975 Entry{PhaseEntry{&this->viscosity_,
976 [this](const unsigned phaseIdx, const Context& ectx)
977 {
978 if (this->extboC_.allocated() && phaseIdx == oilPhaseIdx) {
979 return getValue(ectx.intQuants.oilViscosity());
980 }
981 else if (this->extboC_.allocated() && phaseIdx == gasPhaseIdx) {
982 return getValue(ectx.intQuants.gasViscosity());
983 }
984 else {
985 return getValue(ectx.fs.viscosity(phaseIdx));
986 }
987 }
988 }
989 },
990 Entry{PhaseEntry{&this->residual_,
991 [&modelResid = this->simulator_.model().linearizer().residual()]
992 (const unsigned phaseIdx, const Context& ectx)
993 {
994 const unsigned sIdx = FluidSystem::solventComponentIndex(phaseIdx);
995 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(sIdx);
996 return modelResid[ectx.globalDofIdx][activeCompIdx];
997 }
998 },
1000 },
1001 Entry{ScalarEntry{&this->rockCompPorvMultiplier_,
1002 [&problem = this->simulator_.problem()](const Context& ectx)
1003 {
1004 return problem.template
1006 ectx.globalDofIdx);
1007 }
1008 }
1009 },
1010 Entry{ScalarEntry{&this->rockCompTransMultiplier_,
1011 [&problem = this->simulator_.problem()](const Context& ectx)
1012 {
1013 return problem.
1014 template rockCompTransMultiplier<Scalar>(ectx.intQuants,
1015 ectx.globalDofIdx);
1016 }}
1017 },
1018 Entry{ScalarEntry{&this->minimumOilPressure_,
1019 [&problem = this->simulator_.problem()](const Context& ectx)
1020 {
1021 return std::min(getValue(ectx.fs.pressure(oilPhaseIdx)),
1022 problem.minOilPressure(ectx.globalDofIdx));
1023 }
1024 }
1025 },
1026 Entry{ScalarEntry{&this->bubblePointPressure_,
1027 [&failedCells = this->failedCellsPb_,
1028 &vanguard = this->simulator_.vanguard()](const Context& ectx)
1029 {
1030 try {
1031 return getValue(
1032 FluidSystem::bubblePointPressure(ectx.fs,
1033 ectx.intQuants.pvtRegionIndex())
1034 );
1035 } catch (const NumericalProblem&) {
1036 const auto cartesianIdx = vanguard.cartesianIndex(ectx.globalDofIdx);
1037 failedCells.push_back(cartesianIdx);
1038 return Scalar{0};
1039 }
1040 }
1041 }
1042 },
1043 Entry{ScalarEntry{&this->dewPointPressure_,
1044 [&failedCells = this->failedCellsPd_,
1045 &vanguard = this->simulator_.vanguard()](const Context& ectx)
1046 {
1047 try {
1048 return getValue(
1049 FluidSystem::dewPointPressure(ectx.fs,
1050 ectx.intQuants.pvtRegionIndex())
1051 );
1052 } catch (const NumericalProblem&) {
1053 const auto cartesianIdx = vanguard.cartesianIndex(ectx.globalDofIdx);
1054 failedCells.push_back(cartesianIdx);
1055 return Scalar{0};
1056 }
1057 }
1058 }
1059 },
1060 Entry{ScalarEntry{&this->overburdenPressure_,
1061 [&problem = simulator_.problem()](const Context& ectx)
1062 { return problem.overburdenPressure(ectx.globalDofIdx); }
1063 }
1064 },
1065 Entry{ScalarEntry{&this->temperature_,
1066 [](const Context& ectx)
1067 { return getValue(ectx.fs.temperature(oilPhaseIdx)); }
1068 }
1069 },
1070 Entry{ScalarEntry{&this->sSol_,
1071 [](const Context& ectx)
1072 { return getValue(ectx.intQuants.solventSaturation()); }
1073 }
1074 },
1075 Entry{ScalarEntry{&this->rswSol_,
1076 [](const Context& ectx)
1077 { return getValue(ectx.intQuants.rsSolw()); }
1078 }
1079 },
1080 Entry{ScalarEntry{&this->cPolymer_,
1081 [](const Context& ectx)
1082 { return getValue(ectx.intQuants.polymerConcentration()); }
1083 }
1084 },
1085 Entry{ScalarEntry{&this->cFoam_,
1086 [](const Context& ectx)
1087 { return getValue(ectx.intQuants.foamConcentration()); }
1088 }
1089 },
1090 Entry{ScalarEntry{&this->cSalt_,
1091 [](const Context& ectx)
1092 { return getValue(ectx.fs.saltConcentration()); }
1093 }
1094 },
1095 Entry{ScalarEntry{&this->pSalt_,
1096 [](const Context& ectx)
1097 { return getValue(ectx.fs.saltSaturation()); }
1098 }
1099 },
1100 Entry{ScalarEntry{&this->permFact_,
1101 [](const Context& ectx)
1102 { return getValue(ectx.intQuants.permFactor()); }
1103 }
1104 },
1105 Entry{ScalarEntry{&this->rPorV_,
1106 [&model = this->simulator_.model()](const Context& ectx)
1107 {
1108 const auto totVolume = model.dofTotalVolume(ectx.globalDofIdx);
1109 return totVolume * getValue(ectx.intQuants.porosity());
1110 }
1111 }
1112 },
1113 Entry{ScalarEntry{&this->rs_,
1114 [](const Context& ectx)
1115 { return getValue(ectx.fs.Rs()); }
1116 }
1117 },
1118 Entry{ScalarEntry{&this->rv_,
1119 [](const Context& ectx)
1120 { return getValue(ectx.fs.Rv()); }
1121 }
1122 },
1123 Entry{ScalarEntry{&this->rsw_,
1124 [](const Context& ectx)
1125 { return getValue(ectx.fs.Rsw()); }
1126 }
1127 },
1128 Entry{ScalarEntry{&this->rvw_,
1129 [](const Context& ectx)
1130 { return getValue(ectx.fs.Rvw()); }
1131 }
1132 },
1133 Entry{ScalarEntry{&this->ppcw_,
1134 [&matLawManager = *this->simulator_.problem().materialLawManager()]
1135 (const Context& ectx)
1136 {
1137 return matLawManager.
1138 oilWaterScaledEpsInfoDrainage(ectx.globalDofIdx).maxPcow;
1139 }
1140 }
1141 },
1142 Entry{ScalarEntry{&this->drsdtcon_,
1143 [&problem = this->simulator_.problem()](const Context& ectx)
1144 {
1145 return problem.drsdtcon(ectx.globalDofIdx,
1146 ectx.episodeIndex);
1147 }
1148 }
1149 },
1150 Entry{ScalarEntry{&this->pcgw_,
1151 [](const Context& ectx)
1152 {
1153 return getValue(ectx.fs.pressure(gasPhaseIdx)) -
1154 getValue(ectx.fs.pressure(waterPhaseIdx));
1155 }
1156 }
1157 },
1158 Entry{ScalarEntry{&this->pcow_,
1159 [](const Context& ectx)
1160 {
1161 return getValue(ectx.fs.pressure(oilPhaseIdx)) -
1162 getValue(ectx.fs.pressure(waterPhaseIdx));
1163 }
1164 }
1165 },
1166 Entry{ScalarEntry{&this->pcog_,
1167 [](const Context& ectx)
1168 {
1169 return getValue(ectx.fs.pressure(gasPhaseIdx)) -
1170 getValue(ectx.fs.pressure(oilPhaseIdx));
1171 }
1172 }
1173 },
1174 Entry{ScalarEntry{&this->fluidPressure_,
1175 [](const Context& ectx)
1176 {
1177 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1178 // Output oil pressure as default
1179 return getValue(ectx.fs.pressure(oilPhaseIdx));
1180 }
1181 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1182 // Output gas if oil is not present
1183 return getValue(ectx.fs.pressure(gasPhaseIdx));
1184 }
1185 else {
1186 // Output water if neither oil nor gas is present
1187 return getValue(ectx.fs.pressure(waterPhaseIdx));
1188 }
1189 }
1190 }
1191 },
1192 Entry{ScalarEntry{&this->gasDissolutionFactor_,
1193 [&problem = this->simulator_.problem()](const Context& ectx)
1194 {
1195 const Scalar SoMax = problem.maxOilSaturation(ectx.globalDofIdx);
1196 return FluidSystem::template
1198 oilPhaseIdx,
1199 ectx.pvtRegionIdx,
1200 SoMax);
1201 }
1202 }
1203 },
1204 Entry{ScalarEntry{&this->oilVaporizationFactor_,
1205 [&problem = this->simulator_.problem()](const Context& ectx)
1206 {
1207 const Scalar SoMax = problem.maxOilSaturation(ectx.globalDofIdx);
1208 return FluidSystem::template
1210 gasPhaseIdx,
1211 ectx.pvtRegionIdx,
1212 SoMax);
1213 }
1214 }
1215 },
1216 Entry{ScalarEntry{&this->gasDissolutionFactorInWater_,
1217 [&problem = this->simulator_.problem()](const Context& ectx)
1218 {
1219 const Scalar SwMax = problem.maxWaterSaturation(ectx.globalDofIdx);
1220 return FluidSystem::template
1222 waterPhaseIdx,
1223 ectx.pvtRegionIdx,
1224 SwMax);
1225 }
1226 }
1227 },
1228 Entry{ScalarEntry{&this->waterVaporizationFactor_,
1229 [](const Context& ectx)
1230 {
1231 return FluidSystem::template
1233 gasPhaseIdx,
1234 ectx.pvtRegionIdx);
1235 }
1236 }
1237 },
1238 Entry{ScalarEntry{&this->gasFormationVolumeFactor_,
1239 [](const Context& ectx)
1240 {
1241 return 1.0 / FluidSystem::template
1243 gasPhaseIdx,
1244 ectx.pvtRegionIdx);
1245 }
1246 }
1247 },
1248 Entry{ScalarEntry{&this->saturatedOilFormationVolumeFactor_,
1249 [](const Context& ectx)
1250 {
1251 return 1.0 / FluidSystem::template
1253 oilPhaseIdx,
1254 ectx.pvtRegionIdx);
1255 }
1256 }
1257 },
1258 Entry{ScalarEntry{&this->oilSaturationPressure_,
1259 [](const Context& ectx)
1260 {
1261 return FluidSystem::template
1263 oilPhaseIdx,
1264 ectx.pvtRegionIdx);
1265 }
1266 }
1267 },
1268 Entry{ScalarEntry{&this->soMax_,
1269 [&problem = this->simulator_.problem()](const Context& ectx)
1270 {
1271 return std::max(getValue(ectx.fs.saturation(oilPhaseIdx)),
1272 problem.maxOilSaturation(ectx.globalDofIdx));
1273 }
1274 },
1275 !hysteresisConfig.enableHysteresis()
1276 },
1277 Entry{ScalarEntry{&this->swMax_,
1278 [&problem = this->simulator_.problem()](const Context& ectx)
1279 {
1280 return std::max(getValue(ectx.fs.saturation(waterPhaseIdx)),
1281 problem.maxWaterSaturation(ectx.globalDofIdx));
1282 }
1283 },
1284 !hysteresisConfig.enableHysteresis()
1285 },
1286 Entry{ScalarEntry{&this->soMax_,
1287 [](const Context& ectx)
1288 { return ectx.hParams.somax; }
1289 },
1290 hysteresisConfig.enableHysteresis() &&
1291 hysteresisConfig.enableNonWettingHysteresis() &&
1292 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1293 FluidSystem::phaseIsActive(waterPhaseIdx)
1294 },
1295 Entry{ScalarEntry{&this->swMax_,
1296 [](const Context& ectx)
1297 { return ectx.hParams.swmax; }
1298 },
1299 hysteresisConfig.enableHysteresis() &&
1300 hysteresisConfig.enableWettingHysteresis() &&
1301 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1302 FluidSystem::phaseIsActive(waterPhaseIdx)
1303 },
1304 Entry{ScalarEntry{&this->swmin_,
1305 [](const Context& ectx)
1306 { return ectx.hParams.swmin; }
1307 },
1308 hysteresisConfig.enableHysteresis() &&
1309 hysteresisConfig.enablePCHysteresis() &&
1310 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1311 FluidSystem::phaseIsActive(waterPhaseIdx)
1312 },
1313 Entry{ScalarEntry{&this->sgmax_,
1314 [](const Context& ectx)
1315 { return ectx.hParams.sgmax; }
1316 },
1317 hysteresisConfig.enableHysteresis() &&
1318 hysteresisConfig.enableNonWettingHysteresis() &&
1319 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1320 FluidSystem::phaseIsActive(gasPhaseIdx)
1321 },
1322 Entry{ScalarEntry{&this->shmax_,
1323 [](const Context& ectx)
1324 { return ectx.hParams.shmax; }
1325 },
1326 hysteresisConfig.enableHysteresis() &&
1327 hysteresisConfig.enableWettingHysteresis() &&
1328 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1329 FluidSystem::phaseIsActive(gasPhaseIdx)
1330 },
1331 Entry{ScalarEntry{&this->somin_,
1332 [](const Context& ectx)
1333 { return ectx.hParams.somin; }
1334 },
1335 hysteresisConfig.enableHysteresis() &&
1336 hysteresisConfig.enablePCHysteresis() &&
1337 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1338 FluidSystem::phaseIsActive(gasPhaseIdx)
1339 },
1340 Entry{[&model = this->simulator_.model(), this](const Context& ectx)
1341 {
1342 // Note: We intentionally exclude effects of rock
1343 // compressibility by using referencePorosity() here.
1344 const auto porv = ectx.intQuants.referencePorosity()
1345 * model.dofTotalVolume(ectx.globalDofIdx);
1346
1347 this->aggregateAverageDensityContributions_(ectx.fs, ectx.globalDofIdx,
1348 static_cast<double>(porv));
1349 }, this->regionAvgDensity_.has_value()
1350 },
1351 Entry{[&extboC = this->extboC_](const Context& ectx)
1352 {
1353 extboC.assignVolumes(ectx.globalDofIdx,
1354 ectx.intQuants.xVolume().value(),
1355 ectx.intQuants.yVolume().value());
1356 extboC.assignZFraction(ectx.globalDofIdx,
1357 ectx.intQuants.zFraction().value());
1358
1359 const Scalar stdVolOil = getValue(ectx.fs.saturation(oilPhaseIdx)) *
1360 getValue(ectx.fs.invB(oilPhaseIdx)) +
1361 getValue(ectx.fs.saturation(gasPhaseIdx)) *
1362 getValue(ectx.fs.invB(gasPhaseIdx)) *
1363 getValue(ectx.fs.Rv());
1364 const Scalar stdVolGas = getValue(ectx.fs.saturation(gasPhaseIdx)) *
1365 getValue(ectx.fs.invB(gasPhaseIdx)) *
1366 (1.0 - ectx.intQuants.yVolume().value()) +
1367 getValue(ectx.fs.saturation(oilPhaseIdx)) *
1368 getValue(ectx.fs.invB(oilPhaseIdx)) *
1369 getValue(ectx.fs.Rs()) *
1370 (1.0 - ectx.intQuants.xVolume().value());
1371 const Scalar stdVolCo2 = getValue(ectx.fs.saturation(gasPhaseIdx)) *
1372 getValue(ectx.fs.invB(gasPhaseIdx)) *
1373 ectx.intQuants.yVolume().value() +
1374 getValue(ectx.fs.saturation(oilPhaseIdx)) *
1375 getValue(ectx.fs.invB(oilPhaseIdx)) *
1376 getValue(ectx.fs.Rs()) *
1377 ectx.intQuants.xVolume().value();
1378 const Scalar rhoO = FluidSystem::referenceDensity(gasPhaseIdx, ectx.pvtRegionIdx);
1379 const Scalar rhoG = FluidSystem::referenceDensity(gasPhaseIdx, ectx.pvtRegionIdx);
1380 const Scalar rhoCO2 = ectx.intQuants.zRefDensity();
1381 const Scalar stdMassTotal = 1.0e-10 + stdVolOil * rhoO + stdVolGas * rhoG + stdVolCo2 * rhoCO2;
1382 extboC.assignMassFractions(ectx.globalDofIdx,
1386 }, this->extboC_.allocated()
1387 },
1388 Entry{[&micpC = this->micpC_](const Context& ectx)
1389 {
1390 micpC.assign(ectx.globalDofIdx,
1391 ectx.intQuants.microbialConcentration().value(),
1392 ectx.intQuants.oxygenConcentration().value(),
1393 // Rescaling back the urea concentration (see WellInterface_impl.hpp)
1394 10 * ectx.intQuants.ureaConcentration().value(),
1395 ectx.intQuants.biofilmConcentration().value(),
1396 ectx.intQuants.calciteConcentration().value());
1397 }, this->micpC_.allocated()
1398 },
1399 Entry{[&rftC = this->rftC_,
1400 &vanguard = this->simulator_.vanguard()](const Context& ectx)
1401 {
1402 const auto cartesianIdx = vanguard.cartesianIndex(ectx.globalDofIdx);
1403 rftC.assign(cartesianIdx,
1404 [&fs = ectx.fs]() { return getValue(fs.pressure(oilPhaseIdx)); },
1405 [&fs = ectx.fs]() { return getValue(fs.saturation(waterPhaseIdx)); },
1406 [&fs = ectx.fs]() { return getValue(fs.saturation(gasPhaseIdx)); });
1407 }
1408 },
1409 Entry{[&tC = this->tracerC_,
1410 &tM = this->simulator_.problem().tracerModel()](const Context& ectx)
1411 {
1412 tC.assignFreeConcentrations(ectx.globalDofIdx,
1413 [gIdx = ectx.globalDofIdx, &tM](const unsigned tracerIdx)
1414 { return tM.freeTracerConcentration(tracerIdx, gIdx); });
1415 tC.assignSolConcentrations(ectx.globalDofIdx,
1416 [gIdx = ectx.globalDofIdx, &tM](const unsigned tracerIdx)
1417 { return tM.solTracerConcentration(tracerIdx, gIdx); });
1418 }
1419 },
1420 Entry{[&flowsInf = this->simulator_.problem().model().linearizer().getFlowsInfo(),
1421 &flowsC = this->flowsC_](const Context& ectx)
1422 {
1423 constexpr auto gas_idx = Indices::gasEnabled ?
1424 conti0EqIdx + Indices::canonicalToActiveComponentIndex(gasCompIdx) : -1;
1425 constexpr auto oil_idx = Indices::oilEnabled ?
1426 conti0EqIdx + Indices::canonicalToActiveComponentIndex(oilCompIdx) : -1;
1427 constexpr auto water_idx = Indices::waterEnabled ?
1428 conti0EqIdx + Indices::canonicalToActiveComponentIndex(waterCompIdx) : -1;
1429 const auto& flowsInfos = flowsInf[ectx.globalDofIdx];
1430 for (const auto& flowsInfo : flowsInfos) {
1431 flowsC.assignFlows(ectx.globalDofIdx,
1432 flowsInfo.faceId,
1433 flowsInfo.nncId,
1437 }
1438 }, !this->simulator_.problem().model().linearizer().getFlowsInfo().empty()
1439 },
1440 Entry{[&floresInf = this->simulator_.problem().model().linearizer().getFloresInfo(),
1441 &flowsC = this->flowsC_](const Context& ectx)
1442 {
1443 constexpr auto gas_idx = Indices::gasEnabled ?
1444 conti0EqIdx + Indices::canonicalToActiveComponentIndex(gasCompIdx) : -1;
1445 constexpr auto oil_idx = Indices::oilEnabled ?
1446 conti0EqIdx + Indices::canonicalToActiveComponentIndex(oilCompIdx) : -1;
1447 constexpr auto water_idx = Indices::waterEnabled ?
1448 conti0EqIdx + Indices::canonicalToActiveComponentIndex(waterCompIdx) : -1;
1449 const auto& floresInfos = floresInf[ectx.globalDofIdx];
1450 for (const auto& floresInfo : floresInfos) {
1451 flowsC.assignFlores(ectx.globalDofIdx,
1452 floresInfo.faceId,
1453 floresInfo.nncId,
1457 }
1458 }, !this->simulator_.problem().model().linearizer().getFloresInfo().empty()
1459 },
1460 // hack to make the intial output of rs and rv Ecl compatible.
1461 // For cells with swat == 1 Ecl outputs; rs = rsSat and rv=rvSat, in all but the initial step
1462 // where it outputs rs and rv values calculated by the initialization. To be compatible we overwrite
1463 // rs and rv with the values computed in the initially.
1464 // Volume factors, densities and viscosities need to be recalculated with the updated rs and rv values.
1465 Entry{ScalarEntry{&this->rv_,
1466 [&problem = this->simulator_.problem()](const Context& ectx)
1467 { return problem.initialFluidState(ectx.globalDofIdx).Rv(); }
1468 },
1469 simulator_.episodeIndex() < 0 &&
1470 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1471 FluidSystem::phaseIsActive(gasPhaseIdx)
1472 },
1473 Entry{ScalarEntry{&this->rs_,
1474 [&problem = this->simulator_.problem()](const Context& ectx)
1475 { return problem.initialFluidState(ectx.globalDofIdx).Rs(); }
1476 },
1477 simulator_.episodeIndex() < 0 &&
1478 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1479 FluidSystem::phaseIsActive(gasPhaseIdx)
1480 },
1481 Entry{ScalarEntry{&this->rsw_,
1482 [&problem = this->simulator_.problem()](const Context& ectx)
1483 { return problem.initialFluidState(ectx.globalDofIdx).Rsw(); }
1484 },
1485 simulator_.episodeIndex() < 0 &&
1486 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1487 FluidSystem::phaseIsActive(gasPhaseIdx)
1488 },
1489 Entry{ScalarEntry{&this->rvw_,
1490 [&problem = this->simulator_.problem()](const Context& ectx)
1491 { return problem.initialFluidState(ectx.globalDofIdx).Rvw(); }
1492 },
1493 simulator_.episodeIndex() < 0 &&
1494 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1495 FluidSystem::phaseIsActive(gasPhaseIdx)
1496 },
1497 // re-compute the volume factors, viscosities and densities if asked for
1498 Entry{PhaseEntry{&this->density_,
1499 [&problem = this->simulator_.problem()](const unsigned phase,
1500 const Context& ectx)
1501 {
1502 const auto& fsInitial = problem.initialFluidState(ectx.globalDofIdx);
1503 return FluidSystem::density(fsInitial,
1504 phase,
1505 ectx.intQuants.pvtRegionIndex());
1506 }
1507 },
1508 simulator_.episodeIndex() < 0 &&
1509 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1510 FluidSystem::phaseIsActive(gasPhaseIdx)
1511 },
1512 Entry{PhaseEntry{&this->invB_,
1513 [&problem = this->simulator_.problem()](const unsigned phase,
1514 const Context& ectx)
1515 {
1516 const auto& fsInitial = problem.initialFluidState(ectx.globalDofIdx);
1517 return FluidSystem::inverseFormationVolumeFactor(fsInitial,
1518 phase,
1519 ectx.intQuants.pvtRegionIndex());
1520 }
1521 },
1522 simulator_.episodeIndex() < 0 &&
1523 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1524 FluidSystem::phaseIsActive(gasPhaseIdx)
1525 },
1526 Entry{PhaseEntry{&this->viscosity_,
1527 [&problem = this->simulator_.problem()](const unsigned phase,
1528 const Context& ectx)
1529 {
1530 const auto& fsInitial = problem.initialFluidState(ectx.globalDofIdx);
1531 return FluidSystem::viscosity(fsInitial,
1532 phase,
1533 ectx.intQuants.pvtRegionIndex());
1534 }
1535 },
1536 simulator_.episodeIndex() < 0 &&
1537 FluidSystem::phaseIsActive(oilPhaseIdx) &&
1538 FluidSystem::phaseIsActive(gasPhaseIdx)
1539 },
1540 };
1541
1542 // Setup active extractors
1543 this->extractors_ = Extractor::removeInactive(extractors);
1544
1546 if (this->mech_.allocated()) {
1547 this->extractors_.emplace(
1548 [&mech = this->mech_,
1549 &model = simulator_.problem().geoMechModel()](const Context& ectx)
1550 {
1551 mech.assignDelStress(ectx.globalDofIdx,
1552 model.delstress(ectx.globalDofIdx));
1553
1554 mech.assignDisplacement(ectx.globalDofIdx,
1555 model.disp(ectx.globalDofIdx, /*include_fracture*/true));
1556
1557 // is the tresagii stress which make rock fracture
1558 mech.assignFracStress(ectx.globalDofIdx,
1559 model.fractureStress(ectx.globalDofIdx));
1560
1561 mech.assignLinStress(ectx.globalDofIdx,
1562 model.linstress(ectx.globalDofIdx));
1563
1564 mech.assignPotentialForces(ectx.globalDofIdx,
1565 model.mechPotentialForce(ectx.globalDofIdx),
1566 model.mechPotentialPressForce(ectx.globalDofIdx),
1567 model.mechPotentialTempForce(ectx.globalDofIdx));
1568
1569 mech.assignStrain(ectx.globalDofIdx,
1570 model.strain(ectx.globalDofIdx, /*include_fracture*/true));
1571
1572 // Total stress is not stored but calculated result is Voigt notation
1573 mech.assignStress(ectx.globalDofIdx,
1574 model.stress(ectx.globalDofIdx, /*include_fracture*/true));
1575 }, true
1576 );
1577 }
1578 }
1579 }
1580
1582 void setupBlockExtractors_(const bool isSubStep,
1583 const int reportStepNum)
1584 {
1585 using Entry = typename BlockExtractor::Entry;
1586 using Context = typename BlockExtractor::Context;
1587 using PhaseEntry = typename BlockExtractor::PhaseEntry;
1588 using ScalarEntry = typename BlockExtractor::ScalarEntry;
1589
1590 using namespace std::string_view_literals;
1591
1592 const auto pressure_handler =
1593 Entry{ScalarEntry{std::vector{"BPR"sv, "BPRESSUR"sv},
1594 [](const Context& ectx)
1595 {
1596 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1597 return getValue(ectx.fs.pressure(oilPhaseIdx));
1598 }
1599 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1600 return getValue(ectx.fs.pressure(gasPhaseIdx));
1601 }
1602 else { //if (FluidSystem::phaseIsActive(waterPhaseIdx))
1603 return getValue(ectx.fs.pressure(waterPhaseIdx));
1604 }
1605 }
1606 }
1607 };
1608
1609 const auto handlers = std::array{
1611 Entry{PhaseEntry{std::array{
1612 std::array{"BWSAT"sv, "BOSAT"sv, "BGSAT"sv},
1613 std::array{"BSWAT"sv, "BSOIL"sv, "BSGAS"sv}
1614 },
1615 [](const unsigned phaseIdx, const Context& ectx)
1616 {
1617 return getValue(ectx.fs.saturation(phaseIdx));
1618 }
1619 }
1620 },
1621 Entry{ScalarEntry{"BNSAT",
1622 [](const Context& ectx)
1623 {
1624 return ectx.intQuants.solventSaturation().value();
1625 }
1626 }
1627 },
1628 Entry{ScalarEntry{std::vector{"BTCNFHEA"sv, "BTEMP"sv},
1629 [](const Context& ectx)
1630 {
1631 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1632 return getValue(ectx.fs.temperature(oilPhaseIdx));
1633 }
1634 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1635 return getValue(ectx.fs.temperature(gasPhaseIdx));
1636 }
1637 else { //if (FluidSystem::phaseIsActive(waterPhaseIdx))
1638 return getValue(ectx.fs.temperature(waterPhaseIdx));
1639 }
1640 }
1641 }
1642 },
1643 Entry{PhaseEntry{std::array{
1644 std::array{"BWKR"sv, "BOKR"sv, "BGKR"sv},
1645 std::array{"BKRW"sv, "BKRO"sv, "BKRG"sv}
1646 },
1647 [](const unsigned phaseIdx, const Context& ectx)
1648 {
1649 return getValue(ectx.intQuants.relativePermeability(phaseIdx));
1650 }
1651 }
1652 },
1653 Entry{ScalarEntry{"BKROG",
1654 [&problem = this->simulator_.problem()](const Context& ectx)
1655 {
1656 const auto& materialParams =
1657 problem.materialLawParams(ectx.elemCtx,
1658 ectx.dofIdx,
1659 /* timeIdx = */ 0);
1660 return getValue(MaterialLaw::template
1662 ectx.fs));
1663 }
1664 }
1665 },
1666 Entry{ScalarEntry{"BKROW",
1667 [&problem = this->simulator_.problem()](const Context& ectx)
1668 {
1669 const auto& materialParams = problem.materialLawParams(ectx.elemCtx,
1670 ectx.dofIdx,
1671 /* timeIdx = */ 0);
1672 return getValue(MaterialLaw::template
1674 ectx.fs));
1675 }
1676 }
1677 },
1678 Entry{ScalarEntry{"BWPC",
1679 [](const Context& ectx)
1680 {
1681 return getValue(ectx.fs.pressure(oilPhaseIdx)) -
1682 getValue(ectx.fs.pressure(waterPhaseIdx));
1683 }
1684 }
1685 },
1686 Entry{ScalarEntry{"BGPC",
1687 [](const Context& ectx)
1688 {
1689 return getValue(ectx.fs.pressure(gasPhaseIdx)) -
1690 getValue(ectx.fs.pressure(oilPhaseIdx));
1691 }
1692 }
1693 },
1694 Entry{ScalarEntry{"BWPR",
1695 [](const Context& ectx)
1696 {
1697 return getValue(ectx.fs.pressure(waterPhaseIdx));
1698 }
1699 }
1700 },
1701 Entry{ScalarEntry{"BGPR",
1702 [](const Context& ectx)
1703 {
1704 return getValue(ectx.fs.pressure(gasPhaseIdx));
1705 }
1706 }
1707 },
1708 Entry{PhaseEntry{std::array{
1709 std::array{"BVWAT"sv, "BVOIL"sv, "BVGAS"sv},
1710 std::array{"BWVIS"sv, "BOVIS"sv, "BGVIS"sv}
1711 },
1712 [](const unsigned phaseIdx, const Context& ectx)
1713 {
1714 return getValue(ectx.fs.viscosity(phaseIdx));
1715 }
1716 }
1717 },
1718 Entry{PhaseEntry{std::array{
1719 std::array{"BWDEN"sv, "BODEN"sv, "BGDEN"sv},
1720 std::array{"BDENW"sv, "BDENO"sv, "BDENG"sv}
1721 },
1722 [](const unsigned phaseIdx, const Context& ectx)
1723 {
1724 return getValue(ectx.fs.density(phaseIdx));
1725 }
1726 }
1727 },
1728 Entry{ScalarEntry{"BFLOWI",
1729 [&flowsC = this->flowsC_](const Context& ectx)
1730 {
1731 return flowsC.getFlow(ectx.globalDofIdx, Dir::XPlus, waterCompIdx);
1732 }
1733 }
1734 },
1735 Entry{ScalarEntry{"BFLOWJ",
1736 [&flowsC = this->flowsC_](const Context& ectx)
1737 {
1738 return flowsC.getFlow(ectx.globalDofIdx, Dir::YPlus, waterCompIdx);
1739 }
1740 }
1741 },
1742 Entry{ScalarEntry{"BFLOWK",
1743 [&flowsC = this->flowsC_](const Context& ectx)
1744 {
1745 return flowsC.getFlow(ectx.globalDofIdx, Dir::ZPlus, waterCompIdx);
1746 }
1747 }
1748 },
1749 Entry{ScalarEntry{"BRPV",
1750 [&model = this->simulator_.model()](const Context& ectx)
1751 {
1752 return getValue(ectx.intQuants.porosity()) *
1753 model.dofTotalVolume(ectx.globalDofIdx);
1754 }
1755 }
1756 },
1757 Entry{PhaseEntry{std::array{"BWPV"sv, "BOPV"sv, "BGPV"sv},
1758 [&model = this->simulator_.model()](const unsigned phaseIdx,
1759 const Context& ectx)
1760 {
1761 return getValue(ectx.fs.saturation(phaseIdx)) *
1762 getValue(ectx.intQuants.porosity()) *
1763 model.dofTotalVolume(ectx.globalDofIdx);
1764 }
1765 }
1766 },
1767 Entry{ScalarEntry{"BRS",
1768 [](const Context& ectx)
1769 {
1770 return getValue(ectx.fs.Rs());
1771 }
1772 }
1773 },
1774 Entry{ScalarEntry{"BRV",
1775 [](const Context& ectx)
1776 {
1777 return getValue(ectx.fs.Rv());
1778 }
1779 }
1780 },
1781 Entry{ScalarEntry{"BOIP",
1782 [&model = this->simulator_.model()](const Context& ectx)
1783 {
1784 return (getValue(ectx.fs.invB(oilPhaseIdx)) *
1785 getValue(ectx.fs.saturation(oilPhaseIdx)) +
1786 getValue(ectx.fs.Rv()) *
1787 getValue(ectx.fs.invB(gasPhaseIdx)) *
1788 getValue(ectx.fs.saturation(gasPhaseIdx))) *
1789 model.dofTotalVolume(ectx.globalDofIdx) *
1790 getValue(ectx.intQuants.porosity());
1791 }
1792 }
1793 },
1794 Entry{ScalarEntry{"BGIP",
1795 [&model = this->simulator_.model()](const Context& ectx)
1796 {
1797 Scalar result = getValue(ectx.fs.invB(gasPhaseIdx)) *
1798 getValue(ectx.fs.saturation(gasPhaseIdx));
1799
1800 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1801 result += getValue(ectx.fs.Rs()) *
1802 getValue(ectx.fs.invB(oilPhaseIdx)) *
1803 getValue(ectx.fs.saturation(oilPhaseIdx));
1804 }
1805 else {
1806 result += getValue(ectx.fs.Rsw()) *
1807 getValue(ectx.fs.invB(waterPhaseIdx)) *
1808 getValue(ectx.fs.saturation(waterPhaseIdx));
1809 }
1810
1811 return result *
1812 model.dofTotalVolume(ectx.globalDofIdx) *
1813 getValue(ectx.intQuants.porosity());
1814 }
1815 }
1816 },
1817 Entry{ScalarEntry{"BWIP",
1818 [&model = this->simulator_.model()](const Context& ectx)
1819 {
1820 return getValue(ectx.fs.invB(waterPhaseIdx)) *
1821 getValue(ectx.fs.saturation(waterPhaseIdx)) *
1822 model.dofTotalVolume(ectx.globalDofIdx) *
1823 getValue(ectx.intQuants.porosity());
1824 }
1825 }
1826 },
1827 Entry{ScalarEntry{"BOIPL",
1828 [&model = this->simulator_.model()](const Context& ectx)
1829 {
1830 return getValue(ectx.fs.invB(oilPhaseIdx)) *
1831 getValue(ectx.fs.saturation(oilPhaseIdx)) *
1832 model.dofTotalVolume(ectx.globalDofIdx) *
1833 getValue(ectx.intQuants.porosity());
1834 }
1835 }
1836 },
1837 Entry{ScalarEntry{"BGIPL",
1838 [&model = this->simulator_.model()](const Context& ectx)
1839 {
1840 Scalar result;
1841 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1842 result = getValue(ectx.fs.Rs()) *
1843 getValue(ectx.fs.invB(oilPhaseIdx)) *
1844 getValue(ectx.fs.saturation(oilPhaseIdx));
1845 }
1846 else {
1847 result = getValue(ectx.fs.Rsw()) *
1848 getValue(ectx.fs.invB(waterPhaseIdx)) *
1849 getValue(ectx.fs.saturation(waterPhaseIdx));
1850 }
1851 return result *
1852 model.dofTotalVolume(ectx.globalDofIdx) *
1853 getValue(ectx.intQuants.porosity());
1854 }
1855 }
1856 },
1857 Entry{ScalarEntry{"BGIPG",
1858 [&model = this->simulator_.model()](const Context& ectx)
1859 {
1860 return getValue(ectx.fs.invB(gasPhaseIdx)) *
1861 getValue(ectx.fs.saturation(gasPhaseIdx)) *
1862 model.dofTotalVolume(ectx.globalDofIdx) *
1863 getValue(ectx.intQuants.porosity());
1864 }
1865 }
1866 },
1867 Entry{ScalarEntry{"BOIPG",
1868 [&model = this->simulator_.model()](const Context& ectx)
1869 {
1870 return getValue(ectx.fs.Rv()) *
1871 getValue(ectx.fs.invB(gasPhaseIdx)) *
1872 getValue(ectx.fs.saturation(gasPhaseIdx)) *
1873 model.dofTotalVolume(ectx.globalDofIdx) *
1874 getValue(ectx.intQuants.porosity());
1875 }
1876 }
1877 },
1878 Entry{PhaseEntry{std::array{"BPPW"sv, "BPPO"sv, "BPPG"sv},
1879 [&simConfig = this->eclState_.getSimulationConfig(),
1880 &grav = this->simulator_.problem().gravity(),
1881 &regionAvgDensity = this->regionAvgDensity_,
1882 &problem = this->simulator_.problem(),
1883 &regions = this->regions_](const unsigned phaseIdx, const Context& ectx)
1884 {
1885 auto phase = RegionPhasePoreVolAverage::Phase{};
1886 phase.ix = phaseIdx;
1887
1888 // Note different region handling here. FIPNUM is
1889 // one-based, but we need zero-based lookup in
1890 // DatumDepth. On the other hand, pvtRegionIndex is
1891 // zero-based but we need one-based lookup in
1892 // RegionPhasePoreVolAverage.
1893
1894 // Subtract one to convert FIPNUM to region index.
1895 const auto datum = simConfig.datumDepths()(regions["FIPNUM"][ectx.dofIdx] - 1);
1896
1897 // Add one to convert region index to region ID.
1898 const auto region = RegionPhasePoreVolAverage::Region {
1899 ectx.elemCtx.primaryVars(ectx.dofIdx, /*timeIdx=*/0).pvtRegionIndex() + 1
1900 };
1901
1902 const auto density = regionAvgDensity->value("PVTNUM", phase, region);
1903
1904 const auto press = getValue(ectx.fs.pressure(phase.ix));
1905 const auto dz = problem.dofCenterDepth(ectx.globalDofIdx) - datum;
1906 return press - density*dz*grav[GridView::dimensionworld - 1];
1907 }
1908 }
1909 },
1910 };
1911
1912 this->blockExtractors_ = BlockExtractor::setupExecMap(this->blockData_, handlers);
1913
1914 this->extraBlockData_.clear();
1915 if (reportStepNum > 0 && !isSubStep) {
1916 // check we need extra block pressures for RPTSCHED
1917 const auto& rpt = this->schedule_[reportStepNum - 1].rpt_config.get();
1918 if (rpt.contains("WELLS") && rpt.at("WELLS") > 1) {
1919 this->setupExtraBlockData(reportStepNum,
1920 [&c = this->collectOnIORank_](const int idx)
1921 { return c.isCartIdxOnThisRank(idx); });
1922
1923 const auto extraHandlers = std::array{
1925 };
1926
1927 this->extraBlockExtractors_ = BlockExtractor::setupExecMap(this->extraBlockData_, extraHandlers);
1928 }
1929 }
1930 }
1931
1932 const Simulator& simulator_;
1933 const CollectDataOnIORankType& collectOnIORank_;
1934 std::vector<typename Extractor::Entry> extractors_;
1935 typename BlockExtractor::ExecMap blockExtractors_;
1936 typename BlockExtractor::ExecMap extraBlockExtractors_;
1937};
1938
1939} // namespace Opm
1940
1941#endif // OPM_OUTPUT_BLACK_OIL_MODULE_HPP
Helper class for grid instantiation of ECL file-format using problems.
Output module for the results black oil model writing in ECL binary format.
Output module for the results black oil model writing in ECL binary format.
Declares the properties required by the black oil model.
Definition CollectDataOnIORank.hpp:56
The base class for the element-centered finite-volume discretization scheme.
Definition ecfvdiscretization.hh:147
Definition GenericOutputBlackoilModule.hpp:76
Inter-region flow accumulation maps for all region definition arrays.
Definition InterRegFlows.hpp:179
void compress()
Form CSR adjacency matrix representation of input graph from connections established in previous call...
Definition InterRegFlows.cpp:165
void addConnection(const Cell &source, const Cell &destination, const data::InterRegFlowMap::FlowRates &rates)
Add flow rate connection between regions for all region definitions.
Definition InterRegFlows.cpp:156
void clear()
Clear all internal buffers, but preserve allocated capacity.
Definition InterRegFlows.cpp:172
Output module for the results black oil model writing in ECL binary format.
Definition OutputBlackoilModule.hpp:85
void processElement(const ElementContext &elemCtx)
Modify the internal buffers according to the intensive quanties relevant for an element.
Definition OutputBlackoilModule.hpp:233
void initializeFluxData()
Prepare for capturing connection fluxes, particularly to account for inter-region flows.
Definition OutputBlackoilModule.hpp:389
void setupExtractors(const bool isSubStep, const int reportStepNum)
Setup list of active element-level data extractors.
Definition OutputBlackoilModule.hpp:214
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:192
void processFluxes(const ElementContext &elemCtx, ActiveIndex &&activeIndex, CartesianIndex &&cartesianIndex)
Capture connection fluxes, particularly to account for inter-region flows.
Definition OutputBlackoilModule.hpp:352
void clearExtractors()
Clear list of active element-level data extractors.
Definition OutputBlackoilModule.hpp:222
const InterRegFlowMap & getInterRegFlows() const
Get read-only access to collection of inter-region flows.
Definition OutputBlackoilModule.hpp:407
void finalizeFluxData()
Finalize capturing connection fluxes.
Definition OutputBlackoilModule.hpp:399
Declare the properties used by the infrastructure code of the finite volume discretizations.
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
std::string moduleVersionName()
Return the version name of the module, for example "2015.10" (for a release branch) or "2016....
Definition moduleVersion.cpp:34
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition propertysystem.hh:235
This file provides the infrastructure to retrieve run-time parameters.
The Opm property system, traits with inheritance.
Minimal characteristics of a cell from a simulation grid.
Definition InterRegFlows.hpp:50
Wrapping struct holding types used for block-level data extraction.
Definition OutputExtractor.hpp:193
std::unordered_map< int, std::vector< Exec > > ExecMap
A map of extraction executors, keyed by cartesian cell index.
Definition OutputExtractor.hpp:260
std::variant< ScalarEntry, PhaseEntry > Entry
Descriptor for extractors.
Definition OutputExtractor.hpp:245
static ExecMap setupExecMap(std::map< std::pair< std::string, int >, double > &blockData, const std::array< Entry, size > &handlers)
Setup an extractor executor map from a map of evaluations to perform.
Definition OutputExtractor.hpp:264
static void process(const std::vector< Exec > &blockExtractors, const Context &ectx)
Process a list of block extractors.
Definition OutputExtractor.hpp:367
Context passed to extractor functions.
Definition OutputExtractor.hpp:74
int episodeIndex
Current report step.
Definition OutputExtractor.hpp:77
Struct holding hysteresis parameters.
Definition OutputExtractor.hpp:63
Wrapping struct holding types used for element-level data extraction.
Definition OutputExtractor.hpp:54
static void process(const Context &ectx, const std::vector< Entry > &extractors)
Process the given extractor entries.
Definition OutputExtractor.hpp:158
static std::vector< Entry > removeInactive(std::array< Entry, size > &input)
Obtain vector of active extractors from an array of extractors.
Definition OutputExtractor.hpp:120