27#ifndef OPM_WET_GAS_PVT_HPP
28#define OPM_WET_GAS_PVT_HPP
35#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
36#include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp>
37#include <opm/input/eclipse/Schedule/OilVaporizationProperties.hpp>
41#include <opm/common/OpmLog/OpmLog.hpp>
50template <
class Scalar>
53 using SamplingPoints = std::vector<std::pair<Scalar, Scalar>>;
66 const std::vector<TabulatedTwoDFunction>&
inverseGasB,
68 const std::vector<TabulatedTwoDFunction>&
gasMu,
94 void initFromState(
const EclipseState& eclState,
const Schedule& schedule)
96 const auto& pvtgTables = eclState.getTableManager().getPvtgTables();
97 const auto& densityTable = eclState.getTableManager().getDensityTable();
99 assert(pvtgTables.size() == densityTable.size());
104 for (
unsigned regionIdx = 0; regionIdx <
numRegions; ++ regionIdx) {
105 Scalar rhoRefO = densityTable[regionIdx].oil;
106 Scalar rhoRefG = densityTable[regionIdx].gas;
107 Scalar rhoRefW = densityTable[regionIdx].water;
112 for (
unsigned regionIdx = 0; regionIdx <
numRegions; ++ regionIdx) {
113 const auto& pvtgTable = pvtgTables[regionIdx];
115 const auto& saturatedTable = pvtgTable.getSaturatedTable();
116 assert(saturatedTable.numRows() > 1);
118 auto&
gasMu = gasMu_[regionIdx];
119 auto& invGasB = inverseGasB_[regionIdx];
120 auto& invSatGasB = inverseSaturatedGasB_[regionIdx];
121 auto& invSatGasBMu = inverseSaturatedGasBMu_[regionIdx];
122 auto& oilVaporizationFac = saturatedOilVaporizationFactorTable_[regionIdx];
124 oilVaporizationFac.setXYArrays(saturatedTable.numRows(),
125 saturatedTable.getColumn(
"PG"),
126 saturatedTable.getColumn(
"RV"));
128 std::vector<Scalar> invSatGasBArray;
129 std::vector<Scalar> invSatGasBMuArray;
132 for (
unsigned outerIdx = 0; outerIdx < saturatedTable.numRows(); ++ outerIdx) {
133 Scalar pg = saturatedTable.get(
"PG" , outerIdx);
134 Scalar B = saturatedTable.get(
"BG" , outerIdx);
135 Scalar mu = saturatedTable.get(
"MUG" , outerIdx);
137 invGasB.appendXPos(pg);
138 gasMu.appendXPos(pg);
140 invSatGasBArray.push_back(1.0/B);
141 invSatGasBMuArray.push_back(1.0/(mu*B));
143 assert(invGasB.numX() == outerIdx + 1);
144 assert(
gasMu.numX() == outerIdx + 1);
146 const auto& underSaturatedTable = pvtgTable.getUnderSaturatedTable(outerIdx);
147 size_t numRows = underSaturatedTable.numRows();
148 for (
size_t innerIdx = 0; innerIdx < numRows; ++ innerIdx) {
149 Scalar Rv = underSaturatedTable.get(
"RV" , innerIdx);
150 Scalar Bg = underSaturatedTable.get(
"BG" , innerIdx);
151 Scalar mug = underSaturatedTable.get(
"MUG" , innerIdx);
153 invGasB.appendSamplePoint(outerIdx, Rv, 1.0/Bg);
154 gasMu.appendSamplePoint(outerIdx, Rv, mug);
159 std::vector<double> tmpPressure = saturatedTable.getColumn(
"PG").vectorCopy( );
161 invSatGasB.setXYContainers(tmpPressure, invSatGasBArray);
162 invSatGasBMu.setXYContainers(tmpPressure, invSatGasBMuArray);
166 for (
unsigned xIdx = 0; xIdx < invGasB.numX(); ++xIdx) {
168 assert(invGasB.numY(xIdx) > 0);
172 if (invGasB.numY(xIdx) > 1)
178 size_t masterTableIdx = xIdx + 1;
179 for (; masterTableIdx < saturatedTable.numRows(); ++masterTableIdx)
181 if (pvtgTable.getUnderSaturatedTable(masterTableIdx).numRows() > 1)
185 if (masterTableIdx >= saturatedTable.numRows())
186 throw std::runtime_error(
"PVTG tables are invalid: The last table must exhibit at least one "
187 "entry for undersaturated gas!");
191 extendPvtgTable_(regionIdx,
193 pvtgTable.getUnderSaturatedTable(xIdx),
194 pvtgTable.getUnderSaturatedTable(masterTableIdx));
199 const auto& oilVap = schedule[0].oilvap();
200 if (oilVap.getType() == OilVaporizationProperties::OilVaporization::VAPPARS) {
201 vapPar1_ = oilVap.vap1();
208 void extendPvtgTable_(
unsigned regionIdx,
210 const SimpleTable& curTable,
211 const SimpleTable& masterTable)
213 std::vector<double> RvArray = curTable.getColumn(
"RV").vectorCopy();
214 std::vector<double> gasBArray = curTable.getColumn(
"BG").vectorCopy();
215 std::vector<double> gasMuArray = curTable.getColumn(
"MUG").vectorCopy();
217 auto& invGasB = inverseGasB_[regionIdx];
218 auto&
gasMu = gasMu_[regionIdx];
220 for (
size_t newRowIdx = 1; newRowIdx < masterTable.numRows(); ++ newRowIdx) {
221 const auto& RVColumn = masterTable.getColumn(
"RV");
222 const auto& BGColumn = masterTable.getColumn(
"BG");
223 const auto& viscosityColumn = masterTable.getColumn(
"MUG");
226 Scalar diffRv = RVColumn[newRowIdx] - RVColumn[newRowIdx - 1];
227 Scalar newRv = RvArray.back() + diffRv;
230 Scalar B1 = BGColumn[newRowIdx];
231 Scalar B2 = BGColumn[newRowIdx - 1];
232 Scalar x = (B1 - B2)/( (B1 + B2)/2.0 );
236 Scalar newBg = gasBArray.back()*(1.0 + x/2.0)/(1.0 - x/2.0);
239 Scalar mu1 = viscosityColumn[newRowIdx];
240 Scalar mu2 = viscosityColumn[newRowIdx - 1];
241 Scalar xMu = (mu1 - mu2)/( (mu1 + mu2)/2.0 );
245 Scalar newMug = gasMuArray.back()*(1.0 + xMu/2)/(1.0 - xMu/2.0);
249 RvArray.push_back(newRv);
250 gasBArray.push_back(newBg);
251 gasMuArray.push_back(newMug);
254 invGasB.appendSamplePoint(xIdx, newRv, 1.0/newBg);
255 gasMu.appendSamplePoint(xIdx, newRv, newMug);
271 saturatedOilVaporizationFactorTable_.resize(
numRegions);
283 oilReferenceDensity_[regionIdx] = rhoRefOil;
284 gasReferenceDensity_[regionIdx] = rhoRefGas;
293 { saturatedOilVaporizationFactorTable_[regionIdx].setContainerOfTuples(samplePoints); }
306 auto& invGasB = inverseGasB_[regionIdx];
308 const auto& RvTable = saturatedOilVaporizationFactorTable_[regionIdx];
310 constexpr const Scalar T = 273.15 + 15.56;
312 constexpr const Scalar RvMin = 0.0;
313 Scalar RvMax = RvTable.eval(saturatedOilVaporizationFactorTable_[regionIdx].xMax(),
true);
315 Scalar poMin = samplePoints.front().first;
316 Scalar poMax = samplePoints.back().first;
318 constexpr const size_t nRv = 20;
319 size_t nP = samplePoints.size()*2;
321 Scalar rhooRef = oilReferenceDensity_[regionIdx];
326 updateSaturationPressure_(regionIdx);
332 for (
size_t RvIdx = 0; RvIdx < nRv; ++RvIdx) {
333 Scalar Rv = RvMin + (RvMax - RvMin)*RvIdx/nRv;
335 invGasB.appendXPos(Rv);
337 for (
size_t pIdx = 0; pIdx < nP; ++pIdx) {
338 Scalar pg = poMin + (poMax - poMin)*pIdx/nP;
341 Scalar BgSat = gasFormationVolumeFactor.
eval(poSat,
true);
342 Scalar drhoo_dp = (1.1200 - 1.1189)/((5000 - 4000)*6894.76);
343 Scalar rhoo = rhooRef/BgSat*(1 + drhoo_dp*(pg - poSat));
345 Scalar Bg = rhooRef/rhoo;
347 invGasB.appendSamplePoint(RvIdx, pg, 1.0/Bg);
365 { inverseGasB_[regionIdx] = invBg; }
373 { gasMu_[regionIdx] = mug; }
384 auto& oilVaporizationFac = saturatedOilVaporizationFactorTable_[regionIdx];
386 constexpr const Scalar RvMin = 0.0;
387 Scalar RvMax = oilVaporizationFac.eval(saturatedOilVaporizationFactorTable_[regionIdx].xMax(),
true);
389 Scalar poMin = samplePoints.front().first;
390 Scalar poMax = samplePoints.back().first;
392 constexpr const size_t nRv = 20;
393 size_t nP = samplePoints.size()*2;
400 for (
size_t RvIdx = 0; RvIdx < nRv; ++RvIdx) {
401 Scalar Rv = RvMin + (RvMax - RvMin)*RvIdx/nRv;
403 gasMu_[regionIdx].appendXPos(Rv);
405 for (
size_t pIdx = 0; pIdx < nP; ++pIdx) {
406 Scalar pg = poMin + (poMax - poMin)*pIdx/nP;
407 Scalar mug = mugTable.
eval(pg,
true);
409 gasMu_[regionIdx].appendSamplePoint(RvIdx, pg, mug);
421 for (
unsigned regionIdx = 0; regionIdx <
numRegions; ++ regionIdx) {
424 const auto&
gasMu = gasMu_[regionIdx];
425 const auto& invGasB = inverseGasB_[regionIdx];
426 assert(
gasMu.numX() == invGasB.numX());
428 auto& invGasBMu = inverseGasBMu_[regionIdx];
429 auto& invSatGasB = inverseSaturatedGasB_[regionIdx];
430 auto& invSatGasBMu = inverseSaturatedGasBMu_[regionIdx];
432 std::vector<Scalar> satPressuresArray;
433 std::vector<Scalar> invSatGasBArray;
434 std::vector<Scalar> invSatGasBMuArray;
435 for (
size_t pIdx = 0; pIdx <
gasMu.numX(); ++pIdx) {
436 invGasBMu.appendXPos(
gasMu.xAt(pIdx));
438 assert(
gasMu.numY(pIdx) == invGasB.numY(pIdx));
440 size_t numRv =
gasMu.numY(pIdx);
441 for (
size_t rvIdx = 0; rvIdx < numRv; ++rvIdx)
442 invGasBMu.appendSamplePoint(pIdx,
443 gasMu.yAt(pIdx, rvIdx),
444 invGasB.valueAt(pIdx, rvIdx)
445 /
gasMu.valueAt(pIdx, rvIdx));
450 satPressuresArray.push_back(
gasMu.xAt(pIdx));
451 invSatGasBArray.push_back(invGasB.valueAt(pIdx, numRv - 1));
452 invSatGasBMuArray.push_back(invGasBMu.valueAt(pIdx, numRv - 1));
455 invSatGasB.setXYContainers(satPressuresArray, invSatGasBArray);
456 invSatGasBMu.setXYContainers(satPressuresArray, invSatGasBMuArray);
458 updateSaturationPressure_(regionIdx);
466 {
return gasReferenceDensity_.size(); }
471 template <
class Evaluation>
475 const Evaluation&)
const
477 throw std::runtime_error(
"Requested the enthalpy of gas but the thermal option is not enabled");
483 template <
class Evaluation>
486 const Evaluation& pressure,
487 const Evaluation& Rv,
488 const Evaluation& )
const
490 const Evaluation& invBg = inverseGasB_[regionIdx].eval(pressure, Rv,
true);
491 const Evaluation& invMugBg = inverseGasBMu_[regionIdx].eval(pressure, Rv,
true);
493 return invBg/invMugBg;
499 template <
class Evaluation>
502 const Evaluation& pressure)
const
504 const Evaluation& invBg = inverseSaturatedGasB_[regionIdx].eval(pressure,
true);
505 const Evaluation& invMugBg = inverseSaturatedGasBMu_[regionIdx].eval(pressure,
true);
507 return invBg/invMugBg;
513 template <
class Evaluation>
516 const Evaluation& pressure,
517 const Evaluation& Rv,
518 const Evaluation& )
const
519 {
return inverseGasB_[regionIdx].eval(pressure, Rv,
true); }
524 template <
class Evaluation>
527 const Evaluation& pressure)
const
528 {
return inverseSaturatedGasB_[regionIdx].eval(pressure,
true); }
533 template <
class Evaluation>
536 const Evaluation& )
const
542 template <
class Evaluation = Scalar>
546 const Evaluation& )
const
552 template <
class Evaluation>
555 const Evaluation& pressure)
const
557 return saturatedOilVaporizationFactorTable_[regionIdx].eval(pressure,
true);
567 template <
class Evaluation>
570 const Evaluation& pressure,
571 const Evaluation& oilSaturation,
572 Evaluation maxOilSaturation)
const
575 saturatedOilVaporizationFactorTable_[regionIdx].eval(pressure,
true);
579 maxOilSaturation =
min(maxOilSaturation, Scalar(1.0));
580 if (vapPar1_ > 0.0 && maxOilSaturation > 0.01 && oilSaturation < maxOilSaturation) {
581 constexpr const Scalar eps = 0.001;
582 const Evaluation& So =
max(oilSaturation, eps);
583 tmp *=
max(1e-3,
pow(So/maxOilSaturation, vapPar1_));
599 template <
class Evaluation>
602 const Evaluation& Rv)
const
606 const auto& RvTable = saturatedOilVaporizationFactorTable_[regionIdx];
607 constexpr const Scalar eps = std::numeric_limits<typename Toolbox::Scalar>::epsilon()*1e6;
610 Evaluation pSat = saturationPressure_[regionIdx].eval(Rv,
true);
615 bool onProbation =
false;
616 for (
unsigned i = 0; i < 20; ++i) {
617 const Evaluation& f = RvTable.eval(pSat,
true) - Rv;
618 const Evaluation& fPrime = RvTable.evalDerivative(pSat,
true);
626 const Evaluation& delta = f/fPrime;
644 std::stringstream errlog;
645 errlog <<
"Finding saturation pressure did not converge:"
646 <<
" pSat = " << pSat
649 OpmLog::debug(
"Wet gas saturation pressure", errlog.str());
654 template <
class Evaluation>
659 throw std::runtime_error(
"Not implemented: The PVT model does not provide a diffusionCoefficient()");
663 {
return gasReferenceDensity_[regionIdx]; }
666 {
return oilReferenceDensity_[regionIdx]; }
673 return inverseSaturatedGasB_;
676 const std::vector<TabulatedTwoDFunction>&
gasMu()
const {
681 return inverseGasBMu_;
685 return inverseSaturatedGasBMu_;
689 return saturatedOilVaporizationFactorTable_;
693 return saturationPressure_;
702 return this->gasReferenceDensity_ == data.gasReferenceDensity_ &&
703 this->oilReferenceDensity_ == data.oilReferenceDensity_ &&
715 void updateSaturationPressure_(
unsigned regionIdx)
717 const auto& oilVaporizationFac = saturatedOilVaporizationFactorTable_[regionIdx];
721 size_t n = oilVaporizationFac.numSamples();
722 Scalar delta = (oilVaporizationFac.xMax() - oilVaporizationFac.xMin())/Scalar(n + 1);
724 SamplingPoints pSatSamplePoints;
726 for (
size_t i = 0; i <= n; ++ i) {
727 Scalar pSat = oilVaporizationFac.xMin() + Scalar(i)*delta;
730 pSatSamplePoints.emplace_back(Rv, pSat);
734 auto x_coord_comparator = [](
const auto& a,
const auto& b) {
return a.first == b.first; };
735 auto last = std::unique(pSatSamplePoints.begin(), pSatSamplePoints.end(), x_coord_comparator);
736 if (std::distance(pSatSamplePoints.begin(), last) > 1)
737 pSatSamplePoints.erase(last, pSatSamplePoints.end());
739 saturationPressure_[regionIdx].setContainerOfTuples(pSatSamplePoints);
742 std::vector<Scalar> gasReferenceDensity_;
743 std::vector<Scalar> oilReferenceDensity_;
744 std::vector<TabulatedTwoDFunction> inverseGasB_;
745 std::vector<TabulatedOneDFunction> inverseSaturatedGasB_;
746 std::vector<TabulatedTwoDFunction> gasMu_;
747 std::vector<TabulatedTwoDFunction> inverseGasBMu_;
748 std::vector<TabulatedOneDFunction> inverseSaturatedGasBMu_;
749 std::vector<TabulatedOneDFunction> saturatedOilVaporizationFactorTable_;
750 std::vector<TabulatedOneDFunction> saturationPressure_;
Definition: Exceptions.hpp:46
Implements a linearly interpolated scalar function that depends on one variable.
Definition: Tabulated1DFunction.hpp:48
Evaluation eval(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline at a given position.
Definition: Tabulated1DFunction.hpp:259
void setContainerOfTuples(const XYContainer &points, bool sortInputs=true)
Set the sampling points of the piecewise linear function using a STL-compatible container of tuple-li...
Definition: Tabulated1DFunction.hpp:185
This class represents the Pressure-Volume-Temperature relations of the gas phas with vaporized oil.
Definition: WetGasPvt.hpp:52
const std::vector< TabulatedTwoDFunction > & inverseGasB() const
Definition: WetGasPvt.hpp:668
const std::vector< TabulatedOneDFunction > & inverseSaturatedGasBMu() const
Definition: WetGasPvt.hpp:684
const std::vector< TabulatedTwoDFunction > & gasMu() const
Definition: WetGasPvt.hpp:676
Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &, const Evaluation &pressure) const
Returns the formation volume factor [-] of oil saturated gas at a given pressure.
Definition: WetGasPvt.hpp:525
Evaluation saturatedOilVaporizationFactor(unsigned regionIdx, const Evaluation &, const Evaluation &pressure) const
Returns the oil vaporization factor [m^3/m^3] of the gas phase.
Definition: WetGasPvt.hpp:553
void initEnd()
Finish initializing the gas phase PVT properties.
Definition: WetGasPvt.hpp:417
const std::vector< TabulatedOneDFunction > & inverseSaturatedGasB() const
Definition: WetGasPvt.hpp:672
bool operator==(const WetGasPvt< Scalar > &data) const
Definition: WetGasPvt.hpp:700
const std::vector< TabulatedOneDFunction > & saturatedOilVaporizationFactorTable() const
Definition: WetGasPvt.hpp:688
const std::vector< TabulatedTwoDFunction > & inverseGasBMu() const
Definition: WetGasPvt.hpp:680
Evaluation saturationPressure(unsigned regionIdx, const Evaluation &, const Evaluation &Rv) const
Returns the saturation pressure of the gas phase [Pa] depending on its mass fraction of the oil compo...
Definition: WetGasPvt.hpp:600
WetGasPvt()
Definition: WetGasPvt.hpp:59
Evaluation saturatedOilVaporizationFactor(unsigned regionIdx, const Evaluation &, const Evaluation &pressure, const Evaluation &oilSaturation, Evaluation maxOilSaturation) const
Returns the oil vaporization factor [m^3/m^3] of the gas phase.
Definition: WetGasPvt.hpp:568
void setGasViscosity(unsigned regionIdx, const TabulatedTwoDFunction &mug)
Initialize the viscosity of the gas phase.
Definition: WetGasPvt.hpp:372
Evaluation internalEnergy(unsigned, const Evaluation &, const Evaluation &, const Evaluation &) const
Returns the specific enthalpy [J/kg] of gas given a set of parameters.
Definition: WetGasPvt.hpp:472
const std::vector< TabulatedOneDFunction > & saturationPressure() const
Definition: WetGasPvt.hpp:692
Evaluation saturatedWaterVaporizationFactor(unsigned, const Evaluation &, const Evaluation &) const
Returns the water vaporization factor [m^3/m^3] of the gasphase.
Definition: WetGasPvt.hpp:534
Evaluation saturatedWaterVaporizationFactor(unsigned, const Evaluation &, const Evaluation &, const Evaluation &) const
Returns the water vaporization factor [m^3/m^3] of water saturated gas.
Definition: WetGasPvt.hpp:543
unsigned numRegions() const
Return the number of PVT regions which are considered by this PVT-object.
Definition: WetGasPvt.hpp:465
void setNumRegions(size_t numRegions)
Definition: WetGasPvt.hpp:262
void setSaturatedGasViscosity(unsigned regionIdx, const SamplingPoints &samplePoints)
Initialize the phase viscosity for oil saturated gas.
Definition: WetGasPvt.hpp:382
void setSaturatedGasFormationVolumeFactor(unsigned regionIdx, const SamplingPoints &samplePoints)
Initialize the function for the gas formation volume factor.
Definition: WetGasPvt.hpp:304
Evaluation inverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &, const Evaluation &pressure, const Evaluation &Rv, const Evaluation &) const
Returns the formation volume factor [-] of the fluid phase.
Definition: WetGasPvt.hpp:514
const Scalar oilReferenceDensity(unsigned regionIdx) const
Definition: WetGasPvt.hpp:665
void setSaturatedGasOilVaporizationFactor(unsigned regionIdx, const SamplingPoints &samplePoints)
Initialize the function for the oil vaporization factor .
Definition: WetGasPvt.hpp:292
Evaluation diffusionCoefficient(const Evaluation &, const Evaluation &, unsigned) const
Definition: WetGasPvt.hpp:655
void setInverseGasFormationVolumeFactor(unsigned regionIdx, const TabulatedTwoDFunction &invBg)
Initialize the function for the gas formation volume factor.
Definition: WetGasPvt.hpp:364
const Scalar gasReferenceDensity(unsigned regionIdx) const
Definition: WetGasPvt.hpp:662
WetGasPvt(const std::vector< Scalar > &gasReferenceDensity, const std::vector< Scalar > &oilReferenceDensity, const std::vector< TabulatedTwoDFunction > &inverseGasB, const std::vector< TabulatedOneDFunction > &inverseSaturatedGasB, const std::vector< TabulatedTwoDFunction > &gasMu, const std::vector< TabulatedTwoDFunction > &inverseGasBMu, const std::vector< TabulatedOneDFunction > &inverseSaturatedGasBMu, const std::vector< TabulatedOneDFunction > &saturatedOilVaporizationFactorTable, const std::vector< TabulatedOneDFunction > &saturationPressure, Scalar vapPar1)
Definition: WetGasPvt.hpp:64
Scalar vapPar1() const
Definition: WetGasPvt.hpp:696
Evaluation saturatedViscosity(unsigned regionIdx, const Evaluation &, const Evaluation &pressure) const
Returns the dynamic viscosity [Pa s] of oil saturated gas at a given pressure.
Definition: WetGasPvt.hpp:500
Evaluation viscosity(unsigned regionIdx, const Evaluation &, const Evaluation &pressure, const Evaluation &Rv, const Evaluation &) const
Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
Definition: WetGasPvt.hpp:484
void setReferenceDensities(unsigned regionIdx, Scalar rhoRefOil, Scalar rhoRefGas, Scalar)
Initialize the reference densities of all fluids for a given PVT region.
Definition: WetGasPvt.hpp:278
Definition: Air_Mesitylene.hpp:34
ReturnEval_< Evaluation1, Evaluation2 >::type min(const Evaluation1 &arg1, const Evaluation2 &arg2)
Definition: MathToolbox.hpp:346
ReturnEval_< Evaluation1, Evaluation2 >::type max(const Evaluation1 &arg1, const Evaluation2 &arg2)
Definition: MathToolbox.hpp:341
auto scalarValue(const Evaluation &val) -> decltype(MathToolbox< Evaluation >::scalarValue(val))
Definition: MathToolbox.hpp:335
Evaluation abs(const Evaluation &value)
Definition: MathToolbox.hpp:350
ReturnEval_< Evaluation1, Evaluation2 >::type pow(const Evaluation1 &base, const Evaluation2 &exp)
Definition: MathToolbox.hpp:416