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; }
84 span[0] = comm.min(span[0]);
85 span[1] = comm.max(span[1]);
89void subdivisionCentrePoints(
const Scalar left,
91 const int numIntervals,
92 std::vector<std::pair<Scalar, Scalar>>& subdiv)
94 const auto h = (right - left) / numIntervals;
97 for (
auto i = 0*numIntervals; i < numIntervals; ++i) {
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,
109 const int numIntervals)
111 auto subdiv = std::vector<std::pair<Scalar, Scalar>>{};
112 subdiv.reserve(2 * numIntervals);
114 if (topbot.first > topbot.second) {
115 throw std::out_of_range {
116 "Negative thickness (inverted top/bottom faces) in cell "
117 + std::to_string(cell)
121 subdivisionCentrePoints(topbot.first, topbot.second,
122 2*numIntervals, subdiv);
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)
137 zz += geometry.corner(i)[zCoord];
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();
152 assert(corners == 8);
153 for (
int i=0; i < 4; ++i)
154 bot += geometry.corner(i)[zCoord];
155 for (
int i=4; i < corners; ++i)
156 top += geometry.corner(i)[zCoord];
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();
168 assert(corners == 8);
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 = ®
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;
664 const auto isIncr =
true;
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();
672 sg = this->fromDepthTable(gas_contact,
673 this->gasPos(), isIncr);
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(),
707 this->waterPos(), isIncr);
719 const auto pcow = this->press_.oil - this->press_.water;
721 if (this->swatInit_.empty()) {
722 sw = this->invertCapPress(pcow, this->waterPos(), isIncr);
725 auto [swout, newSwatInit] = this->applySwatInit(pcow);
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()) {
754 auto [swout, newSwatInit] = this->applySwatInit(pcgw, sw);
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;
788 auto sw = watActive? this->sat_.water : 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);
795 const auto& scaledDrainageInfo = this->matLawMgr_
796 .oilWaterScaledEpsInfoDrainage(this->evalPt_.position->cell);
798 const auto thresholdSat = 1.0e-6;
799 if (watActive && ((sw + thresholdSat) > scaledDrainageInfo.Swu)) {
803 this->fluidState_.setSaturation(this->waterPos(), scaledDrainageInfo.Swu);
805 this->fluidState_.setSaturation(this->oilPos(), so + sw - scaledDrainageInfo.Swu);
806 }
else if (gasActive) {
807 this->fluidState_.setSaturation(this->gasPos(), sg + sw - scaledDrainageInfo.Swu);
809 sw = scaledDrainageInfo.Swu;
810 this->computeMaterialLawCapPress();
814 this->press_.oil = this->press_.water + this->materialLawCapPressOilWater();
817 this->press_.gas = this->press_.water + this->materialLawCapPressGasWater();
821 if (gasActive && ((sg + thresholdSat) > scaledDrainageInfo.Sgu)) {
825 this->fluidState_.setSaturation(this->gasPos(), scaledDrainageInfo.Sgu);
827 this->fluidState_.setSaturation(this->oilPos(), so + sg - scaledDrainageInfo.Sgu);
828 }
else if (watActive) {
829 this->fluidState_.setSaturation(this->waterPos(), sw + sg - scaledDrainageInfo.Sgu);
831 sg = scaledDrainageInfo.Sgu;
832 this->computeMaterialLawCapPress();
836 this->press_.oil = this->press_.gas - this->materialLawCapPressGasOil();
839 this->press_.water = this->press_.gas - this->materialLawCapPressGasWater();
843 if (watActive && ((sw - thresholdSat) < scaledDrainageInfo.Swl)) {
847 this->fluidState_.setSaturation(this->waterPos(), scaledDrainageInfo.Swl);
849 this->fluidState_.setSaturation(this->oilPos(), so + sw - scaledDrainageInfo.Swl);
850 }
else if (gasActive) {
851 this->fluidState_.setSaturation(this->gasPos(), sg + sw - scaledDrainageInfo.Swl);
853 sw = scaledDrainageInfo.Swl;
854 this->computeMaterialLawCapPress();
858 this->press_.water = this->press_.oil - this->materialLawCapPressOilWater();
861 this->press_.water = this->press_.gas - this->materialLawCapPressGasWater();
865 if (gasActive && ((sg - thresholdSat) < scaledDrainageInfo.Sgl)) {
869 this->fluidState_.setSaturation(this->gasPos(), scaledDrainageInfo.Sgl);
871 this->fluidState_.setSaturation(this->oilPos(), so + sg - scaledDrainageInfo.Sgl);
872 }
else if (watActive) {
873 this->fluidState_.setSaturation(this->waterPos(), sw + sg - scaledDrainageInfo.Sgl);
875 sg = scaledDrainageInfo.Sgl;
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()
908 const auto& matParams = this->matLawMgr_
909 .materialLawParams(this->evalPt_.position->cell);
911 this->matLawCapPress_.fill(0.0);
912 MaterialLaw::capillaryPressures(this->matLawCapPress_,
913 matParams, this->fluidState_);
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>::
963fromDepthTable(
const Scalar contactdepth,
964 const PhaseIdx phasePos,
965 const bool isincr)
const
968 (this->matLawMgr_, this->evalPt_.position->depth,
969 contactdepth,
static_cast<int>(phasePos),
970 this->evalPt_.position->cell, isincr);
973template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
974typename FluidSystem::Scalar
975PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
976invertCapPress(
const Scalar pc,
977 const PhaseIdx phasePos,
978 const bool isincr)
const
981 (this->matLawMgr_,
static_cast<int>(phasePos),
982 this->evalPt_.position->cell, pc, isincr);
985template<
class Flu
idSystem,
class Region>
988 const int samplePoints)
990 , nsample_(samplePoints)
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);
1049 (this->*equil)(reg, span);
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>
1106equil_WOG(
const Region& reg,
const VSpan& span)
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()
1123 this->makeWatPressure(ic, reg, span);
1126 if (this->oilActive()) {
1128 const auto ic =
typename OPress::InitCond {
1130 this->
water(reg.zwoc()) + reg.pcowWoc()
1133 this->makeOilPressure(ic, reg, span);
1136 if (this->gasActive() && this->oilActive()) {
1138 const auto ic =
typename GPress::InitCond {
1140 this->
oil(reg.zgoc()) + reg.pcgoGoc()
1143 this->makeGasPressure(ic, reg, span);
1144 }
else if (this->gasActive() && !this->oilActive()) {
1146 const auto ic =
typename GPress::InitCond {
1148 this->
water(reg.zwoc()) + reg.pcowWoc()
1150 this->makeGasPressure(ic, reg, span);
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()
1173 this->makeGasPressure(ic, reg, span);
1176 if (this->oilActive()) {
1178 const auto ic =
typename OPress::InitCond {
1180 this->
gas(reg.zgoc()) - reg.pcgoGoc()
1182 this->makeOilPressure(ic, reg, span);
1185 if (this->waterActive() && this->oilActive()) {
1187 const auto ic =
typename WPress::InitCond {
1189 this->
oil(reg.zwoc()) - reg.pcowWoc()
1192 this->makeWatPressure(ic, reg, span);
1193 }
else if (this->waterActive() && !this->oilActive()) {
1195 const auto ic =
typename WPress::InitCond {
1197 this->
gas(reg.zwoc()) - reg.pcowWoc()
1199 this->makeWatPressure(ic, reg, span);
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()
1222 this->makeOilPressure(ic, reg, span);
1225 if (this->waterActive()) {
1227 const auto ic =
typename WPress::InitCond {
1229 this->
oil(reg.zwoc()) - reg.pcowWoc()
1232 this->makeWatPressure(ic, reg, span);
1235 if (this->gasActive()) {
1237 const auto ic =
typename GPress::InitCond {
1239 this->
oil(reg.zgoc()) + reg.pcgoGoc()
1241 this->makeGasPressure(ic, reg, span);
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();
1300 return { equil.begin(), equil.end() };
1303template<
class Gr
idView>
1305equilnum(
const EclipseState& eclipseState,
1306 const GridView& gridview)
1308 std::vector<int> eqlnum(gridview.size(0), 0);
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();
1316 if ( std::any_of(eqlnum.begin(), eqlnum.end(), [num_regions](
int n){return n >= num_regions;}) ) {
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,
1342 const CartesianIndexMapper& cartMapper,
1344 const int num_pressure_points,
1345 const bool applySwatInit)
1346 : temperature_(grid.size(0), eclipseState.getTableManager().rtemp()),
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))),
1356 cartesianIndexMapper_(cartMapper),
1357 num_pressure_points_(num_pressure_points)
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();
1375 updateCellProps_(gridView, num_aquifers);
1378 const std::vector<EquilRecord> rec = getEquil(eclipseState);
1379 const auto& tables = eclipseState.getTableManager();
1381 const RegionMapping<> eqlmap(equilnum(eclipseState, grid));
1382 const int invalidRegion = -1;
1383 regionPvtIdx_.resize(rec.size(), invalidRegion);
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;
1395 output.resize(input.size());
1396 std::copy(input.begin(), input.end(), output.begin());
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];
1408 if (!rec[i].liveOilInitConstantRs()) {
1409 const TableContainer& rsvdTables = tables.getRsvdTables();
1410 const TableContainer& pbvdTables = tables.getPbvdTables();
1411 if (rsvdTables.size() > 0) {
1412 const RsvdTable& rsvdTable = rsvdTables.getTable<RsvdTable>(i);
1413 auto depthColumn = getArray(rsvdTable.getColumn(
"DEPTH").vectorCopy());
1414 auto rsColumn = getArray(rsvdTable.getColumn(
"RS").vectorCopy());
1415 rsFunc_.push_back(std::make_shared<Miscibility::RsVD<FluidSystem>>(pvtIdx,
1416 depthColumn, rsColumn));
1417 }
else if (pbvdTables.size() > 0) {
1418 const PbvdTable& pbvdTable = pbvdTables.getTable<PbvdTable>(i);
1419 auto depthColumn = getArray(pbvdTable.getColumn(
"DEPTH").vectorCopy());
1420 auto pbubColumn = getArray(pbvdTable.getColumn(
"PBUB").vectorCopy());
1421 rsFunc_.push_back(std::make_shared<Miscibility::PBVD<FluidSystem>>(pvtIdx,
1422 depthColumn, pbubColumn));
1425 throw std::runtime_error(
"Cannot initialise: RSVD or PBVD table not available.");
1430 if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) {
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];
1455 if (!rec[i].wetGasInitConstantRv()) {
1456 const TableContainer& rvvdTables = tables.getRvvdTables();
1457 const TableContainer& pdvdTables = tables.getPdvdTables();
1459 if (rvvdTables.size() > 0) {
1460 const RvvdTable& rvvdTable = rvvdTables.getTable<RvvdTable>(i);
1461 auto depthColumn = getArray(rvvdTable.getColumn(
"DEPTH").vectorCopy());
1462 auto rvColumn = getArray(rvvdTable.getColumn(
"RV").vectorCopy());
1463 rvFunc_.push_back(std::make_shared<Miscibility::RvVD<FluidSystem>>(pvtIdx,
1464 depthColumn, rvColumn));
1465 }
else if (pdvdTables.size() > 0) {
1466 const PdvdTable& pdvdTable = pdvdTables.getTable<PdvdTable>(i);
1467 auto depthColumn = getArray(pdvdTable.getColumn(
"DEPTH").vectorCopy());
1468 auto pdewColumn = getArray(pdvdTable.getColumn(
"PDEW").vectorCopy());
1469 rvFunc_.push_back(std::make_shared<Miscibility::PDVD<FluidSystem>>(pvtIdx,
1470 depthColumn, pdewColumn));
1472 throw std::runtime_error(
"Cannot initialise: RVVD or PDCD table not available.");
1476 if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) {
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];
1502 if (!rec[i].humidGasInitConstantRvw()) {
1503 const TableContainer& rvwvdTables = tables.getRvwvdTables();
1505 if (rvwvdTables.size() > 0) {
1506 const RvwvdTable& rvwvdTable = rvwvdTables.getTable<RvwvdTable>(i);
1507 auto depthColumn = getArray(rvwvdTable.getColumn(
"DEPTH").vectorCopy());
1508 auto rvwvdColumn = getArray(rvwvdTable.getColumn(
"RVWVD").vectorCopy());
1509 rvwFunc_.push_back(std::make_shared<Miscibility::RvwVD<FluidSystem>>(pvtIdx,
1510 depthColumn, rvwvdColumn));
1512 throw std::runtime_error(
"Cannot initialise: RVWVD table not available.");
1516 const auto oilActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
1518 if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) {
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));
1535 if (rec[i].waterOilContactDepth() != rec[i].datumDepth()) {
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)
1591 const int numEquilReg = rsFunc_.size();
1592 tempVdTable_.resize(numEquilReg);
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);
1602 const TableContainer& tempvdTables = tables.getRtempvdTables();
1603 for (std::size_t i = 0; i < tempvdTables.size(); ++i) {
1604 const RtempvdTable& tempvdTable = tempvdTables.getTable<RtempvdTable>(i);
1605 tempVdTable_[i].setXYContainers(tempvdTable.getDepthColumn(), tempvdTable.getTemperatureColumn());
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)
1628 const int numEquilReg = rsFunc_.size();
1629 saltVdTable_.resize(numEquilReg);
1630 const auto& tables = eclState.getTableManager();
1631 const TableContainer& saltvdTables = tables.getSaltvdTables();
1634 if (saltvdTables.empty()) {
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) {
1642 const SaltvdTable& saltvdTable = saltvdTables.getTable<SaltvdTable>(i);
1643 saltVdTable_[i].setXYContainers(saltvdTable.getDepthColumn(), saltvdTable.getSaltColumn());
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)
1667 const int numEquilReg = rsFunc_.size();
1668 saltpVdTable_.resize(numEquilReg);
1669 const auto& tables = eclState.getTableManager();
1670 const TableContainer& saltpvdTables = tables.getSaltpvdTables();
1672 for (std::size_t i = 0; i < saltpvdTables.size(); ++i) {
1673 const SaltpvdTable& saltpvdTable = saltpvdTables.getTable<SaltpvdTable>(i);
1674 saltpVdTable_[i].setXYContainers(saltpvdTable.getDepthColumn(), saltpvdTable.getSaltpColumn());
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,
1696 const NumericalAquifers& aquifer)
1698 ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
1699 int numElements = gridView.size(0);
1700 cellCenterDepth_.resize(numElements);
1701 cellZSpan_.resize(numElements);
1702 cellZMinMax_.resize(numElements);
1704 auto elemIt = gridView.template begin<0>();
1705 const auto& elemEndIt = gridView.template end<0>();
1706 const auto num_aqu_cells = aquifer.allAquiferCells();
1707 for (; elemIt != elemEndIt; ++elemIt) {
1708 const Element& element = *elemIt;
1709 const unsigned int elemIdx = elemMapper.index(element);
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);
1714 if (!num_aqu_cells.empty()) {
1715 const auto search = num_aqu_cells.find(cartIx);
1716 if (search != num_aqu_cells.end()) {
1717 const auto* aqu_cell = num_aqu_cells.at(cartIx);
1718 const Scalar depth_change_num_aqu = aqu_cell->depth - cellCenterDepth_[elemIdx];
1719 cellCenterDepth_[elemIdx] += depth_change_num_aqu;
1720 cellZSpan_[elemIdx].first += depth_change_num_aqu;
1721 cellZSpan_[elemIdx].second += depth_change_num_aqu;
1722 cellZMinMax_[elemIdx].first += depth_change_num_aqu;
1723 cellZMinMax_[elemIdx].second += depth_change_num_aqu;
1729template<
class FluidSystem,
1732 class ElementMapper,
1733 class CartesianIndexMapper>
1734void InitialStateComputer<FluidSystem,
1738 CartesianIndexMapper>::
1739applyNumericalAquifers_(
const GridView& gridView,
1740 const NumericalAquifers& aquifer,
1741 const bool co2store_or_h2store)
1743 const auto num_aqu_cells = aquifer.allAquiferCells();
1744 if (num_aqu_cells.empty())
return;
1747 bool oil_as_brine = co2store_or_h2store && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
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;
1758 for (; elemIt != elemEndIt; ++elemIt) {
1759 const Element& element = *elemIt;
1760 const unsigned int elemIdx = elemMapper.index(element);
1761 const auto cartIx = cartesianIndexMapper_.cartesianIndex(elemIdx);
1762 const auto search = num_aqu_cells.find(cartIx);
1763 if (search != num_aqu_cells.end()) {
1765 this->sat_[watPos][elemIdx] = 1.;
1767 if (!co2store_or_h2store && FluidSystem::phaseIsActive(oilPos)) {
1768 this->sat_[oilPos][elemIdx] = 0.;
1771 if (FluidSystem::phaseIsActive(gasPos)) {
1772 this->sat_[gasPos][elemIdx] = 0.;
1774 const auto* aqu_cell = num_aqu_cells.at(cartIx);
1775 const auto msg = fmt::format(
"FOR AQUIFER CELL AT ({}, {}, {}) OF NUMERICAL "
1776 "AQUIFER {}, WATER SATURATION IS SET TO BE UNITY",
1777 aqu_cell->I+1, aqu_cell->J+1, aqu_cell->K+1, aqu_cell->aquifer_id);
1782 if (aqu_cell->init_pressure) {
1783 const Scalar pres = *(aqu_cell->init_pressure);
1784 this->pp_[watPos][elemIdx] = pres;
1785 if (FluidSystem::phaseIsActive(gasPos)) {
1786 this->pp_[gasPos][elemIdx] = pres;
1788 if (FluidSystem::phaseIsActive(oilPos)) {
1789 this->pp_[oilPos][elemIdx] = pres;
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");
1811 for (
const auto& r : reg.activeRegions()) {
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<
1835 MaterialLawManager, FluidSystem, EquilReg<Scalar>,
typename RMap::CellId
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>{};
1842 std::vector<int> regionIsEmpty(rec.size(), 0);
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()) {
1857 regionIsEmpty[r] = 1;
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]
1867 vspan[0] = std::min(vspan[0], std::min(eqreg.zgoc(), eqreg.zwoc()));
1868 vspan[1] = std::max(vspan[1], std::max(eqreg.zgoc(), eqreg.zwoc()));
1870 ptable.equilibrate(eqreg, vspan);
1874 this->equilibrateCellCentres(cells, eqreg, ptable, psat);
1878 this->equilibrateHorizontal(cells, eqreg, -acc,
1887 comm.min(regionIsEmpty.data(),regionIsEmpty.size());
1888 if (comm.rank() == 0) {
1889 for (std::size_t r = 0; r < rec.size(); ++r) {
1890 if (regionIsEmpty[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>::
1908cellLoop(
const CellRange& cells,
1909 EquilibrationMethod&& eqmethod)
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);
1917 const auto watActive = FluidSystem::phaseIsActive(watPos);
1919 auto pressures = Details::PhaseQuantityValue<Scalar>{};
1920 auto saturations = Details::PhaseQuantityValue<Scalar>{};
1925 for (
const auto& cell : cells) {
1926 eqmethod(cell, pressures, saturations, Rs, Rv, Rvw);
1929 this->pp_ [oilPos][cell] = pressures.oil;
1930 this->sat_[oilPos][cell] = saturations.oil;
1934 this->pp_ [gasPos][cell] = pressures.gas;
1935 this->sat_[gasPos][cell] = saturations.gas;
1939 this->pp_ [watPos][cell] = pressures.water;
1940 this->sat_[watPos][cell] = saturations.water;
1943 if (oilActive && gasActive) {
1944 this->rs_[cell] = Rs;
1945 this->rv_[cell] = Rv;
1948 if (watActive && gasActive) {
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,
1966 const EquilReg<Scalar>& eqreg,
1967 const PressTable& ptable,
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,
1976 Details::PhaseQuantityValue<Scalar>& saturations,
1979 Scalar& Rvw) ->
void
1981 const auto pos = CellPos {
1982 cell, cellCenterDepth_[cell]
1985 saturations = psat.deriveSaturations(pos, eqreg, ptable);
1986 pressures = psat.correctedPhasePressures();
1988 const auto temp = this->temperature_[cell];
1990 Rs = eqreg.dissolutionCalculator()
1991 (pos.depth, pressures.oil, temp, saturations.gas);
1993 Rv = eqreg.evaporationCalculator()
1994 (pos.depth, pressures.gas, temp, saturations.oil);
1996 Rvw = eqreg.waterEvaporationCalculator()
1997 (pos.depth, pressures.gas, temp, saturations.water);
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,
2013 const EquilReg<Scalar>& eqreg,
2015 const PressTable& ptable,
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,
2025 Details::PhaseQuantityValue<Scalar>& saturations,
2028 Scalar& Rvw) ->
void
2031 saturations.reset();
2033 Scalar totfrac = 0.0;
2034 for (
const auto& [depth, frac] : Details::horizontalSubdivision(cell, cellZSpan_[cell], acc)) {
2035 const auto pos = CellPos { cell, depth };
2037 saturations.axpy(psat.deriveSaturations(pos, eqreg, ptable), frac);
2038 pressures .axpy(psat.correctedPhasePressures(), frac);
2044 saturations /= totfrac;
2045 pressures /= totfrac;
2048 const auto pos = CellPos {
2049 cell, cellCenterDepth_[cell]
2052 saturations = psat.deriveSaturations(pos, eqreg, ptable);
2053 pressures = psat.correctedPhasePressures();
2056 const auto temp = this->temperature_[cell];
2057 const auto cz = cellCenterDepth_[cell];
2059 Rs = eqreg.dissolutionCalculator()
2060 (cz, pressures.oil, temp, saturations.gas);
2062 Rv = eqreg.evaporationCalculator()
2063 (cz, pressures.gas, temp, saturations.oil);
2065 Rvw = eqreg.waterEvaporationCalculator()
2066 (cz, pressures.gas, temp, saturations.water);
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
FluidSystem::Scalar satFromDepth(const MaterialLawManager &materialLawManager, const typename FluidSystem::Scalar cellDepth, const typename FluidSystem::Scalar contactDepth, const int phase, const int cell, const bool increasing=false)
Compute saturation from depth. Used for constant capillary pressure function.
Definition EquilibrationHelpers_impl.hpp:650
FluidSystem::Scalar satFromPc(const MaterialLawManager &materialLawManager, const int phase, const int cell, const typename FluidSystem::Scalar targetPc, const bool increasing=false)
Compute saturation of some phase corresponding to a given capillary pressure.
Definition EquilibrationHelpers_impl.hpp:586
bool isConstPc(const MaterialLawManager &materialLawManager, const int phase, const int cell)
Return true if capillary pressure function is constant.
Definition EquilibrationHelpers_impl.hpp:670
FluidSystem::Scalar satFromSumOfPcs(const MaterialLawManager &materialLawManager, const int phase1, const int phase2, const int cell, const typename FluidSystem::Scalar targetPc)
Compute saturation of some phase corresponding to a given capillary pressure, where the capillary pre...
Definition EquilibrationHelpers_impl.hpp:619
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
Simple set of per-phase (named by primary component) quantities.
Definition InitStateEquil.hpp:338
Evaluation point within a model geometry.
Definition InitStateEquil.hpp:392