27 #ifndef OPM_LIVE_OIL_PVT_HPP
28 #define OPM_LIVE_OIL_PVT_HPP
37 #include <opm/parser/eclipse/Deck/Deck.hpp>
38 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
39 #include <opm/parser/eclipse/EclipseState/Tables/SimpleTable.hpp>
40 #include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
48 template <
class Scalar>
53 typedef std::vector<std::pair<Scalar, Scalar> > SamplingPoints;
65 void initFromDeck(
const Deck& deck,
const EclipseState& eclState)
67 const auto& pvtoTables = eclState.getTableManager().getPvtoTables();
68 const auto& densityKeyword = deck.getKeyword(
"DENSITY");
70 assert(pvtoTables.size() == densityKeyword.size());
73 setNumRegions(numRegions);
75 for (
unsigned regionIdx = 0; regionIdx <
numRegions; ++ regionIdx) {
76 Scalar rhoRefO = densityKeyword.getRecord(regionIdx).getItem(
"OIL").getSIDouble(0);
77 Scalar rhoRefG = densityKeyword.getRecord(regionIdx).getItem(
"GAS").getSIDouble(0);
78 Scalar rhoRefW = densityKeyword.getRecord(regionIdx).getItem(
"WATER").getSIDouble(0);
84 for (
unsigned regionIdx = 0; regionIdx <
numRegions; ++ regionIdx) {
85 const auto& pvtoTable = pvtoTables[regionIdx];
87 const auto& saturatedTable = pvtoTable.getSaturatedTable();
88 assert(saturatedTable.numRows() > 1);
90 auto& oilMu = oilMuTable_[regionIdx];
91 auto& satOilMu = saturatedOilMuTable_[regionIdx];
92 auto& invOilB = inverseOilBTable_[regionIdx];
93 auto& invSatOilB = inverseSaturatedOilBTable_[regionIdx];
94 auto& gasDissolutionFac = saturatedGasDissolutionFactorTable_[regionIdx];
95 std::vector<Scalar> invSatOilBArray;
96 std::vector<Scalar> satOilMuArray;
99 for (
unsigned outerIdx = 0; outerIdx < saturatedTable.numRows(); ++ outerIdx) {
100 Scalar Rs = saturatedTable.get(
"RS", outerIdx);
101 Scalar BoSat = saturatedTable.get(
"BO", outerIdx);
102 Scalar muoSat = saturatedTable.get(
"MU", outerIdx);
104 satOilMuArray.push_back(muoSat);
105 invSatOilBArray.push_back(1.0/BoSat);
107 invOilB.appendXPos(Rs);
108 oilMu.appendXPos(Rs);
110 assert(invOilB.numX() == outerIdx + 1);
111 assert(oilMu.numX() == outerIdx + 1);
113 const auto& underSaturatedTable = pvtoTable.getUnderSaturatedTable(outerIdx);
114 size_t numRows = underSaturatedTable.numRows();
115 for (
unsigned innerIdx = 0; innerIdx < numRows; ++ innerIdx) {
116 Scalar po = underSaturatedTable.get(
"P", innerIdx);
117 Scalar Bo = underSaturatedTable.get(
"BO", innerIdx);
118 Scalar muo = underSaturatedTable.get(
"MU", innerIdx);
120 invOilB.appendSamplePoint(outerIdx, po, 1.0/Bo);
121 oilMu.appendSamplePoint(outerIdx, po, muo);
128 const auto& tmpPressureColumn = saturatedTable.getColumn(
"P");
129 const auto& tmpGasSolubilityColumn = saturatedTable.getColumn(
"RS");
131 invSatOilB.setXYContainers(tmpPressureColumn, invSatOilBArray);
132 satOilMu.setXYContainers(tmpPressureColumn, satOilMuArray);
133 gasDissolutionFac.setXYContainers(tmpPressureColumn, tmpGasSolubilityColumn);
136 updateSaturationPressure_(regionIdx);
138 for (
unsigned xIdx = 0; xIdx < invOilB.numX(); ++xIdx) {
140 assert(invOilB.numY(xIdx) > 0);
144 if (invOilB.numY(xIdx) > 1)
150 size_t masterTableIdx = xIdx + 1;
151 for (; masterTableIdx < saturatedTable.numRows(); ++masterTableIdx)
153 if (pvtoTable.getUnderSaturatedTable(masterTableIdx).numRows() > 1)
157 if (masterTableIdx >= saturatedTable.numRows())
158 OPM_THROW(std::runtime_error,
159 "PVTO tables are invalid: The last table must exhibit at least one "
160 "entry for undersaturated oil!");
163 extendPvtoTable_(regionIdx,
165 pvtoTable.getUnderSaturatedTable(xIdx),
166 pvtoTable.getUnderSaturatedTable(masterTableIdx));
171 if (deck.hasKeyword(
"VAPPARS")) {
172 const auto& vapParsKeyword = deck.getKeyword(
"VAPPARS");
173 vapPar2_ = vapParsKeyword.getRecord(0).getItem(
"OIL_DENSITY_PROPENSITY").template get<double>(0);
180 void extendPvtoTable_(
unsigned regionIdx,
182 const SimpleTable& curTable,
183 const SimpleTable& masterTable)
185 std::vector<double> pressuresArray = curTable.getColumn(
"P").vectorCopy();
186 std::vector<double> oilBArray = curTable.getColumn(
"BO").vectorCopy();
187 std::vector<double> oilMuArray = curTable.getColumn(
"MU").vectorCopy();
189 auto& invOilB = inverseOilBTable_[regionIdx];
190 auto& oilMu = oilMuTable_[regionIdx];
192 for (
unsigned newRowIdx = 1; newRowIdx < masterTable.numRows(); ++ newRowIdx) {
193 const auto& pressureColumn = masterTable.getColumn(
"P");
194 const auto& BOColumn = masterTable.getColumn(
"BO");
195 const auto& viscosityColumn = masterTable.getColumn(
"MU");
198 Scalar diffPo = pressureColumn[newRowIdx] - pressureColumn[newRowIdx - 1];
199 Scalar newPo = pressuresArray.back() + diffPo;
202 Scalar B1 = BOColumn[newRowIdx];
203 Scalar B2 = BOColumn[newRowIdx - 1];
204 Scalar x = (B1 - B2)/( (B1 + B2)/2.0 );
208 Scalar newBo = oilBArray.back()*(1.0 + x/2.0)/(1.0 - x/2.0);
211 Scalar mu1 = viscosityColumn[newRowIdx];
212 Scalar mu2 = viscosityColumn[newRowIdx - 1];
213 Scalar xMu = (mu1 - mu2)/( (mu1 + mu2)/2.0 );
217 Scalar newMuo = oilMuArray.back()*(1.0 + xMu/2)/(1.0 - xMu/2.0);
221 pressuresArray.push_back(newPo);
222 oilBArray.push_back(newBo);
223 oilMuArray.push_back(newMuo);
226 invOilB.appendSamplePoint(xIdx, newPo, 1.0/newBo);
227 oilMu.appendSamplePoint(xIdx, newPo, newMuo);
232 #endif // HAVE_OPM_PARSER
236 oilReferenceDensity_.resize(numRegions);
237 gasReferenceDensity_.resize(numRegions);
238 inverseOilBTable_.resize(numRegions);
239 inverseOilBMuTable_.resize(numRegions);
240 inverseSaturatedOilBTable_.resize(numRegions);
241 inverseSaturatedOilBMuTable_.resize(numRegions);
242 oilMuTable_.resize(numRegions);
243 saturatedOilMuTable_.resize(numRegions);
244 saturatedGasDissolutionFactorTable_.resize(numRegions);
245 saturationPressure_.resize(numRegions);
256 oilReferenceDensity_[regionIdx] = rhoRefOil;
257 gasReferenceDensity_[regionIdx] = rhoRefGas;
266 { saturatedGasDissolutionFactorTable_[regionIdx].setContainerOfTuples(samplePoints); }
279 Scalar T = 273.15 + 15.56;
280 auto& invOilB = inverseOilBTable_[regionIdx];
282 updateSaturationPressure_(regionIdx);
285 for (
size_t pIdx = 0; pIdx < samplePoints.size(); ++pIdx) {
286 Scalar p1 = std::get<0>(samplePoints[pIdx]);
287 Scalar p2 = p1 * 2.0;
289 Scalar Bo1 = std::get<1>(samplePoints[pIdx]);
290 Scalar drhoo_dp = (1.1200 - 1.1189)/((5000 - 4000)*6894.76);
291 Scalar Bo2 = Bo1/(1.0 + (p2 - p1)*drhoo_dp);
295 invOilB.appendXPos(Rs);
296 invOilB.appendSamplePoint(pIdx, p1, 1.0/Bo1);
297 invOilB.appendSamplePoint(pIdx, p2, 1.0/Bo2);
314 { inverseOilBTable_[regionIdx] = invBo; }
322 { oilMuTable_[regionIdx] = muo; }
333 Scalar T = 273.15 + 15.56;
336 saturatedOilMuTable_[regionIdx].setContainerOfTuples(samplePoints);
340 for (
size_t pIdx = 0; pIdx < samplePoints.size(); ++pIdx) {
341 Scalar p1 = std::get<0>(samplePoints[pIdx]);
342 Scalar p2 = p1 * 2.0;
345 Scalar mu1 = std::get<1>(samplePoints[pIdx]);
350 oilMuTable_[regionIdx].appendXPos(Rs);
351 oilMuTable_[regionIdx].appendSamplePoint(pIdx, p1, mu1);
352 oilMuTable_[regionIdx].appendSamplePoint(pIdx, p2, mu2);
362 size_t numRegions = oilMuTable_.size();
363 for (
unsigned regionIdx = 0; regionIdx <
numRegions; ++ regionIdx) {
366 const auto& oilMu = oilMuTable_[regionIdx];
367 const auto& satOilMu = saturatedOilMuTable_[regionIdx];
368 const auto& invOilB = inverseOilBTable_[regionIdx];
369 assert(oilMu.numX() == invOilB.numX());
371 auto& invOilBMu = inverseOilBMuTable_[regionIdx];
372 auto& invSatOilB = inverseSaturatedOilBTable_[regionIdx];
373 auto& invSatOilBMu = inverseSaturatedOilBMuTable_[regionIdx];
375 std::vector<Scalar> satPressuresArray;
376 std::vector<Scalar> invSatOilBArray;
377 std::vector<Scalar> invSatOilBMuArray;
378 for (
unsigned rsIdx = 0; rsIdx < oilMu.numX(); ++rsIdx) {
379 invOilBMu.appendXPos(oilMu.xAt(rsIdx));
381 assert(oilMu.numY(rsIdx) == invOilB.numY(rsIdx));
383 size_t numPressures = oilMu.numY(rsIdx);
384 for (
unsigned pIdx = 0; pIdx < numPressures; ++pIdx)
385 invOilBMu.appendSamplePoint(rsIdx,
386 oilMu.yAt(rsIdx, pIdx),
387 invOilB.valueAt(rsIdx, pIdx)
388 / oilMu.valueAt(rsIdx, pIdx));
393 satPressuresArray.push_back(oilMu.yAt(rsIdx, 0));
394 invSatOilBArray.push_back(invOilB.valueAt(rsIdx, 0));
395 invSatOilBMuArray.push_back(invSatOilBArray.back()/satOilMu.valueAt(rsIdx));
398 invSatOilB.setXYContainers(satPressuresArray, invSatOilBArray);
399 invSatOilBMu.setXYContainers(satPressuresArray, invSatOilBMuArray);
401 updateSaturationPressure_(regionIdx);
409 {
return inverseOilBMuTable_.size(); }
414 template <
class Evaluation>
417 const Evaluation& pressure,
418 const Evaluation& Rs)
const
421 const Evaluation& invBo = inverseOilBTable_[regionIdx].eval(Rs, pressure,
true);
422 const Evaluation& invMuoBo = inverseOilBMuTable_[regionIdx].eval(Rs, pressure,
true);
424 return invBo/invMuoBo;
430 template <
class Evaluation>
433 const Evaluation& pressure)
const
436 const Evaluation& invBo = inverseSaturatedOilBTable_[regionIdx].eval(pressure,
true);
437 const Evaluation& invMuoBo = inverseSaturatedOilBMuTable_[regionIdx].eval(pressure,
true);
439 return invBo/invMuoBo;
445 template <
class Evaluation>
448 const Evaluation& pressure,
449 const Evaluation& Rs)
const
452 return inverseOilBTable_[regionIdx].eval(Rs, pressure,
true);
458 template <
class Evaluation>
461 const Evaluation& pressure)
const
464 return inverseSaturatedOilBTable_[regionIdx].eval(pressure,
true);
470 template <
class Evaluation>
473 const Evaluation& pressure)
const
474 {
return saturatedGasDissolutionFactorTable_[regionIdx].eval(pressure,
true); }
483 template <
class Evaluation>
486 const Evaluation& pressure,
487 const Evaluation& oilSaturation,
488 Scalar maxOilSaturation)
const
491 saturatedGasDissolutionFactorTable_[regionIdx].eval(pressure,
true);
495 maxOilSaturation = std::min(maxOilSaturation, Scalar(1.0));
496 if (vapPar2_ > 0.0 && maxOilSaturation > 0.01 && oilSaturation < maxOilSaturation) {
497 static const Scalar eps = 0.001;
498 const Evaluation& So = Opm::max(oilSaturation, eps);
499 tmp *= Opm::max(1e-3, Opm::pow(So/maxOilSaturation, vapPar2_));
511 template <
class Evaluation>
513 const Evaluation& temperature OPM_UNUSED,
514 const Evaluation& Rs)
const
518 const auto& RsTable = saturatedGasDissolutionFactorTable_[regionIdx];
519 const Scalar eps = std::numeric_limits<typename Toolbox::Scalar>::epsilon()*1e6;
522 Evaluation pSat = saturationPressure_[regionIdx].eval(Rs,
true);
527 bool onProbation =
false;
528 for (
int i = 0; i < 20; ++i) {
529 const Evaluation& f = RsTable.eval(pSat,
true) - Rs;
530 const Evaluation& fPrime = RsTable.evalDerivative(pSat,
true);
534 if (std::abs(Opm::scalarValue(fPrime)) < 1.0e-30) {
538 const Evaluation& delta = f/fPrime;
552 if (std::abs(Opm::scalarValue(delta)) < std::abs(Opm::scalarValue(pSat))*eps)
556 std::stringstream errlog;
557 errlog <<
"Finding saturation pressure did not converge:"
558 <<
" pSat = " << pSat
560 OpmLog::debug(
"Live oil saturation pressure", errlog.str());
561 OPM_THROW_NOLOG(NumericalProblem, errlog.str());
565 void updateSaturationPressure_(
unsigned regionIdx)
567 typedef std::pair<Scalar, Scalar> Pair;
568 const auto& gasDissolutionFac = saturatedGasDissolutionFactorTable_[regionIdx];
572 size_t n = gasDissolutionFac.numSamples();
573 Scalar delta = (gasDissolutionFac.xMax() - gasDissolutionFac.xMin())/Scalar(n + 1);
575 SamplingPoints pSatSamplePoints;
577 for (
size_t i=0; i <= n; ++ i) {
578 Scalar pSat = gasDissolutionFac.xMin() + Scalar(i)*delta;
584 pSatSamplePoints.push_back(val);
588 auto x_coord_comparator = [](
const Pair& a,
const Pair& b) {
return a.first == b.first; };
589 auto last = std::unique(pSatSamplePoints.begin(), pSatSamplePoints.end(), x_coord_comparator);
590 pSatSamplePoints.erase(last, pSatSamplePoints.end());
592 saturationPressure_[regionIdx].setContainerOfTuples(pSatSamplePoints);
595 std::vector<Scalar> gasReferenceDensity_;
596 std::vector<Scalar> oilReferenceDensity_;
597 std::vector<TabulatedTwoDFunction> inverseOilBTable_;
598 std::vector<TabulatedTwoDFunction> oilMuTable_;
599 std::vector<TabulatedTwoDFunction> inverseOilBMuTable_;
600 std::vector<TabulatedOneDFunction> saturatedOilMuTable_;
601 std::vector<TabulatedOneDFunction> inverseSaturatedOilBTable_;
602 std::vector<TabulatedOneDFunction> inverseSaturatedOilBMuTable_;
603 std::vector<TabulatedOneDFunction> saturatedGasDissolutionFactorTable_;
604 std::vector<TabulatedOneDFunction> saturationPressure_;
Evaluation saturationPressure(unsigned regionIdx, const Evaluation &temperature OPM_UNUSED, const Evaluation &Rs) const
Returns the saturation pressure of the oil phase [Pa] depending on its mass fraction of the gas compo...
Definition: LiveOilPvt.hpp:512
void setOilViscosity(unsigned regionIdx, const TabulatedTwoDFunction &muo)
Initialize the viscosity of the oil phase.
Definition: LiveOilPvt.hpp:321
Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &, const Evaluation &pressure) const
Returns the formation volume factor [-] of the fluid phase.
Definition: LiveOilPvt.hpp:459
void initEnd()
Finish initializing the oil phase PVT properties.
Definition: LiveOilPvt.hpp:359
Evaluation inverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &, const Evaluation &pressure, const Evaluation &Rs) const
Returns the formation volume factor [-] of the fluid phase.
Definition: LiveOilPvt.hpp:446
This class represents the Pressure-Volume-Temperature relations of the oil phas with dissolved gas...
Definition: LiveOilPvt.hpp:49
Evaluation saturatedViscosity(unsigned regionIdx, const Evaluation &, const Evaluation &pressure) const
Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
Definition: LiveOilPvt.hpp:431
void setReferenceDensities(unsigned regionIdx, Scalar rhoRefOil, Scalar rhoRefGas, Scalar)
Initialize the reference densities of all fluids for a given PVT region.
Definition: LiveOilPvt.hpp:251
void setSaturatedOilViscosity(unsigned regionIdx, const SamplingPoints &samplePoints)
Initialize the phase viscosity for gas saturated oil.
Definition: LiveOilPvt.hpp:331
void setSaturatedOilFormationVolumeFactor(unsigned regionIdx, const SamplingPoints &samplePoints)
Initialize the function for the oil formation volume factor.
Definition: LiveOilPvt.hpp:277
void setSaturatedOilGasDissolutionFactor(unsigned regionIdx, const SamplingPoints &samplePoints)
Initialize the function for the gas dissolution factor .
Definition: LiveOilPvt.hpp:265
Evaluation saturatedGasDissolutionFactor(unsigned regionIdx, const Evaluation &, const Evaluation &pressure, const Evaluation &oilSaturation, Scalar maxOilSaturation) const
Returns the gas dissolution factor [m^3/m^3] of the oil phase.
Definition: LiveOilPvt.hpp:484
void setInverseOilFormationVolumeFactor(unsigned regionIdx, const TabulatedTwoDFunction &invBo)
Initialize the function for the oil formation volume factor.
Definition: LiveOilPvt.hpp:313
Evaluation viscosity(unsigned regionIdx, const Evaluation &, const Evaluation &pressure, const Evaluation &Rs) const
Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
Definition: LiveOilPvt.hpp:415
Implements a linearly interpolated scalar function that depends on one variable.
A central place for various physical constants occuring in some equations.
This file provides a wrapper around the "final" C++-2011 statement.
Implements a linearly interpolated scalar function that depends on one variable.
Definition: Tabulated1DFunction.hpp:47
unsigned numRegions() const
Return the number of PVT regions which are considered by this PVT-object.
Definition: LiveOilPvt.hpp:408
Evaluation saturatedGasDissolutionFactor(unsigned regionIdx, const Evaluation &, const Evaluation &pressure) const
Returns the gas dissolution factor [m^3/m^3] of the oil phase.
Definition: LiveOilPvt.hpp:471