27#ifndef OPM_DRY_HUMID_GAS_PVT_HPP
28#define OPM_DRY_HUMID_GAS_PVT_HPP
36#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
37#include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp>
42#include <opm/common/OpmLog/OpmLog.hpp>
51template <
class Scalar>
54 using SamplingPoints = std::vector<std::pair<Scalar, Scalar>>;
67 const std::vector<TabulatedTwoDFunction>&
inverseGasB,
69 const std::vector<TabulatedTwoDFunction>&
gasMu,
95 void initFromState(
const EclipseState& eclState,
const Schedule&)
97 const auto& pvtgwTables = eclState.getTableManager().getPvtgwTables();
98 const auto& densityTable = eclState.getTableManager().getDensityTable();
100 assert(pvtgwTables.size() == densityTable.size());
105 for (
unsigned regionIdx = 0; regionIdx <
numRegions; ++ regionIdx) {
106 Scalar rhoRefO = densityTable[regionIdx].oil;
107 Scalar rhoRefG = densityTable[regionIdx].gas;
108 Scalar rhoRefW = densityTable[regionIdx].water;
113 enableRwgSalt_ = !eclState.getTableManager().getRwgSaltTables().empty();
116 const auto& rwgsaltTables = eclState.getTableManager().getRwgSaltTables();
118 for (
unsigned regionIdx = 0; regionIdx <
numRegions; ++ regionIdx) {
119 const auto& rwgsaltTable = rwgsaltTables[regionIdx];
120 const auto& saturatedTable = rwgsaltTable.getSaturatedTable();
121 assert(saturatedTable.numRows() > 1);
123 auto& waterVaporizationFac = saturatedWaterVaporizationSaltFactorTable_[regionIdx];
124 for (
unsigned outerIdx = 0; outerIdx < saturatedTable.numRows(); ++ outerIdx) {
125 const auto& underSaturatedTable = rwgsaltTable.getUnderSaturatedTable(outerIdx);
126 Scalar pg = saturatedTable.get(
"PG" , outerIdx);
127 waterVaporizationFac.appendXPos(pg);
129 size_t numRows = underSaturatedTable.numRows();
130 for (
size_t innerIdx = 0; innerIdx < numRows; ++ innerIdx) {
131 Scalar saltConcentration = underSaturatedTable.get(
"C_SALT" , innerIdx);
132 Scalar rvwSat= underSaturatedTable.get(
"RVW" , innerIdx);
134 waterVaporizationFac.appendSamplePoint(outerIdx, saltConcentration, rvwSat);
140 for (
unsigned regionIdx = 0; regionIdx <
numRegions; ++ regionIdx) {
141 const auto& pvtgwTable = pvtgwTables[regionIdx];
143 const auto& saturatedTable = pvtgwTable.getSaturatedTable();
144 assert(saturatedTable.numRows() > 1);
146 auto&
gasMu = gasMu_[regionIdx];
147 auto& invGasB = inverseGasB_[regionIdx];
148 auto& invSatGasB = inverseSaturatedGasB_[regionIdx];
149 auto& invSatGasBMu = inverseSaturatedGasBMu_[regionIdx];
150 auto& waterVaporizationFac = saturatedWaterVaporizationFactorTable_[regionIdx];
152 waterVaporizationFac.setXYArrays(saturatedTable.numRows(),
153 saturatedTable.getColumn(
"PG"),
154 saturatedTable.getColumn(
"RW"));
156 std::vector<Scalar> invSatGasBArray;
157 std::vector<Scalar> invSatGasBMuArray;
160 for (
unsigned outerIdx = 0; outerIdx < saturatedTable.numRows(); ++ outerIdx) {
161 Scalar pg = saturatedTable.get(
"PG" , outerIdx);
162 Scalar B = saturatedTable.get(
"BG" , outerIdx);
163 Scalar mu = saturatedTable.get(
"MUG" , outerIdx);
165 invGasB.appendXPos(pg);
166 gasMu.appendXPos(pg);
168 invSatGasBArray.push_back(1.0/B);
169 invSatGasBMuArray.push_back(1.0/(mu*B));
171 assert(invGasB.numX() == outerIdx + 1);
172 assert(
gasMu.numX() == outerIdx + 1);
174 const auto& underSaturatedTable = pvtgwTable.getUnderSaturatedTable(outerIdx);
175 size_t numRows = underSaturatedTable.numRows();
176 for (
size_t innerIdx = 0; innerIdx < numRows; ++ innerIdx) {
177 Scalar Rw = underSaturatedTable.get(
"RW" , innerIdx);
178 Scalar Bg = underSaturatedTable.get(
"BG" , innerIdx);
179 Scalar mug = underSaturatedTable.get(
"MUG" , innerIdx);
181 invGasB.appendSamplePoint(outerIdx, Rw, 1.0/Bg);
182 gasMu.appendSamplePoint(outerIdx, Rw, mug);
187 std::vector<double> tmpPressure = saturatedTable.getColumn(
"PG").vectorCopy( );
189 invSatGasB.setXYContainers(tmpPressure, invSatGasBArray);
190 invSatGasBMu.setXYContainers(tmpPressure, invSatGasBMuArray);
194 for (
unsigned xIdx = 0; xIdx < invGasB.numX(); ++xIdx) {
196 assert(invGasB.numY(xIdx) > 0);
200 if (invGasB.numY(xIdx) > 1)
206 size_t masterTableIdx = xIdx + 1;
207 for (; masterTableIdx < saturatedTable.numRows(); ++masterTableIdx)
209 if (pvtgwTable.getUnderSaturatedTable(masterTableIdx).numRows() > 1)
213 if (masterTableIdx >= saturatedTable.numRows())
214 throw std::runtime_error(
"PVTGW tables are invalid: The last table must exhibit at least one "
215 "entry for undersaturated gas!");
219 extendPvtgwTable_(regionIdx,
221 pvtgwTable.getUnderSaturatedTable(xIdx),
222 pvtgwTable.getUnderSaturatedTable(masterTableIdx));
232 void extendPvtgwTable_(
unsigned regionIdx,
234 const SimpleTable& curTable,
235 const SimpleTable& masterTable)
237 std::vector<double> RwArray = curTable.getColumn(
"RW").vectorCopy();
238 std::vector<double> gasBArray = curTable.getColumn(
"BG").vectorCopy();
239 std::vector<double> gasMuArray = curTable.getColumn(
"MUG").vectorCopy();
241 auto& invGasB = inverseGasB_[regionIdx];
242 auto&
gasMu = gasMu_[regionIdx];
244 for (
size_t newRowIdx = 1; newRowIdx < masterTable.numRows(); ++ newRowIdx) {
245 const auto& RWColumn = masterTable.getColumn(
"RW");
246 const auto& BGColumn = masterTable.getColumn(
"BG");
247 const auto& viscosityColumn = masterTable.getColumn(
"MUG");
250 Scalar diffRw = RWColumn[newRowIdx] - RWColumn[newRowIdx - 1];
251 Scalar newRw = RwArray.back() + diffRw;
254 Scalar B1 = BGColumn[newRowIdx];
255 Scalar B2 = BGColumn[newRowIdx - 1];
256 Scalar x = (B1 - B2)/( (B1 + B2)/2.0 );
260 Scalar newBg = gasBArray.back()*(1.0 + x/2.0)/(1.0 - x/2.0);
263 Scalar mu1 = viscosityColumn[newRowIdx];
264 Scalar mu2 = viscosityColumn[newRowIdx - 1];
265 Scalar xMu = (mu1 - mu2)/( (mu1 + mu2)/2.0 );
269 Scalar newMug = gasMuArray.back()*(1.0 + xMu/2)/(1.0 - xMu/2.0);
273 RwArray.push_back(newRw);
274 gasBArray.push_back(newBg);
275 gasMuArray.push_back(newMug);
278 invGasB.appendSamplePoint(xIdx, newRw, 1.0/newBg);
279 gasMu.appendSamplePoint(xIdx, newRw, newMug);
295 saturatedWaterVaporizationFactorTable_.resize(
numRegions);
307 waterReferenceDensity_[regionIdx] = rhoRefWater;
308 gasReferenceDensity_[regionIdx] = rhoRefGas;
317 { saturatedWaterVaporizationFactorTable_[regionIdx].setContainerOfTuples(samplePoints); }
332 { inverseGasB_[regionIdx] = invBg; }
340 { gasMu_[regionIdx] = mug; }
351 auto& waterVaporizationFac = saturatedWaterVaporizationFactorTable_[regionIdx];
353 constexpr const Scalar RwMin = 0.0;
354 Scalar RwMax = waterVaporizationFac.eval(saturatedWaterVaporizationFactorTable_[regionIdx].xMax(),
true);
356 Scalar poMin = samplePoints.front().first;
357 Scalar poMax = samplePoints.back().first;
359 constexpr const size_t nRw = 20;
360 size_t nP = samplePoints.size()*2;
367 for (
size_t RwIdx = 0; RwIdx < nRw; ++RwIdx) {
368 Scalar Rw = RwMin + (RwMax - RwMin)*RwIdx/nRw;
370 gasMu_[regionIdx].appendXPos(Rw);
372 for (
size_t pIdx = 0; pIdx < nP; ++pIdx) {
373 Scalar pg = poMin + (poMax - poMin)*pIdx/nP;
374 Scalar mug = mugTable.
eval(pg,
true);
376 gasMu_[regionIdx].appendSamplePoint(RwIdx, pg, mug);
388 for (
unsigned regionIdx = 0; regionIdx <
numRegions; ++ regionIdx) {
391 const auto&
gasMu = gasMu_[regionIdx];
392 const auto& invGasB = inverseGasB_[regionIdx];
393 assert(
gasMu.numX() == invGasB.numX());
395 auto& invGasBMu = inverseGasBMu_[regionIdx];
396 auto& invSatGasB = inverseSaturatedGasB_[regionIdx];
397 auto& invSatGasBMu = inverseSaturatedGasBMu_[regionIdx];
399 std::vector<Scalar> satPressuresArray;
400 std::vector<Scalar> invSatGasBArray;
401 std::vector<Scalar> invSatGasBMuArray;
402 for (
size_t pIdx = 0; pIdx <
gasMu.numX(); ++pIdx) {
403 invGasBMu.appendXPos(
gasMu.xAt(pIdx));
405 assert(
gasMu.numY(pIdx) == invGasB.numY(pIdx));
407 size_t numRw =
gasMu.numY(pIdx);
408 for (
size_t RwIdx = 0; RwIdx < numRw; ++RwIdx)
409 invGasBMu.appendSamplePoint(pIdx,
410 gasMu.yAt(pIdx, RwIdx),
411 invGasB.valueAt(pIdx, RwIdx)
412 /
gasMu.valueAt(pIdx, RwIdx));
417 satPressuresArray.push_back(
gasMu.xAt(pIdx));
418 invSatGasBArray.push_back(invGasB.valueAt(pIdx, numRw - 1));
419 invSatGasBMuArray.push_back(invGasBMu.valueAt(pIdx, numRw - 1));
422 invSatGasB.setXYContainers(satPressuresArray, invSatGasBArray);
423 invSatGasBMu.setXYContainers(satPressuresArray, invSatGasBMuArray);
425 updateSaturationPressure_(regionIdx);
433 {
return gasReferenceDensity_.size(); }
438 template <
class Evaluation>
442 const Evaluation&)
const
444 throw std::runtime_error(
"Requested the enthalpy of gas but the thermal option is not enabled");
450 template <
class Evaluation>
453 const Evaluation& pressure,
455 const Evaluation& Rvw)
const
457 const Evaluation& invBg = inverseGasB_[regionIdx].eval(pressure, Rvw,
true);
458 const Evaluation& invMugBg = inverseGasBMu_[regionIdx].eval(pressure, Rvw,
true);
460 return invBg/invMugBg;
466 template <
class Evaluation>
469 const Evaluation& pressure)
const
471 const Evaluation& invBg = inverseSaturatedGasB_[regionIdx].eval(pressure,
true);
472 const Evaluation& invMugBg = inverseSaturatedGasBMu_[regionIdx].eval(pressure,
true);
474 return invBg/invMugBg;
480 template <
class Evaluation>
483 const Evaluation& pressure,
485 const Evaluation& Rvw)
const
486 {
return inverseGasB_[regionIdx].eval(pressure, Rvw,
true); }
491 template <
class Evaluation>
494 const Evaluation& pressure)
const
495 {
return inverseSaturatedGasB_[regionIdx].eval(pressure,
true); }
500 template <
class Evaluation>
503 const Evaluation& pressure)
const
505 return saturatedWaterVaporizationFactorTable_[regionIdx].eval(pressure,
true);
511 template <
class Evaluation>
514 const Evaluation& pressure,
515 const Evaluation& saltConcentration)
const
518 return saturatedWaterVaporizationSaltFactorTable_[regionIdx].eval(pressure, saltConcentration,
true);
520 return saturatedWaterVaporizationFactorTable_[regionIdx].eval(pressure,
true);
527 template <
class Evaluation>
532 const Evaluation& )
const
538 template <
class Evaluation>
541 const Evaluation& )
const
551 template <
class Evaluation>
554 const Evaluation& Rw)
const
558 const auto& RwTable = saturatedWaterVaporizationFactorTable_[regionIdx];
559 const Scalar eps = std::numeric_limits<typename Toolbox::Scalar>::epsilon()*1e6;
562 Evaluation pSat = saturationPressure_[regionIdx].eval(Rw,
true);
567 bool onProbation =
false;
568 for (
unsigned i = 0; i < 20; ++i) {
569 const Evaluation& f = RwTable.eval(pSat,
true) - Rw;
570 const Evaluation& fPrime = RwTable.evalDerivative(pSat,
true);
578 const Evaluation& delta = f/fPrime;
596 std::stringstream errlog;
597 errlog <<
"Finding saturation pressure did not converge:"
598 <<
" pSat = " << pSat
601 OpmLog::debug(
"Wet gas saturation pressure", errlog.str());
606 template <
class Evaluation>
611 throw std::runtime_error(
"Not implemented: The PVT model does not provide a diffusionCoefficient()");
615 {
return gasReferenceDensity_[regionIdx]; }
618 {
return waterReferenceDensity_[regionIdx]; }
625 return inverseSaturatedGasB_;
628 const std::vector<TabulatedTwoDFunction>&
gasMu()
const {
633 return inverseGasBMu_;
637 return inverseSaturatedGasBMu_;
641 return saturatedWaterVaporizationFactorTable_;
645 return saturatedWaterVaporizationSaltFactorTable_;
649 return saturationPressure_;
658 return this->gasReferenceDensity_ == data.gasReferenceDensity_ &&
659 this->waterReferenceDensity_ == data.waterReferenceDensity_ &&
671 void updateSaturationPressure_(
unsigned regionIdx)
673 typedef std::pair<Scalar, Scalar> Pair;
674 const auto& waterVaporizationFac = saturatedWaterVaporizationFactorTable_[regionIdx];
678 size_t n = waterVaporizationFac.numSamples();
679 Scalar delta = (waterVaporizationFac.xMax() - waterVaporizationFac.xMin())/Scalar(n + 1);
681 SamplingPoints pSatSamplePoints;
683 for (
size_t i = 0; i <= n; ++ i) {
684 Scalar pSat = waterVaporizationFac.xMin() + Scalar(i)*delta;
688 pSatSamplePoints.push_back(val);
692 auto x_coord_comparator = [](
const Pair& a,
const Pair& b) {
return a.first == b.first; };
693 auto last = std::unique(pSatSamplePoints.begin(), pSatSamplePoints.end(), x_coord_comparator);
694 if (std::distance(pSatSamplePoints.begin(), last) > 1)
695 pSatSamplePoints.erase(last, pSatSamplePoints.end());
697 saturationPressure_[regionIdx].setContainerOfTuples(pSatSamplePoints);
700 std::vector<Scalar> gasReferenceDensity_;
701 std::vector<Scalar> waterReferenceDensity_;
702 std::vector<TabulatedTwoDFunction> inverseGasB_;
703 std::vector<TabulatedOneDFunction> inverseSaturatedGasB_;
704 std::vector<TabulatedTwoDFunction> gasMu_;
705 std::vector<TabulatedTwoDFunction> inverseGasBMu_;
706 std::vector<TabulatedOneDFunction> inverseSaturatedGasBMu_;
707 std::vector<TabulatedOneDFunction> saturatedWaterVaporizationFactorTable_;
708 std::vector<TabulatedTwoDFunction> saturatedWaterVaporizationSaltFactorTable_;
709 std::vector<TabulatedOneDFunction> saturationPressure_;
This class represents the Pressure-Volume-Temperature relations of the gas phase with vaporized water...
Definition: DryHumidGasPvt.hpp:53
Scalar vapPar1() const
Definition: DryHumidGasPvt.hpp:652
const Scalar gasReferenceDensity(unsigned regionIdx) const
Definition: DryHumidGasPvt.hpp:614
Evaluation inverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &, const Evaluation &pressure, const Evaluation &, const Evaluation &Rvw) const
Returns the formation volume factor [-] of the fluid phase.
Definition: DryHumidGasPvt.hpp:481
Evaluation saturatedOilVaporizationFactor(unsigned, const Evaluation &, const Evaluation &, const Evaluation &, const Evaluation &) const
Returns the oil vaporization factor [m^3/m^3] of the oil phase.
Definition: DryHumidGasPvt.hpp:528
DryHumidGasPvt()
Definition: DryHumidGasPvt.hpp:60
const std::vector< TabulatedOneDFunction > & inverseSaturatedGasB() const
Definition: DryHumidGasPvt.hpp:624
void setGasViscosity(unsigned regionIdx, const TabulatedTwoDFunction &mug)
Initialize the viscosity of the gas phase.
Definition: DryHumidGasPvt.hpp:339
Evaluation saturationPressure(unsigned regionIdx, const Evaluation &, const Evaluation &Rw) const
Returns the saturation pressure of the gas phase [Pa] depending on its mass fraction of the water com...
Definition: DryHumidGasPvt.hpp:552
DryHumidGasPvt(const std::vector< Scalar > &gasReferenceDensity, const std::vector< Scalar > &waterReferenceDensity, 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 > &saturatedWaterVaporizationFactorTable, const std::vector< TabulatedOneDFunction > &saturationPressure, Scalar vapPar1)
Definition: DryHumidGasPvt.hpp:65
bool operator==(const DryHumidGasPvt< Scalar > &data) const
Definition: DryHumidGasPvt.hpp:656
void setInverseGasFormationVolumeFactor(unsigned regionIdx, const TabulatedTwoDFunction &invBg)
Initialize the function for the gas formation volume factor.
Definition: DryHumidGasPvt.hpp:331
const std::vector< TabulatedTwoDFunction > & saturatedWaterVaporizationSaltFactorTable() const
Definition: DryHumidGasPvt.hpp:644
void initEnd()
Finish initializing the gas phase PVT properties.
Definition: DryHumidGasPvt.hpp:384
unsigned numRegions() const
Return the number of PVT regions which are considered by this PVT-object.
Definition: DryHumidGasPvt.hpp:432
void setReferenceDensities(unsigned regionIdx, Scalar, Scalar rhoRefGas, Scalar rhoRefWater)
Initialize the reference densities of all fluids for a given PVT region.
Definition: DryHumidGasPvt.hpp:302
const std::vector< TabulatedTwoDFunction > & inverseGasB() const
Definition: DryHumidGasPvt.hpp:620
const std::vector< TabulatedTwoDFunction > & gasMu() const
Definition: DryHumidGasPvt.hpp:628
Evaluation saturatedWaterVaporizationFactor(unsigned regionIdx, const Evaluation &, const Evaluation &pressure) const
Returns the water vaporization factor [m^3/m^3] of the water phase.
Definition: DryHumidGasPvt.hpp:501
const std::vector< TabulatedOneDFunction > & inverseSaturatedGasBMu() const
Definition: DryHumidGasPvt.hpp:636
void setSaturatedGasWaterVaporizationFactor(unsigned regionIdx, const SamplingPoints &samplePoints)
Initialize the function for the oil vaporization factor .
Definition: DryHumidGasPvt.hpp:316
const std::vector< TabulatedTwoDFunction > & inverseGasBMu() const
Definition: DryHumidGasPvt.hpp:632
Evaluation viscosity(unsigned regionIdx, const Evaluation &, const Evaluation &pressure, const Evaluation &, const Evaluation &Rvw) const
Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
Definition: DryHumidGasPvt.hpp:451
void setSaturatedGasViscosity(unsigned regionIdx, const SamplingPoints &samplePoints)
Initialize the phase viscosity for oil saturated gas.
Definition: DryHumidGasPvt.hpp:349
const Scalar waterReferenceDensity(unsigned regionIdx) const
Definition: DryHumidGasPvt.hpp:617
Evaluation internalEnergy(unsigned, const Evaluation &, const Evaluation &, const Evaluation &) const
Returns the specific enthalpy [J/kg] of gas given a set of parameters.
Definition: DryHumidGasPvt.hpp:439
Evaluation saturatedOilVaporizationFactor(unsigned, const Evaluation &, const Evaluation &) const
Returns the oil vaporization factor [m^3/m^3] of the oil phase.
Definition: DryHumidGasPvt.hpp:539
Evaluation diffusionCoefficient(const Evaluation &, const Evaluation &, unsigned) const
Definition: DryHumidGasPvt.hpp:607
Evaluation saturatedWaterVaporizationFactor(unsigned regionIdx, const Evaluation &, const Evaluation &pressure, const Evaluation &saltConcentration) const
Returns the water vaporization factor [m^3/m^3] of the water phase.
Definition: DryHumidGasPvt.hpp:512
const std::vector< TabulatedOneDFunction > & saturatedWaterVaporizationFactorTable() const
Definition: DryHumidGasPvt.hpp:640
void setNumRegions(size_t numRegions)
Definition: DryHumidGasPvt.hpp:286
const std::vector< TabulatedOneDFunction > & saturationPressure() const
Definition: DryHumidGasPvt.hpp:648
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: DryHumidGasPvt.hpp:467
Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &, const Evaluation &pressure) const
Returns the formation volume factor [-] of water saturated gas at a given pressure.
Definition: DryHumidGasPvt.hpp:492
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
Definition: Air_Mesitylene.hpp:34
auto scalarValue(const Evaluation &val) -> decltype(MathToolbox< Evaluation >::scalarValue(val))
Definition: MathToolbox.hpp:335
Evaluation abs(const Evaluation &value)
Definition: MathToolbox.hpp:350