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