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>
50#include <fmt/format.h>
63template <
typename CellRange,
typename Comm>
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) {
84 span[0] = comm.min(span[0]);
85 span[1] = comm.max(span[1]);
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>>
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 "
121 2*numIntervals, subdiv);
126template <
class 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>
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);
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));
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>
263 const double press)
const
265 return this->density(depth, press) * g_;
268template<
class Flu
idSystem>
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>
297 const double press)
const
299 return this->density(depth, press) * g_;
302template<
class Flu
idSystem,
class RS>
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>
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>
345 const double press)
const
347 return this->density(depth, press) * g_;
350template<
class Flu
idSystem,
class RV,
class RVW>
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>&
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>&
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>
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>
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;});
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");
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()) {
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();
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();
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;
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()) {
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();
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();
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;
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()) {
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();
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;
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;
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>
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);
1662 const auto cartIx = cartesianIndexMapper_.cartesianIndex(elemIdx);
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);
1799 const auto acc = rec[r].initializationTargetAccuracy();
1801 throw std::runtime_error {
1802 "Cannot initialise model: Positive item 9 is not supported "
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])
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();
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);
#define OPM_END_PARALLEL_TRY_CATCH(prefix, comm)
Catch exception and throw in a parallel try-catch clause.
Definition: DeferredLoggingErrorHelpers.hpp:172
#define OPM_BEGIN_PARALLEL_TRY_CATCH()
Macro to setup the try of a parallel try-catch.
Definition: DeferredLoggingErrorHelpers.hpp:138
Auxiliary routines that to solve the ODEs that emerge from the hydrostatic equilibrium problem.
Dune::OwnerOverlapCopyCommunication< int, int > Comm
Definition: FlexibleSolver_impl.hpp:270
Routines that actually solve the ODEs that emerge from the hydrostatic equilibrium problem.
Definition: InitStateEquil.hpp:680
Definition: InitStateEquil.hpp:132
Gas(const TabulatedFunction &tempVdTable, const RV &rv, const RVW &rvw, const int pvtRegionIdx, const double normGrav)
Definition: InitStateEquil_impl.hpp:329
double operator()(const double depth, const double press) const
Definition: InitStateEquil_impl.hpp:344
Definition: InitStateEquil.hpp:109
double operator()(const double depth, const double press) const
Definition: InitStateEquil_impl.hpp:296
Oil(const TabulatedFunction &tempVdTable, const RS &rs, const int pvtRegionIdx, const double normGrav)
Definition: InitStateEquil_impl.hpp:283
Definition: InitStateEquil.hpp:86
Water(const TabulatedFunction &tempVdTable, const TabulatedFunction &saltVdTable, const int pvtRegionIdx, const double normGrav)
Definition: InitStateEquil_impl.hpp:249
double operator()(const double depth, const double press) const
Definition: InitStateEquil_impl.hpp:262
Definition: InitStateEquil.hpp:376
const PhaseQuantityValue & deriveSaturations(const Position &x, const Region ®, const PTable &ptable)
Definition: InitStateEquil_impl.hpp:574
PhaseSaturations(MaterialLawManager &matLawMgr, const std::vector< double > &swatInit)
Definition: InitStateEquil_impl.hpp:550
Definition: InitStateEquil.hpp:159
PressureTable & operator=(const PressureTable &rhs)
Definition: InitStateEquil_impl.hpp:990
double gas(const double depth) const
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)
Definition: InitStateEquil_impl.hpp:960
void equilibrate(const Region ®, const VSpan &span)
Definition: InitStateEquil_impl.hpp:1016
double water(const double depth) const
Definition: InitStateEquil_impl.hpp:1067
std::array< double, 2 > VSpan
Definition: InitStateEquil.hpp:161
double oil(const double depth) const
Definition: InitStateEquil_impl.hpp:1048
Definition: InitStateEquil.hpp:64
double operator()(const double x) const
Definition: InitStateEquil_impl.hpp:215
RK4IVP(const RHS &f, const std::array< double, 2 > &span, const double y0, const int N)
Definition: InitStateEquil_impl.hpp:180
Definition: EquilibrationHelpers.hpp:220
Definition: EquilibrationHelpers.hpp:271
Definition: EquilibrationHelpers.hpp:168
Definition: EquilibrationHelpers.hpp:322
Definition: EquilibrationHelpers.hpp:372
std::vector< EquilRecord > getEquil(const EclipseState &state)
Definition: InitStateEquil_impl.hpp:1261
std::vector< int > equilnum(const EclipseState &eclipseState, const GridView &gridview)
Definition: InitStateEquil_impl.hpp:1275
std::pair< double, double > cellZMinMax(const Element &element)
Definition: InitStateEquil_impl.hpp:161
void verticalExtent(const CellRange &cells, const std::vector< std::pair< double, double > > &cellZMinMax, const Comm &comm, std::array< double, 2 > &span)
Definition: InitStateEquil_impl.hpp:64
double cellCenterDepth(const Element &element)
Definition: InitStateEquil_impl.hpp:127
std::vector< std::pair< double, double > > horizontalSubdivision(const CellID cell, const std::pair< double, double > topbot, const int numIntervals)
Definition: InitStateEquil_impl.hpp:106
std::pair< double, double > cellZSpan(const Element &element)
Definition: InitStateEquil_impl.hpp:142
void subdivisionCentrePoints(const double left, const double right, const int numIntervals, std::vector< std::pair< double, double > > &subdiv)
Definition: InitStateEquil_impl.hpp:88
bool water(const PhaseUsage &pu)
Definition: RegionAttributeHelpers.hpp:308
bool oil(const PhaseUsage &pu)
Definition: RegionAttributeHelpers.hpp:321
bool gas(const PhaseUsage &pu)
Definition: RegionAttributeHelpers.hpp:334
Definition: BlackoilPhases.hpp:27
std::string to_string(const ConvergenceReport::ReservoirFailure::Type t)
Simple set of per-phase (named by primary component) quantities.
Definition: InitStateEquil.hpp:328
Definition: InitStateEquil.hpp:381