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.
waterActive()) { this->deriveWaterSat(); }
582 if (ptable.
gasActive()) { this->deriveGasSat(); }
584 if (this->isOverlappingTransition()) {
585 this->fixUnphysicalTransition();
588 if (ptable.
oilActive()) { this->deriveOilSat(); }
590 this->accountForScaledSaturations();
595template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
599 const PTable& ptable)
601 this->evalPt_.position = &x;
602 this->evalPt_.region = ®
603 this->evalPt_.ptable = &ptable;
606template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
607void PhaseSaturations<MaterialLawManager,FluidSystem,Region,CellID>::
608initializePhaseQuantities()
611 this->press_.reset();
613 const auto depth = this->evalPt_.position->depth;
614 const auto& ptable = *this->evalPt_.ptable;
616 if (ptable.oilActive()) {
617 this->press_.oil = ptable.oil(depth);
620 if (ptable.gasActive()) {
621 this->press_.gas = ptable.gas(depth);
624 if (ptable.waterActive()) {
625 this->press_.water = ptable.water(depth);
629template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
630void PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::deriveOilSat()
632 this->sat_.oil = 1.0 - this->sat_.water - this->sat_.gas;
635template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
636void PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::deriveGasSat()
638 auto& sg = this->sat_.gas;
640 const auto isIncr =
true;
641 const auto oilActive = this->evalPt_.ptable->oilActive();
643 if (this->isConstCapPress(this->gasPos())) {
647 const auto gas_contact = oilActive? this->evalPt_.region->zgoc() : this->evalPt_.region->zwoc();
648 sg = this->fromDepthTable(gas_contact,
649 this->gasPos(), isIncr);
659 const auto pw = oilActive? this->press_.oil : this->press_.water;
660 const auto pcgo = this->press_.gas - pw;
661 sg = this->invertCapPress(pcgo, this->gasPos(), isIncr);
665template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
666void PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::deriveWaterSat()
668 auto& sw = this->sat_.water;
670 const auto isIncr =
false;
672 if (this->isConstCapPress(this->waterPos())) {
676 sw = this->fromDepthTable(this->evalPt_.region->zwoc(),
677 this->waterPos(), isIncr);
689 const auto pcow = this->press_.oil - this->press_.water;
691 if (this->swatInit_.empty()) {
692 sw = this->invertCapPress(pcow, this->waterPos(), isIncr);
695 auto [swout, newSwatInit] = this->applySwatInit(pcow);
697 sw = this->invertCapPress(pcow, this->waterPos(), isIncr);
705template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
706void PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
707fixUnphysicalTransition()
709 auto& sg = this->sat_.gas;
710 auto& sw = this->sat_.water;
718 const auto pcgw = this->press_.gas - this->press_.water;
719 if (! this->swatInit_.empty()) {
723 auto [swout, newSwatInit] = this->applySwatInit(pcgw, sw);
725 const auto isIncr =
false;
726 sw = this->invertCapPress(pcgw, this->waterPos(), isIncr);
733 sw = satFromSumOfPcs<FluidSystem>
734 (this->matLawMgr_, this->waterPos(), this->gasPos(),
735 this->evalPt_.position->cell, pcgw);
738 this->fluidState_.setSaturation(this->oilPos(), 1.0 - sw - sg);
739 this->fluidState_.setSaturation(this->gasPos(), sg);
740 this->fluidState_.setSaturation(this->waterPos(), this->evalPt_
741 .ptable->waterActive() ? sw : 0.0);
744 this->computeMaterialLawCapPress();
745 this->press_.oil = this->press_.gas - this->materialLawCapPressGasOil();
748template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
749void PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
750accountForScaledSaturations()
752 const auto gasActive = this->evalPt_.ptable->gasActive();
753 const auto watActive = this->evalPt_.ptable->waterActive();
754 const auto oilActive = this->evalPt_.ptable->oilActive();
756 auto sg = gasActive? this->sat_.gas : 0.0;
757 auto sw = watActive? this->sat_.water : 0.0;
758 auto so = oilActive? this->sat_.oil : 0.0;
760 this->fluidState_.setSaturation(this->waterPos(), sw);
761 this->fluidState_.setSaturation(this->oilPos(), so);
762 this->fluidState_.setSaturation(this->gasPos(), sg);
764 const auto& scaledDrainageInfo = this->matLawMgr_
765 .oilWaterScaledEpsInfoDrainage(this->evalPt_.position->cell);
767 const auto thresholdSat = 1.0e-6;
768 if (watActive && ((sw + thresholdSat) > scaledDrainageInfo.Swu)) {
772 this->fluidState_.setSaturation(this->waterPos(), scaledDrainageInfo.Swu);
774 this->fluidState_.setSaturation(this->oilPos(), so + sw - scaledDrainageInfo.Swu);
775 }
else if (gasActive) {
776 this->fluidState_.setSaturation(this->gasPos(), sg + sw - scaledDrainageInfo.Swu);
778 sw = scaledDrainageInfo.Swu;
779 this->computeMaterialLawCapPress();
783 this->press_.oil = this->press_.water + this->materialLawCapPressOilWater();
786 this->press_.gas = this->press_.water + this->materialLawCapPressGasWater();
790 if (gasActive && ((sg + thresholdSat) > scaledDrainageInfo.Sgu)) {
794 this->fluidState_.setSaturation(this->gasPos(), scaledDrainageInfo.Sgu);
796 this->fluidState_.setSaturation(this->oilPos(), so + sg - scaledDrainageInfo.Sgu);
797 }
else if (watActive) {
798 this->fluidState_.setSaturation(this->waterPos(), sw + sg - scaledDrainageInfo.Sgu);
800 sg = scaledDrainageInfo.Sgu;
801 this->computeMaterialLawCapPress();
805 this->press_.oil = this->press_.gas - this->materialLawCapPressGasOil();
808 this->press_.water = this->press_.gas - this->materialLawCapPressGasWater();
812 if (watActive && ((sw - thresholdSat) < scaledDrainageInfo.Swl)) {
816 this->fluidState_.setSaturation(this->waterPos(), scaledDrainageInfo.Swl);
818 this->fluidState_.setSaturation(this->oilPos(), so + sw - scaledDrainageInfo.Swl);
819 }
else if (gasActive) {
820 this->fluidState_.setSaturation(this->gasPos(), sg + sw - scaledDrainageInfo.Swl);
822 sw = scaledDrainageInfo.Swl;
823 this->computeMaterialLawCapPress();
827 this->press_.water = this->press_.oil - this->materialLawCapPressOilWater();
830 this->press_.water = this->press_.gas - this->materialLawCapPressGasWater();
834 if (gasActive && ((sg - thresholdSat) < scaledDrainageInfo.Sgl)) {
838 this->fluidState_.setSaturation(this->gasPos(), scaledDrainageInfo.Sgl);
840 this->fluidState_.setSaturation(this->oilPos(), so + sg - scaledDrainageInfo.Sgl);
841 }
else if (watActive) {
842 this->fluidState_.setSaturation(this->waterPos(), sw + sg - scaledDrainageInfo.Sgl);
844 sg = scaledDrainageInfo.Sgl;
845 this->computeMaterialLawCapPress();
849 this->press_.gas = this->press_.oil + this->materialLawCapPressGasOil();
852 this->press_.gas = this->press_.water + this->materialLawCapPressGasWater();
857template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
858std::pair<double, bool>
859PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
860applySwatInit(
const double pcow)
862 return this->applySwatInit(pcow, this->swatInit_[this->evalPt_.position->cell]);
865template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
866std::pair<double, bool>
867PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
868applySwatInit(
const double pcow,
const double sw)
870 return this->matLawMgr_.applySwatinit(this->evalPt_.position->cell, pcow, sw);
873template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
874void PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
875computeMaterialLawCapPress()
877 const auto& matParams = this->matLawMgr_
878 .materialLawParams(this->evalPt_.position->cell);
880 this->matLawCapPress_.fill(0.0);
881 MaterialLaw::capillaryPressures(this->matLawCapPress_,
882 matParams, this->fluidState_);
885template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
886double PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
887materialLawCapPressGasOil()
const
889 return this->matLawCapPress_[this->oilPos()]
890 + this->matLawCapPress_[this->gasPos()];
893template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
894double PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
895materialLawCapPressOilWater()
const
897 return this->matLawCapPress_[this->oilPos()]
898 - this->matLawCapPress_[this->waterPos()];
901template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
902double PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
903materialLawCapPressGasWater()
const
905 return this->matLawCapPress_[this->gasPos()]
906 - this->matLawCapPress_[this->waterPos()];
909template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
910bool PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
911isConstCapPress(
const PhaseIdx phaseIdx)
const
913 return isConstPc<FluidSystem>
914 (this->matLawMgr_, phaseIdx, this->evalPt_.position->cell);
917template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
918bool PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
919isOverlappingTransition()
const
921 return this->evalPt_.ptable->gasActive()
922 && this->evalPt_.ptable->waterActive()
923 && ((this->sat_.gas + this->sat_.water) > 1.0);
926template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
927double PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
928fromDepthTable(
const double contactdepth,
929 const PhaseIdx phasePos,
930 const bool isincr)
const
932 return satFromDepth<FluidSystem>
933 (this->matLawMgr_, this->evalPt_.position->depth,
934 contactdepth,
static_cast<int>(phasePos),
935 this->evalPt_.position->cell, isincr);
938template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
939double PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
940invertCapPress(
const double pc,
941 const PhaseIdx phasePos,
942 const bool isincr)
const
944 return satFromPc<FluidSystem>
945 (this->matLawMgr_,
static_cast<int>(phasePos),
946 this->evalPt_.position->cell, pc, isincr);
949template<
class Flu
idSystem,
class Region>
952 const int samplePoints)
954 , nsample_(samplePoints)
958template <
class Flu
idSystem,
class Region>
961 : gravity_(rhs.gravity_)
962 , nsample_(rhs.nsample_)
964 this->copyInPointers(rhs);
967template <
class Flu
idSystem,
class Region>
970 : gravity_(rhs.gravity_)
971 , nsample_(rhs.nsample_)
972 , oil_ (std::move(rhs.oil_))
973 , gas_ (std::move(rhs.gas_))
974 , wat_ (std::move(rhs.wat_))
978template <
class Flu
idSystem,
class Region>
983 this->gravity_ = rhs.gravity_;
984 this->nsample_ = rhs.nsample_;
985 this->copyInPointers(rhs);
990template <
class Flu
idSystem,
class Region>
995 this->gravity_ = rhs.gravity_;
996 this->nsample_ = rhs.nsample_;
998 this->oil_ = std::move(rhs.oil_);
999 this->gas_ = std::move(rhs.gas_);
1000 this->wat_ = std::move(rhs.wat_);
1005template <
class Flu
idSystem,
class Region>
1011 auto equil = this->selectEquilibrationStrategy(reg);
1013 (this->*equil)(reg, span);
1016template <
class Flu
idSystem,
class Region>
1020 return FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
1023template <
class Flu
idSystem,
class Region>
1027 return FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
1030template <
class Flu
idSystem,
class Region>
1034 return FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx);
1037template <
class Flu
idSystem,
class Region>
1039oil(
const double depth)
const
1041 this->checkPtr(this->oil_.get(),
"OIL");
1043 return this->oil_->value(depth);
1046template <
class Flu
idSystem,
class Region>
1048gas(
const double depth)
const
1050 this->checkPtr(this->gas_.get(),
"GAS");
1052 return this->gas_->value(depth);
1056template <
class Flu
idSystem,
class Region>
1058water(
const double depth)
const
1060 this->checkPtr(this->wat_.get(),
"WATER");
1062 return this->wat_->value(depth);
1065template <
class Flu
idSystem,
class Region>
1067equil_WOG(
const Region& reg,
const VSpan& span)
1072 if (! this->waterActive()) {
1073 throw std::invalid_argument {
1074 "Don't know how to interpret EQUIL datum depth in "
1075 "WATER zone in model without active water phase"
1080 const auto ic =
typename WPress::InitCond {
1081 reg.datum(), reg.pressure()
1084 this->makeWatPressure(ic, reg, span);
1087 if (this->oilActive()) {
1089 const auto ic =
typename OPress::InitCond {
1091 this->
water(reg.zwoc()) + reg.pcowWoc()
1094 this->makeOilPressure(ic, reg, span);
1097 if (this->gasActive() && this->oilActive()) {
1099 const auto ic =
typename GPress::InitCond {
1101 this->
oil(reg.zgoc()) + reg.pcgoGoc()
1104 this->makeGasPressure(ic, reg, span);
1105 }
else if (this->gasActive() && !this->oilActive()) {
1107 const auto ic =
typename GPress::InitCond {
1109 this->
water(reg.zwoc()) + reg.pcowWoc()
1111 this->makeGasPressure(ic, reg, span);
1115template <
class Flu
idSystem,
class Region>
1116void PressureTable<FluidSystem, Region>::
1117equil_GOW(
const Region& reg,
const VSpan& span)
1122 if (! this->gasActive()) {
1123 throw std::invalid_argument {
1124 "Don't know how to interpret EQUIL datum depth in "
1125 "GAS zone in model without active gas phase"
1130 const auto ic =
typename GPress::InitCond {
1131 reg.datum(), reg.pressure()
1134 this->makeGasPressure(ic, reg, span);
1137 if (this->oilActive()) {
1139 const auto ic =
typename OPress::InitCond {
1141 this->
gas(reg.zgoc()) - reg.pcgoGoc()
1143 this->makeOilPressure(ic, reg, span);
1146 if (this->waterActive() && this->oilActive()) {
1148 const auto ic =
typename WPress::InitCond {
1150 this->
oil(reg.zwoc()) - reg.pcowWoc()
1153 this->makeWatPressure(ic, reg, span);
1154 }
else if (this->waterActive() && !this->oilActive()) {
1156 const auto ic =
typename WPress::InitCond {
1158 this->
gas(reg.zwoc()) - reg.pcowWoc()
1160 this->makeWatPressure(ic, reg, span);
1164template <
class Flu
idSystem,
class Region>
1165void PressureTable<FluidSystem, Region>::
1166equil_OWG(
const Region& reg,
const VSpan& span)
1171 if (! this->oilActive()) {
1172 throw std::invalid_argument {
1173 "Don't know how to interpret EQUIL datum depth in "
1174 "OIL zone in model without active oil phase"
1179 const auto ic =
typename OPress::InitCond {
1180 reg.datum(), reg.pressure()
1183 this->makeOilPressure(ic, reg, span);
1186 if (this->waterActive()) {
1188 const auto ic =
typename WPress::InitCond {
1190 this->
oil(reg.zwoc()) - reg.pcowWoc()
1193 this->makeWatPressure(ic, reg, span);
1196 if (this->gasActive()) {
1198 const auto ic =
typename GPress::InitCond {
1200 this->
oil(reg.zgoc()) + reg.pcgoGoc()
1202 this->makeGasPressure(ic, reg, span);
1206template <
class Flu
idSystem,
class Region>
1207void PressureTable<FluidSystem, Region>::
1208makeOilPressure(
const typename OPress::InitCond& ic,
1212 const auto drho = OilPressODE {
1213 reg.tempVdTable(), reg.dissolutionCalculator(),
1214 reg.pvtIdx(), this->gravity_
1217 this->oil_ = std::make_unique<OPress>(drho, ic, this->nsample_, span);
1220template <
class Flu
idSystem,
class Region>
1221void PressureTable<FluidSystem, Region>::
1222makeGasPressure(
const typename GPress::InitCond& ic,
1226 const auto drho = GasPressODE {
1227 reg.tempVdTable(), reg.evaporationCalculator(), reg.waterEvaporationCalculator(),
1228 reg.pvtIdx(), this->gravity_
1231 this->gas_ = std::make_unique<GPress>(drho, ic, this->nsample_, span);
1234template <
class Flu
idSystem,
class Region>
1235void PressureTable<FluidSystem, Region>::
1236makeWatPressure(
const typename WPress::InitCond& ic,
1240 const auto drho = WatPressODE {
1241 reg.tempVdTable(), reg.saltVdTable(), reg.pvtIdx(), this->gravity_
1244 this->wat_ = std::make_unique<WPress>(drho, ic, this->nsample_, span);
1249namespace DeckDependent {
1251std::vector<EquilRecord>
1252getEquil(
const EclipseState& state)
1254 const auto& init = state.getInitConfig();
1256 if(!init.hasEquil()) {
1257 throw std::domain_error(
"Deck does not provide equilibration data.");
1260 const auto& equil = init.getEquil();
1261 return { equil.begin(), equil.end() };
1264template<
class Gr
idView>
1266equilnum(
const EclipseState& eclipseState,
1267 const GridView& gridview)
1269 std::vector<int> eqlnum(gridview.size(0), 0);
1271 if (eclipseState.fieldProps().has_int(
"EQLNUM")) {
1272 const auto& e = eclipseState.fieldProps().get_int(
"EQLNUM");
1273 std::transform(e.begin(), e.end(), eqlnum.begin(), [](
int n){ return n - 1;});
1275 OPM_BEGIN_PARALLEL_TRY_CATCH();
1276 const int num_regions = eclipseState.getTableManager().getEqldims().getNumEquilRegions();
1277 if ( std::any_of(eqlnum.begin(), eqlnum.end(), [num_regions](
int n){return n >= num_regions;}) ) {
1278 throw std::runtime_error(
"Values larger than maximum Equil regions " + std::to_string(num_regions) +
" provided in EQLNUM");
1280 if ( std::any_of(eqlnum.begin(), eqlnum.end(), [](
int n){return n < 0;}) ) {
1281 throw std::runtime_error(
"zero or negative values provided in EQLNUM");
1283 OPM_END_PARALLEL_TRY_CATCH(
"Invalied EQLNUM numbers: ", gridview.comm());
1288template<
class FluidSystem,
1291 class ElementMapper,
1292 class CartesianIndexMapper>
1293template<
class MaterialLawManager>
1294InitialStateComputer<FluidSystem,
1298 CartesianIndexMapper>::
1299InitialStateComputer(MaterialLawManager& materialLawManager,
1300 const EclipseState& eclipseState,
1302 const GridView& gridView,
1303 const CartesianIndexMapper& cartMapper,
1305 const int num_pressure_points,
1306 const bool applySwatInit)
1307 : temperature_(grid.size(0), eclipseState.getTableManager().rtemp()),
1308 saltConcentration_(grid.size(0)),
1309 saltSaturation_(grid.size(0)),
1310 pp_(FluidSystem::numPhases,
1311 std::vector<double>(grid.size(0))),
1312 sat_(FluidSystem::numPhases,
1313 std::vector<double>(grid.size(0))),
1317 cartesianIndexMapper_(cartMapper),
1318 num_pressure_points_(num_pressure_points)
1321 if (applySwatInit) {
1322 if (eclipseState.fieldProps().has_double(
"SWATINIT")) {
1323 swatInit_ = eclipseState.fieldProps().get_double(
"SWATINIT");
1329 const auto& num_aquifers = eclipseState.aquifer().numericalAquifers();
1330 updateCellProps_(gridView, num_aquifers);
1333 const std::vector<EquilRecord> rec = getEquil(eclipseState);
1334 const auto& tables = eclipseState.getTableManager();
1336 const RegionMapping<> eqlmap(equilnum(eclipseState, grid));
1337 const int invalidRegion = -1;
1338 regionPvtIdx_.resize(rec.size(), invalidRegion);
1339 setRegionPvtIdx(eclipseState, eqlmap);
1342 rsFunc_.reserve(rec.size());
1343 if (FluidSystem::enableDissolvedGas()) {
1344 for (std::size_t i = 0; i < rec.size(); ++i) {
1345 if (eqlmap.cells(i).empty()) {
1346 rsFunc_.push_back(std::shared_ptr<Miscibility::RsVD<FluidSystem>>());
1349 const int pvtIdx = regionPvtIdx_[i];
1350 if (!rec[i].liveOilInitConstantRs()) {
1351 const TableContainer& rsvdTables = tables.getRsvdTables();
1352 const TableContainer& pbvdTables = tables.getPbvdTables();
1353 if (rsvdTables.size() > 0) {
1355 const RsvdTable& rsvdTable = rsvdTables.getTable<RsvdTable>(i);
1356 std::vector<double> depthColumn = rsvdTable.getColumn(
"DEPTH").vectorCopy();
1357 std::vector<double> rsColumn = rsvdTable.getColumn(
"RS").vectorCopy();
1358 rsFunc_.push_back(std::make_shared<Miscibility::RsVD<FluidSystem>>(pvtIdx,
1359 depthColumn, rsColumn));
1360 }
else if (pbvdTables.size() > 0) {
1361 const PbvdTable& pbvdTable = pbvdTables.getTable<PbvdTable>(i);
1362 std::vector<double> depthColumn = pbvdTable.getColumn(
"DEPTH").vectorCopy();
1363 std::vector<double> pbubColumn = pbvdTable.getColumn(
"PBUB").vectorCopy();
1364 rsFunc_.push_back(std::make_shared<Miscibility::PBVD<FluidSystem>>(pvtIdx,
1365 depthColumn, pbubColumn));
1368 throw std::runtime_error(
"Cannot initialise: RSVD or PBVD table not available.");
1373 if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) {
1374 throw std::runtime_error(
"Cannot initialise: when no explicit RSVD table is given, \n"
1375 "datum depth must be at the gas-oil-contact. "
1376 "In EQUIL region "+std::to_string(i + 1)+
" (counting from 1), this does not hold.");
1378 const double pContact = rec[i].datumDepthPressure();
1379 const double TContact = 273.15 + 20;
1380 rsFunc_.push_back(std::make_shared<Miscibility::RsSatAtContact<FluidSystem>>(pvtIdx, pContact, TContact));
1385 for (std::size_t i = 0; i < rec.size(); ++i) {
1386 rsFunc_.push_back(std::make_shared<Miscibility::NoMixing>());
1390 rvFunc_.reserve(rec.size());
1391 if (FluidSystem::enableVaporizedOil()) {
1392 for (std::size_t i = 0; i < rec.size(); ++i) {
1393 if (eqlmap.cells(i).empty()) {
1394 rvFunc_.push_back(std::shared_ptr<Miscibility::RvVD<FluidSystem>>());
1397 const int pvtIdx = regionPvtIdx_[i];
1398 if (!rec[i].wetGasInitConstantRv()) {
1399 const TableContainer& rvvdTables = tables.getRvvdTables();
1400 const TableContainer& pdvdTables = tables.getPdvdTables();
1402 if (rvvdTables.size() > 0) {
1403 const RvvdTable& rvvdTable = rvvdTables.getTable<RvvdTable>(i);
1404 std::vector<double> depthColumn = rvvdTable.getColumn(
"DEPTH").vectorCopy();
1405 std::vector<double> rvColumn = rvvdTable.getColumn(
"RV").vectorCopy();
1406 rvFunc_.push_back(std::make_shared<Miscibility::RvVD<FluidSystem>>(pvtIdx,
1407 depthColumn, rvColumn));
1408 }
else if (pdvdTables.size() > 0) {
1409 const PdvdTable& pdvdTable = pdvdTables.getTable<PdvdTable>(i);
1410 std::vector<double> depthColumn = pdvdTable.getColumn(
"DEPTH").vectorCopy();
1411 std::vector<double> pdewColumn = pdvdTable.getColumn(
"PDEW").vectorCopy();
1412 rvFunc_.push_back(std::make_shared<Miscibility::PDVD<FluidSystem>>(pvtIdx,
1413 depthColumn, pdewColumn));
1415 throw std::runtime_error(
"Cannot initialise: RVVD or PDCD table not available.");
1419 if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) {
1420 throw std::runtime_error(
1421 "Cannot initialise: when no explicit RVVD table is given, \n"
1422 "datum depth must be at the gas-oil-contact. "
1423 "In EQUIL region "+std::to_string(i + 1)+
" (counting from 1), this does not hold.");
1425 const double pContact = rec[i].datumDepthPressure() + rec[i].gasOilContactCapillaryPressure();
1426 const double TContact = 273.15 + 20;
1427 rvFunc_.push_back(std::make_shared<Miscibility::RvSatAtContact<FluidSystem>>(pvtIdx,pContact, TContact));
1432 for (std::size_t i = 0; i < rec.size(); ++i) {
1433 rvFunc_.push_back(std::make_shared<Miscibility::NoMixing>());
1437 rvwFunc_.reserve(rec.size());
1438 if (FluidSystem::enableVaporizedWater()) {
1439 for (std::size_t i = 0; i < rec.size(); ++i) {
1440 if (eqlmap.cells(i).empty()) {
1441 rvwFunc_.push_back(std::shared_ptr<Miscibility::RvwVD<FluidSystem>>());
1444 const int pvtIdx = regionPvtIdx_[i];
1445 if (!rec[i].humidGasInitConstantRvw()) {
1446 const TableContainer& rvwvdTables = tables.getRvwvdTables();
1448 if (rvwvdTables.size() > 0) {
1449 const RvwvdTable& rvwvdTable = rvwvdTables.getTable<RvwvdTable>(i);
1450 std::vector<double> depthColumn = rvwvdTable.getColumn(
"DEPTH").vectorCopy();
1451 std::vector<double> rvwvdColumn = rvwvdTable.getColumn(
"RVWVD").vectorCopy();
1452 rvwFunc_.push_back(std::make_shared<Miscibility::RvwVD<FluidSystem>>(pvtIdx,
1453 depthColumn, rvwvdColumn));
1455 throw std::runtime_error(
"Cannot initialise: RVWVD table not available.");
1459 const auto oilActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
1461 if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) {
1462 rvwFunc_.push_back(std::make_shared<Miscibility::NoMixing>());
1463 const auto msg =
"No explicit RVWVD table is given for EQUIL region " + std::to_string(i + 1) +
". \n"
1464 "and datum depth is not at the gas-oil-contact. \n"
1465 "Rvw is set to 0.0 in all cells. \n";
1466 OpmLog::warning(msg);
1470 const double pContact = rec[i].datumDepthPressure() + rec[i].gasOilContactCapillaryPressure();
1471 const double TContact = 273.15 + 20;
1472 rvwFunc_.push_back(std::make_shared<Miscibility::RvwSatAtContact<FluidSystem>>(pvtIdx,pContact, TContact));
1478 if (rec[i].waterOilContactDepth() != rec[i].datumDepth()) {
1479 rvwFunc_.push_back(std::make_shared<Miscibility::NoMixing>());
1480 const auto msg =
"No explicit RVWVD table is given for EQUIL region " + std::to_string(i + 1) +
". \n"
1481 "and datum depth is not at the gas-water-contact. \n"
1482 "Rvw is set to 0.0 in all cells. \n";
1483 OpmLog::warning(msg);
1486 const double pContact = rec[i].datumDepthPressure() + rec[i].waterOilContactCapillaryPressure();
1487 const double TContact = 273.15 + 20;
1488 rvwFunc_.push_back(std::make_shared<Miscibility::RvwSatAtContact<FluidSystem>>(pvtIdx,pContact, TContact));
1495 for (std::size_t i = 0; i < rec.size(); ++i) {
1496 rvwFunc_.push_back(std::make_shared<Miscibility::NoMixing>());
1502 updateInitialTemperature_(eclipseState, eqlmap);
1505 updateInitialSaltConcentration_(eclipseState, eqlmap);
1508 updateInitialSaltSaturation_(eclipseState, eqlmap);
1511 const auto& comm = grid.comm();
1512 calcPressSatRsRv(eqlmap, rec, materialLawManager, comm, grav);
1515 applyNumericalAquifers_(gridView, num_aquifers, eclipseState.runspec().co2Storage() || eclipseState.runspec().h2Storage());
1521template<
class FluidSystem,
1524 class ElementMapper,
1525 class CartesianIndexMapper>
1527void InitialStateComputer<FluidSystem,
1531 CartesianIndexMapper>::
1532updateInitialTemperature_(
const EclipseState& eclState,
const RMap& reg)
1534 const int numEquilReg = rsFunc_.size();
1535 tempVdTable_.resize(numEquilReg);
1536 const auto& tables = eclState.getTableManager();
1537 if (!tables.hasTables(
"RTEMPVD")) {
1538 std::vector<double> x = {0.0,1.0};
1539 std::vector<double> y = {tables.rtemp(),tables.rtemp()};
1540 for (
auto& table : this->tempVdTable_) {
1541 table.setXYContainers(x, y);
1544 const TableContainer& tempvdTables = tables.getRtempvdTables();
1545 for (std::size_t i = 0; i < tempvdTables.size(); ++i) {
1546 const RtempvdTable& tempvdTable = tempvdTables.getTable<RtempvdTable>(i);
1547 tempVdTable_[i].setXYContainers(tempvdTable.getDepthColumn(), tempvdTable.getTemperatureColumn());
1548 const auto& cells = reg.cells(i);
1549 for (
const auto& cell : cells) {
1550 const double depth = cellCenterDepth_[cell];
1551 this->temperature_[cell] = tempVdTable_[i].eval(depth,
true);
1557template<
class FluidSystem,
1560 class ElementMapper,
1561 class CartesianIndexMapper>
1563void InitialStateComputer<FluidSystem,
1567 CartesianIndexMapper>::
1568updateInitialSaltConcentration_(
const EclipseState& eclState,
const RMap& reg)
1570 const int numEquilReg = rsFunc_.size();
1571 saltVdTable_.resize(numEquilReg);
1572 const auto& tables = eclState.getTableManager();
1573 const TableContainer& saltvdTables = tables.getSaltvdTables();
1576 if (saltvdTables.empty()) {
1577 std::vector<double> x = {0.0,1.0};
1578 std::vector<double> y = {0.0,0.0};
1579 for (
auto& table : this->saltVdTable_) {
1580 table.setXYContainers(x, y);
1583 for (std::size_t i = 0; i < saltvdTables.size(); ++i) {
1584 const SaltvdTable& saltvdTable = saltvdTables.getTable<SaltvdTable>(i);
1585 saltVdTable_[i].setXYContainers(saltvdTable.getDepthColumn(), saltvdTable.getSaltColumn());
1587 const auto& cells = reg.cells(i);
1588 for (
const auto& cell : cells) {
1589 const double depth = cellCenterDepth_[cell];
1590 this->saltConcentration_[cell] = saltVdTable_[i].eval(depth,
true);
1596template<
class FluidSystem,
1599 class ElementMapper,
1600 class CartesianIndexMapper>
1602void InitialStateComputer<FluidSystem,
1606 CartesianIndexMapper>::
1607updateInitialSaltSaturation_(
const EclipseState& eclState,
const RMap& reg)
1609 const int numEquilReg = rsFunc_.size();
1610 saltpVdTable_.resize(numEquilReg);
1611 const auto& tables = eclState.getTableManager();
1612 const TableContainer& saltpvdTables = tables.getSaltpvdTables();
1614 for (std::size_t i = 0; i < saltpvdTables.size(); ++i) {
1615 const SaltpvdTable& saltpvdTable = saltpvdTables.getTable<SaltpvdTable>(i);
1616 saltpVdTable_[i].setXYContainers(saltpvdTable.getDepthColumn(), saltpvdTable.getSaltpColumn());
1618 const auto& cells = reg.cells(i);
1619 for (
const auto& cell : cells) {
1620 const double depth = cellCenterDepth_[cell];
1621 this->saltSaturation_[cell] = saltpVdTable_[i].eval(depth,
true);
1627template<
class FluidSystem,
1630 class ElementMapper,
1631 class CartesianIndexMapper>
1632void InitialStateComputer<FluidSystem,
1636 CartesianIndexMapper>::
1637updateCellProps_(
const GridView& gridView,
1638 const NumericalAquifers& aquifer)
1640 ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
1641 int numElements = gridView.size(0);
1642 cellCenterDepth_.resize(numElements);
1643 cellZSpan_.resize(numElements);
1644 cellZMinMax_.resize(numElements);
1646 auto elemIt = gridView.template begin<0>();
1647 const auto& elemEndIt = gridView.template end<0>();
1648 const auto num_aqu_cells = aquifer.allAquiferCells();
1649 for (; elemIt != elemEndIt; ++elemIt) {
1650 const Element& element = *elemIt;
1651 const unsigned int elemIdx = elemMapper.index(element);
1652 cellCenterDepth_[elemIdx] = Details::cellCenterDepth(element);
1653 const auto cartIx = cartesianIndexMapper_.cartesianIndex(elemIdx);
1654 cellZSpan_[elemIdx] = Details::cellZSpan(element);
1655 cellZMinMax_[elemIdx] = Details::cellZMinMax(element);
1656 if (!num_aqu_cells.empty()) {
1657 const auto search = num_aqu_cells.find(cartIx);
1658 if (search != num_aqu_cells.end()) {
1659 const auto* aqu_cell = num_aqu_cells.at(cartIx);
1660 const double depth_change_num_aqu = aqu_cell->depth - cellCenterDepth_[elemIdx];
1661 cellCenterDepth_[elemIdx] += depth_change_num_aqu;
1662 cellZSpan_[elemIdx].first += depth_change_num_aqu;
1663 cellZSpan_[elemIdx].second += depth_change_num_aqu;
1664 cellZMinMax_[elemIdx].first += depth_change_num_aqu;
1665 cellZMinMax_[elemIdx].second += depth_change_num_aqu;
1671template<
class FluidSystem,
1674 class ElementMapper,
1675 class CartesianIndexMapper>
1676void InitialStateComputer<FluidSystem,
1680 CartesianIndexMapper>::
1681applyNumericalAquifers_(
const GridView& gridView,
1682 const NumericalAquifers& aquifer,
1683 const bool co2store_or_h2store)
1685 const auto num_aqu_cells = aquifer.allAquiferCells();
1686 if (num_aqu_cells.empty())
return;
1689 bool oil_as_brine = co2store_or_h2store && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
1690 const auto watPos = oil_as_brine? FluidSystem::oilPhaseIdx : FluidSystem::waterPhaseIdx;
1691 if (!FluidSystem::phaseIsActive(watPos)){
1692 throw std::logic_error {
"Water phase has to be active for numerical aquifer case" };
1695 ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
1696 auto elemIt = gridView.template begin<0>();
1697 const auto& elemEndIt = gridView.template end<0>();
1698 const auto oilPos = FluidSystem::oilPhaseIdx;
1699 const auto gasPos = FluidSystem::gasPhaseIdx;
1700 for (; elemIt != elemEndIt; ++elemIt) {
1701 const Element& element = *elemIt;
1702 const unsigned int elemIdx = elemMapper.index(element);
1703 const auto cartIx = cartesianIndexMapper_.cartesianIndex(elemIdx);
1704 const auto search = num_aqu_cells.find(cartIx);
1705 if (search != num_aqu_cells.end()) {
1707 this->sat_[watPos][elemIdx] = 1.;
1709 if (!co2store_or_h2store && FluidSystem::phaseIsActive(oilPos)) {
1710 this->sat_[oilPos][elemIdx] = 0.;
1713 if (FluidSystem::phaseIsActive(gasPos)) {
1714 this->sat_[gasPos][elemIdx] = 0.;
1716 const auto* aqu_cell = num_aqu_cells.at(cartIx);
1717 const auto msg = fmt::format(
"FOR AQUIFER CELL AT ({}, {}, {}) OF NUMERICAL "
1718 "AQUIFER {}, WATER SATURATION IS SET TO BE UNITY",
1719 aqu_cell->I+1, aqu_cell->J+1, aqu_cell->K+1, aqu_cell->aquifer_id);
1724 if (aqu_cell->init_pressure) {
1725 const double pres = *(aqu_cell->init_pressure);
1726 this->pp_[watPos][elemIdx] = pres;
1727 if (FluidSystem::phaseIsActive(gasPos)) {
1728 this->pp_[gasPos][elemIdx] = pres;
1730 if (FluidSystem::phaseIsActive(oilPos)) {
1731 this->pp_[oilPos][elemIdx] = pres;
1738template<
class FluidSystem,
1741 class ElementMapper,
1742 class CartesianIndexMapper>
1744void InitialStateComputer<FluidSystem,
1748 CartesianIndexMapper>::
1749setRegionPvtIdx(
const EclipseState& eclState,
const RMap& reg)
1751 const auto& pvtnumData = eclState.fieldProps().get_int(
"PVTNUM");
1753 for (
const auto& r : reg.activeRegions()) {
1754 const auto& cells = reg.cells(r);
1755 regionPvtIdx_[r] = pvtnumData[*cells.begin()] - 1;
1759template<
class FluidSystem,
1762 class ElementMapper,
1763 class CartesianIndexMapper>
1764template<
class RMap,
class MaterialLawManager,
class Comm>
1765void InitialStateComputer<FluidSystem,
1769 CartesianIndexMapper>::
1770calcPressSatRsRv(
const RMap& reg,
1771 const std::vector<EquilRecord>& rec,
1772 MaterialLawManager& materialLawManager,
1776 using PhaseSat = Details::PhaseSaturations<
1777 MaterialLawManager, FluidSystem, EquilReg,
typename RMap::CellId
1780 auto ptable = Details::PressureTable<FluidSystem, EquilReg>{ grav, this->num_pressure_points_ };
1781 auto psat = PhaseSat { materialLawManager, this->swatInit_ };
1782 auto vspan = std::array<double, 2>{};
1784 std::vector<int> regionIsEmpty(rec.size(), 0);
1785 for (std::size_t r = 0; r < rec.size(); ++r) {
1786 const auto& cells = reg.cells(r);
1788 Details::verticalExtent(cells, cellZMinMax_, comm, vspan);
1790 const auto acc = rec[r].initializationTargetAccuracy();
1792 throw std::runtime_error {
1793 "Cannot initialise model: Positive item 9 is not supported "
1794 "in EQUIL keyword, record " + std::to_string(r + 1)
1798 if (cells.empty()) {
1799 regionIsEmpty[r] = 1;
1803 const auto eqreg = EquilReg {
1804 rec[r], this->rsFunc_[r], this->rvFunc_[r], this->rvwFunc_[r], this->tempVdTable_[r], this->saltVdTable_[r], this->regionPvtIdx_[r]
1809 vspan[0] = std::min(vspan[0], std::min(eqreg.zgoc(), eqreg.zwoc()));
1810 vspan[1] = std::max(vspan[1], std::max(eqreg.zgoc(), eqreg.zwoc()));
1812 ptable.equilibrate(eqreg, vspan);
1816 this->equilibrateCellCentres(cells, eqreg, ptable, psat);
1820 this->equilibrateHorizontal(cells, eqreg, -acc,
1829 comm.min(regionIsEmpty.data(),regionIsEmpty.size());
1830 if (comm.rank() == 0) {
1831 for (std::size_t r = 0; r < rec.size(); ++r) {
1832 if (regionIsEmpty[r])
1833 OpmLog::warning(
"Equilibration region " + std::to_string(r + 1)
1834 +
" has no active cells");
1839template<
class FluidSystem,
1842 class ElementMapper,
1843 class CartesianIndexMapper>
1844template<
class CellRange,
class EquilibrationMethod>
1845void InitialStateComputer<FluidSystem,
1849 CartesianIndexMapper>::
1850cellLoop(
const CellRange& cells,
1851 EquilibrationMethod&& eqmethod)
1853 const auto oilPos = FluidSystem::oilPhaseIdx;
1854 const auto gasPos = FluidSystem::gasPhaseIdx;
1855 const auto watPos = FluidSystem::waterPhaseIdx;
1857 const auto oilActive = FluidSystem::phaseIsActive(oilPos);
1858 const auto gasActive = FluidSystem::phaseIsActive(gasPos);
1859 const auto watActive = FluidSystem::phaseIsActive(watPos);
1861 auto pressures = Details::PhaseQuantityValue{};
1862 auto saturations = Details::PhaseQuantityValue{};
1867 for (
const auto& cell : cells) {
1868 eqmethod(cell, pressures, saturations, Rs, Rv, Rvw);
1871 this->pp_ [oilPos][cell] = pressures.oil;
1872 this->sat_[oilPos][cell] = saturations.oil;
1876 this->pp_ [gasPos][cell] = pressures.gas;
1877 this->sat_[gasPos][cell] = saturations.gas;
1881 this->pp_ [watPos][cell] = pressures.water;
1882 this->sat_[watPos][cell] = saturations.water;
1885 if (oilActive && gasActive) {
1886 this->rs_[cell] = Rs;
1887 this->rv_[cell] = Rv;
1890 if (watActive && gasActive) {
1891 this->rvw_[cell] = Rvw;
1896template<
class FluidSystem,
1899 class ElementMapper,
1900 class CartesianIndexMapper>
1901template<
class CellRange,
class PressTable,
class PhaseSat>
1902void InitialStateComputer<FluidSystem,
1906 CartesianIndexMapper>::
1907equilibrateCellCentres(
const CellRange& cells,
1908 const EquilReg& eqreg,
1909 const PressTable& ptable,
1912 using CellPos =
typename PhaseSat::Position;
1913 using CellID = std::remove_cv_t<std::remove_reference_t<
1914 decltype(std::declval<CellPos>().cell)>>;
1915 this->cellLoop(cells, [
this, &eqreg, &ptable, &psat]
1917 Details::PhaseQuantityValue& pressures,
1918 Details::PhaseQuantityValue& saturations,
1921 double& Rvw) ->
void
1923 const auto pos = CellPos {
1924 cell, cellCenterDepth_[cell]
1927 saturations = psat.deriveSaturations(pos, eqreg, ptable);
1928 pressures = psat.correctedPhasePressures();
1930 const auto temp = this->temperature_[cell];
1932 Rs = eqreg.dissolutionCalculator()
1933 (pos.depth, pressures.oil, temp, saturations.gas);
1935 Rv = eqreg.evaporationCalculator()
1936 (pos.depth, pressures.gas, temp, saturations.oil);
1938 Rvw = eqreg.waterEvaporationCalculator()
1939 (pos.depth, pressures.gas, temp, saturations.water);
1943template<
class FluidSystem,
1946 class ElementMapper,
1947 class CartesianIndexMapper>
1948template<
class CellRange,
class PressTable,
class PhaseSat>
1949void InitialStateComputer<FluidSystem,
1953 CartesianIndexMapper>::
1954equilibrateHorizontal(
const CellRange& cells,
1955 const EquilReg& eqreg,
1957 const PressTable& ptable,
1960 using CellPos =
typename PhaseSat::Position;
1961 using CellID = std::remove_cv_t<std::remove_reference_t<
1962 decltype(std::declval<CellPos>().cell)>>;
1964 this->cellLoop(cells, [
this, acc, &eqreg, &ptable, &psat]
1966 Details::PhaseQuantityValue& pressures,
1967 Details::PhaseQuantityValue& saturations,
1970 double& Rvw) ->
void
1973 saturations.reset();
1976 for (
const auto& [depth, frac] : Details::horizontalSubdivision(cell, cellZSpan_[cell], acc)) {
1977 const auto pos = CellPos { cell, depth };
1979 saturations.axpy(psat.deriveSaturations(pos, eqreg, ptable), frac);
1980 pressures .axpy(psat.correctedPhasePressures(), frac);
1986 saturations /= totfrac;
1987 pressures /= totfrac;
1990 const auto pos = CellPos {
1991 cell, cellCenterDepth_[cell]
1994 saturations = psat.deriveSaturations(pos, eqreg, ptable);
1995 pressures = psat.correctedPhasePressures();
1998 const auto temp = this->temperature_[cell];
1999 const auto cz = cellCenterDepth_[cell];
2001 Rs = eqreg.dissolutionCalculator()
2002 (cz, pressures.oil, temp, saturations.gas);
2004 Rv = eqreg.evaporationCalculator()
2005 (cz, pressures.gas, temp, saturations.oil);
2007 Rvw = eqreg.waterEvaporationCalculator()
2008 (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:981
double gas(const double depth) const
Evaluate gas phase pressure at specified depth.
Definition InitStateEquil_impl.hpp:1048
bool waterActive() const
Predicate for whether or not water is an active phase.
Definition InitStateEquil_impl.hpp:1032
bool gasActive() const
Predicate for whether or not gas is an active phase.
Definition InitStateEquil_impl.hpp:1025
bool oilActive() const
Predicate for whether or not oil is an active phase.
Definition InitStateEquil_impl.hpp:1018
PressureTable(const double gravity, const int samplePoints=2000)
Constructor.
Definition InitStateEquil_impl.hpp:951
double water(const double depth) const
Evaluate water phase pressure at specified depth.
Definition InitStateEquil_impl.hpp:1058
double oil(const double depth) const
Evaluate oil phase pressure at specified depth.
Definition InitStateEquil_impl.hpp:1039
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