23#ifndef OPM_INIT_STATE_EQUIL_IMPL_HPP
24#define OPM_INIT_STATE_EQUIL_IMPL_HPP
26#include <dune/grid/common/mcmgmapper.hh>
28#include <opm/common/OpmLog/OpmLog.hpp>
30#include <opm/grid/utility/RegionMapping.hpp>
32#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
33#include <opm/input/eclipse/EclipseState/Tables/RsvdTable.hpp>
34#include <opm/input/eclipse/EclipseState/Tables/RvvdTable.hpp>
35#include <opm/input/eclipse/EclipseState/Tables/RvwvdTable.hpp>
36#include <opm/input/eclipse/EclipseState/Tables/PbvdTable.hpp>
37#include <opm/input/eclipse/EclipseState/Tables/PdvdTable.hpp>
38#include <opm/input/eclipse/EclipseState/Tables/SaltvdTable.hpp>
39#include <opm/input/eclipse/EclipseState/Tables/RtempvdTable.hpp>
41#include <opm/input/eclipse/EclipseState/Tables/SaltpvdTable.hpp>
43#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
44#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
48#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
50#include <fmt/format.h>
63template <
typename CellRange,
class Scalar>
64void verticalExtent(
const CellRange& cells,
65 const std::vector<std::pair<Scalar, Scalar>>& cellZMinMax,
66 const Parallel::Communication& comm,
67 std::array<Scalar,2>&
span)
69 span[0] = std::numeric_limits<Scalar>::max();
70 span[1] = std::numeric_limits<Scalar>::lowest();
80 for (
const auto& cell : cells) {
81 if (cellZMinMax[cell].
first <
span[0]) {
span[0] = cellZMinMax[cell].first; }
82 if (cellZMinMax[cell].
second >
span[1]) {
span[1] = cellZMinMax[cell].second; }
89void subdivisionCentrePoints(
const Scalar
left,
92 std::vector<std::pair<Scalar, Scalar>>&
subdiv)
98 const auto start = end;
99 end =
left + (i + 1)*
h;
101 subdiv.emplace_back((start + end) / 2,
h);
105template <
typename CellID,
typename Scalar>
106std::vector<std::pair<Scalar, Scalar>>
107horizontalSubdivision(
const CellID cell,
108 const std::pair<Scalar, Scalar>
topbot,
111 auto subdiv = std::vector<std::pair<Scalar, Scalar>>{};
115 throw std::out_of_range {
116 "Negative thickness (inverted top/bottom faces) in cell "
117 + std::to_string(cell)
127template <
class Scalar,
class Element>
128Scalar cellCenterDepth(
const Element& element)
130 typedef typename Element::Geometry Geometry;
131 static constexpr int zCoord = Element::dimension - 1;
134 const Geometry& geometry = element.geometry();
135 const int corners = geometry.corners();
136 for (
int i=0; i <
corners; ++i)
142template <
class Scalar,
class Element>
143std::pair<Scalar,Scalar> cellZSpan(
const Element& element)
145 typedef typename Element::Geometry Geometry;
146 static constexpr int zCoord = Element::dimension - 1;
150 const Geometry& geometry = element.geometry();
151 const int corners = geometry.corners();
153 for (
int i=0; i < 4; ++i)
155 for (
int i=4; i <
corners; ++i)
158 return std::make_pair(
bot/4,
top/4);
161template <
class Scalar,
class Element>
162std::pair<Scalar,Scalar> cellZMinMax(
const Element& element)
164 typedef typename Element::Geometry Geometry;
165 static constexpr int zCoord = Element::dimension - 1;
166 const Geometry& geometry = element.geometry();
167 const int corners = geometry.corners();
169 auto min = std::numeric_limits<Scalar>::max();
170 auto max = std::numeric_limits<Scalar>::lowest();
173 for (
int i=0; i <
corners; ++i) {
174 min = std::min(
min,
static_cast<Scalar
>(geometry.corner(i)[
zCoord]));
175 max = std::max(
max,
static_cast<Scalar
>(geometry.corner(i)[
zCoord]));
177 return std::make_pair(
min,
max);
180template<
class Scalar,
class RHS>
181RK4IVP<Scalar,RHS>::RK4IVP(
const RHS&
f,
182 const std::array<Scalar,2>&
span,
188 const Scalar h = stepsize();
189 const Scalar h2 = h / 2;
190 const Scalar h6 = h / 6;
196 f_.push_back(f(span_[0], y0));
198 for (
int i = 0; i < N; ++i) {
199 const Scalar x = span_[0] + i*h;
200 const Scalar y = y_.back();
202 const Scalar k1 = f_[i];
203 const Scalar k2 = f(x + h2, y + h2*k1);
204 const Scalar k3 = f(x + h2, y + h2*k2);
205 const Scalar k4 = f(x + h, y + h*k3);
207 y_.push_back(y + h6*(k1 + 2*(k2 + k3) + k4));
208 f_.push_back(f(x + h, y_.back()));
211 assert (y_.size() ==
typename std::vector<Scalar>::size_type(N + 1));
214template<
class Scalar,
class RHS>
215Scalar RK4IVP<Scalar,RHS>::
216operator()(
const Scalar x)
const
220 const Scalar h = stepsize();
221 int i = (x - span_[0]) / h;
222 const Scalar t = (x - (span_[0] + i*h)) / h;
225 if (i < 0) { i = 0; }
226 if (N_ <= i) { i = N_ - 1; }
228 const Scalar y0 = y_[i], y1 = y_[i + 1];
229 const Scalar f0 = f_[i], f1 = f_[i + 1];
231 Scalar u = (1 - 2*t) * (y1 - y0);
232 u += h * ((t - 1)*f0 + t*f1);
234 u += (1 - t)*y0 + t*y1;
239template<
class Scalar,
class RHS>
240Scalar RK4IVP<Scalar,RHS>::
243 return (span_[1] - span_[0]) / N_;
246namespace PhasePressODE {
248template<
class Flu
idSystem>
250Water(
const TabulatedFunction& tempVdTable,
251 const TabulatedFunction& saltVdTable,
252 const int pvtRegionIdx,
253 const Scalar normGrav)
254 : tempVdTable_(tempVdTable)
255 , saltVdTable_(saltVdTable)
256 , pvtRegionIdx_(pvtRegionIdx)
261template<
class Flu
idSystem>
262typename Water<FluidSystem>::Scalar
264operator()(
const Scalar depth,
265 const Scalar press)
const
267 return this->density(depth, press) * g_;
270template<
class Flu
idSystem>
271typename Water<FluidSystem>::Scalar
273density(
const Scalar depth,
274 const Scalar press)
const
277 Scalar saltConcentration = saltVdTable_.eval(depth,
true);
278 Scalar temp = tempVdTable_.eval(depth,
true);
279 Scalar rho = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx_,
284 rho *= FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_);
288template<
class Flu
idSystem,
class RS>
290Oil(
const TabulatedFunction& tempVdTable,
292 const int pvtRegionIdx,
293 const Scalar normGrav)
294 : tempVdTable_(tempVdTable)
296 , pvtRegionIdx_(pvtRegionIdx)
301template<
class Flu
idSystem,
class RS>
302typename Oil<FluidSystem,RS>::Scalar
304operator()(
const Scalar depth,
305 const Scalar press)
const
307 return this->density(depth, press) * g_;
310template<
class Flu
idSystem,
class RS>
311typename Oil<FluidSystem,RS>::Scalar
313density(
const Scalar depth,
314 const Scalar press)
const
316 const Scalar temp = tempVdTable_.eval(depth,
true);
318 if (FluidSystem::enableDissolvedGas())
319 rs = rs_(depth, press, temp);
322 if (rs >= FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press)) {
323 bOil = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp, press);
326 bOil = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, press, rs);
328 Scalar rho = bOil * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_);
329 if (FluidSystem::enableDissolvedGas()) {
330 rho += rs * bOil * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
336template<
class Flu
idSystem,
class RV,
class RVW>
337Gas<FluidSystem,RV,RVW>::
338Gas(
const TabulatedFunction& tempVdTable,
341 const int pvtRegionIdx,
342 const Scalar normGrav)
343 : tempVdTable_(tempVdTable)
346 , pvtRegionIdx_(pvtRegionIdx)
351template<
class Flu
idSystem,
class RV,
class RVW>
352typename Gas<FluidSystem,RV,RVW>::Scalar
353Gas<FluidSystem,RV,RVW>::
354operator()(
const Scalar depth,
355 const Scalar press)
const
357 return this->density(depth, press) * g_;
360template<
class Flu
idSystem,
class RV,
class RVW>
361typename Gas<FluidSystem,RV,RVW>::Scalar
362Gas<FluidSystem,RV,RVW>::
363density(
const Scalar depth,
364 const Scalar press)
const
366 const Scalar temp = tempVdTable_.eval(depth,
true);
368 if (FluidSystem::enableVaporizedOil())
369 rv = rv_(depth, press, temp);
372 if (FluidSystem::enableVaporizedWater())
373 rvw = rvw_(depth, press, temp);
377 if (FluidSystem::enableVaporizedOil() && FluidSystem::enableVaporizedWater()) {
378 if (rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press)
379 && rvw >= FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp, press))
381 bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp, press);
383 bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, press, rv, rvw);
385 Scalar rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
386 rho += rv * bGas * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_)
387 + rvw * bGas * FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_);
391 if (FluidSystem::enableVaporizedOil()){
392 if (rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press)) {
393 bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp, press);
395 bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_,
401 Scalar rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
402 rho += rv * bGas * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_);
406 if (FluidSystem::enableVaporizedWater()){
407 if (rvw >= FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp, press)) {
408 bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp, press);
411 bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_,
417 Scalar rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
418 rho += rvw * bGas * FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_);
423 bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp,
427 Scalar rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
434template<
class Flu
idSystem,
class Region>
436PressureTable<FluidSystem,Region>::
437PressureFunction<ODE>::PressureFunction(
const ODE& ode,
443 this->value_[Direction::Up] = std::make_unique<Distribution>
444 (ode, VSpan {{ ic.depth, span[0] }}, ic.pressure, nsample);
446 this->value_[Direction::Down] = std::make_unique<Distribution>
447 (ode, VSpan {{ ic.depth, span[1] }}, ic.pressure, nsample);
450template<
class Flu
idSystem,
class Region>
452PressureTable<FluidSystem,Region>::
453PressureFunction<ODE>::PressureFunction(
const PressureFunction& rhs)
454 : initial_(rhs.initial_)
456 this->value_[Direction::Up] =
457 std::make_unique<Distribution>(*rhs.value_[Direction::Up]);
459 this->value_[Direction::Down] =
460 std::make_unique<Distribution>(*rhs.value_[Direction::Down]);
463template<
class Flu
idSystem,
class Region>
465typename PressureTable<FluidSystem,Region>::template PressureFunction<ODE>&
466PressureTable<FluidSystem,Region>::
467PressureFunction<ODE>::
468operator=(
const PressureFunction& rhs)
470 this->initial_ = rhs.initial_;
472 this->value_[Direction::Up] =
473 std::make_unique<Distribution>(*rhs.value_[Direction::Up]);
475 this->value_[Direction::Down] =
476 std::make_unique<Distribution>(*rhs.value_[Direction::Down]);
481template<
class Flu
idSystem,
class Region>
483typename PressureTable<FluidSystem,Region>::template PressureFunction<ODE>&
484PressureTable<FluidSystem,Region>::
485PressureFunction<ODE>::
486operator=(PressureFunction&& rhs)
488 this->initial_ = rhs.initial_;
489 this->value_ = std::move(rhs.value_);
494template<
class Flu
idSystem,
class Region>
496typename PressureTable<FluidSystem,Region>::Scalar
497PressureTable<FluidSystem,Region>::
498PressureFunction<ODE>::
499value(
const Scalar depth)
const
501 if (depth < this->initial_.depth) {
503 return (*this->value_[Direction::Up])(depth);
505 else if (depth > this->initial_.depth) {
507 return (*this->value_[Direction::Down])(depth);
511 return this->initial_.pressure;
516template<
class Flu
idSystem,
class Region>
517template<
typename PressFunc>
518void PressureTable<FluidSystem,Region>::
519checkPtr(
const PressFunc* phasePress,
520 const std::string& phaseName)
const
522 if (phasePress !=
nullptr) {
return; }
524 throw std::invalid_argument {
525 "Phase pressure function for \"" + phaseName
526 +
"\" most not be null"
530template<
class Flu
idSystem,
class Region>
531typename PressureTable<FluidSystem,Region>::Strategy
532PressureTable<FluidSystem,Region>::
533selectEquilibrationStrategy(
const Region& reg)
const
535 if (!this->oilActive()) {
536 if (reg.datum() > reg.zwoc()) {
537 return &PressureTable::equil_WOG;
539 return &PressureTable::equil_GOW;
542 if (reg.datum() > reg.zwoc()) {
543 return &PressureTable::equil_WOG;
545 else if (reg.datum() < reg.zgoc()) {
546 return &PressureTable::equil_GOW;
549 return &PressureTable::equil_OWG;
553template<
class Flu
idSystem,
class Region>
554void PressureTable<FluidSystem,Region>::
555copyInPointers(
const PressureTable& rhs)
557 if (rhs.oil_ !=
nullptr) {
558 this->oil_ = std::make_unique<OPress>(*rhs.oil_);
561 if (rhs.gas_ !=
nullptr) {
562 this->gas_ = std::make_unique<GPress>(*rhs.gas_);
565 if (rhs.wat_ !=
nullptr) {
566 this->wat_ = std::make_unique<WPress>(*rhs.wat_);
570template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
573 const std::vector<Scalar>& swatInit)
574 : matLawMgr_(matLawMgr)
575 , swatInit_ (swatInit)
579template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
582 : matLawMgr_(rhs.matLawMgr_)
583 , swatInit_ (rhs.swatInit_)
585 , press_ (rhs.press_)
588 this->setEvaluationPoint(*rhs.evalPt_.position,
590 *rhs.evalPt_.ptable);
593template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
600 this->setEvaluationPoint(x,
reg, ptable);
601 this->initializePhaseQuantities();
603 if (ptable.
gasActive()) { this->deriveGasSat(); }
605 if (ptable.
waterActive()) { this->deriveWaterSat(); }
608 if (this->isOverlappingTransition()) {
609 this->fixUnphysicalTransition();
612 if (ptable.
oilActive()) { this->deriveOilSat(); }
614 this->accountForScaledSaturations();
619template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
623 const PTable& ptable)
625 this->evalPt_.position = &x;
626 this->evalPt_.region = &
reg;
627 this->evalPt_.ptable = &ptable;
630template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
631void PhaseSaturations<MaterialLawManager,FluidSystem,Region,CellID>::
632initializePhaseQuantities()
635 this->press_.reset();
637 const auto depth = this->evalPt_.position->depth;
638 const auto& ptable = *this->evalPt_.ptable;
640 if (ptable.oilActive()) {
641 this->press_.oil = ptable.oil(depth);
644 if (ptable.gasActive()) {
645 this->press_.gas = ptable.gas(depth);
648 if (ptable.waterActive()) {
649 this->press_.water = ptable.water(depth);
653template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
654void PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::deriveOilSat()
656 this->sat_.oil = 1.0 - this->sat_.water - this->sat_.gas;
659template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
660void PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::deriveGasSat()
662 auto&
sg = this->sat_.gas;
665 const auto oilActive = this->evalPt_.ptable->oilActive();
667 if (this->isConstCapPress(this->gasPos())) {
671 const auto gas_contact = oilActive? this->evalPt_.region->zgoc() : this->evalPt_.region->zwoc();
683 const auto pw = oilActive? this->press_.oil : this->press_.water;
684 const auto pcgo = this->press_.gas - pw;
685 sg = this->invertCapPress(
pcgo, this->gasPos(),
isIncr);
689template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
690void PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::deriveWaterSat()
692 auto&
sw = this->sat_.water;
694 const auto oilActive = this->evalPt_.ptable->oilActive();
697 sw = 1.0 - this->sat_.gas;
700 const auto isIncr =
false;
702 if (this->isConstCapPress(this->waterPos())) {
706 sw = this->fromDepthTable(this->evalPt_.region->zwoc(),
719 const auto pcow = this->press_.oil - this->press_.water;
721 if (this->swatInit_.empty()) {
722 sw = this->invertCapPress(
pcow, this->waterPos(),
isIncr);
727 sw = this->invertCapPress(
pcow, this->waterPos(),
isIncr);
736template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
737void PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
738fixUnphysicalTransition()
740 auto&
sg = this->sat_.gas;
741 auto&
sw = this->sat_.water;
749 const auto pcgw = this->press_.gas - this->press_.water;
750 if (! this->swatInit_.empty()) {
756 const auto isIncr =
false;
757 sw = this->invertCapPress(
pcgw, this->waterPos(),
isIncr);
765 (this->matLawMgr_, this->waterPos(), this->gasPos(),
766 this->evalPt_.position->cell,
pcgw);
769 this->fluidState_.setSaturation(this->oilPos(), 1.0 -
sw -
sg);
770 this->fluidState_.setSaturation(this->gasPos(),
sg);
771 this->fluidState_.setSaturation(this->waterPos(), this->evalPt_
772 .ptable->waterActive() ?
sw : 0.0);
775 this->computeMaterialLawCapPress();
776 this->press_.oil = this->press_.gas - this->materialLawCapPressGasOil();
779template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
780void PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
781accountForScaledSaturations()
783 const auto gasActive = this->evalPt_.ptable->gasActive();
784 const auto watActive = this->evalPt_.ptable->waterActive();
785 const auto oilActive = this->evalPt_.ptable->oilActive();
787 auto sg = gasActive? this->sat_.gas : 0.0;
789 auto so = oilActive? this->sat_.oil : 0.0;
791 this->fluidState_.setSaturation(this->waterPos(),
sw);
792 this->fluidState_.setSaturation(this->oilPos(),
so);
793 this->fluidState_.setSaturation(this->gasPos(),
sg);
796 .oilWaterScaledEpsInfoDrainage(this->evalPt_.position->cell);
806 }
else if (gasActive) {
810 this->computeMaterialLawCapPress();
814 this->press_.oil = this->press_.water + this->materialLawCapPressOilWater();
817 this->press_.gas = this->press_.water + this->materialLawCapPressGasWater();
832 this->computeMaterialLawCapPress();
836 this->press_.oil = this->press_.gas - this->materialLawCapPressGasOil();
839 this->press_.water = this->press_.gas - this->materialLawCapPressGasWater();
850 }
else if (gasActive) {
854 this->computeMaterialLawCapPress();
858 this->press_.water = this->press_.oil - this->materialLawCapPressOilWater();
861 this->press_.water = this->press_.gas - this->materialLawCapPressGasWater();
876 this->computeMaterialLawCapPress();
880 this->press_.gas = this->press_.oil + this->materialLawCapPressGasOil();
883 this->press_.gas = this->press_.water + this->materialLawCapPressGasWater();
888template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
889std::pair<typename FluidSystem::Scalar, bool>
890PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
891applySwatInit(
const Scalar
pcow)
893 return this->applySwatInit(
pcow, this->swatInit_[this->evalPt_.position->cell]);
896template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
897std::pair<typename FluidSystem::Scalar, bool>
898PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
899applySwatInit(
const Scalar
pcow,
const Scalar
sw)
901 return this->matLawMgr_.applySwatinit(this->evalPt_.position->cell,
pcow,
sw);
904template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
905void PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
906computeMaterialLawCapPress()
909 .materialLawParams(this->evalPt_.position->cell);
911 this->matLawCapPress_.fill(0.0);
912 MaterialLaw::capillaryPressures(this->matLawCapPress_,
916template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
917typename FluidSystem::Scalar
918PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
919materialLawCapPressGasOil()
const
921 return this->matLawCapPress_[this->oilPos()]
922 + this->matLawCapPress_[this->gasPos()];
925template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
926typename FluidSystem::Scalar
927PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
928materialLawCapPressOilWater()
const
930 return this->matLawCapPress_[this->oilPos()]
931 - this->matLawCapPress_[this->waterPos()];
934template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
935typename FluidSystem::Scalar
936PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
937materialLawCapPressGasWater()
const
939 return this->matLawCapPress_[this->gasPos()]
940 - this->matLawCapPress_[this->waterPos()];
943template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
944bool PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
945isConstCapPress(
const PhaseIdx
phaseIdx)
const
948 (this->matLawMgr_,
phaseIdx, this->evalPt_.position->cell);
951template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
952bool PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
953isOverlappingTransition()
const
955 return this->evalPt_.ptable->gasActive()
956 && this->evalPt_.ptable->waterActive()
957 && ((this->sat_.gas + this->sat_.water) > 1.0);
960template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
961typename FluidSystem::Scalar
962PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
968 (this->matLawMgr_, this->evalPt_.position->depth,
973template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
974typename FluidSystem::Scalar
975PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
976invertCapPress(
const Scalar
pc,
981 (this->matLawMgr_,
static_cast<int>(
phasePos),
982 this->evalPt_.position->cell,
pc,
isincr);
985template<
class Flu
idSystem,
class Region>
994template <
class Flu
idSystem,
class Region>
997 : gravity_(rhs.gravity_)
998 , nsample_(rhs.nsample_)
1000 this->copyInPointers(rhs);
1003template <
class Flu
idSystem,
class Region>
1006 : gravity_(rhs.gravity_)
1007 , nsample_(rhs.nsample_)
1008 , oil_ (std::
move(rhs.oil_))
1009 , gas_ (std::
move(rhs.gas_))
1010 , wat_ (std::
move(rhs.wat_))
1014template <
class Flu
idSystem,
class Region>
1019 this->gravity_ = rhs.gravity_;
1020 this->nsample_ = rhs.nsample_;
1021 this->copyInPointers(rhs);
1026template <
class Flu
idSystem,
class Region>
1031 this->gravity_ = rhs.gravity_;
1032 this->nsample_ = rhs.nsample_;
1034 this->oil_ = std::move(rhs.oil_);
1035 this->gas_ = std::move(rhs.gas_);
1036 this->wat_ = std::move(rhs.wat_);
1041template <
class Flu
idSystem,
class Region>
1047 auto equil = this->selectEquilibrationStrategy(
reg);
1052template <
class Flu
idSystem,
class Region>
1056 return FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
1059template <
class Flu
idSystem,
class Region>
1063 return FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
1066template <
class Flu
idSystem,
class Region>
1070 return FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx);
1073template <
class Flu
idSystem,
class Region>
1074typename FluidSystem::Scalar
1076oil(
const Scalar depth)
const
1078 this->checkPtr(this->oil_.get(),
"OIL");
1080 return this->oil_->value(depth);
1083template <
class Flu
idSystem,
class Region>
1084typename FluidSystem::Scalar
1086gas(
const Scalar depth)
const
1088 this->checkPtr(this->gas_.get(),
"GAS");
1090 return this->gas_->value(depth);
1094template <
class Flu
idSystem,
class Region>
1095typename FluidSystem::Scalar
1097water(
const Scalar depth)
const
1099 this->checkPtr(this->wat_.get(),
"WATER");
1101 return this->wat_->value(depth);
1104template <
class Flu
idSystem,
class Region>
1111 if (! this->waterActive()) {
1112 throw std::invalid_argument {
1113 "Don't know how to interpret EQUIL datum depth in "
1114 "WATER zone in model without active water phase"
1119 const auto ic =
typename WPress::InitCond {
1120 reg.datum(),
reg.pressure()
1126 if (this->oilActive()) {
1128 const auto ic =
typename OPress::InitCond {
1136 if (this->gasActive() && this->oilActive()) {
1138 const auto ic =
typename GPress::InitCond {
1144 }
else if (this->gasActive() && !this->oilActive()) {
1146 const auto ic =
typename GPress::InitCond {
1154template <
class Flu
idSystem,
class Region>
1155void PressureTable<FluidSystem, Region>::
1156equil_GOW(
const Region&
reg,
const VSpan&
span)
1161 if (! this->gasActive()) {
1162 throw std::invalid_argument {
1163 "Don't know how to interpret EQUIL datum depth in "
1164 "GAS zone in model without active gas phase"
1169 const auto ic =
typename GPress::InitCond {
1170 reg.datum(),
reg.pressure()
1176 if (this->oilActive()) {
1178 const auto ic =
typename OPress::InitCond {
1185 if (this->waterActive() && this->oilActive()) {
1187 const auto ic =
typename WPress::InitCond {
1193 }
else if (this->waterActive() && !this->oilActive()) {
1195 const auto ic =
typename WPress::InitCond {
1203template <
class Flu
idSystem,
class Region>
1204void PressureTable<FluidSystem, Region>::
1205equil_OWG(
const Region&
reg,
const VSpan&
span)
1210 if (! this->oilActive()) {
1211 throw std::invalid_argument {
1212 "Don't know how to interpret EQUIL datum depth in "
1213 "OIL zone in model without active oil phase"
1218 const auto ic =
typename OPress::InitCond {
1219 reg.datum(),
reg.pressure()
1225 if (this->waterActive()) {
1227 const auto ic =
typename WPress::InitCond {
1235 if (this->gasActive()) {
1237 const auto ic =
typename GPress::InitCond {
1245template <
class Flu
idSystem,
class Region>
1246void PressureTable<FluidSystem, Region>::
1247makeOilPressure(
const typename OPress::InitCond&
ic,
1251 const auto drho = OilPressODE {
1252 reg.tempVdTable(),
reg.dissolutionCalculator(),
1253 reg.pvtIdx(), this->gravity_
1256 this->oil_ = std::make_unique<OPress>(
drho,
ic, this->nsample_,
span);
1259template <
class Flu
idSystem,
class Region>
1260void PressureTable<FluidSystem, Region>::
1261makeGasPressure(
const typename GPress::InitCond&
ic,
1265 const auto drho = GasPressODE {
1266 reg.tempVdTable(),
reg.evaporationCalculator(),
reg.waterEvaporationCalculator(),
1267 reg.pvtIdx(), this->gravity_
1270 this->gas_ = std::make_unique<GPress>(
drho,
ic, this->nsample_,
span);
1273template <
class Flu
idSystem,
class Region>
1274void PressureTable<FluidSystem, Region>::
1275makeWatPressure(
const typename WPress::InitCond&
ic,
1279 const auto drho = WatPressODE {
1280 reg.tempVdTable(),
reg.saltVdTable(),
reg.pvtIdx(), this->gravity_
1283 this->wat_ = std::make_unique<WPress>(
drho,
ic, this->nsample_,
span);
1288namespace DeckDependent {
1290std::vector<EquilRecord>
1291getEquil(
const EclipseState&
state)
1293 const auto& init =
state.getInitConfig();
1295 if(!init.hasEquil()) {
1296 throw std::domain_error(
"Deck does not provide equilibration data.");
1299 const auto&
equil = init.getEquil();
1303template<
class Gr
idView>
1305equilnum(
const EclipseState& eclipseState,
1310 if (eclipseState.fieldProps().has_int(
"EQLNUM")) {
1311 const auto&
e = eclipseState.fieldProps().get_int(
"EQLNUM");
1312 std::transform(
e.begin(),
e.end(),
eqlnum.begin(), [](
int n){ return n - 1;});
1314 OPM_BEGIN_PARALLEL_TRY_CATCH();
1315 const int num_regions = eclipseState.getTableManager().getEqldims().getNumEquilRegions();
1317 throw std::runtime_error(
"Values larger than maximum Equil regions " + std::to_string(
num_regions) +
" provided in EQLNUM");
1319 if ( std::any_of(
eqlnum.begin(),
eqlnum.end(), [](
int n){return n < 0;}) ) {
1320 throw std::runtime_error(
"zero or negative values provided in EQLNUM");
1322 OPM_END_PARALLEL_TRY_CATCH(
"Invalied EQLNUM numbers: ",
gridview.comm());
1327template<
class FluidSystem,
1330 class ElementMapper,
1331 class CartesianIndexMapper>
1332template<
class MaterialLawManager>
1333InitialStateComputer<FluidSystem,
1337 CartesianIndexMapper>::
1338InitialStateComputer(MaterialLawManager& materialLawManager,
1339 const EclipseState& eclipseState,
1341 const GridView& gridView,
1345 const bool applySwatInit)
1347 saltConcentration_(grid.size(0)),
1348 saltSaturation_(grid.size(0)),
1349 pp_(FluidSystem::numPhases,
1350 std::
vector<Scalar>(grid.size(0))),
1351 sat_(FluidSystem::numPhases,
1352 std::
vector<Scalar>(grid.size(0))),
1360 if (applySwatInit) {
1361 if (eclipseState.fieldProps().has_double(
"SWATINIT")) {
1362 if constexpr (std::is_same_v<Scalar,double>) {
1363 swatInit_ = eclipseState.fieldProps().get_double(
"SWATINIT");
1365 const auto&
input = eclipseState.fieldProps().get_double(
"SWATINIT");
1366 swatInit_.resize(
input.size());
1367 std::copy(
input.begin(),
input.end(), swatInit_.begin());
1374 const auto&
num_aquifers = eclipseState.aquifer().numericalAquifers();
1378 const std::vector<EquilRecord>
rec = getEquil(eclipseState);
1379 const auto&
tables = eclipseState.getTableManager();
1384 setRegionPvtIdx(eclipseState,
eqlmap);
1387 rsFunc_.reserve(
rec.size());
1389 auto getArray = [](
const std::vector<double>&
input)
1391 if constexpr (std::is_same_v<Scalar,double>) {
1394 std::vector<Scalar>
output;
1401 if (FluidSystem::enableDissolvedGas()) {
1402 for (std::size_t i = 0; i <
rec.size(); ++i) {
1403 if (
eqlmap.cells(i).empty()) {
1404 rsFunc_.push_back(std::shared_ptr<Miscibility::RsVD<FluidSystem>>());
1407 const int pvtIdx = regionPvtIdx_[i];
1415 rsFunc_.push_back(std::make_shared<Miscibility::RsVD<FluidSystem>>(pvtIdx,
1421 rsFunc_.push_back(std::make_shared<Miscibility::PBVD<FluidSystem>>(pvtIdx,
1425 throw std::runtime_error(
"Cannot initialise: RSVD or PBVD table not available.");
1431 throw std::runtime_error(
"Cannot initialise: when no explicit RSVD table is given, \n"
1432 "datum depth must be at the gas-oil-contact. "
1433 "In EQUIL region "+std::to_string(i + 1)+
" (counting from 1), this does not hold.");
1435 const Scalar
pContact =
rec[i].datumDepthPressure();
1436 const Scalar
TContact = 273.15 + 20;
1437 rsFunc_.push_back(std::make_shared<Miscibility::RsSatAtContact<FluidSystem>>(pvtIdx,
pContact,
TContact));
1442 for (std::size_t i = 0; i <
rec.size(); ++i) {
1443 rsFunc_.push_back(std::make_shared<Miscibility::NoMixing<Scalar>>());
1447 rvFunc_.reserve(
rec.size());
1448 if (FluidSystem::enableVaporizedOil()) {
1449 for (std::size_t i = 0; i <
rec.size(); ++i) {
1450 if (
eqlmap.cells(i).empty()) {
1451 rvFunc_.push_back(std::shared_ptr<Miscibility::RvVD<FluidSystem>>());
1454 const int pvtIdx = regionPvtIdx_[i];
1463 rvFunc_.push_back(std::make_shared<Miscibility::RvVD<FluidSystem>>(pvtIdx,
1469 rvFunc_.push_back(std::make_shared<Miscibility::PDVD<FluidSystem>>(pvtIdx,
1472 throw std::runtime_error(
"Cannot initialise: RVVD or PDCD table not available.");
1477 throw std::runtime_error(
1478 "Cannot initialise: when no explicit RVVD table is given, \n"
1479 "datum depth must be at the gas-oil-contact. "
1480 "In EQUIL region "+std::to_string(i + 1)+
" (counting from 1), this does not hold.");
1482 const Scalar
pContact =
rec[i].datumDepthPressure() +
rec[i].gasOilContactCapillaryPressure();
1483 const Scalar
TContact = 273.15 + 20;
1484 rvFunc_.push_back(std::make_shared<Miscibility::RvSatAtContact<FluidSystem>>(pvtIdx,
pContact,
TContact));
1489 for (std::size_t i = 0; i <
rec.size(); ++i) {
1490 rvFunc_.push_back(std::make_shared<Miscibility::NoMixing<Scalar>>());
1494 rvwFunc_.reserve(
rec.size());
1495 if (FluidSystem::enableVaporizedWater()) {
1496 for (std::size_t i = 0; i <
rec.size(); ++i) {
1497 if (
eqlmap.cells(i).empty()) {
1498 rvwFunc_.push_back(std::shared_ptr<Miscibility::RvwVD<FluidSystem>>());
1501 const int pvtIdx = regionPvtIdx_[i];
1509 rvwFunc_.push_back(std::make_shared<Miscibility::RvwVD<FluidSystem>>(pvtIdx,
1512 throw std::runtime_error(
"Cannot initialise: RVWVD table not available.");
1516 const auto oilActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
1519 rvwFunc_.push_back(std::make_shared<Miscibility::NoMixing<Scalar>>());
1520 const auto msg =
"No explicit RVWVD table is given for EQUIL region " + std::to_string(i + 1) +
". \n"
1521 "and datum depth is not at the gas-oil-contact. \n"
1522 "Rvw is set to 0.0 in all cells. \n";
1523 OpmLog::warning(
msg);
1527 const Scalar
pContact =
rec[i].datumDepthPressure() +
rec[i].gasOilContactCapillaryPressure();
1528 const Scalar
TContact = 273.15 + 20;
1529 rvwFunc_.push_back(std::make_shared<Miscibility::RvwSatAtContact<FluidSystem>>(pvtIdx,
pContact,
TContact));
1536 rvwFunc_.push_back(std::make_shared<Miscibility::NoMixing<Scalar>>());
1537 const auto msg =
"No explicit RVWVD table is given for EQUIL region " + std::to_string(i + 1) +
". \n"
1538 "and datum depth is not at the gas-water-contact. \n"
1539 "Rvw is set to 0.0 in all cells. \n";
1540 OpmLog::warning(
msg);
1543 const Scalar
pContact =
rec[i].datumDepthPressure() +
rec[i].waterOilContactCapillaryPressure();
1544 const Scalar
TContact = 273.15 + 20;
1545 rvwFunc_.push_back(std::make_shared<Miscibility::RvwSatAtContact<FluidSystem>>(pvtIdx,
pContact,
TContact));
1552 for (std::size_t i = 0; i <
rec.size(); ++i) {
1553 rvwFunc_.push_back(std::make_shared<Miscibility::NoMixing<Scalar>>());
1559 updateInitialTemperature_(eclipseState,
eqlmap);
1562 updateInitialSaltConcentration_(eclipseState,
eqlmap);
1565 updateInitialSaltSaturation_(eclipseState,
eqlmap);
1568 const auto& comm = grid.comm();
1569 calcPressSatRsRv(
eqlmap,
rec, materialLawManager, comm,
grav);
1572 applyNumericalAquifers_(gridView,
num_aquifers, eclipseState.runspec().co2Storage() || eclipseState.runspec().h2Storage());
1578template<
class FluidSystem,
1581 class ElementMapper,
1582 class CartesianIndexMapper>
1584void InitialStateComputer<FluidSystem,
1588 CartesianIndexMapper>::
1589updateInitialTemperature_(
const EclipseState& eclState,
const RMap&
reg)
1593 const auto&
tables = eclState.getTableManager();
1594 if (!
tables.hasTables(
"RTEMPVD")) {
1595 std::vector<Scalar> x = {0.0,1.0};
1596 std::vector<Scalar>
y = {
static_cast<Scalar
>(
tables.rtemp()),
1597 static_cast<Scalar
>(
tables.rtemp())};
1598 for (
auto&
table :
this->tempVdTable_) {
1599 table.setXYContainers(x,
y);
1603 for (std::size_t i = 0; i <
tempvdTables.size(); ++i) {
1606 const auto& cells =
reg.cells(i);
1607 for (
const auto& cell : cells) {
1608 const Scalar depth = cellCenterDepth_[cell];
1609 this->temperature_[cell] = tempVdTable_[i].eval(depth,
true);
1615template<
class FluidSystem,
1618 class ElementMapper,
1619 class CartesianIndexMapper>
1621void InitialStateComputer<FluidSystem,
1625 CartesianIndexMapper>::
1626updateInitialSaltConcentration_(
const EclipseState& eclState,
const RMap&
reg)
1630 const auto&
tables = eclState.getTableManager();
1635 std::vector<Scalar> x = {0.0,1.0};
1636 std::vector<Scalar>
y = {0.0,0.0};
1637 for (
auto&
table :
this->saltVdTable_) {
1638 table.setXYContainers(x,
y);
1641 for (std::size_t i = 0; i <
saltvdTables.size(); ++i) {
1645 const auto& cells =
reg.cells(i);
1646 for (
const auto& cell : cells) {
1647 const Scalar depth = cellCenterDepth_[cell];
1648 this->saltConcentration_[cell] = saltVdTable_[i].eval(depth,
true);
1654template<
class FluidSystem,
1657 class ElementMapper,
1658 class CartesianIndexMapper>
1660void InitialStateComputer<FluidSystem,
1664 CartesianIndexMapper>::
1665updateInitialSaltSaturation_(
const EclipseState& eclState,
const RMap&
reg)
1669 const auto&
tables = eclState.getTableManager();
1676 const auto& cells =
reg.cells(i);
1677 for (
const auto& cell : cells) {
1678 const Scalar depth = cellCenterDepth_[cell];
1679 this->saltSaturation_[cell] = saltpVdTable_[i].eval(depth,
true);
1685template<
class FluidSystem,
1688 class ElementMapper,
1689 class CartesianIndexMapper>
1690void InitialStateComputer<FluidSystem,
1694 CartesianIndexMapper>::
1695updateCellProps_(
const GridView& gridView,
1698 ElementMapper
elemMapper(gridView, Dune::mcmgElementLayout());
1704 auto elemIt = gridView.template begin<0>();
1705 const auto&
elemEndIt = gridView.template end<0>();
1708 const Element& element = *
elemIt;
1710 cellCenterDepth_[
elemIdx] = Details::cellCenterDepth<Scalar>(element);
1711 const auto cartIx = cartesianIndexMapper_.cartesianIndex(
elemIdx);
1712 cellZSpan_[
elemIdx] = Details::cellZSpan<Scalar>(element);
1713 cellZMinMax_[
elemIdx] = Details::cellZMinMax<Scalar>(element);
1729template<
class FluidSystem,
1732 class ElementMapper,
1733 class CartesianIndexMapper>
1734void InitialStateComputer<FluidSystem,
1738 CartesianIndexMapper>::
1739applyNumericalAquifers_(
const GridView& gridView,
1748 const auto watPos =
oil_as_brine? FluidSystem::oilPhaseIdx : FluidSystem::waterPhaseIdx;
1749 if (!FluidSystem::phaseIsActive(
watPos)){
1750 throw std::logic_error {
"Water phase has to be active for numerical aquifer case" };
1753 ElementMapper
elemMapper(gridView, Dune::mcmgElementLayout());
1754 auto elemIt = gridView.template begin<0>();
1755 const auto&
elemEndIt = gridView.template end<0>();
1756 const auto oilPos = FluidSystem::oilPhaseIdx;
1757 const auto gasPos = FluidSystem::gasPhaseIdx;
1759 const Element& element = *
elemIt;
1761 const auto cartIx = cartesianIndexMapper_.cartesianIndex(
elemIdx);
1768 this->sat_[oilPos][
elemIdx] = 0.;
1771 if (FluidSystem::phaseIsActive(gasPos)) {
1772 this->sat_[gasPos][
elemIdx] = 0.;
1775 const auto msg = fmt::format(
"FOR AQUIFER CELL AT ({}, {}, {}) OF NUMERICAL "
1776 "AQUIFER {}, WATER SATURATION IS SET TO BE UNITY",
1785 if (FluidSystem::phaseIsActive(gasPos)) {
1788 if (FluidSystem::phaseIsActive(oilPos)) {
1796template<
class FluidSystem,
1799 class ElementMapper,
1800 class CartesianIndexMapper>
1802void InitialStateComputer<FluidSystem,
1806 CartesianIndexMapper>::
1807setRegionPvtIdx(
const EclipseState& eclState,
const RMap&
reg)
1809 const auto&
pvtnumData = eclState.fieldProps().get_int(
"PVTNUM");
1812 const auto& cells =
reg.cells(
r);
1813 regionPvtIdx_[
r] =
pvtnumData[*cells.begin()] - 1;
1817template<
class FluidSystem,
1820 class ElementMapper,
1821 class CartesianIndexMapper>
1822template<
class RMap,
class MaterialLawManager,
class Comm>
1823void InitialStateComputer<FluidSystem,
1827 CartesianIndexMapper>::
1828calcPressSatRsRv(
const RMap&
reg,
1829 const std::vector<EquilRecord>&
rec,
1830 MaterialLawManager& materialLawManager,
1834 using PhaseSat = Details::PhaseSaturations<
1838 auto ptable = Details::PressureTable<FluidSystem, EquilReg<Scalar>>{
grav, this->num_pressure_points_ };
1839 auto psat =
PhaseSat { materialLawManager, this->swatInit_ };
1840 auto vspan = std::array<Scalar, 2>{};
1843 for (std::size_t
r = 0;
r <
rec.size(); ++
r) {
1844 const auto& cells =
reg.cells(
r);
1846 Details::verticalExtent(cells, cellZMinMax_, comm,
vspan);
1848 const auto acc =
rec[
r].initializationTargetAccuracy();
1850 throw std::runtime_error {
1851 "Cannot initialise model: Positive item 9 is not supported "
1852 "in EQUIL keyword, record " + std::to_string(
r + 1)
1856 if (cells.empty()) {
1861 const auto eqreg = EquilReg {
1862 rec[
r], this->rsFunc_[
r], this->rvFunc_[
r], this->rvwFunc_[
r], this->tempVdTable_[
r], this->saltVdTable_[
r], this->regionPvtIdx_[
r]
1874 this->equilibrateCellCentres(cells,
eqreg, ptable,
psat);
1878 this->equilibrateHorizontal(cells,
eqreg, -
acc,
1888 if (comm.rank() == 0) {
1889 for (std::size_t
r = 0;
r <
rec.size(); ++
r) {
1891 OpmLog::warning(
"Equilibration region " + std::to_string(
r + 1)
1892 +
" has no active cells");
1897template<
class FluidSystem,
1900 class ElementMapper,
1901 class CartesianIndexMapper>
1902template<
class CellRange,
class EquilibrationMethod>
1903void InitialStateComputer<FluidSystem,
1907 CartesianIndexMapper>::
1911 const auto oilPos = FluidSystem::oilPhaseIdx;
1912 const auto gasPos = FluidSystem::gasPhaseIdx;
1913 const auto watPos = FluidSystem::waterPhaseIdx;
1915 const auto oilActive = FluidSystem::phaseIsActive(oilPos);
1916 const auto gasActive = FluidSystem::phaseIsActive(gasPos);
1919 auto pressures = Details::PhaseQuantityValue<Scalar>{};
1920 auto saturations = Details::PhaseQuantityValue<Scalar>{};
1925 for (
const auto& cell : cells) {
1929 this->pp_ [oilPos][cell] =
pressures.oil;
1934 this->pp_ [gasPos][cell] =
pressures.gas;
1943 if (oilActive && gasActive) {
1944 this->rs_[cell] = Rs;
1945 this->rv_[cell] = Rv;
1949 this->rvw_[cell] = Rvw;
1954template<
class FluidSystem,
1957 class ElementMapper,
1958 class CartesianIndexMapper>
1959template<
class CellRange,
class PressTable,
class PhaseSat>
1960void InitialStateComputer<FluidSystem,
1964 CartesianIndexMapper>::
1965equilibrateCellCentres(
const CellRange& cells,
1970 using CellPos =
typename PhaseSat::Position;
1971 using CellID = std::remove_cv_t<std::remove_reference_t<
1972 decltype(std::declval<CellPos>().cell)>>;
1973 this->cellLoop(cells, [
this, &
eqreg, &ptable, &
psat]
1975 Details::PhaseQuantityValue<Scalar>&
pressures,
1979 Scalar& Rvw) ->
void
1982 cell, cellCenterDepth_[cell]
1988 const auto temp = this->temperature_[cell];
1990 Rs =
eqreg.dissolutionCalculator()
1993 Rv =
eqreg.evaporationCalculator()
1996 Rvw =
eqreg.waterEvaporationCalculator()
2001template<
class FluidSystem,
2004 class ElementMapper,
2005 class CartesianIndexMapper>
2006template<
class CellRange,
class PressTable,
class PhaseSat>
2007void InitialStateComputer<FluidSystem,
2011 CartesianIndexMapper>::
2012equilibrateHorizontal(
const CellRange& cells,
2018 using CellPos =
typename PhaseSat::Position;
2019 using CellID = std::remove_cv_t<std::remove_reference_t<
2020 decltype(std::declval<CellPos>().cell)>>;
2022 this->cellLoop(cells, [
this,
acc, &
eqreg, &ptable, &
psat]
2024 Details::PhaseQuantityValue<Scalar>&
pressures,
2028 Scalar& Rvw) ->
void
2034 for (
const auto& [depth,
frac] : Details::horizontalSubdivision(cell, cellZSpan_[cell],
acc)) {
2035 const auto pos =
CellPos { cell, depth };
2049 cell, cellCenterDepth_[cell]
2056 const auto temp = this->temperature_[cell];
2057 const auto cz = cellCenterDepth_[cell];
2059 Rs =
eqreg.dissolutionCalculator()
2062 Rv =
eqreg.evaporationCalculator()
2065 Rvw =
eqreg.waterEvaporationCalculator()
Auxiliary routines that to solve the ODEs that emerge from the hydrostatic equilibrium problem.
Routines that actually solve the ODEs that emerge from the hydrostatic equilibrium problem.
Calculator for phase saturations.
Definition InitStateEquil.hpp:386
const PhaseQuantityValue< Scalar > & deriveSaturations(const Position &x, const Region ®, const PTable &ptable)
Calculate phase saturations at particular point of the simulation model geometry.
Definition InitStateEquil_impl.hpp:596
PhaseSaturations(MaterialLawManager &matLawMgr, const std::vector< Scalar > &swatInit)
Constructor.
Definition InitStateEquil_impl.hpp:572
Definition InitStateEquil.hpp:167
PressureTable & operator=(const PressureTable &rhs)
Assignment operator.
Definition InitStateEquil_impl.hpp:1017
Scalar water(const Scalar depth) const
Evaluate water phase pressure at specified depth.
Definition InitStateEquil_impl.hpp:1097
Scalar gas(const Scalar depth) const
Evaluate gas phase pressure at specified depth.
Definition InitStateEquil_impl.hpp:1086
bool waterActive() const
Predicate for whether or not water is an active phase.
Definition InitStateEquil_impl.hpp:1068
bool gasActive() const
Predicate for whether or not gas is an active phase.
Definition InitStateEquil_impl.hpp:1061
Scalar oil(const Scalar depth) const
Evaluate oil phase pressure at specified depth.
Definition InitStateEquil_impl.hpp:1076
bool oilActive() const
Predicate for whether or not oil is an active phase.
Definition InitStateEquil_impl.hpp:1054
PressureTable(const Scalar gravity, const int samplePoints=2000)
Constructor.
Definition InitStateEquil_impl.hpp:987
bool water(const PhaseUsage &pu)
Active water predicate.
Definition RegionAttributeHelpers.hpp:309
bool oil(const PhaseUsage &pu)
Active oil predicate.
Definition RegionAttributeHelpers.hpp:322
bool gas(const PhaseUsage &pu)
Active gas predicate.
Definition RegionAttributeHelpers.hpp:335
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
Simple set of per-phase (named by primary component) quantities.
Definition InitStateEquil.hpp:338
Evaluation point within a model geometry.
Definition InitStateEquil.hpp:392