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,
typename Comm>
64void verticalExtent(
const CellRange& cells,
65 const std::vector<std::pair<double, double>>& cellZMinMax,
67 std::array<double,2>& span)
69 span[0] = std::numeric_limits<double>::max();
70 span[1] = std::numeric_limits<double>::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]);
88void subdivisionCentrePoints(
const double left,
90 const int numIntervals,
91 std::vector<std::pair<double, double>>& subdiv)
93 const auto h = (right - left) / numIntervals;
96 for (
auto i = 0*numIntervals; i < numIntervals; ++i) {
97 const auto start = end;
98 end = left + (i + 1)*h;
100 subdiv.emplace_back((start + end) / 2, h);
104template <
typename CellID>
105std::vector<std::pair<double, double>>
106horizontalSubdivision(
const CellID cell,
107 const std::pair<double, double> topbot,
108 const int numIntervals)
110 auto subdiv = std::vector<std::pair<double, double>>{};
111 subdiv.reserve(2 * numIntervals);
113 if (topbot.first > topbot.second) {
114 throw std::out_of_range {
115 "Negative thickness (inverted top/bottom faces) in cell "
116 + std::to_string(cell)
120 subdivisionCentrePoints(topbot.first, topbot.second,
121 2*numIntervals, subdiv);
126template <
class Element>
127double cellCenterDepth(
const Element& element)
129 typedef typename Element::Geometry Geometry;
130 static constexpr int zCoord = Element::dimension - 1;
133 const Geometry& geometry = element.geometry();
134 const int corners = geometry.corners();
135 for (
int i=0; i < corners; ++i)
136 zz += geometry.corner(i)[zCoord];
141template <
class Element>
142std::pair<double,double> cellZSpan(
const Element& element)
144 typedef typename Element::Geometry Geometry;
145 static constexpr int zCoord = Element::dimension - 1;
149 const Geometry& geometry = element.geometry();
150 const int corners = geometry.corners();
151 assert(corners == 8);
152 for (
int i=0; i < 4; ++i)
153 bot += geometry.corner(i)[zCoord];
154 for (
int i=4; i < corners; ++i)
155 top += geometry.corner(i)[zCoord];
157 return std::make_pair(bot/4, top/4);
160template <
class Element>
161std::pair<double,double> cellZMinMax(
const Element& element)
163 typedef typename Element::Geometry Geometry;
164 static constexpr int zCoord = Element::dimension - 1;
165 const Geometry& geometry = element.geometry();
166 const int corners = geometry.corners();
167 assert(corners == 8);
168 auto min = std::numeric_limits<double>::max();
169 auto max = std::numeric_limits<double>::lowest();
172 for (
int i=0; i < corners; ++i) {
173 min = std::min(min, geometry.corner(i)[zCoord]);
174 max = std::max(max, geometry.corner(i)[zCoord]);
176 return std::make_pair(min, max);
180RK4IVP<RHS>::RK4IVP(
const RHS& f,
181 const std::array<double,2>& span,
187 const double h = stepsize();
188 const double h2 = h / 2;
189 const double h6 = h / 6;
195 f_.push_back(f(span_[0], y0));
197 for (
int i = 0; i < N; ++i) {
198 const double x = span_[0] + i*h;
199 const double y = y_.back();
201 const double k1 = f_[i];
202 const double k2 = f(x + h2, y + h2*k1);
203 const double k3 = f(x + h2, y + h2*k2);
204 const double k4 = f(x + h, y + h*k3);
206 y_.push_back(y + h6*(k1 + 2*(k2 + k3) + k4));
207 f_.push_back(f(x + h, y_.back()));
210 assert (y_.size() == std::vector<double>::size_type(N + 1));
215operator()(
const double x)
const
219 const double h = stepsize();
220 int i = (x - span_[0]) / h;
221 const double t = (x - (span_[0] + i*h)) / h;
224 if (i < 0) { i = 0; }
225 if (N_ <= i) { i = N_ - 1; }
227 const double y0 = y_[i], y1 = y_[i + 1];
228 const double f0 = f_[i], f1 = f_[i + 1];
230 double u = (1 - 2*t) * (y1 - y0);
231 u += h * ((t - 1)*f0 + t*f1);
233 u += (1 - t)*y0 + t*y1;
242 return (span_[1] - span_[0]) / N_;
245namespace PhasePressODE {
247template<
class Flu
idSystem>
249Water(
const TabulatedFunction& tempVdTable,
250 const TabulatedFunction& saltVdTable,
251 const int pvtRegionIdx,
252 const double normGrav)
253 : tempVdTable_(tempVdTable)
254 , saltVdTable_(saltVdTable)
255 , pvtRegionIdx_(pvtRegionIdx)
260template<
class Flu
idSystem>
261double Water<FluidSystem>::
262operator()(
const double depth,
263 const double press)
const
265 return this->density(depth, press) * g_;
268template<
class Flu
idSystem>
269double Water<FluidSystem>::
270density(
const double depth,
271 const double press)
const
274 double saltConcentration = saltVdTable_.eval(depth,
true);
275 double temp = tempVdTable_.eval(depth,
true);
276 double rho = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, press, 0.0 , saltConcentration);
277 rho *= FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_);
281template<
class Flu
idSystem,
class RS>
283Oil(
const TabulatedFunction& tempVdTable,
285 const int pvtRegionIdx,
286 const double normGrav)
287 : tempVdTable_(tempVdTable)
289 , pvtRegionIdx_(pvtRegionIdx)
294template<
class Flu
idSystem,
class RS>
295double Oil<FluidSystem,RS>::
296operator()(
const double depth,
297 const double press)
const
299 return this->density(depth, press) * g_;
302template<
class Flu
idSystem,
class RS>
303double Oil<FluidSystem,RS>::
304density(
const double depth,
305 const double press)
const
307 const double temp = tempVdTable_.eval(depth,
true);
309 if(FluidSystem::enableDissolvedGas())
310 rs = rs_(depth, press, temp);
313 if (rs >= FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press)) {
314 bOil = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp, press);
317 bOil = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, press, rs);
319 double rho = bOil * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_);
320 if (FluidSystem::enableDissolvedGas()) {
321 rho += rs * bOil * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
327template<
class Flu
idSystem,
class RV,
class RVW>
328Gas<FluidSystem,RV,RVW>::
329Gas(
const TabulatedFunction& tempVdTable,
332 const int pvtRegionIdx,
333 const double normGrav)
334 : tempVdTable_(tempVdTable)
337 , pvtRegionIdx_(pvtRegionIdx)
342template<
class Flu
idSystem,
class RV,
class RVW>
343double Gas<FluidSystem,RV,RVW>::
344operator()(
const double depth,
345 const double press)
const
347 return this->density(depth, press) * g_;
350template<
class Flu
idSystem,
class RV,
class RVW>
351double Gas<FluidSystem,RV,RVW>::
352density(
const double depth,
353 const double press)
const
355 const double temp = tempVdTable_.eval(depth,
true);
357 if (FluidSystem::enableVaporizedOil())
358 rv = rv_(depth, press, temp);
361 if (FluidSystem::enableVaporizedWater())
362 rvw = rvw_(depth, press, temp);
366 if (FluidSystem::enableVaporizedOil() && FluidSystem::enableVaporizedWater()) {
367 if (rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press)
368 && rvw >= FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp, press))
370 bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp, press);
372 bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, press, rv, rvw);
374 double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
375 rho += rv * bGas * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_)
376 + rvw * bGas * FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_);
380 if (FluidSystem::enableVaporizedOil()){
381 if (rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press)) {
382 bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp, press);
384 bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, press, rv, 0.0);
386 double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
387 rho += rv * bGas * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_);
391 if (FluidSystem::enableVaporizedWater()){
392 if (rvw >= FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp, press)) {
393 bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp, press);
396 bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, press, 0.0, rvw);
398 double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
399 rho += rvw * bGas * FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_);
404 bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, press, 0.0, 0.0);
405 double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
412template<
class Flu
idSystem,
class Region>
414PressureTable<FluidSystem,Region>::
415PressureFunction<ODE>::PressureFunction(
const ODE& ode,
421 this->value_[Direction::Up] = std::make_unique<Distribution>
422 (ode, VSpan {{ ic.depth, span[0] }}, ic.pressure, nsample);
424 this->value_[Direction::Down] = std::make_unique<Distribution>
425 (ode, VSpan {{ ic.depth, span[1] }}, ic.pressure, nsample);
428template<
class Flu
idSystem,
class Region>
430PressureTable<FluidSystem,Region>::
431PressureFunction<ODE>::PressureFunction(
const PressureFunction& rhs)
432 : initial_(rhs.initial_)
434 this->value_[Direction::Up] =
435 std::make_unique<Distribution>(*rhs.value_[Direction::Up]);
437 this->value_[Direction::Down] =
438 std::make_unique<Distribution>(*rhs.value_[Direction::Down]);
441template<
class Flu
idSystem,
class Region>
443typename PressureTable<FluidSystem,Region>::template PressureFunction<ODE>&
444PressureTable<FluidSystem,Region>::
445PressureFunction<ODE>::
446operator=(
const PressureFunction& rhs)
448 this->initial_ = rhs.initial_;
450 this->value_[Direction::Up] =
451 std::make_unique<Distribution>(*rhs.value_[Direction::Up]);
453 this->value_[Direction::Down] =
454 std::make_unique<Distribution>(*rhs.value_[Direction::Down]);
459template<
class Flu
idSystem,
class Region>
461typename PressureTable<FluidSystem,Region>::template PressureFunction<ODE>&
462PressureTable<FluidSystem,Region>::
463PressureFunction<ODE>::
464operator=(PressureFunction&& rhs)
466 this->initial_ = rhs.initial_;
467 this->value_ = std::move(rhs.value_);
472template<
class Flu
idSystem,
class Region>
475PressureTable<FluidSystem,Region>::
476PressureFunction<ODE>::
477value(
const double depth)
const
479 if (depth < this->initial_.depth) {
481 return (*this->value_[Direction::Up])(depth);
483 else if (depth > this->initial_.depth) {
485 return (*this->value_[Direction::Down])(depth);
489 return this->initial_.pressure;
494template<
class Flu
idSystem,
class Region>
495template<
typename PressFunc>
496void PressureTable<FluidSystem,Region>::
497checkPtr(
const PressFunc* phasePress,
498 const std::string& phaseName)
const
500 if (phasePress !=
nullptr) {
return; }
502 throw std::invalid_argument {
503 "Phase pressure function for \"" + phaseName
504 +
"\" most not be null"
508template<
class Flu
idSystem,
class Region>
509typename PressureTable<FluidSystem,Region>::Strategy
510PressureTable<FluidSystem,Region>::
511selectEquilibrationStrategy(
const Region& reg)
const
513 if (!this->oilActive()) {
514 if (reg.datum() > reg.zwoc()) {
515 return &PressureTable::equil_WOG;
517 return &PressureTable::equil_GOW;
520 if (reg.datum() > reg.zwoc()) {
521 return &PressureTable::equil_WOG;
523 else if (reg.datum() < reg.zgoc()) {
524 return &PressureTable::equil_GOW;
527 return &PressureTable::equil_OWG;
531template<
class Flu
idSystem,
class Region>
532void PressureTable<FluidSystem,Region>::
533copyInPointers(
const PressureTable& rhs)
535 if (rhs.oil_ !=
nullptr) {
536 this->oil_ = std::make_unique<OPress>(*rhs.oil_);
539 if (rhs.gas_ !=
nullptr) {
540 this->gas_ = std::make_unique<GPress>(*rhs.gas_);
543 if (rhs.wat_ !=
nullptr) {
544 this->wat_ = std::make_unique<WPress>(*rhs.wat_);
548template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
551 const std::vector<double>& swatInit)
552 : matLawMgr_(matLawMgr)
553 , swatInit_ (swatInit)
557template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
560 : matLawMgr_(rhs.matLawMgr_)
561 , swatInit_ (rhs.swatInit_)
563 , press_ (rhs.press_)
566 this->setEvaluationPoint(*rhs.evalPt_.position,
568 *rhs.evalPt_.ptable);
571template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
578 this->setEvaluationPoint(x, reg, ptable);
579 this->initializePhaseQuantities();
581 if (ptable.
gasActive()) { this->deriveGasSat(); }
583 if (ptable.
waterActive()) { this->deriveWaterSat(); }
586 if (this->isOverlappingTransition()) {
587 this->fixUnphysicalTransition();
590 if (ptable.
oilActive()) { this->deriveOilSat(); }
592 this->accountForScaledSaturations();
597template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
601 const PTable& ptable)
603 this->evalPt_.position = &x;
604 this->evalPt_.region = ®
605 this->evalPt_.ptable = &ptable;
608template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
609void PhaseSaturations<MaterialLawManager,FluidSystem,Region,CellID>::
610initializePhaseQuantities()
613 this->press_.reset();
615 const auto depth = this->evalPt_.position->depth;
616 const auto& ptable = *this->evalPt_.ptable;
618 if (ptable.oilActive()) {
619 this->press_.oil = ptable.oil(depth);
622 if (ptable.gasActive()) {
623 this->press_.gas = ptable.gas(depth);
626 if (ptable.waterActive()) {
627 this->press_.water = ptable.water(depth);
631template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
632void PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::deriveOilSat()
634 this->sat_.oil = 1.0 - this->sat_.water - this->sat_.gas;
637template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
638void PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::deriveGasSat()
640 auto& sg = this->sat_.gas;
642 const auto isIncr =
true;
643 const auto oilActive = this->evalPt_.ptable->oilActive();
645 if (this->isConstCapPress(this->gasPos())) {
649 const auto gas_contact = oilActive? this->evalPt_.region->zgoc() : this->evalPt_.region->zwoc();
650 sg = this->fromDepthTable(gas_contact,
651 this->gasPos(), isIncr);
661 const auto pw = oilActive? this->press_.oil : this->press_.water;
662 const auto pcgo = this->press_.gas - pw;
663 sg = this->invertCapPress(pcgo, this->gasPos(), isIncr);
667template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
668void PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::deriveWaterSat()
670 auto& sw = this->sat_.water;
672 const auto oilActive = this->evalPt_.ptable->oilActive();
675 sw = 1.0 - this->sat_.gas;
678 const auto isIncr =
false;
680 if (this->isConstCapPress(this->waterPos())) {
684 sw = this->fromDepthTable(this->evalPt_.region->zwoc(),
685 this->waterPos(), isIncr);
697 const auto pcow = this->press_.oil - this->press_.water;
699 if (this->swatInit_.empty()) {
700 sw = this->invertCapPress(pcow, this->waterPos(), isIncr);
703 auto [swout, newSwatInit] = this->applySwatInit(pcow);
705 sw = this->invertCapPress(pcow, this->waterPos(), isIncr);
714template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
715void PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
716fixUnphysicalTransition()
718 auto& sg = this->sat_.gas;
719 auto& sw = this->sat_.water;
727 const auto pcgw = this->press_.gas - this->press_.water;
728 if (! this->swatInit_.empty()) {
732 auto [swout, newSwatInit] = this->applySwatInit(pcgw, sw);
734 const auto isIncr =
false;
735 sw = this->invertCapPress(pcgw, this->waterPos(), isIncr);
742 sw = satFromSumOfPcs<FluidSystem>
743 (this->matLawMgr_, this->waterPos(), this->gasPos(),
744 this->evalPt_.position->cell, pcgw);
747 this->fluidState_.setSaturation(this->oilPos(), 1.0 - sw - sg);
748 this->fluidState_.setSaturation(this->gasPos(), sg);
749 this->fluidState_.setSaturation(this->waterPos(), this->evalPt_
750 .ptable->waterActive() ? sw : 0.0);
753 this->computeMaterialLawCapPress();
754 this->press_.oil = this->press_.gas - this->materialLawCapPressGasOil();
757template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
758void PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
759accountForScaledSaturations()
761 const auto gasActive = this->evalPt_.ptable->gasActive();
762 const auto watActive = this->evalPt_.ptable->waterActive();
763 const auto oilActive = this->evalPt_.ptable->oilActive();
765 auto sg = gasActive? this->sat_.gas : 0.0;
766 auto sw = watActive? this->sat_.water : 0.0;
767 auto so = oilActive? this->sat_.oil : 0.0;
769 this->fluidState_.setSaturation(this->waterPos(), sw);
770 this->fluidState_.setSaturation(this->oilPos(), so);
771 this->fluidState_.setSaturation(this->gasPos(), sg);
773 const auto& scaledDrainageInfo = this->matLawMgr_
774 .oilWaterScaledEpsInfoDrainage(this->evalPt_.position->cell);
776 const auto thresholdSat = 1.0e-6;
777 if (watActive && ((sw + thresholdSat) > scaledDrainageInfo.Swu)) {
781 this->fluidState_.setSaturation(this->waterPos(), scaledDrainageInfo.Swu);
783 this->fluidState_.setSaturation(this->oilPos(), so + sw - scaledDrainageInfo.Swu);
784 }
else if (gasActive) {
785 this->fluidState_.setSaturation(this->gasPos(), sg + sw - scaledDrainageInfo.Swu);
787 sw = scaledDrainageInfo.Swu;
788 this->computeMaterialLawCapPress();
792 this->press_.oil = this->press_.water + this->materialLawCapPressOilWater();
795 this->press_.gas = this->press_.water + this->materialLawCapPressGasWater();
799 if (gasActive && ((sg + thresholdSat) > scaledDrainageInfo.Sgu)) {
803 this->fluidState_.setSaturation(this->gasPos(), scaledDrainageInfo.Sgu);
805 this->fluidState_.setSaturation(this->oilPos(), so + sg - scaledDrainageInfo.Sgu);
806 }
else if (watActive) {
807 this->fluidState_.setSaturation(this->waterPos(), sw + sg - scaledDrainageInfo.Sgu);
809 sg = scaledDrainageInfo.Sgu;
810 this->computeMaterialLawCapPress();
814 this->press_.oil = this->press_.gas - this->materialLawCapPressGasOil();
817 this->press_.water = this->press_.gas - this->materialLawCapPressGasWater();
821 if (watActive && ((sw - thresholdSat) < scaledDrainageInfo.Swl)) {
825 this->fluidState_.setSaturation(this->waterPos(), scaledDrainageInfo.Swl);
827 this->fluidState_.setSaturation(this->oilPos(), so + sw - scaledDrainageInfo.Swl);
828 }
else if (gasActive) {
829 this->fluidState_.setSaturation(this->gasPos(), sg + sw - scaledDrainageInfo.Swl);
831 sw = scaledDrainageInfo.Swl;
832 this->computeMaterialLawCapPress();
836 this->press_.water = this->press_.oil - this->materialLawCapPressOilWater();
839 this->press_.water = this->press_.gas - this->materialLawCapPressGasWater();
843 if (gasActive && ((sg - thresholdSat) < scaledDrainageInfo.Sgl)) {
847 this->fluidState_.setSaturation(this->gasPos(), scaledDrainageInfo.Sgl);
849 this->fluidState_.setSaturation(this->oilPos(), so + sg - scaledDrainageInfo.Sgl);
850 }
else if (watActive) {
851 this->fluidState_.setSaturation(this->waterPos(), sw + sg - scaledDrainageInfo.Sgl);
853 sg = scaledDrainageInfo.Sgl;
854 this->computeMaterialLawCapPress();
858 this->press_.gas = this->press_.oil + this->materialLawCapPressGasOil();
861 this->press_.gas = this->press_.water + this->materialLawCapPressGasWater();
866template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
867std::pair<double, bool>
868PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
869applySwatInit(
const double pcow)
871 return this->applySwatInit(pcow, this->swatInit_[this->evalPt_.position->cell]);
874template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
875std::pair<double, bool>
876PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
877applySwatInit(
const double pcow,
const double sw)
879 return this->matLawMgr_.applySwatinit(this->evalPt_.position->cell, pcow, sw);
882template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
883void PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
884computeMaterialLawCapPress()
886 const auto& matParams = this->matLawMgr_
887 .materialLawParams(this->evalPt_.position->cell);
889 this->matLawCapPress_.fill(0.0);
890 MaterialLaw::capillaryPressures(this->matLawCapPress_,
891 matParams, this->fluidState_);
894template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
895double PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
896materialLawCapPressGasOil()
const
898 return this->matLawCapPress_[this->oilPos()]
899 + this->matLawCapPress_[this->gasPos()];
902template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
903double PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
904materialLawCapPressOilWater()
const
906 return this->matLawCapPress_[this->oilPos()]
907 - this->matLawCapPress_[this->waterPos()];
910template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
911double PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
912materialLawCapPressGasWater()
const
914 return this->matLawCapPress_[this->gasPos()]
915 - this->matLawCapPress_[this->waterPos()];
918template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
919bool PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
920isConstCapPress(
const PhaseIdx phaseIdx)
const
922 return isConstPc<FluidSystem>
923 (this->matLawMgr_, phaseIdx, this->evalPt_.position->cell);
926template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
927bool PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
928isOverlappingTransition()
const
930 return this->evalPt_.ptable->gasActive()
931 && this->evalPt_.ptable->waterActive()
932 && ((this->sat_.gas + this->sat_.water) > 1.0);
935template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
936double PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
937fromDepthTable(
const double contactdepth,
938 const PhaseIdx phasePos,
939 const bool isincr)
const
941 return satFromDepth<FluidSystem>
942 (this->matLawMgr_, this->evalPt_.position->depth,
943 contactdepth,
static_cast<int>(phasePos),
944 this->evalPt_.position->cell, isincr);
947template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
948double PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
949invertCapPress(
const double pc,
950 const PhaseIdx phasePos,
951 const bool isincr)
const
953 return satFromPc<FluidSystem>
954 (this->matLawMgr_,
static_cast<int>(phasePos),
955 this->evalPt_.position->cell, pc, isincr);
958template<
class Flu
idSystem,
class Region>
961 const int samplePoints)
963 , nsample_(samplePoints)
967template <
class Flu
idSystem,
class Region>
970 : gravity_(rhs.gravity_)
971 , nsample_(rhs.nsample_)
973 this->copyInPointers(rhs);
976template <
class Flu
idSystem,
class Region>
979 : gravity_(rhs.gravity_)
980 , nsample_(rhs.nsample_)
981 , oil_ (std::move(rhs.oil_))
982 , gas_ (std::move(rhs.gas_))
983 , wat_ (std::move(rhs.wat_))
987template <
class Flu
idSystem,
class Region>
992 this->gravity_ = rhs.gravity_;
993 this->nsample_ = rhs.nsample_;
994 this->copyInPointers(rhs);
999template <
class Flu
idSystem,
class Region>
1004 this->gravity_ = rhs.gravity_;
1005 this->nsample_ = rhs.nsample_;
1007 this->oil_ = std::move(rhs.oil_);
1008 this->gas_ = std::move(rhs.gas_);
1009 this->wat_ = std::move(rhs.wat_);
1014template <
class Flu
idSystem,
class Region>
1020 auto equil = this->selectEquilibrationStrategy(reg);
1022 (this->*equil)(reg, span);
1025template <
class Flu
idSystem,
class Region>
1029 return FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
1032template <
class Flu
idSystem,
class Region>
1036 return FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
1039template <
class Flu
idSystem,
class Region>
1043 return FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx);
1046template <
class Flu
idSystem,
class Region>
1048oil(
const double depth)
const
1050 this->checkPtr(this->oil_.get(),
"OIL");
1052 return this->oil_->value(depth);
1055template <
class Flu
idSystem,
class Region>
1057gas(
const double depth)
const
1059 this->checkPtr(this->gas_.get(),
"GAS");
1061 return this->gas_->value(depth);
1065template <
class Flu
idSystem,
class Region>
1067water(
const double depth)
const
1069 this->checkPtr(this->wat_.get(),
"WATER");
1071 return this->wat_->value(depth);
1074template <
class Flu
idSystem,
class Region>
1076equil_WOG(
const Region& reg,
const VSpan& span)
1081 if (! this->waterActive()) {
1082 throw std::invalid_argument {
1083 "Don't know how to interpret EQUIL datum depth in "
1084 "WATER zone in model without active water phase"
1089 const auto ic =
typename WPress::InitCond {
1090 reg.datum(), reg.pressure()
1093 this->makeWatPressure(ic, reg, span);
1096 if (this->oilActive()) {
1098 const auto ic =
typename OPress::InitCond {
1100 this->
water(reg.zwoc()) + reg.pcowWoc()
1103 this->makeOilPressure(ic, reg, span);
1106 if (this->gasActive() && this->oilActive()) {
1108 const auto ic =
typename GPress::InitCond {
1110 this->
oil(reg.zgoc()) + reg.pcgoGoc()
1113 this->makeGasPressure(ic, reg, span);
1114 }
else if (this->gasActive() && !this->oilActive()) {
1116 const auto ic =
typename GPress::InitCond {
1118 this->
water(reg.zwoc()) + reg.pcowWoc()
1120 this->makeGasPressure(ic, reg, span);
1124template <
class Flu
idSystem,
class Region>
1125void PressureTable<FluidSystem, Region>::
1126equil_GOW(
const Region& reg,
const VSpan& span)
1131 if (! this->gasActive()) {
1132 throw std::invalid_argument {
1133 "Don't know how to interpret EQUIL datum depth in "
1134 "GAS zone in model without active gas phase"
1139 const auto ic =
typename GPress::InitCond {
1140 reg.datum(), reg.pressure()
1143 this->makeGasPressure(ic, reg, span);
1146 if (this->oilActive()) {
1148 const auto ic =
typename OPress::InitCond {
1150 this->
gas(reg.zgoc()) - reg.pcgoGoc()
1152 this->makeOilPressure(ic, reg, span);
1155 if (this->waterActive() && this->oilActive()) {
1157 const auto ic =
typename WPress::InitCond {
1159 this->
oil(reg.zwoc()) - reg.pcowWoc()
1162 this->makeWatPressure(ic, reg, span);
1163 }
else if (this->waterActive() && !this->oilActive()) {
1165 const auto ic =
typename WPress::InitCond {
1167 this->
gas(reg.zwoc()) - reg.pcowWoc()
1169 this->makeWatPressure(ic, reg, span);
1173template <
class Flu
idSystem,
class Region>
1174void PressureTable<FluidSystem, Region>::
1175equil_OWG(
const Region& reg,
const VSpan& span)
1180 if (! this->oilActive()) {
1181 throw std::invalid_argument {
1182 "Don't know how to interpret EQUIL datum depth in "
1183 "OIL zone in model without active oil phase"
1188 const auto ic =
typename OPress::InitCond {
1189 reg.datum(), reg.pressure()
1192 this->makeOilPressure(ic, reg, span);
1195 if (this->waterActive()) {
1197 const auto ic =
typename WPress::InitCond {
1199 this->
oil(reg.zwoc()) - reg.pcowWoc()
1202 this->makeWatPressure(ic, reg, span);
1205 if (this->gasActive()) {
1207 const auto ic =
typename GPress::InitCond {
1209 this->
oil(reg.zgoc()) + reg.pcgoGoc()
1211 this->makeGasPressure(ic, reg, span);
1215template <
class Flu
idSystem,
class Region>
1216void PressureTable<FluidSystem, Region>::
1217makeOilPressure(
const typename OPress::InitCond& ic,
1221 const auto drho = OilPressODE {
1222 reg.tempVdTable(), reg.dissolutionCalculator(),
1223 reg.pvtIdx(), this->gravity_
1226 this->oil_ = std::make_unique<OPress>(drho, ic, this->nsample_, span);
1229template <
class Flu
idSystem,
class Region>
1230void PressureTable<FluidSystem, Region>::
1231makeGasPressure(
const typename GPress::InitCond& ic,
1235 const auto drho = GasPressODE {
1236 reg.tempVdTable(), reg.evaporationCalculator(), reg.waterEvaporationCalculator(),
1237 reg.pvtIdx(), this->gravity_
1240 this->gas_ = std::make_unique<GPress>(drho, ic, this->nsample_, span);
1243template <
class Flu
idSystem,
class Region>
1244void PressureTable<FluidSystem, Region>::
1245makeWatPressure(
const typename WPress::InitCond& ic,
1249 const auto drho = WatPressODE {
1250 reg.tempVdTable(), reg.saltVdTable(), reg.pvtIdx(), this->gravity_
1253 this->wat_ = std::make_unique<WPress>(drho, ic, this->nsample_, span);
1258namespace DeckDependent {
1260std::vector<EquilRecord>
1261getEquil(
const EclipseState& state)
1263 const auto& init = state.getInitConfig();
1265 if(!init.hasEquil()) {
1266 throw std::domain_error(
"Deck does not provide equilibration data.");
1269 const auto& equil = init.getEquil();
1270 return { equil.begin(), equil.end() };
1273template<
class Gr
idView>
1275equilnum(
const EclipseState& eclipseState,
1276 const GridView& gridview)
1278 std::vector<int> eqlnum(gridview.size(0), 0);
1280 if (eclipseState.fieldProps().has_int(
"EQLNUM")) {
1281 const auto& e = eclipseState.fieldProps().get_int(
"EQLNUM");
1282 std::transform(e.begin(), e.end(), eqlnum.begin(), [](
int n){ return n - 1;});
1284 OPM_BEGIN_PARALLEL_TRY_CATCH();
1285 const int num_regions = eclipseState.getTableManager().getEqldims().getNumEquilRegions();
1286 if ( std::any_of(eqlnum.begin(), eqlnum.end(), [num_regions](
int n){return n >= num_regions;}) ) {
1287 throw std::runtime_error(
"Values larger than maximum Equil regions " + std::to_string(num_regions) +
" provided in EQLNUM");
1289 if ( std::any_of(eqlnum.begin(), eqlnum.end(), [](
int n){return n < 0;}) ) {
1290 throw std::runtime_error(
"zero or negative values provided in EQLNUM");
1292 OPM_END_PARALLEL_TRY_CATCH(
"Invalied EQLNUM numbers: ", gridview.comm());
1297template<
class FluidSystem,
1300 class ElementMapper,
1301 class CartesianIndexMapper>
1302template<
class MaterialLawManager>
1303InitialStateComputer<FluidSystem,
1307 CartesianIndexMapper>::
1308InitialStateComputer(MaterialLawManager& materialLawManager,
1309 const EclipseState& eclipseState,
1311 const GridView& gridView,
1312 const CartesianIndexMapper& cartMapper,
1314 const int num_pressure_points,
1315 const bool applySwatInit)
1316 : temperature_(grid.size(0), eclipseState.getTableManager().rtemp()),
1317 saltConcentration_(grid.size(0)),
1318 saltSaturation_(grid.size(0)),
1319 pp_(FluidSystem::numPhases,
1320 std::vector<double>(grid.size(0))),
1321 sat_(FluidSystem::numPhases,
1322 std::vector<double>(grid.size(0))),
1326 cartesianIndexMapper_(cartMapper),
1327 num_pressure_points_(num_pressure_points)
1330 if (applySwatInit) {
1331 if (eclipseState.fieldProps().has_double(
"SWATINIT")) {
1332 swatInit_ = eclipseState.fieldProps().get_double(
"SWATINIT");
1338 const auto& num_aquifers = eclipseState.aquifer().numericalAquifers();
1339 updateCellProps_(gridView, num_aquifers);
1342 const std::vector<EquilRecord> rec = getEquil(eclipseState);
1343 const auto& tables = eclipseState.getTableManager();
1345 const RegionMapping<> eqlmap(equilnum(eclipseState, grid));
1346 const int invalidRegion = -1;
1347 regionPvtIdx_.resize(rec.size(), invalidRegion);
1348 setRegionPvtIdx(eclipseState, eqlmap);
1351 rsFunc_.reserve(rec.size());
1352 if (FluidSystem::enableDissolvedGas()) {
1353 for (std::size_t i = 0; i < rec.size(); ++i) {
1354 if (eqlmap.cells(i).empty()) {
1355 rsFunc_.push_back(std::shared_ptr<Miscibility::RsVD<FluidSystem>>());
1358 const int pvtIdx = regionPvtIdx_[i];
1359 if (!rec[i].liveOilInitConstantRs()) {
1360 const TableContainer& rsvdTables = tables.getRsvdTables();
1361 const TableContainer& pbvdTables = tables.getPbvdTables();
1362 if (rsvdTables.size() > 0) {
1364 const RsvdTable& rsvdTable = rsvdTables.getTable<RsvdTable>(i);
1365 std::vector<double> depthColumn = rsvdTable.getColumn(
"DEPTH").vectorCopy();
1366 std::vector<double> rsColumn = rsvdTable.getColumn(
"RS").vectorCopy();
1367 rsFunc_.push_back(std::make_shared<Miscibility::RsVD<FluidSystem>>(pvtIdx,
1368 depthColumn, rsColumn));
1369 }
else if (pbvdTables.size() > 0) {
1370 const PbvdTable& pbvdTable = pbvdTables.getTable<PbvdTable>(i);
1371 std::vector<double> depthColumn = pbvdTable.getColumn(
"DEPTH").vectorCopy();
1372 std::vector<double> pbubColumn = pbvdTable.getColumn(
"PBUB").vectorCopy();
1373 rsFunc_.push_back(std::make_shared<Miscibility::PBVD<FluidSystem>>(pvtIdx,
1374 depthColumn, pbubColumn));
1377 throw std::runtime_error(
"Cannot initialise: RSVD or PBVD table not available.");
1382 if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) {
1383 throw std::runtime_error(
"Cannot initialise: when no explicit RSVD table is given, \n"
1384 "datum depth must be at the gas-oil-contact. "
1385 "In EQUIL region "+std::to_string(i + 1)+
" (counting from 1), this does not hold.");
1387 const double pContact = rec[i].datumDepthPressure();
1388 const double TContact = 273.15 + 20;
1389 rsFunc_.push_back(std::make_shared<Miscibility::RsSatAtContact<FluidSystem>>(pvtIdx, pContact, TContact));
1394 for (std::size_t i = 0; i < rec.size(); ++i) {
1395 rsFunc_.push_back(std::make_shared<Miscibility::NoMixing>());
1399 rvFunc_.reserve(rec.size());
1400 if (FluidSystem::enableVaporizedOil()) {
1401 for (std::size_t i = 0; i < rec.size(); ++i) {
1402 if (eqlmap.cells(i).empty()) {
1403 rvFunc_.push_back(std::shared_ptr<Miscibility::RvVD<FluidSystem>>());
1406 const int pvtIdx = regionPvtIdx_[i];
1407 if (!rec[i].wetGasInitConstantRv()) {
1408 const TableContainer& rvvdTables = tables.getRvvdTables();
1409 const TableContainer& pdvdTables = tables.getPdvdTables();
1411 if (rvvdTables.size() > 0) {
1412 const RvvdTable& rvvdTable = rvvdTables.getTable<RvvdTable>(i);
1413 std::vector<double> depthColumn = rvvdTable.getColumn(
"DEPTH").vectorCopy();
1414 std::vector<double> rvColumn = rvvdTable.getColumn(
"RV").vectorCopy();
1415 rvFunc_.push_back(std::make_shared<Miscibility::RvVD<FluidSystem>>(pvtIdx,
1416 depthColumn, rvColumn));
1417 }
else if (pdvdTables.size() > 0) {
1418 const PdvdTable& pdvdTable = pdvdTables.getTable<PdvdTable>(i);
1419 std::vector<double> depthColumn = pdvdTable.getColumn(
"DEPTH").vectorCopy();
1420 std::vector<double> pdewColumn = pdvdTable.getColumn(
"PDEW").vectorCopy();
1421 rvFunc_.push_back(std::make_shared<Miscibility::PDVD<FluidSystem>>(pvtIdx,
1422 depthColumn, pdewColumn));
1424 throw std::runtime_error(
"Cannot initialise: RVVD or PDCD table not available.");
1428 if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) {
1429 throw std::runtime_error(
1430 "Cannot initialise: when no explicit RVVD table is given, \n"
1431 "datum depth must be at the gas-oil-contact. "
1432 "In EQUIL region "+std::to_string(i + 1)+
" (counting from 1), this does not hold.");
1434 const double pContact = rec[i].datumDepthPressure() + rec[i].gasOilContactCapillaryPressure();
1435 const double TContact = 273.15 + 20;
1436 rvFunc_.push_back(std::make_shared<Miscibility::RvSatAtContact<FluidSystem>>(pvtIdx,pContact, TContact));
1441 for (std::size_t i = 0; i < rec.size(); ++i) {
1442 rvFunc_.push_back(std::make_shared<Miscibility::NoMixing>());
1446 rvwFunc_.reserve(rec.size());
1447 if (FluidSystem::enableVaporizedWater()) {
1448 for (std::size_t i = 0; i < rec.size(); ++i) {
1449 if (eqlmap.cells(i).empty()) {
1450 rvwFunc_.push_back(std::shared_ptr<Miscibility::RvwVD<FluidSystem>>());
1453 const int pvtIdx = regionPvtIdx_[i];
1454 if (!rec[i].humidGasInitConstantRvw()) {
1455 const TableContainer& rvwvdTables = tables.getRvwvdTables();
1457 if (rvwvdTables.size() > 0) {
1458 const RvwvdTable& rvwvdTable = rvwvdTables.getTable<RvwvdTable>(i);
1459 std::vector<double> depthColumn = rvwvdTable.getColumn(
"DEPTH").vectorCopy();
1460 std::vector<double> rvwvdColumn = rvwvdTable.getColumn(
"RVWVD").vectorCopy();
1461 rvwFunc_.push_back(std::make_shared<Miscibility::RvwVD<FluidSystem>>(pvtIdx,
1462 depthColumn, rvwvdColumn));
1464 throw std::runtime_error(
"Cannot initialise: RVWVD table not available.");
1468 const auto oilActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
1470 if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) {
1471 rvwFunc_.push_back(std::make_shared<Miscibility::NoMixing>());
1472 const auto msg =
"No explicit RVWVD table is given for EQUIL region " + std::to_string(i + 1) +
". \n"
1473 "and datum depth is not at the gas-oil-contact. \n"
1474 "Rvw is set to 0.0 in all cells. \n";
1475 OpmLog::warning(msg);
1479 const double pContact = rec[i].datumDepthPressure() + rec[i].gasOilContactCapillaryPressure();
1480 const double TContact = 273.15 + 20;
1481 rvwFunc_.push_back(std::make_shared<Miscibility::RvwSatAtContact<FluidSystem>>(pvtIdx,pContact, TContact));
1487 if (rec[i].waterOilContactDepth() != rec[i].datumDepth()) {
1488 rvwFunc_.push_back(std::make_shared<Miscibility::NoMixing>());
1489 const auto msg =
"No explicit RVWVD table is given for EQUIL region " + std::to_string(i + 1) +
". \n"
1490 "and datum depth is not at the gas-water-contact. \n"
1491 "Rvw is set to 0.0 in all cells. \n";
1492 OpmLog::warning(msg);
1495 const double pContact = rec[i].datumDepthPressure() + rec[i].waterOilContactCapillaryPressure();
1496 const double TContact = 273.15 + 20;
1497 rvwFunc_.push_back(std::make_shared<Miscibility::RvwSatAtContact<FluidSystem>>(pvtIdx,pContact, TContact));
1504 for (std::size_t i = 0; i < rec.size(); ++i) {
1505 rvwFunc_.push_back(std::make_shared<Miscibility::NoMixing>());
1511 updateInitialTemperature_(eclipseState, eqlmap);
1514 updateInitialSaltConcentration_(eclipseState, eqlmap);
1517 updateInitialSaltSaturation_(eclipseState, eqlmap);
1520 const auto& comm = grid.comm();
1521 calcPressSatRsRv(eqlmap, rec, materialLawManager, comm, grav);
1524 applyNumericalAquifers_(gridView, num_aquifers, eclipseState.runspec().co2Storage() || eclipseState.runspec().h2Storage());
1530template<
class FluidSystem,
1533 class ElementMapper,
1534 class CartesianIndexMapper>
1536void InitialStateComputer<FluidSystem,
1540 CartesianIndexMapper>::
1541updateInitialTemperature_(
const EclipseState& eclState,
const RMap& reg)
1543 const int numEquilReg = rsFunc_.size();
1544 tempVdTable_.resize(numEquilReg);
1545 const auto& tables = eclState.getTableManager();
1546 if (!tables.hasTables(
"RTEMPVD")) {
1547 std::vector<double> x = {0.0,1.0};
1548 std::vector<double> y = {tables.rtemp(),tables.rtemp()};
1549 for (
auto& table : this->tempVdTable_) {
1550 table.setXYContainers(x, y);
1553 const TableContainer& tempvdTables = tables.getRtempvdTables();
1554 for (std::size_t i = 0; i < tempvdTables.size(); ++i) {
1555 const RtempvdTable& tempvdTable = tempvdTables.getTable<RtempvdTable>(i);
1556 tempVdTable_[i].setXYContainers(tempvdTable.getDepthColumn(), tempvdTable.getTemperatureColumn());
1557 const auto& cells = reg.cells(i);
1558 for (
const auto& cell : cells) {
1559 const double depth = cellCenterDepth_[cell];
1560 this->temperature_[cell] = tempVdTable_[i].eval(depth,
true);
1566template<
class FluidSystem,
1569 class ElementMapper,
1570 class CartesianIndexMapper>
1572void InitialStateComputer<FluidSystem,
1576 CartesianIndexMapper>::
1577updateInitialSaltConcentration_(
const EclipseState& eclState,
const RMap& reg)
1579 const int numEquilReg = rsFunc_.size();
1580 saltVdTable_.resize(numEquilReg);
1581 const auto& tables = eclState.getTableManager();
1582 const TableContainer& saltvdTables = tables.getSaltvdTables();
1585 if (saltvdTables.empty()) {
1586 std::vector<double> x = {0.0,1.0};
1587 std::vector<double> y = {0.0,0.0};
1588 for (
auto& table : this->saltVdTable_) {
1589 table.setXYContainers(x, y);
1592 for (std::size_t i = 0; i < saltvdTables.size(); ++i) {
1593 const SaltvdTable& saltvdTable = saltvdTables.getTable<SaltvdTable>(i);
1594 saltVdTable_[i].setXYContainers(saltvdTable.getDepthColumn(), saltvdTable.getSaltColumn());
1596 const auto& cells = reg.cells(i);
1597 for (
const auto& cell : cells) {
1598 const double depth = cellCenterDepth_[cell];
1599 this->saltConcentration_[cell] = saltVdTable_[i].eval(depth,
true);
1605template<
class FluidSystem,
1608 class ElementMapper,
1609 class CartesianIndexMapper>
1611void InitialStateComputer<FluidSystem,
1615 CartesianIndexMapper>::
1616updateInitialSaltSaturation_(
const EclipseState& eclState,
const RMap& reg)
1618 const int numEquilReg = rsFunc_.size();
1619 saltpVdTable_.resize(numEquilReg);
1620 const auto& tables = eclState.getTableManager();
1621 const TableContainer& saltpvdTables = tables.getSaltpvdTables();
1623 for (std::size_t i = 0; i < saltpvdTables.size(); ++i) {
1624 const SaltpvdTable& saltpvdTable = saltpvdTables.getTable<SaltpvdTable>(i);
1625 saltpVdTable_[i].setXYContainers(saltpvdTable.getDepthColumn(), saltpvdTable.getSaltpColumn());
1627 const auto& cells = reg.cells(i);
1628 for (
const auto& cell : cells) {
1629 const double depth = cellCenterDepth_[cell];
1630 this->saltSaturation_[cell] = saltpVdTable_[i].eval(depth,
true);
1636template<
class FluidSystem,
1639 class ElementMapper,
1640 class CartesianIndexMapper>
1641void InitialStateComputer<FluidSystem,
1645 CartesianIndexMapper>::
1646updateCellProps_(
const GridView& gridView,
1647 const NumericalAquifers& aquifer)
1649 ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
1650 int numElements = gridView.size(0);
1651 cellCenterDepth_.resize(numElements);
1652 cellZSpan_.resize(numElements);
1653 cellZMinMax_.resize(numElements);
1655 auto elemIt = gridView.template begin<0>();
1656 const auto& elemEndIt = gridView.template end<0>();
1657 const auto num_aqu_cells = aquifer.allAquiferCells();
1658 for (; elemIt != elemEndIt; ++elemIt) {
1659 const Element& element = *elemIt;
1660 const unsigned int elemIdx = elemMapper.index(element);
1661 cellCenterDepth_[elemIdx] = Details::cellCenterDepth(element);
1662 const auto cartIx = cartesianIndexMapper_.cartesianIndex(elemIdx);
1663 cellZSpan_[elemIdx] = Details::cellZSpan(element);
1664 cellZMinMax_[elemIdx] = Details::cellZMinMax(element);
1665 if (!num_aqu_cells.empty()) {
1666 const auto search = num_aqu_cells.find(cartIx);
1667 if (search != num_aqu_cells.end()) {
1668 const auto* aqu_cell = num_aqu_cells.at(cartIx);
1669 const double depth_change_num_aqu = aqu_cell->depth - cellCenterDepth_[elemIdx];
1670 cellCenterDepth_[elemIdx] += depth_change_num_aqu;
1671 cellZSpan_[elemIdx].first += depth_change_num_aqu;
1672 cellZSpan_[elemIdx].second += depth_change_num_aqu;
1673 cellZMinMax_[elemIdx].first += depth_change_num_aqu;
1674 cellZMinMax_[elemIdx].second += depth_change_num_aqu;
1680template<
class FluidSystem,
1683 class ElementMapper,
1684 class CartesianIndexMapper>
1685void InitialStateComputer<FluidSystem,
1689 CartesianIndexMapper>::
1690applyNumericalAquifers_(
const GridView& gridView,
1691 const NumericalAquifers& aquifer,
1692 const bool co2store_or_h2store)
1694 const auto num_aqu_cells = aquifer.allAquiferCells();
1695 if (num_aqu_cells.empty())
return;
1698 bool oil_as_brine = co2store_or_h2store && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
1699 const auto watPos = oil_as_brine? FluidSystem::oilPhaseIdx : FluidSystem::waterPhaseIdx;
1700 if (!FluidSystem::phaseIsActive(watPos)){
1701 throw std::logic_error {
"Water phase has to be active for numerical aquifer case" };
1704 ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
1705 auto elemIt = gridView.template begin<0>();
1706 const auto& elemEndIt = gridView.template end<0>();
1707 const auto oilPos = FluidSystem::oilPhaseIdx;
1708 const auto gasPos = FluidSystem::gasPhaseIdx;
1709 for (; elemIt != elemEndIt; ++elemIt) {
1710 const Element& element = *elemIt;
1711 const unsigned int elemIdx = elemMapper.index(element);
1712 const auto cartIx = cartesianIndexMapper_.cartesianIndex(elemIdx);
1713 const auto search = num_aqu_cells.find(cartIx);
1714 if (search != num_aqu_cells.end()) {
1716 this->sat_[watPos][elemIdx] = 1.;
1718 if (!co2store_or_h2store && FluidSystem::phaseIsActive(oilPos)) {
1719 this->sat_[oilPos][elemIdx] = 0.;
1722 if (FluidSystem::phaseIsActive(gasPos)) {
1723 this->sat_[gasPos][elemIdx] = 0.;
1725 const auto* aqu_cell = num_aqu_cells.at(cartIx);
1726 const auto msg = fmt::format(
"FOR AQUIFER CELL AT ({}, {}, {}) OF NUMERICAL "
1727 "AQUIFER {}, WATER SATURATION IS SET TO BE UNITY",
1728 aqu_cell->I+1, aqu_cell->J+1, aqu_cell->K+1, aqu_cell->aquifer_id);
1733 if (aqu_cell->init_pressure) {
1734 const double pres = *(aqu_cell->init_pressure);
1735 this->pp_[watPos][elemIdx] = pres;
1736 if (FluidSystem::phaseIsActive(gasPos)) {
1737 this->pp_[gasPos][elemIdx] = pres;
1739 if (FluidSystem::phaseIsActive(oilPos)) {
1740 this->pp_[oilPos][elemIdx] = pres;
1747template<
class FluidSystem,
1750 class ElementMapper,
1751 class CartesianIndexMapper>
1753void InitialStateComputer<FluidSystem,
1757 CartesianIndexMapper>::
1758setRegionPvtIdx(
const EclipseState& eclState,
const RMap& reg)
1760 const auto& pvtnumData = eclState.fieldProps().get_int(
"PVTNUM");
1762 for (
const auto& r : reg.activeRegions()) {
1763 const auto& cells = reg.cells(r);
1764 regionPvtIdx_[r] = pvtnumData[*cells.begin()] - 1;
1768template<
class FluidSystem,
1771 class ElementMapper,
1772 class CartesianIndexMapper>
1773template<
class RMap,
class MaterialLawManager,
class Comm>
1774void InitialStateComputer<FluidSystem,
1778 CartesianIndexMapper>::
1779calcPressSatRsRv(
const RMap& reg,
1780 const std::vector<EquilRecord>& rec,
1781 MaterialLawManager& materialLawManager,
1785 using PhaseSat = Details::PhaseSaturations<
1786 MaterialLawManager, FluidSystem, EquilReg,
typename RMap::CellId
1789 auto ptable = Details::PressureTable<FluidSystem, EquilReg>{ grav, this->num_pressure_points_ };
1790 auto psat = PhaseSat { materialLawManager, this->swatInit_ };
1791 auto vspan = std::array<double, 2>{};
1793 std::vector<int> regionIsEmpty(rec.size(), 0);
1794 for (std::size_t r = 0; r < rec.size(); ++r) {
1795 const auto& cells = reg.cells(r);
1797 Details::verticalExtent(cells, cellZMinMax_, comm, vspan);
1799 const auto acc = rec[r].initializationTargetAccuracy();
1801 throw std::runtime_error {
1802 "Cannot initialise model: Positive item 9 is not supported "
1803 "in EQUIL keyword, record " + std::to_string(r + 1)
1807 if (cells.empty()) {
1808 regionIsEmpty[r] = 1;
1812 const auto eqreg = EquilReg {
1813 rec[r], this->rsFunc_[r], this->rvFunc_[r], this->rvwFunc_[r], this->tempVdTable_[r], this->saltVdTable_[r], this->regionPvtIdx_[r]
1818 vspan[0] = std::min(vspan[0], std::min(eqreg.zgoc(), eqreg.zwoc()));
1819 vspan[1] = std::max(vspan[1], std::max(eqreg.zgoc(), eqreg.zwoc()));
1821 ptable.equilibrate(eqreg, vspan);
1825 this->equilibrateCellCentres(cells, eqreg, ptable, psat);
1829 this->equilibrateHorizontal(cells, eqreg, -acc,
1838 comm.min(regionIsEmpty.data(),regionIsEmpty.size());
1839 if (comm.rank() == 0) {
1840 for (std::size_t r = 0; r < rec.size(); ++r) {
1841 if (regionIsEmpty[r])
1842 OpmLog::warning(
"Equilibration region " + std::to_string(r + 1)
1843 +
" has no active cells");
1848template<
class FluidSystem,
1851 class ElementMapper,
1852 class CartesianIndexMapper>
1853template<
class CellRange,
class EquilibrationMethod>
1854void InitialStateComputer<FluidSystem,
1858 CartesianIndexMapper>::
1859cellLoop(
const CellRange& cells,
1860 EquilibrationMethod&& eqmethod)
1862 const auto oilPos = FluidSystem::oilPhaseIdx;
1863 const auto gasPos = FluidSystem::gasPhaseIdx;
1864 const auto watPos = FluidSystem::waterPhaseIdx;
1866 const auto oilActive = FluidSystem::phaseIsActive(oilPos);
1867 const auto gasActive = FluidSystem::phaseIsActive(gasPos);
1868 const auto watActive = FluidSystem::phaseIsActive(watPos);
1870 auto pressures = Details::PhaseQuantityValue{};
1871 auto saturations = Details::PhaseQuantityValue{};
1876 for (
const auto& cell : cells) {
1877 eqmethod(cell, pressures, saturations, Rs, Rv, Rvw);
1880 this->pp_ [oilPos][cell] = pressures.oil;
1881 this->sat_[oilPos][cell] = saturations.oil;
1885 this->pp_ [gasPos][cell] = pressures.gas;
1886 this->sat_[gasPos][cell] = saturations.gas;
1890 this->pp_ [watPos][cell] = pressures.water;
1891 this->sat_[watPos][cell] = saturations.water;
1894 if (oilActive && gasActive) {
1895 this->rs_[cell] = Rs;
1896 this->rv_[cell] = Rv;
1899 if (watActive && gasActive) {
1900 this->rvw_[cell] = Rvw;
1905template<
class FluidSystem,
1908 class ElementMapper,
1909 class CartesianIndexMapper>
1910template<
class CellRange,
class PressTable,
class PhaseSat>
1911void InitialStateComputer<FluidSystem,
1915 CartesianIndexMapper>::
1916equilibrateCellCentres(
const CellRange& cells,
1917 const EquilReg& eqreg,
1918 const PressTable& ptable,
1921 using CellPos =
typename PhaseSat::Position;
1922 using CellID = std::remove_cv_t<std::remove_reference_t<
1923 decltype(std::declval<CellPos>().cell)>>;
1924 this->cellLoop(cells, [
this, &eqreg, &ptable, &psat]
1926 Details::PhaseQuantityValue& pressures,
1927 Details::PhaseQuantityValue& saturations,
1930 double& Rvw) ->
void
1932 const auto pos = CellPos {
1933 cell, cellCenterDepth_[cell]
1936 saturations = psat.deriveSaturations(pos, eqreg, ptable);
1937 pressures = psat.correctedPhasePressures();
1939 const auto temp = this->temperature_[cell];
1941 Rs = eqreg.dissolutionCalculator()
1942 (pos.depth, pressures.oil, temp, saturations.gas);
1944 Rv = eqreg.evaporationCalculator()
1945 (pos.depth, pressures.gas, temp, saturations.oil);
1947 Rvw = eqreg.waterEvaporationCalculator()
1948 (pos.depth, pressures.gas, temp, saturations.water);
1952template<
class FluidSystem,
1955 class ElementMapper,
1956 class CartesianIndexMapper>
1957template<
class CellRange,
class PressTable,
class PhaseSat>
1958void InitialStateComputer<FluidSystem,
1962 CartesianIndexMapper>::
1963equilibrateHorizontal(
const CellRange& cells,
1964 const EquilReg& eqreg,
1966 const PressTable& ptable,
1969 using CellPos =
typename PhaseSat::Position;
1970 using CellID = std::remove_cv_t<std::remove_reference_t<
1971 decltype(std::declval<CellPos>().cell)>>;
1973 this->cellLoop(cells, [
this, acc, &eqreg, &ptable, &psat]
1975 Details::PhaseQuantityValue& pressures,
1976 Details::PhaseQuantityValue& saturations,
1979 double& Rvw) ->
void
1982 saturations.reset();
1985 for (
const auto& [depth, frac] : Details::horizontalSubdivision(cell, cellZSpan_[cell], acc)) {
1986 const auto pos = CellPos { cell, depth };
1988 saturations.axpy(psat.deriveSaturations(pos, eqreg, ptable), frac);
1989 pressures .axpy(psat.correctedPhasePressures(), frac);
1995 saturations /= totfrac;
1996 pressures /= totfrac;
1999 const auto pos = CellPos {
2000 cell, cellCenterDepth_[cell]
2003 saturations = psat.deriveSaturations(pos, eqreg, ptable);
2004 pressures = psat.correctedPhasePressures();
2007 const auto temp = this->temperature_[cell];
2008 const auto cz = cellCenterDepth_[cell];
2010 Rs = eqreg.dissolutionCalculator()
2011 (cz, pressures.oil, temp, saturations.gas);
2013 Rv = eqreg.evaporationCalculator()
2014 (cz, pressures.gas, temp, saturations.oil);
2016 Rvw = eqreg.waterEvaporationCalculator()
2017 (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:376
const PhaseQuantityValue & deriveSaturations(const Position &x, const Region ®, const PTable &ptable)
Calculate phase saturations at particular point of the simulation model geometry.
Definition InitStateEquil_impl.hpp:574
PhaseSaturations(MaterialLawManager &matLawMgr, const std::vector< double > &swatInit)
Constructor.
Definition InitStateEquil_impl.hpp:550
Definition InitStateEquil.hpp:159
PressureTable & operator=(const PressureTable &rhs)
Assignment operator.
Definition InitStateEquil_impl.hpp:990
double gas(const double depth) const
Evaluate gas phase pressure at specified depth.
Definition InitStateEquil_impl.hpp:1057
bool waterActive() const
Predicate for whether or not water is an active phase.
Definition InitStateEquil_impl.hpp:1041
bool gasActive() const
Predicate for whether or not gas is an active phase.
Definition InitStateEquil_impl.hpp:1034
bool oilActive() const
Predicate for whether or not oil is an active phase.
Definition InitStateEquil_impl.hpp:1027
PressureTable(const double gravity, const int samplePoints=2000)
Constructor.
Definition InitStateEquil_impl.hpp:960
double water(const double depth) const
Evaluate water phase pressure at specified depth.
Definition InitStateEquil_impl.hpp:1067
double oil(const double depth) const
Evaluate oil phase pressure at specified depth.
Definition InitStateEquil_impl.hpp:1048
bool water(const PhaseUsage &pu)
Active water predicate.
Definition RegionAttributeHelpers.hpp:308
bool oil(const PhaseUsage &pu)
Active oil predicate.
Definition RegionAttributeHelpers.hpp:321
bool gas(const PhaseUsage &pu)
Active gas predicate.
Definition RegionAttributeHelpers.hpp:334
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition BlackoilPhases.hpp:27
Simple set of per-phase (named by primary component) quantities.
Definition InitStateEquil.hpp:328
Evaluation point within a model geometry.
Definition InitStateEquil.hpp:381