opm-simulators
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 
46 #include <opm/models/blackoil/blackoilenergymodules.hh>
47 #include <opm/models/blackoil/blackoilproperties.hh>
48 #include <opm/models/discretization/common/fvbaseproperties.hh>
49 #include <opm/models/utils/parametersystem.hpp>
50 #include <opm/models/utils/propertysystem.hh>
51 
52 #include <opm/output/data/Cells.hpp>
53 #include <opm/output/eclipse/EclipseIO.hpp>
54 #include <opm/output/eclipse/Inplace.hpp>
55 
56 #include <opm/simulators/flow/CollectDataOnIORank.hpp>
57 #include <opm/simulators/flow/FlowBaseVanguard.hpp>
58 #include <opm/simulators/flow/GenericOutputBlackoilModule.hpp>
59 #include <opm/simulators/flow/OutputExtractor.hpp>
60 
61 #include <algorithm>
62 #include <array>
63 #include <cassert>
64 #include <cstddef>
65 #include <functional>
66 #include <stdexcept>
67 #include <string>
68 #include <type_traits>
69 #include <utility>
70 #include <vector>
71 
72 namespace Opm {
73 
74 // forward declaration
75 template <class TypeTag>
76 class EcfvDiscretization;
77 
84 template <class TypeTag>
85 class OutputBlackOilModule : public GenericOutputBlackoilModule<GetPropType<TypeTag, Properties::FluidSystem>>
86 {
87  using Simulator = GetPropType<TypeTag, Properties::Simulator>;
88  using Discretization = GetPropType<TypeTag, Properties::Discretization>;
89  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
90  using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
91  using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
92  using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
93  using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
94  using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
95  using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
96  using FluidState = typename IntensiveQuantities::FluidState;
97  using GridView = GetPropType<TypeTag, Properties::GridView>;
98  using Element = typename GridView::template Codim<0>::Entity;
99  using ElementIterator = typename GridView::template Codim<0>::Iterator;
100  using BaseType = GenericOutputBlackoilModule<FluidSystem>;
101  using Indices = GetPropType<TypeTag, Properties::Indices>;
102  using Dir = FaceDir::DirEnum;
103  using BlockExtractor = detail::BlockExtractor<TypeTag>;
104  using Extractor = detail::Extractor<TypeTag>;
105  using Grid = GetPropType<TypeTag, Properties::Grid>;
106  using EquilGrid = GetPropType<TypeTag, Properties::EquilGrid>;
107  using CollectDataOnIORankType = CollectDataOnIORank<Grid,EquilGrid,GridView>;
108 
109  static constexpr int conti0EqIdx = Indices::conti0EqIdx;
110  static constexpr int numPhases = FluidSystem::numPhases;
111  static constexpr int oilPhaseIdx = FluidSystem::oilPhaseIdx;
112  static constexpr int gasPhaseIdx = FluidSystem::gasPhaseIdx;
113  static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
114  static constexpr int gasCompIdx = FluidSystem::gasCompIdx;
115  static constexpr int oilCompIdx = FluidSystem::oilCompIdx;
116  static constexpr int waterCompIdx = FluidSystem::waterCompIdx;
117  static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
118  enum { enableBioeffects = getPropValue<TypeTag, Properties::EnableBioeffects>() };
119  enum { enableMICP = Indices::enableMICP };
120  enum { enableVapwat = getPropValue<TypeTag, Properties::EnableVapwat>() };
121  enum { enableDisgasInWater = getPropValue<TypeTag, Properties::EnableDisgasInWater>() };
122  enum { enableDissolvedGas = Indices::compositionSwitchIdx >= 0 };
123 
124  template<class VectorType>
125  static Scalar value_or_zero(int idx, const VectorType& v)
126  {
127  if (idx == -1) {
128  return 0.0;
129  }
130  return v.empty() ? 0.0 : v[idx];
131  }
132 
133 public:
134  OutputBlackOilModule(const Simulator& simulator,
135  const SummaryConfig& smryCfg,
136  const CollectDataOnIORankType& collectOnIORank)
137  : BaseType(simulator.vanguard().eclState(),
138  simulator.vanguard().schedule(),
139  smryCfg,
140  simulator.vanguard().summaryState(),
141  moduleVersionName(),
142  [this](const int idx)
143  { return simulator_.problem().eclWriter().collectOnIORank().localIdxToGlobalIdx(idx); },
144  [&collectOnIORank](const int idx)
145  { return collectOnIORank.isCartIdxOnThisRank(idx); },
146  simulator.vanguard().grid().comm(),
147  energyModuleType == EnergyModules::FullyImplicitThermal ||
148  energyModuleType == EnergyModules::SequentialImplicitThermal,
149  energyModuleType == EnergyModules::ConstantTemperature,
150  getPropValue<TypeTag, Properties::EnableMech>(),
151  getPropValue<TypeTag, Properties::EnableSolvent>(),
152  getPropValue<TypeTag, Properties::EnablePolymer>(),
153  getPropValue<TypeTag, Properties::EnableFoam>(),
154  getPropValue<TypeTag, Properties::EnableBrine>(),
155  getPropValue<TypeTag, Properties::EnableSaltPrecipitation>(),
156  getPropValue<TypeTag, Properties::EnableExtbo>(),
157  getPropValue<TypeTag, Properties::EnableBioeffects>(),
158  getPropValue<TypeTag, Properties::EnableGeochemistry>())
159  , simulator_(simulator)
160  , collectOnIORank_(collectOnIORank)
161  {
162  for (auto& region_pair : this->regions_) {
163  this->createLocalRegion_(region_pair.second);
164  }
165 
166  auto isCartIdxOnThisRank = [&collectOnIORank](const int idx) {
167  return collectOnIORank.isCartIdxOnThisRank(idx);
168  };
169 
170  this->setupBlockData(isCartIdxOnThisRank);
171 
172  if (! Parameters::Get<Parameters::OwnerCellsFirst>()) {
173  const std::string msg = "The output code does not support --owner-cells-first=false.";
174  if (collectOnIORank.isIORank()) {
175  OpmLog::error(msg);
176  }
177  OPM_THROW_NOLOG(std::runtime_error, msg);
178  }
179 
180  if (smryCfg.match("[FB]PP[OGW]") || smryCfg.match("RPP[OGW]*")) {
181  auto rset = this->eclState_.fieldProps().fip_regions();
182  rset.push_back("PVTNUM");
183 
184  // Note: We explicitly use decltype(auto) here because the
185  // default scheme (-> auto) will deduce an undesirable type. We
186  // need the "reference to vector" semantics in this instance.
187  this->regionAvgDensity_
188  .emplace(this->simulator_.gridView().comm(),
189  FluidSystem::numPhases, rset,
190  [fp = std::cref(this->eclState_.fieldProps())]
191  (const std::string& rsetName) -> decltype(auto)
192  { return fp.get().get_int(rsetName); });
193  }
194  }
195 
200  void
201  allocBuffers(const unsigned bufferSize,
202  const unsigned reportStepNum,
203  const bool substep,
204  const bool log,
205  const bool isRestart)
206  {
207  if (! std::is_same<Discretization, EcfvDiscretization<TypeTag>>::value) {
208  return;
209  }
210 
211  const auto& problem = this->simulator_.problem();
212 
213  this->doAllocBuffers(bufferSize,
214  reportStepNum,
215  substep,
216  log,
217  isRestart,
218  &problem.materialLawManager()->hysteresisConfig(),
219  problem.eclWriter().getOutputNnc().front().size());
220  }
221 
223  void setupExtractors(const bool isSubStep,
224  const int reportStepNum)
225  {
226  this->setupElementExtractors_();
227  this->setupBlockExtractors_(isSubStep, reportStepNum);
228  }
229 
231  void clearExtractors()
232  {
233  this->extractors_.clear();
234  this->blockExtractors_.clear();
235  this->extraBlockExtractors_.clear();
236  }
237 
242  void processElement(const ElementContext& elemCtx)
243  {
244  OPM_TIMEBLOCK_LOCAL(processElement, Subsystem::Output);
245  if (!std::is_same<Discretization, EcfvDiscretization<TypeTag>>::value) {
246  return;
247  }
248 
249  if (this->extractors_.empty()) {
250  assert(0);
251  }
252 
253  const auto& matLawManager = simulator_.problem().materialLawManager();
254 
255  typename Extractor::HysteresisParams hysterParams;
256  for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
257  const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
258  const auto& fs = intQuants.fluidState();
259 
260  const typename Extractor::Context ectx{
261  elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0),
262  elemCtx.primaryVars(dofIdx, /*timeIdx=*/0).pvtRegionIndex(),
263  elemCtx.simulator().episodeIndex(),
264  fs,
265  intQuants,
266  hysterParams
267  };
268 
269  if (matLawManager->enableHysteresis()) {
270  if (FluidSystem::phaseIsActive(oilPhaseIdx) && FluidSystem::phaseIsActive(waterPhaseIdx)) {
271  matLawManager->oilWaterHysteresisParams(hysterParams.somax,
272  hysterParams.swmax,
273  hysterParams.swmin,
274  ectx.globalDofIdx);
275  }
276  if (FluidSystem::phaseIsActive(oilPhaseIdx) && FluidSystem::phaseIsActive(gasPhaseIdx)) {
277  matLawManager->gasOilHysteresisParams(hysterParams.sgmax,
278  hysterParams.shmax,
279  hysterParams.somin,
280  ectx.globalDofIdx);
281  }
282  }
283 
284  Extractor::process(ectx, extractors_);
285  }
286  }
287 
288  void processElementBlockData(const ElementContext& elemCtx)
289  {
290  OPM_TIMEBLOCK_LOCAL(processElementBlockData, Subsystem::Output);
291  if (!std::is_same<Discretization, EcfvDiscretization<TypeTag>>::value) {
292  return;
293  }
294 
295  if (this->blockExtractors_.empty() && this->extraBlockExtractors_.empty()) {
296  return;
297  }
298 
299  for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
300  // Adding block data
301  const auto globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
302  const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx);
303 
304  const auto be_it = this->blockExtractors_.find(cartesianIdx);
305  const auto bee_it = this->extraBlockExtractors_.find(cartesianIdx);
306  if (be_it == this->blockExtractors_.end() &&
307  bee_it == this->extraBlockExtractors_.end())
308  {
309  continue;
310  }
311 
312  const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
313  const auto& fs = intQuants.fluidState();
314 
315  const typename BlockExtractor::Context ectx{
316  globalDofIdx,
317  dofIdx,
318  fs,
319  intQuants,
320  elemCtx,
321  };
322 
323  if (be_it != this->blockExtractors_.end()) {
324  BlockExtractor::process(be_it->second, ectx);
325  }
326  if (bee_it != this->extraBlockExtractors_.end()) {
327  BlockExtractor::process(bee_it->second, ectx);
328  }
329  }
330  }
331 
332  void outputFipAndResvLog(const Inplace& inplace,
333  const std::size_t reportStepNum,
334  double elapsed,
335  boost::posix_time::ptime currentDate,
336  const bool substep,
337  const Parallel::Communication& comm)
338  {
339 
340  if (comm.rank() != 0) {
341  return;
342  }
343 
344  // For report step 0 we use the RPTSOL config, else derive from RPTSCHED
345  std::unique_ptr<FIPConfig> fipSched;
346  if (reportStepNum > 0) {
347  const auto& rpt = this->schedule_[reportStepNum-1].rpt_config.get();
348  fipSched = std::make_unique<FIPConfig>(rpt);
349  }
350  const FIPConfig& fipc = reportStepNum == 0 ? this->eclState_.getEclipseConfig().fip()
351  : *fipSched;
352 
353  if (!substep && !this->forceDisableFipOutput_ && fipc.output(FIPConfig::OutputField::FIELD)) {
354 
355  this->logOutput_.timeStamp("BALANCE", elapsed, reportStepNum, currentDate);
356 
357  const auto& initial_inplace = *this->initialInplace();
358  this->logOutput_.fip(inplace, initial_inplace, "");
359 
360  if (fipc.output(FIPConfig::OutputField::FIPNUM)) {
361  this->logOutput_.fip(inplace, initial_inplace, "FIPNUM");
362 
363  if (fipc.output(FIPConfig::OutputField::RESV))
364  this->logOutput_.fipResv(inplace, "FIPNUM");
365  }
366 
367  if (fipc.output(FIPConfig::OutputField::FIP)) {
368  for (const auto& reg : this->regions_) {
369  if (reg.first != "FIPNUM") {
370  std::ostringstream ss;
371  ss << "BAL" << reg.first.substr(3);
372  this->logOutput_.timeStamp(ss.str(), elapsed, reportStepNum, currentDate);
373  this->logOutput_.fip(inplace, initial_inplace, reg.first);
374 
375  if (fipc.output(FIPConfig::OutputField::RESV))
376  this->logOutput_.fipResv(inplace, reg.first);
377  }
378  }
379  }
380  }
381  }
382 
383  void outputFipAndResvLogToCSV(const std::size_t reportStepNum,
384  const bool substep,
385  const Parallel::Communication& comm)
386  {
387  if (comm.rank() != 0) {
388  return;
389  }
390 
391  if ((reportStepNum == 0) && (!substep) &&
392  (this->schedule_.initialReportConfiguration().has_value()) &&
393  (this->schedule_.initialReportConfiguration()->contains("CSVFIP"))) {
394 
395  std::ostringstream csv_stream;
396 
397  this->logOutput_.csv_header(csv_stream);
398 
399  const auto& initial_inplace = *this->initialInplace();
400 
401  this->logOutput_.fip_csv(csv_stream, initial_inplace, "FIPNUM");
402 
403  for (const auto& reg : this->regions_) {
404  if (reg.first != "FIPNUM") {
405  this->logOutput_.fip_csv(csv_stream, initial_inplace, reg.first);
406  }
407  }
408 
409  const IOConfig& io = this->eclState_.getIOConfig();
410  auto csv_fname = io.getOutputDir() + "/" + io.getBaseName() + ".CSV";
411 
412  std::ofstream outputFile(csv_fname);
413 
414  outputFile << csv_stream.str();
415 
416  outputFile.close();
417  }
418  }
419 
448  template <class ActiveIndex, class CartesianIndex>
449  void processFluxes(const ElementContext& elemCtx,
450  ActiveIndex&& activeIndex,
451  CartesianIndex&& cartesianIndex)
452  {
453  OPM_TIMEBLOCK_LOCAL(processFluxes, Subsystem::Output);
454  const auto identifyCell = [&activeIndex, &cartesianIndex](const Element& elem)
455  -> InterRegFlowMap::Cell
456  {
457  const auto cellIndex = activeIndex(elem);
458 
459  return {
460  static_cast<int>(cellIndex),
461  cartesianIndex(cellIndex),
462  elem.partitionType() == Dune::InteriorEntity
463  };
464  };
465 
466  const auto timeIdx = 0u;
467  const auto& stencil = elemCtx.stencil(timeIdx);
468  const auto numInteriorFaces = elemCtx.numInteriorFaces(timeIdx);
469 
470  for (auto scvfIdx = 0 * numInteriorFaces; scvfIdx < numInteriorFaces; ++scvfIdx) {
471  const auto& face = stencil.interiorFace(scvfIdx);
472  const auto left = identifyCell(stencil.element(face.interiorIndex()));
473  const auto right = identifyCell(stencil.element(face.exteriorIndex()));
474 
475  const auto rates = this->
476  getComponentSurfaceRates(elemCtx, face.area(), scvfIdx, timeIdx);
477 
478  this->interRegionFlows_.addConnection(left, right, rates);
479  }
480  }
481 
486  void initializeFluxData()
487  {
488  // Inter-region flow rates. Note: ".clear()" prepares to accumulate
489  // contributions per bulk connection between FIP regions.
490  this->interRegionFlows_.clear();
491  }
492 
496  void finalizeFluxData()
497  {
498  this->interRegionFlows_.compress();
499  }
500 
504  const InterRegFlowMap& getInterRegFlows() const
505  {
506  return this->interRegionFlows_;
507  }
508 
509  template <class FluidState>
510  void assignToFluidState(FluidState& fs, unsigned elemIdx) const
511  {
512  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
513  if (this->saturation_[phaseIdx].empty())
514  continue;
515 
516  fs.setSaturation(phaseIdx, this->saturation_[phaseIdx][elemIdx]);
517  }
518 
519  if (!this->fluidPressure_.empty()) {
520  // this assumes that capillary pressures only depend on the phase saturations
521  // and possibly on temperature. (this is always the case for ECL problems.)
522  std::array<Scalar, numPhases> pc = {0};
523  const MaterialLawParams& matParams = simulator_.problem().materialLawParams(elemIdx);
524  MaterialLaw::capillaryPressures(pc, matParams, fs);
525  Valgrind::CheckDefined(this->fluidPressure_[elemIdx]);
526  Valgrind::CheckDefined(pc);
527  const auto& pressure = this->fluidPressure_[elemIdx];
528  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
529  if (!FluidSystem::phaseIsActive(phaseIdx))
530  continue;
531 
532  if (Indices::oilEnabled)
533  fs.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx]));
534  else if (Indices::gasEnabled)
535  fs.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx]));
536  else if (Indices::waterEnabled)
537  //single (water) phase
538  fs.setPressure(phaseIdx, pressure);
539  }
540  }
541 
542  if constexpr (energyModuleType != EnergyModules::NoTemperature) {
543  if (!this->temperature_.empty())
544  fs.setTemperature(this->temperature_[elemIdx]);
545  }
546  if constexpr (enableDissolvedGas) {
547  if (!this->rs_.empty())
548  fs.setRs(this->rs_[elemIdx]);
549  if (!this->rv_.empty())
550  fs.setRv(this->rv_[elemIdx]);
551  }
552  if constexpr (enableDisgasInWater) {
553  if (!this->rsw_.empty())
554  fs.setRsw(this->rsw_[elemIdx]);
555  }
556  if constexpr (enableVapwat) {
557  if (!this->rvw_.empty())
558  fs.setRvw(this->rvw_[elemIdx]);
559  }
560  }
561 
562  void initHysteresisParams(Simulator& simulator, unsigned elemIdx) const
563  {
564  if (!this->soMax_.empty())
565  simulator.problem().setMaxOilSaturation(elemIdx, this->soMax_[elemIdx]);
566 
567  if (simulator.problem().materialLawManager()->enableHysteresis()) {
568  auto matLawManager = simulator.problem().materialLawManager();
569 
570  if (FluidSystem::phaseIsActive(oilPhaseIdx)
571  && FluidSystem::phaseIsActive(waterPhaseIdx)) {
572  Scalar somax = 2.0;
573  Scalar swmax = -2.0;
574  Scalar swmin = 2.0;
575 
576  if (matLawManager->enableNonWettingHysteresis()) {
577  if (!this->soMax_.empty()) {
578  somax = this->soMax_[elemIdx];
579  }
580  }
581  if (matLawManager->enableWettingHysteresis()) {
582  if (!this->swMax_.empty()) {
583  swmax = this->swMax_[elemIdx];
584  }
585  }
586  if (matLawManager->enablePCHysteresis()) {
587  if (!this->swmin_.empty()) {
588  swmin = this->swmin_[elemIdx];
589  }
590  }
591  matLawManager->setOilWaterHysteresisParams(
592  somax, swmax, swmin, elemIdx);
593  }
594  if (FluidSystem::phaseIsActive(oilPhaseIdx)
595  && FluidSystem::phaseIsActive(gasPhaseIdx)) {
596  Scalar sgmax = 2.0;
597  Scalar shmax = -2.0;
598  Scalar somin = 2.0;
599 
600  if (matLawManager->enableNonWettingHysteresis()) {
601  if (!this->sgmax_.empty()) {
602  sgmax = this->sgmax_[elemIdx];
603  }
604  }
605  if (matLawManager->enableWettingHysteresis()) {
606  if (!this->shmax_.empty()) {
607  shmax = this->shmax_[elemIdx];
608  }
609  }
610  if (matLawManager->enablePCHysteresis()) {
611  if (!this->somin_.empty()) {
612  somin = this->somin_[elemIdx];
613  }
614  }
615  matLawManager->setGasOilHysteresisParams(
616  sgmax, shmax, somin, elemIdx);
617  }
618 
619  }
620 
621  if (simulator_.vanguard().eclState().fieldProps().has_double("SWATINIT")) {
622  simulator.problem().materialLawManager()
623  ->applyRestartSwatInit(elemIdx, this->ppcw_[elemIdx]);
624  }
625  }
626 
627  void updateFluidInPlace(const ElementContext& elemCtx)
628  {
629  for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
630  updateFluidInPlace_(elemCtx, dofIdx);
631  }
632  }
633 
634  void updateFluidInPlace(const unsigned globalDofIdx,
635  const IntensiveQuantities& intQuants,
636  const double totVolume)
637  {
638  this->updateFluidInPlace_(globalDofIdx, intQuants, totVolume);
639  }
640 
641 private:
642  template <typename T>
643  using RemoveCVR = std::remove_cv_t<std::remove_reference_t<T>>;
644 
645  template <typename, class = void>
646  struct HasGeoMech : public std::false_type {};
647 
648  template <typename Problem>
649  struct HasGeoMech<
650  Problem, std::void_t<decltype(std::declval<Problem>().geoMechModel())>
651  > : public std::true_type {};
652 
653  template <typename, class = void>
654  struct HasGeochemistry : public std::false_type {};
655 
656  template <typename Problem>
657  struct HasGeochemistry<
658  Problem, std::void_t<decltype(std::declval<Problem>().geochemistryModel())>
659  > : public std::true_type {};
660 
661  bool isDefunctParallelWell(const std::string& wname) const override
662  {
663  if (simulator_.gridView().comm().size() == 1)
664  return false;
665  const auto& parallelWells = simulator_.vanguard().parallelWells();
666  std::pair<std::string, bool> value {wname, true};
667  auto candidate = std::lower_bound(parallelWells.begin(), parallelWells.end(), value);
668  return candidate == parallelWells.end() || *candidate != value;
669  }
670 
671  bool isOwnedByCurrentRank(const std::string& wname) const override
672  {
673  return this->simulator_.problem().wellModel().isOwner(wname);
674  }
675 
676  bool isOnCurrentRank(const std::string& wname) const override
677  {
678  return this->simulator_.problem().wellModel().hasLocalCells(wname);
679  }
680 
681  void updateFluidInPlace_(const ElementContext& elemCtx, const unsigned dofIdx)
682  {
683  const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
684  const unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
685  const auto totVolume = elemCtx.simulator().model().dofTotalVolume(globalDofIdx);
686 
687  this->updateFluidInPlace_(globalDofIdx, intQuants, totVolume);
688  }
689 
690  void updateFluidInPlace_(const unsigned globalDofIdx,
691  const IntensiveQuantities& intQuants,
692  const double totVolume)
693  {
694  OPM_TIMEBLOCK_LOCAL(updateFluidInPlace, Subsystem::Output);
695 
696  this->updateTotalVolumesAndPressures_(globalDofIdx, intQuants, totVolume);
697 
698  if (this->computeFip_) {
699  this->updatePhaseInplaceVolumes_(globalDofIdx, intQuants, totVolume);
700  }
701  }
702 
703  void createLocalRegion_(std::vector<int>& region)
704  {
705  // For CpGrid with LGRs, where level zero grid has been distributed,
706  // resize region is needed, since in this case the total amount of
707  // element - per process - in level zero grid and leaf grid do not
708  // coincide, in general.
709  region.resize(simulator_.gridView().size(0));
710  std::size_t elemIdx = 0;
711  for (const auto& elem : elements(simulator_.gridView())) {
712  if (elem.partitionType() != Dune::InteriorEntity) {
713  region[elemIdx] = 0;
714  }
715 
716  ++elemIdx;
717  }
718  }
719 
720  template <typename FluidState>
721  void aggregateAverageDensityContributions_(const FluidState& fs,
722  const unsigned int globalDofIdx,
723  const double porv)
724  {
725  auto pvCellValue = RegionPhasePoreVolAverage::CellValue{};
726  pvCellValue.porv = porv;
727 
728  for (auto phaseIdx = 0*FluidSystem::numPhases;
729  phaseIdx < FluidSystem::numPhases; ++phaseIdx)
730  {
731  if (! FluidSystem::phaseIsActive(phaseIdx)) {
732  continue;
733  }
734 
735  pvCellValue.value = getValue(fs.density(phaseIdx));
736  pvCellValue.sat = getValue(fs.saturation(phaseIdx));
737 
738  this->regionAvgDensity_
739  ->addCell(globalDofIdx,
740  RegionPhasePoreVolAverage::Phase { phaseIdx },
741  pvCellValue);
742  }
743  }
744 
760  data::InterRegFlowMap::FlowRates
761  getComponentSurfaceRates(const ElementContext& elemCtx,
762  const Scalar faceArea,
763  const std::size_t scvfIdx,
764  const std::size_t timeIdx) const
765  {
766  using Component = data::InterRegFlowMap::Component;
767 
768  auto rates = data::InterRegFlowMap::FlowRates {};
769 
770  const auto& extQuant = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
771 
772  const auto alpha = getValue(extQuant.extrusionFactor()) * faceArea;
773 
774  if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
775  const auto& up = elemCtx
776  .intensiveQuantities(extQuant.upstreamIndex(oilPhaseIdx), timeIdx);
777 
778  const auto pvtReg = up.pvtRegionIndex();
779 
780  const auto bO = getValue(getInvB_<FluidSystem, FluidState, Scalar>
781  (up.fluidState(), oilPhaseIdx, pvtReg));
782 
783  const auto qO = alpha * bO * getValue(extQuant.volumeFlux(oilPhaseIdx));
784 
785  rates[Component::Oil] += qO;
786 
787  if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
788  const auto Rs = getValue(
789  BlackOil::getRs_<FluidSystem, FluidState, Scalar>
790  (up.fluidState(), pvtReg));
791 
792  rates[Component::Gas] += qO * Rs;
793  rates[Component::Disgas] += qO * Rs;
794  }
795  }
796 
797  if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
798  const auto& up = elemCtx
799  .intensiveQuantities(extQuant.upstreamIndex(gasPhaseIdx), timeIdx);
800 
801  const auto pvtReg = up.pvtRegionIndex();
802 
803  const auto bG = getValue(getInvB_<FluidSystem, FluidState, Scalar>
804  (up.fluidState(), gasPhaseIdx, pvtReg));
805 
806  const auto qG = alpha * bG * getValue(extQuant.volumeFlux(gasPhaseIdx));
807 
808  rates[Component::Gas] += qG;
809 
810  if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
811  const auto Rv = getValue(
812  BlackOil::getRv_<FluidSystem, FluidState, Scalar>
813  (up.fluidState(), pvtReg));
814 
815  rates[Component::Oil] += qG * Rv;
816  rates[Component::Vapoil] += qG * Rv;
817  }
818  }
819 
820  if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
821  const auto& up = elemCtx
822  .intensiveQuantities(extQuant.upstreamIndex(waterPhaseIdx), timeIdx);
823 
824  const auto pvtReg = up.pvtRegionIndex();
825 
826  const auto bW = getValue(getInvB_<FluidSystem, FluidState, Scalar>
827  (up.fluidState(), waterPhaseIdx, pvtReg));
828 
829  rates[Component::Water] +=
830  alpha * bW * getValue(extQuant.volumeFlux(waterPhaseIdx));
831  }
832 
833  return rates;
834  }
835 
836  template <typename FluidState>
837  Scalar hydroCarbonFraction(const FluidState& fs) const
838  {
839  if (this->eclState_.runspec().co2Storage()) {
840  // CO2 storage: Hydrocarbon volume is full pore-volume.
841  return 1.0;
842  }
843 
844  // Common case. Hydrocarbon volume is fraction occupied by actual
845  // hydrocarbons.
846  auto hydrocarbon = Scalar {0};
847  if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
848  hydrocarbon += getValue(fs.saturation(oilPhaseIdx));
849  }
850 
851  if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
852  hydrocarbon += getValue(fs.saturation(gasPhaseIdx));
853  }
854 
855  return hydrocarbon;
856  }
857 
858  void updateTotalVolumesAndPressures_(const unsigned globalDofIdx,
859  const IntensiveQuantities& intQuants,
860  const double totVolume)
861  {
862  const auto& fs = intQuants.fluidState();
863 
864  const double pv = totVolume * intQuants.porosity().value();
865  const auto hydrocarbon = this->hydroCarbonFraction(fs);
866 
867  if (! this->hydrocarbonPoreVolume_.empty()) {
868  this->fipC_.assignPoreVolume(globalDofIdx,
869  totVolume * intQuants.referencePorosity());
870 
871  this->dynamicPoreVolume_[globalDofIdx] = pv;
872  this->hydrocarbonPoreVolume_[globalDofIdx] = pv * hydrocarbon;
873  }
874 
875  if (!this->pressureTimesHydrocarbonVolume_.empty() &&
876  !this->pressureTimesPoreVolume_.empty())
877  {
878  assert(this->hydrocarbonPoreVolume_.size() == this->pressureTimesHydrocarbonVolume_.size());
879  assert(this->fipC_.get(Inplace::Phase::PoreVolume).size() == this->pressureTimesPoreVolume_.size());
880 
881  if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
882  this->pressureTimesPoreVolume_[globalDofIdx] =
883  getValue(fs.pressure(oilPhaseIdx)) * pv;
884 
885  this->pressureTimesHydrocarbonVolume_[globalDofIdx] =
886  this->pressureTimesPoreVolume_[globalDofIdx] * hydrocarbon;
887  }
888  else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
889  this->pressureTimesPoreVolume_[globalDofIdx] =
890  getValue(fs.pressure(gasPhaseIdx)) * pv;
891 
892  this->pressureTimesHydrocarbonVolume_[globalDofIdx] =
893  this->pressureTimesPoreVolume_[globalDofIdx] * hydrocarbon;
894  }
895  else if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
896  this->pressureTimesPoreVolume_[globalDofIdx] =
897  getValue(fs.pressure(waterPhaseIdx)) * pv;
898  }
899  }
900  }
901 
902  void updatePhaseInplaceVolumes_(const unsigned globalDofIdx,
903  const IntensiveQuantities& intQuants,
904  const double totVolume)
905  {
906  std::array<Scalar, FluidSystem::numPhases> fip {};
907  std::array<Scalar, FluidSystem::numPhases> fipr{}; // at reservoir condition
908 
909  const auto& fs = intQuants.fluidState();
910  const auto pv = totVolume * intQuants.porosity().value();
911 
912  for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
913  if (!FluidSystem::phaseIsActive(phaseIdx)) {
914  continue;
915  }
916 
917  const auto b = getValue(fs.invB(phaseIdx));
918  const auto s = getValue(fs.saturation(phaseIdx));
919 
920  fipr[phaseIdx] = s * pv;
921  fip [phaseIdx] = b * fipr[phaseIdx];
922  }
923 
924  this->fipC_.assignVolumesSurface(globalDofIdx, fip);
925  this->fipC_.assignVolumesReservoir(globalDofIdx,
926  fs.saltConcentration().value(),
927  fipr);
928 
929  if (FluidSystem::phaseIsActive(oilPhaseIdx) &&
930  FluidSystem::phaseIsActive(gasPhaseIdx))
931  {
932  this->updateOilGasDistribution(globalDofIdx, fs, fip);
933  }
934 
935  if (FluidSystem::phaseIsActive(waterPhaseIdx) &&
936  FluidSystem::phaseIsActive(gasPhaseIdx))
937  {
938  this->updateGasWaterDistribution(globalDofIdx, fs, fip);
939  }
940 
941  if (FluidSystem::phaseIsActive(gasPhaseIdx) &&
942  this->fipC_.hasCo2InGas())
943  {
944  this->updateCO2InGas(globalDofIdx, pv, intQuants);
945  }
946 
947  if (this->fipC_.hasCo2InWater() &&
948  (FluidSystem::phaseIsActive(waterPhaseIdx) ||
949  FluidSystem::phaseIsActive(oilPhaseIdx)))
950  {
951  this->updateCO2InWater(globalDofIdx, pv, fs);
952  }
953 
954  if constexpr(enableBioeffects) {
955  const auto surfVolWat = pv * getValue(fs.saturation(waterPhaseIdx)) *
956  getValue(fs.invB(waterPhaseIdx));
957  if (this->fipC_.hasMicrobialMass()) {
958  this->updateMicrobialMass(globalDofIdx, intQuants, surfVolWat);
959  }
960  if (this->fipC_.hasBiofilmMass()) {
961  this->updateBiofilmMass(globalDofIdx, intQuants, totVolume);
962  }
963  if constexpr(enableMICP) {
964  if (this->fipC_.hasOxygenMass()) {
965  this->updateOxygenMass(globalDofIdx, intQuants, surfVolWat);
966  }
967  if (this->fipC_.hasUreaMass()) {
968  this->updateUreaMass(globalDofIdx, intQuants, surfVolWat);
969  }
970  if (this->fipC_.hasCalciteMass()) {
971  this->updateCalciteMass(globalDofIdx, intQuants, totVolume);
972  }
973  }
974  }
975 
976  if (this->fipC_.hasWaterMass() && FluidSystem::phaseIsActive(waterPhaseIdx))
977  {
978  this->updateWaterMass(globalDofIdx, fs, fip);
979  }
980  }
981 
982  template <typename FluidState, typename FIPArray>
983  void updateOilGasDistribution(const unsigned globalDofIdx,
984  const FluidState& fs,
985  const FIPArray& fip)
986  {
987  // Gas dissolved in oil and vaporized oil
988  const auto gasInPlaceLiquid = getValue(fs.Rs()) * fip[oilPhaseIdx];
989  const auto oilInPlaceGas = getValue(fs.Rv()) * fip[gasPhaseIdx];
990 
991  this->fipC_.assignOilGasDistribution(globalDofIdx, gasInPlaceLiquid, oilInPlaceGas);
992  }
993 
994  template <typename FluidState, typename FIPArray>
995  void updateGasWaterDistribution(const unsigned globalDofIdx,
996  const FluidState& fs,
997  const FIPArray& fip)
998  {
999  // Gas dissolved in water and vaporized water
1000  const auto gasInPlaceWater = getValue(fs.Rsw()) * fip[waterPhaseIdx];
1001  const auto waterInPlaceGas = getValue(fs.Rvw()) * fip[gasPhaseIdx];
1002 
1003  this->fipC_.assignGasWater(globalDofIdx, fip, gasInPlaceWater, waterInPlaceGas);
1004  }
1005 
1006  template <typename IntensiveQuantities>
1007  void updateCO2InGas(const unsigned globalDofIdx,
1008  const double pv,
1009  const IntensiveQuantities& intQuants)
1010  {
1011  const auto& scaledDrainageInfo = this->simulator_.problem().materialLawManager()
1012  ->oilWaterScaledEpsInfoDrainage(globalDofIdx);
1013 
1014  const auto& fs = intQuants.fluidState();
1015  Scalar sgcr = scaledDrainageInfo.Sgcr;
1016  if (this->simulator_.problem().materialLawManager()->enableHysteresis()) {
1017  const auto& matParams = simulator_.problem().materialLawParams(globalDofIdx);
1018  sgcr = MaterialLaw::trappedGasSaturation(matParams, /*maximumTrapping*/false);
1019  }
1020 
1021  Scalar trappedGasSaturation = scaledDrainageInfo.Sgcr;
1022  if (this->fipC_.has(Inplace::Phase::CO2MassInGasPhaseMaximumTrapped) ||
1023  this->fipC_.has(Inplace::Phase::CO2MassInGasPhaseMaximumUnTrapped))
1024  {
1025  if (this->simulator_.problem().materialLawManager()->enableHysteresis()) {
1026  const auto& matParams = simulator_.problem().materialLawParams(globalDofIdx);
1027  // Get the maximum trapped gas saturation
1028  trappedGasSaturation = MaterialLaw::trappedGasSaturation(matParams, /*maximumTrapping*/true);
1029  }
1030  }
1031 
1032  const Scalar sg = getValue(fs.saturation(gasPhaseIdx));
1033  Scalar strandedGasSaturation = scaledDrainageInfo.Sgcr;
1034  if (this->fipC_.has(Inplace::Phase::CO2MassInGasPhaseEffectiveTrapped) ||
1035  this->fipC_.has(Inplace::Phase::CO2MassInGasPhaseEffectiveUnTrapped))
1036  {
1037  if (this->simulator_.problem().materialLawManager()->enableHysteresis()) {
1038  const auto& matParams = simulator_.problem().materialLawParams(globalDofIdx);
1039  const double krg = getValue(intQuants.relativePermeability(gasPhaseIdx));
1040  strandedGasSaturation = MaterialLaw::strandedGasSaturation(matParams, sg, krg);
1041  }
1042  }
1043 
1044  const typename FIPContainer<FluidSystem>::Co2InGasInput v{
1045  pv,
1046  sg,
1047  sgcr,
1048  getValue(fs.density(gasPhaseIdx)),
1049  FluidSystem::phaseIsActive(waterPhaseIdx)
1050  ? FluidSystem::convertRvwToXgW(getValue(fs.Rvw()), fs.pvtRegionIndex())
1051  : FluidSystem::convertRvToXgO(getValue(fs.Rv()), fs.pvtRegionIndex()),
1052  FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex()),
1053  trappedGasSaturation,
1054  strandedGasSaturation,
1055  };
1056 
1057  this->fipC_.assignCo2InGas(globalDofIdx, v);
1058  }
1059 
1060  template <typename FluidState>
1061  void updateCO2InWater(const unsigned globalDofIdx,
1062  const double pv,
1063  const FluidState& fs)
1064  {
1065  const auto co2InWater = FluidSystem::phaseIsActive(oilPhaseIdx)
1066  ? this->co2InWaterFromOil(fs, pv)
1067  : this->co2InWaterFromWater(fs, pv);
1068 
1069  const Scalar mM = FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex());
1070 
1071  this->fipC_.assignCo2InWater(globalDofIdx, co2InWater, mM);
1072  }
1073 
1074  template <typename FluidState>
1075  Scalar co2InWaterFromWater(const FluidState& fs, const double pv) const
1076  {
1077  const double rhow = getValue(fs.density(waterPhaseIdx));
1078  const double sw = getValue(fs.saturation(waterPhaseIdx));
1079  const double xwG = FluidSystem::convertRswToXwG(getValue(fs.Rsw()), fs.pvtRegionIndex());
1080 
1081  const Scalar mM = FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex());
1082 
1083  return xwG * pv * rhow * sw / mM;
1084  }
1085 
1086  template <typename FluidState>
1087  Scalar co2InWaterFromOil(const FluidState& fs, const double pv) const
1088  {
1089  const double rhoo = getValue(fs.density(oilPhaseIdx));
1090  const double so = getValue(fs.saturation(oilPhaseIdx));
1091  const double xoG = FluidSystem::convertRsToXoG(getValue(fs.Rs()), fs.pvtRegionIndex());
1092 
1093  const Scalar mM = FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex());
1094 
1095  return xoG * pv * rhoo * so / mM;
1096  }
1097 
1098  template <typename FluidState, typename FIPArray>
1099  void updateWaterMass(const unsigned globalDofIdx,
1100  const FluidState& fs,
1101  const FIPArray& fip
1102  )
1103  {
1104  const Scalar rhoW = FluidSystem::referenceDensity(waterPhaseIdx, fs.pvtRegionIndex());
1105 
1106  this->fipC_.assignWaterMass(globalDofIdx, fip, rhoW);
1107  }
1108 
1109  template <typename IntensiveQuantities>
1110  void updateMicrobialMass(const unsigned globalDofIdx,
1111  const IntensiveQuantities& intQuants,
1112  const double surfVolWat)
1113  {
1114  const Scalar mass = surfVolWat * intQuants.microbialConcentration().value();
1115 
1116  this->fipC_.assignMicrobialMass(globalDofIdx, mass);
1117  }
1118 
1119  template <typename IntensiveQuantities>
1120  void updateOxygenMass(const unsigned globalDofIdx,
1121  const IntensiveQuantities& intQuants,
1122  const double surfVolWat)
1123  {
1124  const Scalar mass = surfVolWat * intQuants.oxygenConcentration().value();
1125 
1126  this->fipC_.assignOxygenMass(globalDofIdx, mass);
1127  }
1128 
1129  template <typename IntensiveQuantities>
1130  void updateUreaMass(const unsigned globalDofIdx,
1131  const IntensiveQuantities& intQuants,
1132  const double surfVolWat)
1133  {
1134  const Scalar mass = surfVolWat * intQuants.ureaConcentration().value();
1135 
1136  this->fipC_.assignUreaMass(globalDofIdx, mass);
1137  }
1138 
1139  template <typename IntensiveQuantities>
1140  void updateBiofilmMass(const unsigned globalDofIdx,
1141  const IntensiveQuantities& intQuants,
1142  const double totVolume)
1143  {
1144  const Scalar mass = totVolume * intQuants.biofilmMass().value();
1145 
1146  this->fipC_.assignBiofilmMass(globalDofIdx, mass);
1147  }
1148 
1149  template <typename IntensiveQuantities>
1150  void updateCalciteMass(const unsigned globalDofIdx,
1151  const IntensiveQuantities& intQuants,
1152  const double totVolume)
1153  {
1154  const Scalar mass = totVolume * intQuants.calciteMass().value();
1155 
1156  this->fipC_.assignCalciteMass(globalDofIdx, mass);
1157  }
1158 
1160  void setupElementExtractors_()
1161  {
1162  using Entry = typename Extractor::Entry;
1163  using Context = typename Extractor::Context;
1164  using ScalarEntry = typename Extractor::ScalarEntry;
1165  using PhaseEntry = typename Extractor::PhaseEntry;
1166 
1167  const bool hasResidual = simulator_.model().linearizer().residual().size() > 0;
1168  const auto& hysteresisConfig = simulator_.problem().materialLawManager()->hysteresisConfig();
1169 
1170  auto extractors = std::array{
1171  Entry{PhaseEntry{&this->saturation_,
1172  [](const unsigned phase, const Context& ectx)
1173  { return getValue(ectx.fs.saturation(phase)); }
1174  }
1175  },
1176  Entry{PhaseEntry{&this->invB_,
1177  [](const unsigned phase, const Context& ectx)
1178  { return getValue(ectx.fs.invB(phase)); }
1179  }
1180  },
1181  Entry{PhaseEntry{&this->density_,
1182  [](const unsigned phase, const Context& ectx)
1183  { return getValue(ectx.fs.density(phase)); }
1184  }
1185  },
1186  Entry{PhaseEntry{&this->relativePermeability_,
1187  [](const unsigned phase, const Context& ectx)
1188  { return getValue(ectx.intQuants.relativePermeability(phase)); }
1189  }
1190  },
1191  Entry{PhaseEntry{&this->viscosity_,
1192  [this](const unsigned phaseIdx, const Context& ectx)
1193  {
1194  if (this->extboC_.allocated() && phaseIdx == oilPhaseIdx) {
1195  return getValue(ectx.intQuants.oilViscosity());
1196  }
1197  else if (this->extboC_.allocated() && phaseIdx == gasPhaseIdx) {
1198  return getValue(ectx.intQuants.gasViscosity());
1199  }
1200  else {
1201  return getValue(ectx.fs.viscosity(phaseIdx));
1202  }
1203  }
1204  }
1205  },
1206  Entry{PhaseEntry{&this->residual_,
1207  [&modelResid = this->simulator_.model().linearizer().residual()]
1208  (const unsigned phaseIdx, const Context& ectx)
1209  {
1210  const unsigned sIdx = FluidSystem::solventComponentIndex(phaseIdx);
1211  const unsigned activeCompIdx = FluidSystem::canonicalToActiveCompIdx(sIdx);
1212  return modelResid[ectx.globalDofIdx][activeCompIdx];
1213  }
1214  },
1215  hasResidual
1216  },
1217  Entry{ScalarEntry{&this->rockCompPorvMultiplier_,
1218  [&problem = this->simulator_.problem()](const Context& ectx)
1219  {
1220  return problem.template
1221  rockCompPoroMultiplier<Scalar>(ectx.intQuants,
1222  ectx.globalDofIdx);
1223  }
1224  }
1225  },
1226  Entry{ScalarEntry{&this->rockCompTransMultiplier_,
1227  [&problem = this->simulator_.problem()](const Context& ectx)
1228  {
1229  return problem.
1230  template rockCompTransMultiplier<Scalar>(ectx.intQuants,
1231  ectx.globalDofIdx);
1232  }}
1233  },
1234  Entry{ScalarEntry{&this->minimumOilPressure_,
1235  [&problem = this->simulator_.problem()](const Context& ectx)
1236  {
1237  return std::min(getValue(ectx.fs.pressure(oilPhaseIdx)),
1238  problem.minOilPressure(ectx.globalDofIdx));
1239  }
1240  }
1241  },
1242  Entry{ScalarEntry{&this->bubblePointPressure_,
1243  [&failedCells = this->failedCellsPb_,
1244  &vanguard = this->simulator_.vanguard()](const Context& ectx)
1245  {
1246  try {
1247  return getValue(
1248  FluidSystem::bubblePointPressure(ectx.fs,
1249  ectx.intQuants.pvtRegionIndex())
1250  );
1251  } catch (const NumericalProblem&) {
1252  const auto cartesianIdx = vanguard.cartesianIndex(ectx.globalDofIdx);
1253  failedCells.push_back(cartesianIdx);
1254  return Scalar{0};
1255  }
1256  }
1257  }
1258  },
1259  Entry{ScalarEntry{&this->dewPointPressure_,
1260  [&failedCells = this->failedCellsPd_,
1261  &vanguard = this->simulator_.vanguard()](const Context& ectx)
1262  {
1263  try {
1264  return getValue(
1265  FluidSystem::dewPointPressure(ectx.fs,
1266  ectx.intQuants.pvtRegionIndex())
1267  );
1268  } catch (const NumericalProblem&) {
1269  const auto cartesianIdx = vanguard.cartesianIndex(ectx.globalDofIdx);
1270  failedCells.push_back(cartesianIdx);
1271  return Scalar{0};
1272  }
1273  }
1274  }
1275  },
1276  Entry{ScalarEntry{&this->overburdenPressure_,
1277  [&problem = simulator_.problem()](const Context& ectx)
1278  { return problem.overburdenPressure(ectx.globalDofIdx); }
1279  }
1280  },
1281  Entry{ScalarEntry{&this->temperature_,
1282  [](const Context& ectx)
1283  { return getValue(ectx.fs.temperature(oilPhaseIdx)); }
1284  }
1285  },
1286  Entry{ScalarEntry{&this->sSol_,
1287  [](const Context& ectx)
1288  { return getValue(ectx.intQuants.solventSaturation()); }
1289  }
1290  },
1291  Entry{ScalarEntry{&this->rswSol_,
1292  [](const Context& ectx)
1293  { return getValue(ectx.intQuants.rsSolw()); }
1294  }
1295  },
1296  Entry{ScalarEntry{&this->cPolymer_,
1297  [](const Context& ectx)
1298  { return getValue(ectx.intQuants.polymerConcentration()); }
1299  }
1300  },
1301  Entry{ScalarEntry{&this->cFoam_,
1302  [](const Context& ectx)
1303  { return getValue(ectx.intQuants.foamConcentration()); }
1304  }
1305  },
1306  Entry{ScalarEntry{&this->cSalt_,
1307  [](const Context& ectx)
1308  { return getValue(ectx.fs.saltConcentration()); }
1309  }
1310  },
1311  Entry{ScalarEntry{&this->pSalt_,
1312  [](const Context& ectx)
1313  { return getValue(ectx.fs.saltSaturation()); }
1314  }
1315  },
1316  Entry{ScalarEntry{&this->permFact_,
1317  [](const Context& ectx)
1318  { return getValue(ectx.intQuants.permFactor()); }
1319  }
1320  },
1321  Entry{ScalarEntry{&this->rPorV_,
1322  [&model = this->simulator_.model()](const Context& ectx)
1323  {
1324  const auto totVolume = model.dofTotalVolume(ectx.globalDofIdx);
1325  return totVolume * getValue(ectx.intQuants.porosity());
1326  }
1327  }
1328  },
1329  Entry{ScalarEntry{&this->rs_,
1330  [](const Context& ectx)
1331  { return getValue(ectx.fs.Rs()); }
1332  }
1333  },
1334  Entry{ScalarEntry{&this->rv_,
1335  [](const Context& ectx)
1336  { return getValue(ectx.fs.Rv()); }
1337  }
1338  },
1339  Entry{ScalarEntry{&this->rsw_,
1340  [](const Context& ectx)
1341  { return getValue(ectx.fs.Rsw()); }
1342  }
1343  },
1344  Entry{ScalarEntry{&this->rvw_,
1345  [](const Context& ectx)
1346  { return getValue(ectx.fs.Rvw()); }
1347  }
1348  },
1349  Entry{ScalarEntry{&this->ppcw_,
1350  [&matLawManager = *this->simulator_.problem().materialLawManager()]
1351  (const Context& ectx)
1352  {
1353  return matLawManager.
1354  oilWaterScaledEpsInfoDrainage(ectx.globalDofIdx).maxPcow;
1355  }
1356  }
1357  },
1358  Entry{ScalarEntry{&this->drsdtcon_,
1359  [&problem = this->simulator_.problem()](const Context& ectx)
1360  {
1361  return problem.drsdtcon(ectx.globalDofIdx,
1362  ectx.episodeIndex);
1363  }
1364  }
1365  },
1366  Entry{ScalarEntry{&this->pcgw_,
1367  [](const Context& ectx)
1368  {
1369  return getValue(ectx.fs.pressure(gasPhaseIdx)) -
1370  getValue(ectx.fs.pressure(waterPhaseIdx));
1371  }
1372  }
1373  },
1374  Entry{ScalarEntry{&this->pcow_,
1375  [](const Context& ectx)
1376  {
1377  return getValue(ectx.fs.pressure(oilPhaseIdx)) -
1378  getValue(ectx.fs.pressure(waterPhaseIdx));
1379  }
1380  }
1381  },
1382  Entry{ScalarEntry{&this->pcog_,
1383  [](const Context& ectx)
1384  {
1385  return getValue(ectx.fs.pressure(gasPhaseIdx)) -
1386  getValue(ectx.fs.pressure(oilPhaseIdx));
1387  }
1388  }
1389  },
1390  Entry{ScalarEntry{&this->fluidPressure_,
1391  [](const Context& ectx)
1392  {
1393  if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1394  // Output oil pressure as default
1395  return getValue(ectx.fs.pressure(oilPhaseIdx));
1396  }
1397  else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1398  // Output gas if oil is not present
1399  return getValue(ectx.fs.pressure(gasPhaseIdx));
1400  }
1401  else {
1402  // Output water if neither oil nor gas is present
1403  return getValue(ectx.fs.pressure(waterPhaseIdx));
1404  }
1405  }
1406  }
1407  },
1408  Entry{ScalarEntry{&this->gasDissolutionFactor_,
1409  [&problem = this->simulator_.problem()](const Context& ectx)
1410  {
1411  const Scalar SoMax = problem.maxOilSaturation(ectx.globalDofIdx);
1412  return FluidSystem::template
1413  saturatedDissolutionFactor<FluidState, Scalar>(ectx.fs,
1414  oilPhaseIdx,
1415  ectx.pvtRegionIdx,
1416  SoMax);
1417  }
1418  }
1419  },
1420  Entry{ScalarEntry{&this->oilVaporizationFactor_,
1421  [&problem = this->simulator_.problem()](const Context& ectx)
1422  {
1423  const Scalar SoMax = problem.maxOilSaturation(ectx.globalDofIdx);
1424  return FluidSystem::template
1425  saturatedDissolutionFactor<FluidState, Scalar>(ectx.fs,
1426  gasPhaseIdx,
1427  ectx.pvtRegionIdx,
1428  SoMax);
1429  }
1430  }
1431  },
1432  Entry{ScalarEntry{&this->gasDissolutionFactorInWater_,
1433  [&problem = this->simulator_.problem()](const Context& ectx)
1434  {
1435  const Scalar SwMax = problem.maxWaterSaturation(ectx.globalDofIdx);
1436  return FluidSystem::template
1437  saturatedDissolutionFactor<FluidState, Scalar>(ectx.fs,
1438  waterPhaseIdx,
1439  ectx.pvtRegionIdx,
1440  SwMax);
1441  }
1442  }
1443  },
1444  Entry{ScalarEntry{&this->waterVaporizationFactor_,
1445  [](const Context& ectx)
1446  {
1447  return FluidSystem::template
1448  saturatedVaporizationFactor<FluidState, Scalar>(ectx.fs,
1449  gasPhaseIdx,
1450  ectx.pvtRegionIdx);
1451  }
1452  }
1453  },
1454  Entry{ScalarEntry{&this->gasFormationVolumeFactor_,
1455  [](const Context& ectx)
1456  {
1457  return 1.0 / FluidSystem::template
1458  inverseFormationVolumeFactor<FluidState, Scalar>(ectx.fs,
1459  gasPhaseIdx,
1460  ectx.pvtRegionIdx);
1461  }
1462  }
1463  },
1464  Entry{ScalarEntry{&this->saturatedOilFormationVolumeFactor_,
1465  [](const Context& ectx)
1466  {
1467  return 1.0 / FluidSystem::template
1468  saturatedInverseFormationVolumeFactor<FluidState, Scalar>(ectx.fs,
1469  oilPhaseIdx,
1470  ectx.pvtRegionIdx);
1471  }
1472  }
1473  },
1474  Entry{ScalarEntry{&this->oilSaturationPressure_,
1475  [](const Context& ectx)
1476  {
1477  return FluidSystem::template
1478  saturationPressure<FluidState, Scalar>(ectx.fs,
1479  oilPhaseIdx,
1480  ectx.pvtRegionIdx);
1481  }
1482  }
1483  },
1484  Entry{ScalarEntry{&this->soMax_,
1485  [&problem = this->simulator_.problem()](const Context& ectx)
1486  {
1487  return std::max(getValue(ectx.fs.saturation(oilPhaseIdx)),
1488  problem.maxOilSaturation(ectx.globalDofIdx));
1489  }
1490  },
1491  !hysteresisConfig.enableHysteresis()
1492  },
1493  Entry{ScalarEntry{&this->swMax_,
1494  [&problem = this->simulator_.problem()](const Context& ectx)
1495  {
1496  return std::max(getValue(ectx.fs.saturation(waterPhaseIdx)),
1497  problem.maxWaterSaturation(ectx.globalDofIdx));
1498  }
1499  },
1500  !hysteresisConfig.enableHysteresis()
1501  },
1502  Entry{ScalarEntry{&this->soMax_,
1503  [](const Context& ectx)
1504  { return ectx.hParams.somax; }
1505  },
1506  hysteresisConfig.enableHysteresis() &&
1507  hysteresisConfig.enableNonWettingHysteresis() &&
1508  FluidSystem::phaseIsActive(oilPhaseIdx) &&
1509  FluidSystem::phaseIsActive(waterPhaseIdx)
1510  },
1511  Entry{ScalarEntry{&this->swMax_,
1512  [](const Context& ectx)
1513  { return ectx.hParams.swmax; }
1514  },
1515  hysteresisConfig.enableHysteresis() &&
1516  hysteresisConfig.enableWettingHysteresis() &&
1517  FluidSystem::phaseIsActive(oilPhaseIdx) &&
1518  FluidSystem::phaseIsActive(waterPhaseIdx)
1519  },
1520  Entry{ScalarEntry{&this->swmin_,
1521  [](const Context& ectx)
1522  { return ectx.hParams.swmin; }
1523  },
1524  hysteresisConfig.enableHysteresis() &&
1525  hysteresisConfig.enablePCHysteresis() &&
1526  FluidSystem::phaseIsActive(oilPhaseIdx) &&
1527  FluidSystem::phaseIsActive(waterPhaseIdx)
1528  },
1529  Entry{ScalarEntry{&this->sgmax_,
1530  [](const Context& ectx)
1531  { return ectx.hParams.sgmax; }
1532  },
1533  hysteresisConfig.enableHysteresis() &&
1534  hysteresisConfig.enableNonWettingHysteresis() &&
1535  FluidSystem::phaseIsActive(oilPhaseIdx) &&
1536  FluidSystem::phaseIsActive(gasPhaseIdx)
1537  },
1538  Entry{ScalarEntry{&this->shmax_,
1539  [](const Context& ectx)
1540  { return ectx.hParams.shmax; }
1541  },
1542  hysteresisConfig.enableHysteresis() &&
1543  hysteresisConfig.enableWettingHysteresis() &&
1544  FluidSystem::phaseIsActive(oilPhaseIdx) &&
1545  FluidSystem::phaseIsActive(gasPhaseIdx)
1546  },
1547  Entry{ScalarEntry{&this->somin_,
1548  [](const Context& ectx)
1549  { return ectx.hParams.somin; }
1550  },
1551  hysteresisConfig.enableHysteresis() &&
1552  hysteresisConfig.enablePCHysteresis() &&
1553  FluidSystem::phaseIsActive(oilPhaseIdx) &&
1554  FluidSystem::phaseIsActive(gasPhaseIdx)
1555  },
1556  Entry{[&model = this->simulator_.model(), this](const Context& ectx)
1557  {
1558  // Note: We intentionally exclude effects of rock
1559  // compressibility by using referencePorosity() here.
1560  const auto porv = ectx.intQuants.referencePorosity()
1561  * model.dofTotalVolume(ectx.globalDofIdx);
1562 
1563  this->aggregateAverageDensityContributions_(ectx.fs, ectx.globalDofIdx,
1564  static_cast<double>(porv));
1565  }, this->regionAvgDensity_.has_value()
1566  },
1567  Entry{[&extboC = this->extboC_](const Context& ectx)
1568  {
1569  extboC.assignVolumes(ectx.globalDofIdx,
1570  ectx.intQuants.xVolume().value(),
1571  ectx.intQuants.yVolume().value());
1572  extboC.assignZFraction(ectx.globalDofIdx,
1573  ectx.intQuants.zFraction().value());
1574 
1575  const Scalar stdVolOil = getValue(ectx.fs.saturation(oilPhaseIdx)) *
1576  getValue(ectx.fs.invB(oilPhaseIdx)) +
1577  getValue(ectx.fs.saturation(gasPhaseIdx)) *
1578  getValue(ectx.fs.invB(gasPhaseIdx)) *
1579  getValue(ectx.fs.Rv());
1580  const Scalar stdVolGas = getValue(ectx.fs.saturation(gasPhaseIdx)) *
1581  getValue(ectx.fs.invB(gasPhaseIdx)) *
1582  (1.0 - ectx.intQuants.yVolume().value()) +
1583  getValue(ectx.fs.saturation(oilPhaseIdx)) *
1584  getValue(ectx.fs.invB(oilPhaseIdx)) *
1585  getValue(ectx.fs.Rs()) *
1586  (1.0 - ectx.intQuants.xVolume().value());
1587  const Scalar stdVolCo2 = getValue(ectx.fs.saturation(gasPhaseIdx)) *
1588  getValue(ectx.fs.invB(gasPhaseIdx)) *
1589  ectx.intQuants.yVolume().value() +
1590  getValue(ectx.fs.saturation(oilPhaseIdx)) *
1591  getValue(ectx.fs.invB(oilPhaseIdx)) *
1592  getValue(ectx.fs.Rs()) *
1593  ectx.intQuants.xVolume().value();
1594  const Scalar rhoO = FluidSystem::referenceDensity(gasPhaseIdx, ectx.pvtRegionIdx);
1595  const Scalar rhoG = FluidSystem::referenceDensity(gasPhaseIdx, ectx.pvtRegionIdx);
1596  const Scalar rhoCO2 = ectx.intQuants.zRefDensity();
1597  const Scalar stdMassTotal = 1.0e-10 + stdVolOil * rhoO + stdVolGas * rhoG + stdVolCo2 * rhoCO2;
1598  extboC.assignMassFractions(ectx.globalDofIdx,
1599  stdVolGas * rhoG / stdMassTotal,
1600  stdVolOil * rhoO / stdMassTotal,
1601  stdVolCo2 * rhoCO2 / stdMassTotal);
1602  }, this->extboC_.allocated()
1603  },
1604  Entry{[&bioeffectsC = this->bioeffectsC_](const Context& ectx)
1605  {
1606  bioeffectsC.assign(ectx.globalDofIdx,
1607  ectx.intQuants.microbialConcentration().value(),
1608  ectx.intQuants.biofilmVolumeFraction().value());
1609  if (Indices::enableMICP) {
1610  bioeffectsC.assign(ectx.globalDofIdx,
1611  ectx.intQuants.oxygenConcentration().value(),
1612  ectx.intQuants.ureaConcentration().value(),
1613  ectx.intQuants.calciteVolumeFraction().value());
1614  }
1615  }, this->bioeffectsC_.allocated()
1616  },
1617  Entry{[&runspec = this->eclState_.runspec(),
1618  &CO2H2C = this->CO2H2C_](const Context& ectx)
1619  {
1620  const auto xwg = FluidSystem::convertRswToXwG(getValue(ectx.fs.Rsw()), ectx.pvtRegionIdx);
1621  const auto xgw = FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.pvtRegionIdx);
1622  CO2H2C.assign(ectx.globalDofIdx,
1623  FluidSystem::convertXwGToxwG(xwg, ectx.pvtRegionIdx),
1624  FluidSystem::convertXgWToxgW(xgw, ectx.pvtRegionIdx),
1625  runspec.co2Storage());
1626  }, this->CO2H2C_.allocated()
1627  },
1628  Entry{[&rftC = this->rftC_,
1629  &vanguard = this->simulator_.vanguard()](const Context& ectx)
1630  {
1631  const auto cartesianIdx = vanguard.cartesianIndex(ectx.globalDofIdx);
1632  rftC.assign(cartesianIdx,
1633  [&fs = ectx.fs]() { return getValue(fs.pressure(oilPhaseIdx)); },
1634  [&fs = ectx.fs]() { return getValue(fs.saturation(waterPhaseIdx)); },
1635  [&fs = ectx.fs]() { return getValue(fs.saturation(gasPhaseIdx)); });
1636  }
1637  },
1638  Entry{[&tC = this->tracerC_,
1639  &tM = this->simulator_.problem().tracerModel()](const Context& ectx)
1640  {
1641  tC.assignFreeConcentrations(ectx.globalDofIdx,
1642  [gIdx = ectx.globalDofIdx, &tM](const unsigned tracerIdx)
1643  { return tM.freeTracerConcentration(tracerIdx, gIdx); });
1644  tC.assignSolConcentrations(ectx.globalDofIdx,
1645  [gIdx = ectx.globalDofIdx, &tM](const unsigned tracerIdx)
1646  { return tM.solTracerConcentration(tracerIdx, gIdx); });
1647  }
1648  },
1649  Entry{[&flowsInf = this->simulator_.problem().model().linearizer().getFlowsInfo(),
1650  &flowsC = this->flowsC_,
1651  &vanguard = this->simulator_.vanguard()](const Context& ectx)
1652  {
1653  const auto gas_idx = Indices::gasEnabled ?
1654  conti0EqIdx + FluidSystem::canonicalToActiveCompIdx(gasCompIdx) : -1;
1655  const auto oil_idx = Indices::oilEnabled ?
1656  conti0EqIdx + FluidSystem::canonicalToActiveCompIdx(oilCompIdx) : -1;
1657  const auto water_idx = Indices::waterEnabled ?
1658  conti0EqIdx + FluidSystem::canonicalToActiveCompIdx(waterCompIdx) : -1;
1659  const auto& flowsInfos = flowsInf[ectx.globalDofIdx];
1660  if (!flowsC.blockFlows().empty()) {
1661  const std::vector<int>& blockIdxs = flowsC.blockFlows();
1662  const unsigned cartesianIdx = vanguard.cartesianIndex(ectx.globalDofIdx);
1663  if (std::ranges::binary_search(blockIdxs, cartesianIdx)) {
1664  const auto compIdxs = std::array{ gasCompIdx, oilCompIdx, waterCompIdx };
1665  const auto compEnabled = std::array{ Indices::gasEnabled, Indices::oilEnabled, Indices::waterEnabled };
1666  for (const auto& flowsInfo : flowsInfos) {
1667  if (flowsInfo.faceId < 0) {
1668  continue;
1669  }
1670  for (unsigned ii = 0; ii < compIdxs.size(); ++ii) {
1671  if (!compEnabled[ii]) {
1672  continue;
1673  }
1674  if (flowsC.hasBlockFlowValue(cartesianIdx, flowsInfo.faceId, compIdxs[ii])) {
1675  flowsC.assignBlockFlows(flowsC.blockFlowsIds(cartesianIdx, flowsInfo.faceId, compIdxs[ii]),
1676  flowsInfo.faceId,
1677  compIdxs[ii],
1678  flowsInfo.flow[conti0EqIdx
1679  + FluidSystem::canonicalToActiveCompIdx(compIdxs[ii])]);
1680  }
1681  }
1682  }
1683  }
1684  }
1685  else {
1686  for (const auto& flowsInfo : flowsInfos) {
1687  flowsC.assignFlows(ectx.globalDofIdx,
1688  flowsInfo.faceId,
1689  flowsInfo.nncId,
1690  value_or_zero(gas_idx, flowsInfo.flow),
1691  value_or_zero(oil_idx, flowsInfo.flow),
1692  value_or_zero(water_idx, flowsInfo.flow));
1693  }
1694  }
1695  }, !this->simulator_.problem().model().linearizer().getFlowsInfo().empty()
1696  },
1697  Entry{[&floresInf = this->simulator_.problem().model().linearizer().getFloresInfo(),
1698  &flowsC = this->flowsC_](const Context& ectx)
1699  {
1700  const auto gas_idx = Indices::gasEnabled ?
1701  conti0EqIdx + FluidSystem::canonicalToActiveCompIdx(gasCompIdx) : -1;
1702  const auto oil_idx = Indices::oilEnabled ?
1703  conti0EqIdx + FluidSystem::canonicalToActiveCompIdx(oilCompIdx) : -1;
1704  const auto water_idx = Indices::waterEnabled ?
1705  conti0EqIdx + FluidSystem::canonicalToActiveCompIdx(waterCompIdx) : -1;
1706  const auto& floresInfos = floresInf[ectx.globalDofIdx];
1707  for (const auto& floresInfo : floresInfos) {
1708  flowsC.assignFlores(ectx.globalDofIdx,
1709  floresInfo.faceId,
1710  floresInfo.nncId,
1711  value_or_zero(gas_idx, floresInfo.flow),
1712  value_or_zero(oil_idx, floresInfo.flow),
1713  value_or_zero(water_idx, floresInfo.flow));
1714  }
1715  }, !this->simulator_.problem().model().linearizer().getFloresInfo().empty()
1716  },
1717  Entry{[&velocityInf = this->simulator_.problem().model().linearizer().getVelocityInfo(),
1718  &flowsC = this->flowsC_,
1719  &vanguard = this->simulator_.vanguard()](const Context& ectx)
1720  {
1721  const auto& velocityInfos = velocityInf[ectx.globalDofIdx];
1722  const std::vector<int>& blockIdxs = flowsC.blockVelocity();
1723  const unsigned cartesianIdx = vanguard.cartesianIndex(ectx.globalDofIdx);
1724  if (std::ranges::binary_search(blockIdxs, cartesianIdx)) {
1725  const auto compIdxs = std::array{ gasCompIdx, oilCompIdx, waterCompIdx };
1726  const auto compEnabled = std::array{ Indices::gasEnabled, Indices::oilEnabled, Indices::waterEnabled };
1727  for (const auto& velocityInfo : velocityInfos) {
1728  if (velocityInfo.faceId < 0) {
1729  continue;
1730  }
1731  for (unsigned ii = 0; ii < compIdxs.size(); ++ii) {
1732  if (!compEnabled[ii]) {
1733  continue;
1734  }
1735  if (flowsC.hasBlockVelocityValue(cartesianIdx, velocityInfo.faceId, compIdxs[ii])) {
1736  flowsC.assignBlockVelocity(flowsC.blockVelocityIds(cartesianIdx, velocityInfo.faceId, compIdxs[ii]),
1737  velocityInfo.faceId,
1738  compIdxs[ii],
1739  velocityInfo.velocity[conti0EqIdx
1740  + FluidSystem::canonicalToActiveCompIdx(compIdxs[ii])]);
1741  }
1742  }
1743  }
1744  }
1745  }, !this->flowsC_.blockVelocity().empty() &&
1746  !this->simulator_.problem().model().linearizer().getVelocityInfo().empty()
1747  },
1748  // hack to make the intial output of rs and rv Ecl compatible.
1749  // For cells with swat == 1 Ecl outputs; rs = rsSat and rv=rvSat, in all but the initial step
1750  // where it outputs rs and rv values calculated by the initialization. To be compatible we overwrite
1751  // rs and rv with the values computed in the initially.
1752  // Volume factors, densities and viscosities need to be recalculated with the updated rs and rv values.
1753  Entry{ScalarEntry{&this->rv_,
1754  [&problem = this->simulator_.problem()](const Context& ectx)
1755  { return problem.initialFluidState(ectx.globalDofIdx).Rv(); }
1756  },
1757  simulator_.episodeIndex() < 0 &&
1758  FluidSystem::phaseIsActive(oilPhaseIdx) &&
1759  FluidSystem::phaseIsActive(gasPhaseIdx)
1760  },
1761  Entry{ScalarEntry{&this->rs_,
1762  [&problem = this->simulator_.problem()](const Context& ectx)
1763  { return problem.initialFluidState(ectx.globalDofIdx).Rs(); }
1764  },
1765  simulator_.episodeIndex() < 0 &&
1766  FluidSystem::phaseIsActive(oilPhaseIdx) &&
1767  FluidSystem::phaseIsActive(gasPhaseIdx)
1768  },
1769  Entry{ScalarEntry{&this->rsw_,
1770  [&problem = this->simulator_.problem()](const Context& ectx)
1771  { return problem.initialFluidState(ectx.globalDofIdx).Rsw(); }
1772  },
1773  simulator_.episodeIndex() < 0 &&
1774  FluidSystem::phaseIsActive(oilPhaseIdx) &&
1775  FluidSystem::phaseIsActive(gasPhaseIdx)
1776  },
1777  Entry{ScalarEntry{&this->rvw_,
1778  [&problem = this->simulator_.problem()](const Context& ectx)
1779  { return problem.initialFluidState(ectx.globalDofIdx).Rvw(); }
1780  },
1781  simulator_.episodeIndex() < 0 &&
1782  FluidSystem::phaseIsActive(oilPhaseIdx) &&
1783  FluidSystem::phaseIsActive(gasPhaseIdx)
1784  },
1785  // re-compute the volume factors, viscosities and densities if asked for
1786  Entry{PhaseEntry{&this->density_,
1787  [&problem = this->simulator_.problem()](const unsigned phase,
1788  const Context& ectx)
1789  {
1790  const auto& fsInitial = problem.initialFluidState(ectx.globalDofIdx);
1791  return FluidSystem::density(fsInitial,
1792  phase,
1793  ectx.intQuants.pvtRegionIndex());
1794  }
1795  },
1796  simulator_.episodeIndex() < 0 &&
1797  FluidSystem::phaseIsActive(oilPhaseIdx) &&
1798  FluidSystem::phaseIsActive(gasPhaseIdx)
1799  },
1800  Entry{PhaseEntry{&this->invB_,
1801  [&problem = this->simulator_.problem()](const unsigned phase,
1802  const Context& ectx)
1803  {
1804  const auto& fsInitial = problem.initialFluidState(ectx.globalDofIdx);
1805  return FluidSystem::inverseFormationVolumeFactor(fsInitial,
1806  phase,
1807  ectx.intQuants.pvtRegionIndex());
1808  }
1809  },
1810  simulator_.episodeIndex() < 0 &&
1811  FluidSystem::phaseIsActive(oilPhaseIdx) &&
1812  FluidSystem::phaseIsActive(gasPhaseIdx)
1813  },
1814  Entry{PhaseEntry{&this->viscosity_,
1815  [&problem = this->simulator_.problem()](const unsigned phase,
1816  const Context& ectx)
1817  {
1818  const auto& fsInitial = problem.initialFluidState(ectx.globalDofIdx);
1819  return FluidSystem::viscosity(fsInitial,
1820  phase,
1821  ectx.intQuants.pvtRegionIndex());
1822  }
1823  },
1824  simulator_.episodeIndex() < 0 &&
1825  FluidSystem::phaseIsActive(oilPhaseIdx) &&
1826  FluidSystem::phaseIsActive(gasPhaseIdx)
1827  },
1828  };
1829 
1830  // Setup active extractors
1831  this->extractors_ = Extractor::removeInactive(extractors);
1832 
1833  // Geochemistry
1834  if constexpr (getPropValue<TypeTag, Properties::EnableGeochemistry>()) {
1835  if (this->geochemC_.allocated()) {
1836  this->extractors_.push_back(
1837  Entry{
1838  [&gC = this->geochemC_,
1839  &gM = this->simulator_.problem().geochemistryModel()](const Context& ectx)
1840  {
1841  gC.assignSpeciesConcentrations(
1842  ectx.globalDofIdx,
1843  [gIdx = ectx.globalDofIdx, &gM](const unsigned speciesIdx)
1844  { return gM.speciesConcentration(speciesIdx, gIdx); }
1845  );
1846  gC.assignMineralConcentrations(
1847  ectx.globalDofIdx,
1848  [gIdx = ectx.globalDofIdx, &gM](const unsigned mineralIdx)
1849  { return gM.mineralConcentration(mineralIdx, gIdx); }
1850  );
1851  gC.assignPH(ectx.globalDofIdx, gM.PH(ectx.globalDofIdx));
1852  }
1853  }
1854  );
1855  }
1856  }
1857 
1858  // Geomechanics
1859  if constexpr (getPropValue<TypeTag, Properties::EnableMech>()) {
1860  if (this->mech_.allocated()) {
1861  this->extractors_.push_back(
1862  Entry{[&mech = this->mech_,
1863  &model = simulator_.problem().geoMechModel()](const Context& ectx)
1864  {
1865  mech.assignDelStress(ectx.globalDofIdx,
1866  model.delstress(ectx.globalDofIdx));
1867 
1868  mech.assignDisplacement(ectx.globalDofIdx,
1869  model.disp(ectx.globalDofIdx, /*include_fracture*/true));
1870 
1871  // is the tresagii stress which make rock fracture
1872  mech.assignFracStress(ectx.globalDofIdx,
1873  model.fractureStress(ectx.globalDofIdx));
1874 
1875  mech.assignLinStress(ectx.globalDofIdx,
1876  model.linstress(ectx.globalDofIdx));
1877 
1878  mech.assignPotentialForces(ectx.globalDofIdx,
1879  model.mechPotentialForce(ectx.globalDofIdx),
1880  model.mechPotentialPressForce(ectx.globalDofIdx),
1881  model.mechPotentialTempForce(ectx.globalDofIdx));
1882 
1883  mech.assignStrain(ectx.globalDofIdx,
1884  model.strain(ectx.globalDofIdx, /*include_fracture*/true));
1885 
1886  // Total stress is not stored but calculated result is Voigt notation
1887  mech.assignStress(ectx.globalDofIdx,
1888  model.stress(ectx.globalDofIdx, /*include_fracture*/true));
1889  }
1890  }
1891  );
1892  }
1893  }
1894  }
1895 
1897  void setupBlockExtractors_(const bool isSubStep,
1898  const int reportStepNum)
1899  {
1900  using Entry = typename BlockExtractor::Entry;
1901  using Context = typename BlockExtractor::Context;
1902  using PhaseEntry = typename BlockExtractor::PhaseEntry;
1903  using ScalarEntry = typename BlockExtractor::ScalarEntry;
1904 
1905  using namespace std::string_view_literals;
1906 
1907  const auto pressure_handler =
1908  Entry{ScalarEntry{std::vector{"BPR"sv, "BPRESSUR"sv},
1909  [](const Context& ectx)
1910  {
1911  if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1912  return getValue(ectx.fs.pressure(oilPhaseIdx));
1913  }
1914  else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1915  return getValue(ectx.fs.pressure(gasPhaseIdx));
1916  }
1917  else { //if (FluidSystem::phaseIsActive(waterPhaseIdx))
1918  return getValue(ectx.fs.pressure(waterPhaseIdx));
1919  }
1920  }
1921  }
1922  };
1923 
1924  const auto handlers = std::array{
1925  pressure_handler,
1926  Entry{PhaseEntry{std::array{
1927  std::array{"BWSAT"sv, "BOSAT"sv, "BGSAT"sv},
1928  std::array{"BSWAT"sv, "BSOIL"sv, "BSGAS"sv}
1929  },
1930  [](const unsigned phaseIdx, const Context& ectx)
1931  {
1932  return getValue(ectx.fs.saturation(phaseIdx));
1933  }
1934  }
1935  },
1936  Entry{ScalarEntry{"BNSAT",
1937  [](const Context& ectx)
1938  {
1939  return ectx.intQuants.solventSaturation().value();
1940  }
1941  }
1942  },
1943  Entry{ScalarEntry{std::vector{"BTCNFHEA"sv, "BTEMP"sv},
1944  [](const Context& ectx)
1945  {
1946  if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1947  return getValue(ectx.fs.temperature(oilPhaseIdx));
1948  }
1949  else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1950  return getValue(ectx.fs.temperature(gasPhaseIdx));
1951  }
1952  else { //if (FluidSystem::phaseIsActive(waterPhaseIdx))
1953  return getValue(ectx.fs.temperature(waterPhaseIdx));
1954  }
1955  }
1956  }
1957  },
1958  Entry{PhaseEntry{std::array{
1959  std::array{"BWKR"sv, "BOKR"sv, "BGKR"sv},
1960  std::array{"BKRW"sv, "BKRO"sv, "BKRG"sv}
1961  },
1962  [](const unsigned phaseIdx, const Context& ectx)
1963  {
1964  return getValue(ectx.intQuants.relativePermeability(phaseIdx));
1965  }
1966  }
1967  },
1968  Entry{ScalarEntry{"BKROG",
1969  [&problem = this->simulator_.problem()](const Context& ectx)
1970  {
1971  const auto& materialParams =
1972  problem.materialLawParams(ectx.elemCtx,
1973  ectx.dofIdx,
1974  /* timeIdx = */ 0);
1975  return getValue(MaterialLaw::template
1976  relpermOilInOilGasSystem<Evaluation>(materialParams,
1977  ectx.fs));
1978  }
1979  }
1980  },
1981  Entry{ScalarEntry{"BKROW",
1982  [&problem = this->simulator_.problem()](const Context& ectx)
1983  {
1984  const auto& materialParams = problem.materialLawParams(ectx.elemCtx,
1985  ectx.dofIdx,
1986  /* timeIdx = */ 0);
1987  return getValue(MaterialLaw::template
1988  relpermOilInOilWaterSystem<Evaluation>(materialParams,
1989  ectx.fs));
1990  }
1991  }
1992  },
1993  Entry{ScalarEntry{"BWPC",
1994  [](const Context& ectx)
1995  {
1996  if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1997  return getValue(ectx.fs.pressure(oilPhaseIdx)) -
1998  getValue(ectx.fs.pressure(waterPhaseIdx));
1999  }
2000  else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
2001  return getValue(ectx.fs.pressure(gasPhaseIdx)) -
2002  getValue(ectx.fs.pressure(waterPhaseIdx));
2003  }
2004  else {
2005  return Scalar(0.0);
2006  }
2007  }
2008  }
2009  },
2010  Entry{ScalarEntry{"BGPC",
2011  [](const Context& ectx)
2012  {
2013  if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
2014  return getValue(ectx.fs.pressure(gasPhaseIdx)) -
2015  getValue(ectx.fs.pressure(oilPhaseIdx));
2016  }
2017  else if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
2018  return getValue(ectx.fs.pressure(gasPhaseIdx)) -
2019  getValue(ectx.fs.pressure(waterPhaseIdx));
2020  }
2021  else {
2022  return Scalar(0.0);
2023  }
2024  }
2025  }
2026  },
2027  Entry{ScalarEntry{"BWPR",
2028  [](const Context& ectx)
2029  {
2030  return getValue(ectx.fs.pressure(waterPhaseIdx));
2031  }
2032  }
2033  },
2034  Entry{ScalarEntry{"BGPR",
2035  [](const Context& ectx)
2036  {
2037  return getValue(ectx.fs.pressure(gasPhaseIdx));
2038  }
2039  }
2040  },
2041  Entry{PhaseEntry{std::array{
2042  std::array{"BVWAT"sv, "BVOIL"sv, "BVGAS"sv},
2043  std::array{"BWVIS"sv, "BOVIS"sv, "BGVIS"sv}
2044  },
2045  [](const unsigned phaseIdx, const Context& ectx)
2046  {
2047  return getValue(ectx.fs.viscosity(phaseIdx));
2048  }
2049  }
2050  },
2051  Entry{PhaseEntry{std::array{
2052  std::array{"BWDEN"sv, "BODEN"sv, "BGDEN"sv},
2053  std::array{"BDENW"sv, "BDENO"sv, "BDENG"sv}
2054  },
2055  [](const unsigned phaseIdx, const Context& ectx)
2056  {
2057  return getValue(ectx.fs.density(phaseIdx));
2058  }
2059  }
2060  },
2061  Entry{ScalarEntry{"BFLOGI",
2062  [&flowsC = this->flowsC_,
2063  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2064  {
2065  const unsigned index = !flowsC.blockFlows().empty() ?
2066  flowsC.blockFlowsIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2067  FaceDir::ToIntersectionIndex(Dir::XPlus), gasCompIdx) : ectx.globalDofIdx;
2068  return flowsC.getFlow(index, Dir::XPlus, gasCompIdx);
2069  }
2070  }
2071  },
2072  Entry{ScalarEntry{"BFLOGI-",
2073  [&flowsC = this->flowsC_,
2074  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2075  {
2076  const unsigned index = !flowsC.blockFlows().empty() ?
2077  flowsC.blockFlowsIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2078  FaceDir::ToIntersectionIndex(Dir::XMinus), gasCompIdx) : ectx.globalDofIdx;
2079  return flowsC.getFlow(index, Dir::XMinus, gasCompIdx);
2080  }
2081  }
2082  },
2083  Entry{ScalarEntry{"BFLOGJ",
2084  [&flowsC = this->flowsC_,
2085  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2086  {
2087  const unsigned index = !flowsC.blockFlows().empty() ?
2088  flowsC.blockFlowsIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2089  FaceDir::ToIntersectionIndex(Dir::YPlus), gasCompIdx) : ectx.globalDofIdx;
2090  return flowsC.getFlow(index, Dir::YPlus, gasCompIdx);
2091  }
2092  }
2093  },
2094  Entry{ScalarEntry{"BFLOGJ-",
2095  [&flowsC = this->flowsC_,
2096  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2097  {
2098  const unsigned index = !flowsC.blockFlows().empty() ?
2099  flowsC.blockFlowsIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2100  FaceDir::ToIntersectionIndex(Dir::YMinus), gasCompIdx) : ectx.globalDofIdx;
2101  return flowsC.getFlow(index, Dir::YMinus, gasCompIdx);
2102  }
2103  }
2104  },
2105  Entry{ScalarEntry{"BFLOGK",
2106  [&flowsC = this->flowsC_,
2107  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2108  {
2109  const unsigned index = !flowsC.blockFlows().empty() ?
2110  flowsC.blockFlowsIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2111  FaceDir::ToIntersectionIndex(Dir::ZPlus), gasCompIdx) : ectx.globalDofIdx;
2112  return flowsC.getFlow(index, Dir::ZPlus, gasCompIdx);
2113  }
2114  }
2115  },
2116  Entry{ScalarEntry{"BFLOGK-",
2117  [&flowsC = this->flowsC_,
2118  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2119  {
2120  const unsigned index = !flowsC.blockFlows().empty() ?
2121  flowsC.blockFlowsIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2122  FaceDir::ToIntersectionIndex(Dir::ZMinus), gasCompIdx) : ectx.globalDofIdx;
2123  return flowsC.getFlow(index, Dir::ZMinus, gasCompIdx);
2124  }
2125  }
2126  },
2127  Entry{ScalarEntry{"BFLOOI",
2128  [&flowsC = this->flowsC_,
2129  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2130  {
2131  const unsigned index = !flowsC.blockFlows().empty() ?
2132  flowsC.blockFlowsIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2133  FaceDir::ToIntersectionIndex(Dir::XPlus), oilCompIdx) : ectx.globalDofIdx;
2134  return flowsC.getFlow(index, Dir::XPlus, oilCompIdx);
2135  }
2136  }
2137  },
2138  Entry{ScalarEntry{"BFLOOI-",
2139  [&flowsC = this->flowsC_,
2140  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2141  {
2142  const unsigned index = !flowsC.blockFlows().empty() ?
2143  flowsC.blockFlowsIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2144  FaceDir::ToIntersectionIndex(Dir::XMinus), oilCompIdx) : ectx.globalDofIdx;
2145  return flowsC.getFlow(index, Dir::XMinus, oilCompIdx);
2146  }
2147  }
2148  },
2149  Entry{ScalarEntry{"BFLOOJ",
2150  [&flowsC = this->flowsC_,
2151  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2152  {
2153  const unsigned index = !flowsC.blockFlows().empty() ?
2154  flowsC.blockFlowsIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2155  FaceDir::ToIntersectionIndex(Dir::YPlus), oilCompIdx) : ectx.globalDofIdx;
2156  return flowsC.getFlow(index, Dir::YPlus, oilCompIdx);
2157  }
2158  }
2159  },
2160  Entry{ScalarEntry{"BFLOOJ-",
2161  [&flowsC = this->flowsC_,
2162  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2163  {
2164  const unsigned index = !flowsC.blockFlows().empty() ?
2165  flowsC.blockFlowsIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2166  FaceDir::ToIntersectionIndex(Dir::YMinus), oilCompIdx) : ectx.globalDofIdx;
2167  return flowsC.getFlow(index, Dir::YMinus, oilCompIdx);
2168  }
2169  }
2170  },
2171  Entry{ScalarEntry{"BFLOOK",
2172  [&flowsC = this->flowsC_,
2173  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2174  {
2175  const unsigned index = !flowsC.blockFlows().empty() ?
2176  flowsC.blockFlowsIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2177  FaceDir::ToIntersectionIndex(Dir::ZPlus), oilCompIdx) : ectx.globalDofIdx;
2178  return flowsC.getFlow(index, Dir::ZPlus, oilCompIdx);
2179  }
2180  }
2181  },
2182  Entry{ScalarEntry{"BFLOOK-",
2183  [&flowsC = this->flowsC_,
2184  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2185  {
2186  const unsigned index = !flowsC.blockFlows().empty() ?
2187  flowsC.blockFlowsIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2188  FaceDir::ToIntersectionIndex(Dir::ZMinus), oilCompIdx) : ectx.globalDofIdx;
2189  return flowsC.getFlow(index, Dir::ZMinus, oilCompIdx);
2190  }
2191  }
2192  },
2193  Entry{ScalarEntry{"BFLOWI",
2194  [&flowsC = this->flowsC_,
2195  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2196  {
2197  const unsigned index = !flowsC.blockFlows().empty() ?
2198  flowsC.blockFlowsIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2199  FaceDir::ToIntersectionIndex(Dir::XPlus), waterCompIdx) : ectx.globalDofIdx;
2200  return flowsC.getFlow(index, Dir::XPlus, waterCompIdx);
2201  }
2202  }
2203  },
2204  Entry{ScalarEntry{"BFLOWI-",
2205  [&flowsC = this->flowsC_,
2206  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2207  {
2208  const unsigned index = !flowsC.blockFlows().empty() ?
2209  flowsC.blockFlowsIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2210  FaceDir::ToIntersectionIndex(Dir::XMinus), waterCompIdx) : ectx.globalDofIdx;
2211  return flowsC.getFlow(index, Dir::XMinus, waterCompIdx);
2212  }
2213  }
2214  },
2215  Entry{ScalarEntry{"BFLOWJ",
2216  [&flowsC = this->flowsC_,
2217  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2218  {
2219  const unsigned index = !flowsC.blockFlows().empty() ?
2220  flowsC.blockFlowsIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2221  FaceDir::ToIntersectionIndex(Dir::YPlus), waterCompIdx) : ectx.globalDofIdx;
2222  return flowsC.getFlow(index, Dir::YPlus, waterCompIdx);
2223  }
2224  }
2225  },
2226  Entry{ScalarEntry{"BFLOWJ-",
2227  [&flowsC = this->flowsC_,
2228  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2229  {
2230  const unsigned index = !flowsC.blockFlows().empty() ?
2231  flowsC.blockFlowsIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2232  FaceDir::ToIntersectionIndex(Dir::YMinus), waterCompIdx) : ectx.globalDofIdx;
2233  return flowsC.getFlow(index, Dir::YMinus, waterCompIdx);
2234  }
2235  }
2236  },
2237  Entry{ScalarEntry{"BFLOWK",
2238  [&flowsC = this->flowsC_,
2239  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2240  {
2241  const unsigned index = !flowsC.blockFlows().empty() ?
2242  flowsC.blockFlowsIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2243  FaceDir::ToIntersectionIndex(Dir::ZPlus), waterCompIdx) : ectx.globalDofIdx;
2244  return flowsC.getFlow(index, Dir::ZPlus, waterCompIdx);
2245  }
2246  }
2247  },
2248  Entry{ScalarEntry{"BFLOWK-",
2249  [&flowsC = this->flowsC_,
2250  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2251  {
2252  const unsigned index = !flowsC.blockFlows().empty() ?
2253  flowsC.blockFlowsIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2254  FaceDir::ToIntersectionIndex(Dir::ZMinus), waterCompIdx) : ectx.globalDofIdx;
2255  return flowsC.getFlow(index, Dir::ZMinus, waterCompIdx);
2256  }
2257  }
2258  },
2259  Entry{ScalarEntry{"BVELGI",
2260  [&flowsC = this->flowsC_,
2261  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2262  {
2263  const unsigned index = flowsC.blockVelocityIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2264  FaceDir::ToIntersectionIndex(Dir::XPlus), gasCompIdx);
2265  return flowsC.getVelocity(index, Dir::XPlus, gasCompIdx);
2266  }
2267  }
2268  },
2269  Entry{ScalarEntry{"BVELGI-",
2270  [&flowsC = this->flowsC_,
2271  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2272  {
2273  const unsigned index = flowsC.blockVelocityIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2274  FaceDir::ToIntersectionIndex(Dir::XMinus), gasCompIdx);
2275  return flowsC.getVelocity(index, Dir::XMinus, gasCompIdx);
2276  }
2277  }
2278  },
2279  Entry{ScalarEntry{"BVELGJ",
2280  [&flowsC = this->flowsC_,
2281  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2282  {
2283  const unsigned index = flowsC.blockVelocityIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2284  FaceDir::ToIntersectionIndex(Dir::YPlus), gasCompIdx);
2285  return flowsC.getVelocity(index, Dir::YPlus, gasCompIdx);
2286  }
2287  }
2288  },
2289  Entry{ScalarEntry{"BVELGJ-",
2290  [&flowsC = this->flowsC_,
2291  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2292  {
2293  const unsigned index = flowsC.blockVelocityIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2294  FaceDir::ToIntersectionIndex(Dir::YMinus), gasCompIdx);
2295  return flowsC.getVelocity(index, Dir::YMinus, gasCompIdx);
2296  }
2297  }
2298  },
2299  Entry{ScalarEntry{"BVELGK",
2300  [&flowsC = this->flowsC_,
2301  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2302  {
2303  const unsigned index = flowsC.blockVelocityIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2304  FaceDir::ToIntersectionIndex(Dir::ZPlus), gasCompIdx);
2305  return flowsC.getVelocity(index, Dir::ZPlus, gasCompIdx);
2306  }
2307  }
2308  },
2309  Entry{ScalarEntry{"BVELGK-",
2310  [&flowsC = this->flowsC_,
2311  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2312  {
2313  const unsigned index = flowsC.blockVelocityIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2314  FaceDir::ToIntersectionIndex(Dir::ZMinus), gasCompIdx);
2315  return flowsC.getVelocity(index, Dir::ZMinus, gasCompIdx);
2316  }
2317  }
2318  },
2319  Entry{ScalarEntry{"BVELOI",
2320  [&flowsC = this->flowsC_,
2321  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2322  {
2323  const unsigned index = flowsC.blockVelocityIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2324  FaceDir::ToIntersectionIndex(Dir::XPlus), oilCompIdx);
2325  return flowsC.getVelocity(index, Dir::XPlus, oilCompIdx);
2326  }
2327  }
2328  },
2329  Entry{ScalarEntry{"BVELOI-",
2330  [&flowsC = this->flowsC_,
2331  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2332  {
2333  const unsigned index = flowsC.blockVelocityIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2334  FaceDir::ToIntersectionIndex(Dir::XMinus), oilCompIdx);
2335  return flowsC.getVelocity(index, Dir::XMinus, oilCompIdx);
2336  }
2337  }
2338  },
2339  Entry{ScalarEntry{"BVELOJ",
2340  [&flowsC = this->flowsC_,
2341  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2342  {
2343  const unsigned index = flowsC.blockVelocityIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2344  FaceDir::ToIntersectionIndex(Dir::YPlus), oilCompIdx);
2345  return flowsC.getVelocity(index, Dir::YPlus, oilCompIdx);
2346  }
2347  }
2348  },
2349  Entry{ScalarEntry{"BVELOJ-",
2350  [&flowsC = this->flowsC_,
2351  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2352  {
2353  const unsigned index = flowsC.blockVelocityIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2354  FaceDir::ToIntersectionIndex(Dir::YMinus), oilCompIdx);
2355  return flowsC.getVelocity(index, Dir::YMinus, oilCompIdx);
2356  }
2357  }
2358  },
2359  Entry{ScalarEntry{"BVELOK",
2360  [&flowsC = this->flowsC_,
2361  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2362  {
2363  const unsigned index = flowsC.blockVelocityIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2364  FaceDir::ToIntersectionIndex(Dir::ZPlus), oilCompIdx);
2365  return flowsC.getVelocity(index, Dir::ZPlus, oilCompIdx);
2366  }
2367  }
2368  },
2369  Entry{ScalarEntry{"BVELOK-",
2370  [&flowsC = this->flowsC_,
2371  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2372  {
2373  const unsigned index = flowsC.blockVelocityIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2374  FaceDir::ToIntersectionIndex(Dir::ZMinus), oilCompIdx);
2375  return flowsC.getVelocity(index, Dir::ZMinus, oilCompIdx);
2376  }
2377  }
2378  },
2379  Entry{ScalarEntry{"BVELWI",
2380  [&flowsC = this->flowsC_,
2381  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2382  {
2383  const unsigned index = flowsC.blockVelocityIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2384  FaceDir::ToIntersectionIndex(Dir::XPlus), waterCompIdx);
2385  return flowsC.getVelocity(index, Dir::XPlus, waterCompIdx);
2386  }
2387  }
2388  },
2389  Entry{ScalarEntry{"BVELWI-",
2390  [&flowsC = this->flowsC_,
2391  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2392  {
2393  const unsigned index = flowsC.blockVelocityIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2394  FaceDir::ToIntersectionIndex(Dir::XMinus), waterCompIdx);
2395  return flowsC.getVelocity(index, Dir::XMinus, waterCompIdx);
2396  }
2397  }
2398  },
2399  Entry{ScalarEntry{"BVELWJ",
2400  [&flowsC = this->flowsC_,
2401  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2402  {
2403  const unsigned index = flowsC.blockVelocityIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2404  FaceDir::ToIntersectionIndex(Dir::YPlus), waterCompIdx);
2405  return flowsC.getVelocity(index, Dir::YPlus, waterCompIdx);
2406  }
2407  }
2408  },
2409  Entry{ScalarEntry{"BVELWJ-",
2410  [&flowsC = this->flowsC_,
2411  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2412  {
2413  const unsigned index = flowsC.blockVelocityIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2414  FaceDir::ToIntersectionIndex(Dir::YMinus), waterCompIdx);
2415  return flowsC.getVelocity(index, Dir::YMinus, waterCompIdx);
2416  }
2417  }
2418  },
2419  Entry{ScalarEntry{"BVELWK",
2420  [&flowsC = this->flowsC_,
2421  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2422  {
2423  const unsigned index = flowsC.blockVelocityIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2424  FaceDir::ToIntersectionIndex(Dir::ZPlus), waterCompIdx);
2425  return flowsC.getVelocity(index, Dir::ZPlus, waterCompIdx);
2426  }
2427  }
2428  },
2429  Entry{ScalarEntry{"BVELWK-",
2430  [&flowsC = this->flowsC_,
2431  &vanguard = this->simulator_.vanguard()](const Context& ectx)
2432  {
2433  const unsigned index = flowsC.blockVelocityIds(vanguard.cartesianIndex(ectx.globalDofIdx),
2434  FaceDir::ToIntersectionIndex(Dir::ZMinus), waterCompIdx);
2435  return flowsC.getVelocity(index, Dir::ZMinus, waterCompIdx);
2436  }
2437  }
2438  },
2439  Entry{ScalarEntry{"BRPV",
2440  [&model = this->simulator_.model()](const Context& ectx)
2441  {
2442  return getValue(ectx.intQuants.porosity()) *
2443  model.dofTotalVolume(ectx.globalDofIdx);
2444  }
2445  }
2446  },
2447  Entry{PhaseEntry{std::array{"BWPV"sv, "BOPV"sv, "BGPV"sv},
2448  [&model = this->simulator_.model()](const unsigned phaseIdx,
2449  const Context& ectx)
2450  {
2451  return getValue(ectx.fs.saturation(phaseIdx)) *
2452  getValue(ectx.intQuants.porosity()) *
2453  model.dofTotalVolume(ectx.globalDofIdx);
2454  }
2455  }
2456  },
2457  Entry{ScalarEntry{"BRS",
2458  [](const Context& ectx)
2459  {
2460  return getValue(ectx.fs.Rs());
2461  }
2462  }
2463  },
2464  Entry{ScalarEntry{"BRV",
2465  [](const Context& ectx)
2466  {
2467  return getValue(ectx.fs.Rv());
2468  }
2469  }
2470  },
2471  Entry{ScalarEntry{"BOIP",
2472  [&model = this->simulator_.model()](const Context& ectx)
2473  {
2474  return (getValue(ectx.fs.invB(oilPhaseIdx)) *
2475  getValue(ectx.fs.saturation(oilPhaseIdx)) +
2476  getValue(ectx.fs.Rv()) *
2477  getValue(ectx.fs.invB(gasPhaseIdx)) *
2478  getValue(ectx.fs.saturation(gasPhaseIdx))) *
2479  model.dofTotalVolume(ectx.globalDofIdx) *
2480  getValue(ectx.intQuants.porosity());
2481  }
2482  }
2483  },
2484  Entry{ScalarEntry{"BGIP",
2485  [&model = this->simulator_.model()](const Context& ectx)
2486  {
2487  Scalar result = getValue(ectx.fs.invB(gasPhaseIdx)) *
2488  getValue(ectx.fs.saturation(gasPhaseIdx));
2489 
2490  if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
2491  result += getValue(ectx.fs.Rs()) *
2492  getValue(ectx.fs.invB(oilPhaseIdx)) *
2493  getValue(ectx.fs.saturation(oilPhaseIdx));
2494  }
2495  else {
2496  result += getValue(ectx.fs.Rsw()) *
2497  getValue(ectx.fs.invB(waterPhaseIdx)) *
2498  getValue(ectx.fs.saturation(waterPhaseIdx));
2499  }
2500 
2501  return result *
2502  model.dofTotalVolume(ectx.globalDofIdx) *
2503  getValue(ectx.intQuants.porosity());
2504  }
2505  }
2506  },
2507  Entry{ScalarEntry{"BWIP",
2508  [&model = this->simulator_.model()](const Context& ectx)
2509  {
2510  return getValue(ectx.fs.invB(waterPhaseIdx)) *
2511  getValue(ectx.fs.saturation(waterPhaseIdx)) *
2512  model.dofTotalVolume(ectx.globalDofIdx) *
2513  getValue(ectx.intQuants.porosity());
2514  }
2515  }
2516  },
2517  Entry{ScalarEntry{"BOIPL",
2518  [&model = this->simulator_.model()](const Context& ectx)
2519  {
2520  return getValue(ectx.fs.invB(oilPhaseIdx)) *
2521  getValue(ectx.fs.saturation(oilPhaseIdx)) *
2522  model.dofTotalVolume(ectx.globalDofIdx) *
2523  getValue(ectx.intQuants.porosity());
2524  }
2525  }
2526  },
2527  Entry{ScalarEntry{"BGIPL",
2528  [&model = this->simulator_.model()](const Context& ectx)
2529  {
2530  Scalar result;
2531  if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
2532  result = getValue(ectx.fs.Rs()) *
2533  getValue(ectx.fs.invB(oilPhaseIdx)) *
2534  getValue(ectx.fs.saturation(oilPhaseIdx));
2535  }
2536  else {
2537  result = getValue(ectx.fs.Rsw()) *
2538  getValue(ectx.fs.invB(waterPhaseIdx)) *
2539  getValue(ectx.fs.saturation(waterPhaseIdx));
2540  }
2541  return result *
2542  model.dofTotalVolume(ectx.globalDofIdx) *
2543  getValue(ectx.intQuants.porosity());
2544  }
2545  }
2546  },
2547  Entry{ScalarEntry{"BGIPG",
2548  [&model = this->simulator_.model()](const Context& ectx)
2549  {
2550  return getValue(ectx.fs.invB(gasPhaseIdx)) *
2551  getValue(ectx.fs.saturation(gasPhaseIdx)) *
2552  model.dofTotalVolume(ectx.globalDofIdx) *
2553  getValue(ectx.intQuants.porosity());
2554  }
2555  }
2556  },
2557  Entry{ScalarEntry{"BOIPG",
2558  [&model = this->simulator_.model()](const Context& ectx)
2559  {
2560  return getValue(ectx.fs.Rv()) *
2561  getValue(ectx.fs.invB(gasPhaseIdx)) *
2562  getValue(ectx.fs.saturation(gasPhaseIdx)) *
2563  model.dofTotalVolume(ectx.globalDofIdx) *
2564  getValue(ectx.intQuants.porosity());
2565  }
2566  }
2567  },
2568  Entry{PhaseEntry{std::array{"BPPW"sv, "BPPO"sv, "BPPG"sv},
2569  [&simConfig = this->eclState_.getSimulationConfig(),
2570  &grav = this->simulator_.problem().gravity(),
2571  &regionAvgDensity = this->regionAvgDensity_,
2572  &problem = this->simulator_.problem(),
2573  &regions = this->regions_](const unsigned phaseIdx, const Context& ectx)
2574  {
2575  auto phase = RegionPhasePoreVolAverage::Phase{};
2576  phase.ix = phaseIdx;
2577 
2578  // Note different region handling here. FIPNUM is
2579  // one-based, but we need zero-based lookup in
2580  // DatumDepth. On the other hand, pvtRegionIndex is
2581  // zero-based but we need one-based lookup in
2582  // RegionPhasePoreVolAverage.
2583 
2584  // Subtract one to convert FIPNUM to region index.
2585  const auto datum = simConfig.datumDepths()(regions["FIPNUM"][ectx.dofIdx] - 1);
2586 
2587  // Add one to convert region index to region ID.
2588  const auto region = RegionPhasePoreVolAverage::Region {
2589  ectx.elemCtx.primaryVars(ectx.dofIdx, /*timeIdx=*/0).pvtRegionIndex() + 1
2590  };
2591 
2592  const auto density = regionAvgDensity->value("PVTNUM", phase, region);
2593 
2594  const auto press = getValue(ectx.fs.pressure(phase.ix));
2595  const auto dz = problem.dofCenterDepth(ectx.globalDofIdx) - datum;
2596  return press - density*dz*grav[GridView::dimensionworld - 1];
2597  }
2598  }
2599  },
2600  Entry{ScalarEntry{"BAMIP",
2601  [&model = this->simulator_.model()](const Context& ectx)
2602  {
2603  const Scalar rhoW = FluidSystem::referenceDensity(waterPhaseIdx,
2604  ectx.intQuants.pvtRegionIndex());
2605  return getValue(ectx.fs.invB(waterPhaseIdx)) *
2606  getValue(ectx.fs.saturation(waterPhaseIdx)) *
2607  rhoW *
2608  model.dofTotalVolume(ectx.globalDofIdx) *
2609  getValue(ectx.intQuants.porosity());
2610  }
2611  }
2612  },
2613  Entry{ScalarEntry{"BMMIP",
2614  [&model = this->simulator_.model()](const Context& ectx)
2615  {
2616  return getValue(ectx.intQuants.microbialConcentration()) *
2617  getValue(ectx.fs.saturation(waterPhaseIdx)) *
2618  getValue(ectx.intQuants.porosity()) *
2619  model.dofTotalVolume(ectx.globalDofIdx);
2620  }
2621  }
2622  },
2623  Entry{ScalarEntry{"BMOIP",
2624  [&model = this->simulator_.model()](const Context& ectx)
2625  {
2626  return getValue(ectx.intQuants.oxygenConcentration()) *
2627  getValue(ectx.intQuants.porosity()) *
2628  model.dofTotalVolume(ectx.globalDofIdx);
2629  }
2630  }
2631  },
2632  Entry{ScalarEntry{"BMUIP",
2633  [&model = this->simulator_.model()](const Context& ectx)
2634  {
2635  return getValue(ectx.intQuants.ureaConcentration()) *
2636  getValue(ectx.intQuants.porosity()) *
2637  model.dofTotalVolume(ectx.globalDofIdx) * 1;
2638  }
2639  }
2640  },
2641  Entry{ScalarEntry{"BMBIP",
2642  [&model = this->simulator_.model()](const Context& ectx)
2643  {
2644  return model.dofTotalVolume(ectx.globalDofIdx) *
2645  getValue(ectx.intQuants.biofilmMass());
2646  }
2647  }
2648  },
2649  Entry{ScalarEntry{"BMCIP",
2650  [&model = this->simulator_.model()](const Context& ectx)
2651  {
2652  return model.dofTotalVolume(ectx.globalDofIdx) *
2653  getValue(ectx.intQuants.calciteMass());
2654  }
2655  }
2656  },
2657  Entry{ScalarEntry{"BGMIP",
2658  [&model = this->simulator_.model()](const Context& ectx)
2659  {
2660  Scalar result = getValue(ectx.fs.invB(gasPhaseIdx)) *
2661  getValue(ectx.fs.saturation(gasPhaseIdx));
2662 
2663  if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
2664  result += getValue(ectx.fs.Rs()) *
2665  getValue(ectx.fs.invB(oilPhaseIdx)) *
2666  getValue(ectx.fs.saturation(oilPhaseIdx));
2667  }
2668  else {
2669  result += getValue(ectx.fs.Rsw()) *
2670  getValue(ectx.fs.invB(waterPhaseIdx)) *
2671  getValue(ectx.fs.saturation(waterPhaseIdx));
2672  }
2673  const Scalar rhoG = FluidSystem::referenceDensity(gasPhaseIdx,
2674  ectx.intQuants.pvtRegionIndex());
2675  return result *
2676  model.dofTotalVolume(ectx.globalDofIdx) *
2677  getValue(ectx.intQuants.porosity()) *
2678  rhoG;
2679  }
2680  }
2681  },
2682  Entry{ScalarEntry{"BGMGP",
2683  [&model = this->simulator_.model()](const Context& ectx)
2684  {
2685  const Scalar rhoG = FluidSystem::referenceDensity(gasPhaseIdx,
2686  ectx.intQuants.pvtRegionIndex());
2687  return getValue(ectx.fs.invB(gasPhaseIdx)) *
2688  getValue(ectx.fs.saturation(gasPhaseIdx)) *
2689  model.dofTotalVolume(ectx.globalDofIdx) *
2690  getValue(ectx.intQuants.porosity()) *
2691  rhoG;
2692  }
2693  }
2694  },
2695  Entry{ScalarEntry{"BGMDS",
2696  [&model = this->simulator_.model()](const Context& ectx)
2697  {
2698  Scalar result;
2699  if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
2700  result = getValue(ectx.fs.Rs()) *
2701  getValue(ectx.fs.invB(oilPhaseIdx)) *
2702  getValue(ectx.fs.saturation(oilPhaseIdx));
2703  }
2704  else {
2705  result = getValue(ectx.fs.Rsw()) *
2706  getValue(ectx.fs.invB(waterPhaseIdx)) *
2707  getValue(ectx.fs.saturation(waterPhaseIdx));
2708  }
2709  const Scalar rhoG = FluidSystem::referenceDensity(gasPhaseIdx,
2710  ectx.intQuants.pvtRegionIndex());
2711  return result *
2712  model.dofTotalVolume(ectx.globalDofIdx) *
2713  getValue(ectx.intQuants.porosity()) *
2714  rhoG;
2715  }
2716  }
2717  },
2718  Entry{ScalarEntry{"BGMST",
2719  [&model = this->simulator_.model(),
2720  &problem = this->simulator_.problem()](const Context& ectx)
2721  {
2722  const auto& scaledDrainageInfo = problem.materialLawManager()
2723  ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2724  const Scalar sg = getValue(ectx.fs.saturation(gasPhaseIdx));
2725  Scalar strandedGas = scaledDrainageInfo.Sgcr;
2726  if (problem.materialLawManager()->enableHysteresis()) {
2727  const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2728  const Scalar krg = getValue(ectx.intQuants.relativePermeability(gasPhaseIdx));
2729  strandedGas = MaterialLaw::strandedGasSaturation(matParams, sg, krg);
2730  }
2731  const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2732  FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2733  : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2734  return (1.0 - xgW) *
2735  model.dofTotalVolume(ectx.globalDofIdx) *
2736  getValue(ectx.intQuants.porosity()) *
2737  getValue(ectx.fs.density(gasPhaseIdx)) *
2738  std::min(strandedGas, sg);
2739  }
2740  }
2741  },
2742  Entry{ScalarEntry{"BGMUS",
2743  [&model = this->simulator_.model(),
2744  &problem = this->simulator_.problem()](const Context& ectx)
2745  {
2746  const auto& scaledDrainageInfo = problem.materialLawManager()
2747  ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2748  const Scalar sg = getValue(ectx.fs.saturation(gasPhaseIdx));
2749  Scalar strandedGas = scaledDrainageInfo.Sgcr;
2750  if (problem.materialLawManager()->enableHysteresis()) {
2751  const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2752  const Scalar krg = getValue(ectx.intQuants.relativePermeability(gasPhaseIdx));
2753  strandedGas = MaterialLaw::strandedGasSaturation(matParams, sg, krg);
2754  }
2755  const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2756  FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2757  : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2758  return (1.0 - xgW) *
2759  model.dofTotalVolume(ectx.globalDofIdx) *
2760  getValue(ectx.intQuants.porosity()) *
2761  getValue(ectx.fs.density(gasPhaseIdx)) *
2762  std::max(Scalar{0.0}, sg - strandedGas);
2763  }
2764  }
2765  },
2766  Entry{ScalarEntry{"BGMTR",
2767  [&model = this->simulator_.model(),
2768  &problem = this->simulator_.problem()](const Context& ectx)
2769  {
2770  const auto& scaledDrainageInfo = problem.materialLawManager()
2771  ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2772  Scalar trappedGas = scaledDrainageInfo.Sgcr;
2773  if (problem.materialLawManager()->enableHysteresis()) {
2774  const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2775  trappedGas = MaterialLaw::trappedGasSaturation(matParams, /*maxTrapping*/true);
2776  }
2777  const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2778  FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2779  : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2780  return (1.0 - xgW) *
2781  model.dofTotalVolume(ectx.globalDofIdx) *
2782  getValue(ectx.intQuants.porosity()) *
2783  getValue(ectx.fs.density(gasPhaseIdx)) *
2784  std::min(trappedGas, getValue(ectx.fs.saturation(gasPhaseIdx)));
2785  }
2786  }
2787  },
2788  Entry{ScalarEntry{"BGMMO",
2789  [&model = this->simulator_.model(),
2790  &problem = this->simulator_.problem()](const Context& ectx)
2791  {
2792  const auto& scaledDrainageInfo = problem.materialLawManager()
2793  ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2794  Scalar trappedGas = scaledDrainageInfo.Sgcr;
2795  if (problem.materialLawManager()->enableHysteresis()) {
2796  const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2797  trappedGas = MaterialLaw::trappedGasSaturation(matParams, /*maxTrapping*/true);
2798  }
2799  const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2800  FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2801  : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2802  return (1.0 - xgW) *
2803  model.dofTotalVolume(ectx.globalDofIdx) *
2804  getValue(ectx.intQuants.porosity()) *
2805  getValue(ectx.fs.density(gasPhaseIdx)) *
2806  std::max(Scalar{0.0}, getValue(ectx.fs.saturation(gasPhaseIdx)) - trappedGas);
2807  }
2808  }
2809  },
2810  Entry{ScalarEntry{"BGKTR",
2811  [&model = this->simulator_.model(),
2812  &problem = this->simulator_.problem()](const Context& ectx)
2813  {
2814  const auto& scaledDrainageInfo = problem.materialLawManager()
2815  ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2816  const Scalar sg = getValue(ectx.fs.saturation(gasPhaseIdx));
2817  Scalar sgcr = scaledDrainageInfo.Sgcr;
2818  if (problem.materialLawManager()->enableHysteresis()) {
2819  const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2820  sgcr = MaterialLaw::trappedGasSaturation(matParams, /*maxTrapping*/false);
2821  }
2822  if (sg > sgcr) {
2823  return 0.0;
2824  }
2825  else {
2826  const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2827  FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2828  : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2829  return (1.0 - xgW) *
2830  model.dofTotalVolume(ectx.globalDofIdx) *
2831  getValue(ectx.intQuants.porosity()) *
2832  getValue(ectx.fs.density(gasPhaseIdx)) *
2833  getValue(ectx.fs.saturation(gasPhaseIdx));
2834  }
2835  }
2836  }
2837  },
2838  Entry{ScalarEntry{"BGKMO",
2839  [&model = this->simulator_.model(),
2840  &problem = this->simulator_.problem()](const Context& ectx)
2841  {
2842  const auto& scaledDrainageInfo = problem.materialLawManager()
2843  ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2844  const Scalar sg = getValue(ectx.fs.saturation(gasPhaseIdx));
2845  Scalar sgcr = scaledDrainageInfo.Sgcr;
2846  if (problem.materialLawManager()->enableHysteresis()) {
2847  const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2848  sgcr = MaterialLaw::trappedGasSaturation(matParams, /*maxTrapping*/false);
2849  }
2850  if (sgcr >= sg) {
2851  return 0.0;
2852  }
2853  else {
2854  const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2855  FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2856  : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2857  return (1.0 - xgW) *
2858  model.dofTotalVolume(ectx.globalDofIdx) *
2859  getValue(ectx.intQuants.porosity()) *
2860  getValue(ectx.fs.density(gasPhaseIdx)) *
2861  getValue(ectx.fs.saturation(gasPhaseIdx));
2862  }
2863  }
2864  }
2865  },
2866  Entry{ScalarEntry{"BGCDI",
2867  [&model = this->simulator_.model(),
2868  &problem = this->simulator_.problem()](const Context& ectx)
2869  {
2870  const auto& scaledDrainageInfo = problem.materialLawManager()
2871  ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2872  Scalar sgcr = scaledDrainageInfo.Sgcr;
2873  if (problem.materialLawManager()->enableHysteresis()) {
2874  const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2875  sgcr = MaterialLaw::trappedGasSaturation(matParams, /*maxTrapping*/false);
2876  }
2877  const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2878  FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2879  : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2880  return (1.0 - xgW) *
2881  model.dofTotalVolume(ectx.globalDofIdx) *
2882  getValue(ectx.intQuants.porosity()) *
2883  getValue(ectx.fs.density(gasPhaseIdx)) *
2884  std::min(sgcr, getValue(ectx.fs.saturation(gasPhaseIdx))) /
2885  FluidSystem::molarMass(gasCompIdx, ectx.intQuants.pvtRegionIndex());
2886  }
2887  }
2888  },
2889  Entry{ScalarEntry{"BGCDM",
2890  [&model = this->simulator_.model(),
2891  &problem = this->simulator_.problem()](const Context& ectx)
2892  {
2893  const auto& scaledDrainageInfo = problem.materialLawManager()
2894  ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2895  Scalar sgcr = scaledDrainageInfo.Sgcr;
2896  if (problem.materialLawManager()->enableHysteresis()) {
2897  const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2898  sgcr = MaterialLaw::trappedGasSaturation(matParams, /*maxTrapping*/false);
2899  }
2900  const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2901  FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2902  : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2903  return (1.0 - xgW) *
2904  model.dofTotalVolume(ectx.globalDofIdx) *
2905  getValue(ectx.intQuants.porosity()) *
2906  getValue(ectx.fs.density(gasPhaseIdx)) *
2907  std::max(Scalar{0.0}, getValue(ectx.fs.saturation(gasPhaseIdx)) - sgcr) /
2908  FluidSystem::molarMass(gasCompIdx, ectx.intQuants.pvtRegionIndex());
2909  }
2910  }
2911  },
2912  Entry{ScalarEntry{"BGKDI",
2913  [&model = this->simulator_.model(),
2914  &problem = this->simulator_.problem()](const Context& ectx)
2915  {
2916  const auto& scaledDrainageInfo = problem.materialLawManager()
2917  ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2918  const Scalar sg = getValue(ectx.fs.saturation(gasPhaseIdx));
2919  Scalar sgcr = scaledDrainageInfo.Sgcr;
2920  if (problem.materialLawManager()->enableHysteresis()) {
2921  const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2922  sgcr = MaterialLaw::trappedGasSaturation(matParams, /*maxTrapping*/false);
2923  }
2924  if (sg > sgcr) {
2925  return 0.0;
2926  }
2927  else {
2928  const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2929  FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2930  : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2931  return (1.0 - xgW) *
2932  model.dofTotalVolume(ectx.globalDofIdx) *
2933  getValue(ectx.intQuants.porosity()) *
2934  getValue(ectx.fs.density(gasPhaseIdx)) *
2935  getValue(ectx.fs.saturation(gasPhaseIdx)) /
2936  FluidSystem::molarMass(gasCompIdx, ectx.intQuants.pvtRegionIndex());
2937  }
2938  }
2939  }
2940  },
2941  Entry{ScalarEntry{"BGKDM",
2942  [&model = this->simulator_.model(),
2943  &problem = this->simulator_.problem()](const Context& ectx)
2944  {
2945  const auto& scaledDrainageInfo = problem.materialLawManager()
2946  ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
2947  const Scalar sg = getValue(ectx.fs.saturation(gasPhaseIdx));
2948  Scalar sgcr = scaledDrainageInfo.Sgcr;
2949  if (problem.materialLawManager()->enableHysteresis()) {
2950  const auto& matParams = problem.materialLawParams(ectx.dofIdx);
2951  sgcr = MaterialLaw::trappedGasSaturation(matParams, /*maxTrapping*/false);
2952  }
2953  if (sgcr >= sg) {
2954  return 0.0;
2955  }
2956  else {
2957  const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
2958  FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
2959  : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
2960  return (1.0 - xgW) *
2961  model.dofTotalVolume(ectx.globalDofIdx) *
2962  getValue(ectx.intQuants.porosity()) *
2963  getValue(ectx.fs.density(gasPhaseIdx)) *
2964  getValue(ectx.fs.saturation(gasPhaseIdx)) /
2965  FluidSystem::molarMass(gasCompIdx, ectx.intQuants.pvtRegionIndex());
2966  }
2967  }
2968  }
2969  },
2970  Entry{ScalarEntry{"BWCD",
2971  [&model = this->simulator_.model()](const Context& ectx)
2972  {
2973  Scalar result;
2974  if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
2975  result = getValue(ectx.fs.Rs()) *
2976  getValue(ectx.fs.invB(oilPhaseIdx)) *
2977  getValue(ectx.fs.saturation(oilPhaseIdx));
2978  }
2979  else {
2980  result = getValue(ectx.fs.Rsw()) *
2981  getValue(ectx.fs.invB(waterPhaseIdx)) *
2982  getValue(ectx.fs.saturation(waterPhaseIdx));
2983  }
2984  const Scalar rhoG = FluidSystem::referenceDensity(gasPhaseIdx,
2985  ectx.intQuants.pvtRegionIndex());
2986  return result *
2987  model.dofTotalVolume(ectx.globalDofIdx) *
2988  getValue(ectx.intQuants.porosity()) *
2989  rhoG /
2990  FluidSystem::molarMass(gasCompIdx, ectx.intQuants.pvtRegionIndex());
2991  }
2992  }
2993  },
2994  Entry{ScalarEntry{"BWIPG",
2995  [&model = this->simulator_.model()](const Context& ectx)
2996  {
2997  Scalar result = 0.0;
2998  if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
2999  result = getValue(ectx.fs.Rvw()) *
3000  getValue(ectx.fs.invB(gasPhaseIdx)) *
3001  getValue(ectx.fs.saturation(gasPhaseIdx));
3002  }
3003  return result *
3004  model.dofTotalVolume(ectx.globalDofIdx) *
3005  getValue(ectx.intQuants.porosity());
3006  }
3007  }
3008  },
3009  Entry{ScalarEntry{"BWIPL",
3010  [&model = this->simulator_.model()](const Context& ectx)
3011  {
3012  return getValue(ectx.fs.invB(waterPhaseIdx)) *
3013  getValue(ectx.fs.saturation(waterPhaseIdx)) *
3014  model.dofTotalVolume(ectx.globalDofIdx) *
3015  getValue(ectx.intQuants.porosity());
3016  }
3017  }
3018  },
3019  };
3020 
3021  this->blockExtractors_ = BlockExtractor::setupExecMap(this->blockData_, handlers);
3022 
3023  this->extraBlockData_.clear();
3024  if (reportStepNum > 0 && !isSubStep) {
3025  // check we need extra block pressures for RPTSCHED
3026  const auto& rpt = this->schedule_[reportStepNum - 1].rpt_config.get();
3027  if (rpt.contains("WELLS") && rpt.at("WELLS") > 1) {
3028  this->setupExtraBlockData(reportStepNum,
3029  [&c = this->collectOnIORank_](const int idx)
3030  { return c.isCartIdxOnThisRank(idx); });
3031 
3032  const auto extraHandlers = std::array{
3033  pressure_handler,
3034  };
3035 
3036  this->extraBlockExtractors_ = BlockExtractor::setupExecMap(this->extraBlockData_, extraHandlers);
3037  }
3038  }
3039  }
3040 
3041  const Simulator& simulator_;
3042  const CollectDataOnIORankType& collectOnIORank_;
3043  std::vector<typename Extractor::Entry> extractors_;
3044  typename BlockExtractor::ExecMap blockExtractors_;
3045  typename BlockExtractor::ExecMap extraBlockExtractors_;
3046 };
3047 
3048 } // namespace Opm
3049 
3050 #endif // OPM_OUTPUT_BLACK_OIL_MODULE_HPP
Opm::GetPropType
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
Opm::EcfvDiscretization
The base class for the element-centered finite-volume discretization scheme.
Definition: fvbasegradientcalculator.hh:42
Opm::moduleVersionName
std::string moduleVersionName()
Return the version name of the module, for example "2015.10" (for a release branch) or "2016...
Definition: moduleVersion.cpp:34
Opm::CollectDataOnIORank
Definition: CollectDataOnIORank.hpp:55
FlowBaseVanguard.hpp
Helper class for grid instantiation of ECL file-format using problems.
parametersystem.hpp
This file provides the infrastructure to retrieve run-time parameters.
Opm::detail::Extractor
Wrapping struct holding types used for element-level data extraction.
Definition: OutputExtractor.hpp:53
GenericOutputBlackoilModule.hpp
Output module for the results black oil model writing in ECL binary format.
Opm::detail::BlockExtractor::process
static void process(const std::vector< Exec > &blockExtractors, const Context &ectx)
Process a list of block extractors.
Definition: OutputExtractor.hpp:364
Opm::InterRegFlowMap::compress
void compress()
Form CSR adjacency matrix representation of input graph from connections established in previous call...
Definition: InterRegFlows.cpp:164
Opm::OutputBlackOilModule::allocBuffers
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:201
Opm::detail::Extractor::HysteresisParams::sgmax
Scalar sgmax
Max gas saturation.
Definition: OutputExtractor.hpp:67
Opm
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: blackoilbioeffectsmodules.hh:45
Opm::detail::Extractor::HysteresisParams
Struct holding hysteresis parameters.
Definition: OutputExtractor.hpp:62
Opm::detail::BlockExtractor::setupExecMap
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:263
Opm::InterRegFlowMapSingleFIP::Cell
Minimal characteristics of a cell from a simulation grid.
Definition: InterRegFlows.hpp:50
blackoilproperties.hh
Declares the properties required by the black oil model.
Opm::GenericOutputBlackoilModule
Definition: GenericOutputBlackoilModule.hpp:78
fvbaseproperties.hh
Declare the properties used by the infrastructure code of the finite volume discretizations.
Opm::detail::Extractor::HysteresisParams::shmax
Scalar shmax
Max something.
Definition: OutputExtractor.hpp:68
Opm::OutputBlackOilModule::getInterRegFlows
const InterRegFlowMap & getInterRegFlows() const
Get read-only access to collection of inter-region flows.
Definition: OutputBlackoilModule.hpp:504
Opm::InterRegFlowMap
Inter-region flow accumulation maps for all region definition arrays.
Definition: InterRegFlows.hpp:178
OutputExtractor.hpp
Output module for the results black oil model writing in ECL binary format.
Opm::OutputBlackOilModule::initializeFluxData
void initializeFluxData()
Prepare for capturing connection fluxes, particularly to account for inter-region flows...
Definition: OutputBlackoilModule.hpp:486
Opm::detail::BlockExtractor::Entry
std::variant< ScalarEntry, PhaseEntry > Entry
Descriptor for extractors.
Definition: OutputExtractor.hpp:244
Opm::OutputBlackOilModule
Output module for the results black oil model writing in ECL binary format.
Definition: OutputBlackoilModule.hpp:85
Opm::detail::BlockExtractor::ExecMap
std::unordered_map< int, std::vector< Exec > > ExecMap
A map of extraction executors, keyed by cartesian cell index.
Definition: OutputExtractor.hpp:259
Opm::InterRegFlowMap::addConnection
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:155
Opm::OutputBlackOilModule::clearExtractors
void clearExtractors()
Clear list of active element-level data extractors.
Definition: OutputBlackoilModule.hpp:231
Opm::OutputBlackOilModule::finalizeFluxData
void finalizeFluxData()
Finalize capturing connection fluxes.
Definition: OutputBlackoilModule.hpp:496
Opm::OutputBlackOilModule::processElement
void processElement(const ElementContext &elemCtx)
Modify the internal buffers according to the intensive quanties relevant for an element.
Definition: OutputBlackoilModule.hpp:242
Opm::detail::Extractor::HysteresisParams::swmax
Scalar swmax
Max water saturation.
Definition: OutputExtractor.hpp:65
Opm::detail::Extractor::process
static void process(const Context &ectx, const std::vector< Entry > &extractors)
Process the given extractor entries.
Definition: OutputExtractor.hpp:157
Opm::InterRegFlowMap::clear
void clear()
Clear all internal buffers, but preserve allocated capacity.
Definition: InterRegFlows.cpp:171
Opm::detail::Extractor::HysteresisParams::somin
Scalar somin
Min oil saturation.
Definition: OutputExtractor.hpp:69
Opm::detail::BlockExtractor
Wrapping struct holding types used for block-level data extraction.
Definition: OutputExtractor.hpp:191
propertysystem.hh
The Opm property system, traits with inheritance.
Opm::detail::Extractor::removeInactive
static std::vector< Entry > removeInactive(std::array< Entry, size > &input)
Obtain vector of active extractors from an array of extractors.
Definition: OutputExtractor.hpp:120
Opm::detail::Extractor::Context
Context passed to extractor functions.
Definition: OutputExtractor.hpp:73
Opm::detail::Extractor::HysteresisParams::swmin
Scalar swmin
Min water saturation.
Definition: OutputExtractor.hpp:66
blackoilenergymodules.hh
Contains the classes required to extend the black-oil model by energy.
Opm::detail::Extractor::HysteresisParams::somax
Scalar somax
Max oil saturation.
Definition: OutputExtractor.hpp:64
Opm::detail::Extractor::Context::episodeIndex
int episodeIndex
Current report step.
Definition: OutputExtractor.hpp:77
Opm::OutputBlackOilModule::processFluxes
void processFluxes(const ElementContext &elemCtx, ActiveIndex &&activeIndex, CartesianIndex &&cartesianIndex)
Capture connection fluxes, particularly to account for inter-region flows.
Definition: OutputBlackoilModule.hpp:449
Opm::OutputBlackOilModule::setupExtractors
void setupExtractors(const bool isSubStep, const int reportStepNum)
Setup list of active element-level data extractors.
Definition: OutputBlackoilModule.hpp:223