My Project
LiveOilPvt.hpp
Go to the documentation of this file.
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 /*
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 2 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 
19  Consult the COPYING file in the top-level source directory of this
20  module for the precise wording of the license and the list of
21  copyright holders.
22 */
27 #ifndef OPM_LIVE_OIL_PVT_HPP
28 #define OPM_LIVE_OIL_PVT_HPP
29 
33 
34 #if HAVE_ECL_INPUT
35 #include <opm/input/eclipse/EclipseState/EclipseState.hpp>
36 #include <opm/input/eclipse/Schedule/OilVaporizationProperties.hpp>
37 #include <opm/input/eclipse/Schedule/Schedule.hpp>
38 #include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp>
39 #endif
40 
41 #if HAVE_OPM_COMMON
42 #include <opm/common/OpmLog/OpmLog.hpp>
43 #endif
44 
45 namespace Opm {
50 template <class Scalar>
52 {
53  using SamplingPoints = std::vector<std::pair<Scalar, Scalar>>;
54 
55 public:
58 
59  LiveOilPvt()
60  {
61  vapPar2_ = 0.0;
62  }
63 
64  LiveOilPvt(const std::vector<Scalar>& gasReferenceDensity,
65  const std::vector<Scalar>& oilReferenceDensity,
66  const std::vector<TabulatedTwoDFunction>& inverseOilBTable,
67  const std::vector<TabulatedTwoDFunction>& oilMuTable,
68  const std::vector<TabulatedTwoDFunction>& inverseOilBMuTable,
69  const std::vector<TabulatedOneDFunction>& saturatedOilMuTable,
70  const std::vector<TabulatedOneDFunction>& inverseSaturatedOilBTable,
71  const std::vector<TabulatedOneDFunction>& inverseSaturatedOilBMuTable,
72  const std::vector<TabulatedOneDFunction>& saturatedGasDissolutionFactorTable,
73  const std::vector<TabulatedOneDFunction>& saturationPressure,
74  Scalar vapPar2)
75  : gasReferenceDensity_(gasReferenceDensity)
76  , oilReferenceDensity_(oilReferenceDensity)
77  , inverseOilBTable_(inverseOilBTable)
78  , oilMuTable_(oilMuTable)
79  , inverseOilBMuTable_(inverseOilBMuTable)
80  , saturatedOilMuTable_(saturatedOilMuTable)
81  , inverseSaturatedOilBTable_(inverseSaturatedOilBTable)
82  , inverseSaturatedOilBMuTable_(inverseSaturatedOilBMuTable)
83  , saturatedGasDissolutionFactorTable_(saturatedGasDissolutionFactorTable)
84  , saturationPressure_(saturationPressure)
85  , vapPar2_(vapPar2)
86  { }
87 
88 #if HAVE_ECL_INPUT
92  void initFromState(const EclipseState& eclState, const Schedule& schedule)
93  {
94  const auto& pvtoTables = eclState.getTableManager().getPvtoTables();
95  const auto& densityTable = eclState.getTableManager().getDensityTable();
96 
97  assert(pvtoTables.size() == densityTable.size());
98 
99  size_t numRegions = pvtoTables.size();
100  setNumRegions(numRegions);
101 
102  for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
103  Scalar rhoRefO = densityTable[regionIdx].oil;
104  Scalar rhoRefG = densityTable[regionIdx].gas;
105  Scalar rhoRefW = densityTable[regionIdx].water;
106 
107  setReferenceDensities(regionIdx, rhoRefO, rhoRefG, rhoRefW);
108  }
109 
110  // initialize the internal table objects
111  for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
112  const auto& pvtoTable = pvtoTables[regionIdx];
113 
114  const auto& saturatedTable = pvtoTable.getSaturatedTable();
115  assert(saturatedTable.numRows() > 1);
116 
117  auto& oilMu = oilMuTable_[regionIdx];
118  auto& satOilMu = saturatedOilMuTable_[regionIdx];
119  auto& invOilB = inverseOilBTable_[regionIdx];
120  auto& invSatOilB = inverseSaturatedOilBTable_[regionIdx];
121  auto& gasDissolutionFac = saturatedGasDissolutionFactorTable_[regionIdx];
122  std::vector<Scalar> invSatOilBArray;
123  std::vector<Scalar> satOilMuArray;
124 
125  // extract the table for the gas dissolution and the oil formation volume factors
126  for (unsigned outerIdx = 0; outerIdx < saturatedTable.numRows(); ++ outerIdx) {
127  Scalar Rs = saturatedTable.get("RS", outerIdx);
128  Scalar BoSat = saturatedTable.get("BO", outerIdx);
129  Scalar muoSat = saturatedTable.get("MU", outerIdx);
130 
131  satOilMuArray.push_back(muoSat);
132  invSatOilBArray.push_back(1.0/BoSat);
133 
134  invOilB.appendXPos(Rs);
135  oilMu.appendXPos(Rs);
136 
137  assert(invOilB.numX() == outerIdx + 1);
138  assert(oilMu.numX() == outerIdx + 1);
139 
140  const auto& underSaturatedTable = pvtoTable.getUnderSaturatedTable(outerIdx);
141  size_t numRows = underSaturatedTable.numRows();
142  for (unsigned innerIdx = 0; innerIdx < numRows; ++ innerIdx) {
143  Scalar po = underSaturatedTable.get("P", innerIdx);
144  Scalar Bo = underSaturatedTable.get("BO", innerIdx);
145  Scalar muo = underSaturatedTable.get("MU", innerIdx);
146 
147  invOilB.appendSamplePoint(outerIdx, po, 1.0/Bo);
148  oilMu.appendSamplePoint(outerIdx, po, muo);
149  }
150  }
151 
152  // update the tables for the formation volume factor and for the gas
153  // dissolution factor of saturated oil
154  {
155  const auto& tmpPressureColumn = saturatedTable.getColumn("P");
156  const auto& tmpGasSolubilityColumn = saturatedTable.getColumn("RS");
157 
158  invSatOilB.setXYContainers(tmpPressureColumn, invSatOilBArray);
159  satOilMu.setXYContainers(tmpPressureColumn, satOilMuArray);
160  gasDissolutionFac.setXYContainers(tmpPressureColumn, tmpGasSolubilityColumn);
161  }
162 
163  updateSaturationPressure_(regionIdx);
164  // make sure to have at least two sample points per Rs value
165  for (unsigned xIdx = 0; xIdx < invOilB.numX(); ++xIdx) {
166  // a single sample point is definitely needed
167  assert(invOilB.numY(xIdx) > 0);
168 
169  // everything is fine if the current table has two or more sampling points
170  // for a given mole fraction
171  if (invOilB.numY(xIdx) > 1)
172  continue;
173 
174  // find the master table which will be used as a template to extend the
175  // current line. We define master table as the first table which has values
176  // for undersaturated oil...
177  size_t masterTableIdx = xIdx + 1;
178  for (; masterTableIdx < saturatedTable.numRows(); ++masterTableIdx)
179  {
180  if (pvtoTable.getUnderSaturatedTable(masterTableIdx).numRows() > 1)
181  break;
182  }
183 
184  if (masterTableIdx >= saturatedTable.numRows())
185  throw std::runtime_error("PVTO tables are invalid: The last table must exhibit at least one "
186  "entry for undersaturated oil!");
187 
188  // extend the current table using the master table.
189  extendPvtoTable_(regionIdx,
190  xIdx,
191  pvtoTable.getUnderSaturatedTable(xIdx),
192  pvtoTable.getUnderSaturatedTable(masterTableIdx));
193  }
194  }
195 
196  vapPar2_ = 0.0;
197  const auto& oilVap = schedule[0].oilvap();
198  if (oilVap.getType() == OilVaporizationProperties::OilVaporization::VAPPARS) {
199  vapPar2_ = oilVap.vap2();
200  }
201 
202  initEnd();
203  }
204 
205 private:
206  void extendPvtoTable_(unsigned regionIdx,
207  unsigned xIdx,
208  const SimpleTable& curTable,
209  const SimpleTable& masterTable)
210  {
211  std::vector<double> pressuresArray = curTable.getColumn("P").vectorCopy();
212  std::vector<double> oilBArray = curTable.getColumn("BO").vectorCopy();
213  std::vector<double> oilMuArray = curTable.getColumn("MU").vectorCopy();
214 
215  auto& invOilB = inverseOilBTable_[regionIdx];
216  auto& oilMu = oilMuTable_[regionIdx];
217 
218  for (unsigned newRowIdx = 1; newRowIdx < masterTable.numRows(); ++ newRowIdx) {
219  const auto& pressureColumn = masterTable.getColumn("P");
220  const auto& BOColumn = masterTable.getColumn("BO");
221  const auto& viscosityColumn = masterTable.getColumn("MU");
222 
223  // compute the oil pressure for the new entry
224  Scalar diffPo = pressureColumn[newRowIdx] - pressureColumn[newRowIdx - 1];
225  Scalar newPo = pressuresArray.back() + diffPo;
226 
227  // calculate the compressibility of the master table
228  Scalar B1 = BOColumn[newRowIdx];
229  Scalar B2 = BOColumn[newRowIdx - 1];
230  Scalar x = (B1 - B2)/( (B1 + B2)/2.0 );
231 
232  // calculate the oil formation volume factor which exhibits the same
233  // compressibility for the new pressure
234  Scalar newBo = oilBArray.back()*(1.0 + x/2.0)/(1.0 - x/2.0);
235 
236  // calculate the "viscosibility" of the master table
237  Scalar mu1 = viscosityColumn[newRowIdx];
238  Scalar mu2 = viscosityColumn[newRowIdx - 1];
239  Scalar xMu = (mu1 - mu2)/( (mu1 + mu2)/2.0 );
240 
241  // calculate the oil formation volume factor which exhibits the same
242  // compressibility for the new pressure
243  Scalar newMuo = oilMuArray.back()*(1.0 + xMu/2)/(1.0 - xMu/2.0);
244 
245  // append the new values to the arrays which we use to compute the additional
246  // values ...
247  pressuresArray.push_back(newPo);
248  oilBArray.push_back(newBo);
249  oilMuArray.push_back(newMuo);
250 
251  // ... and register them with the internal table objects
252  invOilB.appendSamplePoint(xIdx, newPo, 1.0/newBo);
253  oilMu.appendSamplePoint(xIdx, newPo, newMuo);
254  }
255  }
256 
257 public:
258 #endif // HAVE_ECL_INPUT
259 
260  void setNumRegions(size_t numRegions)
261  {
262  oilReferenceDensity_.resize(numRegions);
263  gasReferenceDensity_.resize(numRegions);
264  inverseOilBTable_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::LeftExtreme});
265  inverseOilBMuTable_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::LeftExtreme});
266  inverseSaturatedOilBTable_.resize(numRegions);
267  inverseSaturatedOilBMuTable_.resize(numRegions);
268  oilMuTable_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::LeftExtreme});
269  saturatedOilMuTable_.resize(numRegions);
270  saturatedGasDissolutionFactorTable_.resize(numRegions);
271  saturationPressure_.resize(numRegions);
272  }
273 
277  void setReferenceDensities(unsigned regionIdx,
278  Scalar rhoRefOil,
279  Scalar rhoRefGas,
280  Scalar /*rhoRefWater*/)
281  {
282  oilReferenceDensity_[regionIdx] = rhoRefOil;
283  gasReferenceDensity_[regionIdx] = rhoRefGas;
284  }
285 
291  void setSaturatedOilGasDissolutionFactor(unsigned regionIdx, const SamplingPoints& samplePoints)
292  { saturatedGasDissolutionFactorTable_[regionIdx].setContainerOfTuples(samplePoints); }
293 
303  void setSaturatedOilFormationVolumeFactor(unsigned regionIdx, const SamplingPoints& samplePoints)
304  {
305  constexpr const Scalar T = 273.15 + 15.56; // [K]
306  auto& invOilB = inverseOilBTable_[regionIdx];
307 
308  updateSaturationPressure_(regionIdx);
309 
310  // calculate a table of estimated densities of undersatured gas
311  for (size_t pIdx = 0; pIdx < samplePoints.size(); ++pIdx) {
312  Scalar p1 = std::get<0>(samplePoints[pIdx]);
313  Scalar p2 = p1 * 2.0;
314 
315  Scalar Bo1 = std::get<1>(samplePoints[pIdx]);
316  Scalar drhoo_dp = (1.1200 - 1.1189)/((5000 - 4000)*6894.76);
317  Scalar Bo2 = Bo1/(1.0 + (p2 - p1)*drhoo_dp);
318 
319  Scalar Rs = saturatedGasDissolutionFactor(regionIdx, T, p1);
320 
321  invOilB.appendXPos(Rs);
322  invOilB.appendSamplePoint(pIdx, p1, 1.0/Bo1);
323  invOilB.appendSamplePoint(pIdx, p2, 1.0/Bo2);
324  }
325  }
326 
339  void setInverseOilFormationVolumeFactor(unsigned regionIdx, const TabulatedTwoDFunction& invBo)
340  { inverseOilBTable_[regionIdx] = invBo; }
341 
347  void setOilViscosity(unsigned regionIdx, const TabulatedTwoDFunction& muo)
348  { oilMuTable_[regionIdx] = muo; }
349 
357  void setSaturatedOilViscosity(unsigned regionIdx, const SamplingPoints& samplePoints)
358  {
359  constexpr const Scalar T = 273.15 + 15.56; // [K]
360 
361  // update the table for the saturated oil
362  saturatedOilMuTable_[regionIdx].setContainerOfTuples(samplePoints);
363 
364  // calculate a table of estimated viscosities depending on pressure and gas mass
365  // fraction for untersaturated oil to make the other code happy
366  for (size_t pIdx = 0; pIdx < samplePoints.size(); ++pIdx) {
367  Scalar p1 = std::get<0>(samplePoints[pIdx]);
368  Scalar p2 = p1 * 2.0;
369 
370  // no pressure dependence of the viscosity
371  Scalar mu1 = std::get<1>(samplePoints[pIdx]);
372  Scalar mu2 = mu1;
373 
374  Scalar Rs = saturatedGasDissolutionFactor(regionIdx, T, p1);
375 
376  oilMuTable_[regionIdx].appendXPos(Rs);
377  oilMuTable_[regionIdx].appendSamplePoint(pIdx, p1, mu1);
378  oilMuTable_[regionIdx].appendSamplePoint(pIdx, p2, mu2);
379  }
380  }
381 
385  void initEnd()
386  {
387  // calculate the final 2D functions which are used for interpolation.
388  size_t numRegions = oilMuTable_.size();
389  for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
390  // calculate the table which stores the inverse of the product of the oil
391  // formation volume factor and the oil viscosity
392  const auto& oilMu = oilMuTable_[regionIdx];
393  const auto& satOilMu = saturatedOilMuTable_[regionIdx];
394  const auto& invOilB = inverseOilBTable_[regionIdx];
395  assert(oilMu.numX() == invOilB.numX());
396 
397  auto& invOilBMu = inverseOilBMuTable_[regionIdx];
398  auto& invSatOilB = inverseSaturatedOilBTable_[regionIdx];
399  auto& invSatOilBMu = inverseSaturatedOilBMuTable_[regionIdx];
400 
401  std::vector<Scalar> satPressuresArray;
402  std::vector<Scalar> invSatOilBArray;
403  std::vector<Scalar> invSatOilBMuArray;
404  for (unsigned rsIdx = 0; rsIdx < oilMu.numX(); ++rsIdx) {
405  invOilBMu.appendXPos(oilMu.xAt(rsIdx));
406 
407  assert(oilMu.numY(rsIdx) == invOilB.numY(rsIdx));
408 
409  size_t numPressures = oilMu.numY(rsIdx);
410  for (unsigned pIdx = 0; pIdx < numPressures; ++pIdx)
411  invOilBMu.appendSamplePoint(rsIdx,
412  oilMu.yAt(rsIdx, pIdx),
413  invOilB.valueAt(rsIdx, pIdx)
414  / oilMu.valueAt(rsIdx, pIdx));
415 
416  // the sampling points in UniformXTabulated2DFunction are always sorted
417  // in ascending order. Thus, the value for saturated oil is the first one
418  // (i.e., the one for the lowest pressure value)
419  satPressuresArray.push_back(oilMu.yAt(rsIdx, 0));
420  invSatOilBArray.push_back(invOilB.valueAt(rsIdx, 0));
421  invSatOilBMuArray.push_back(invSatOilBArray.back()/satOilMu.valueAt(rsIdx));
422  }
423 
424  invSatOilB.setXYContainers(satPressuresArray, invSatOilBArray);
425  invSatOilBMu.setXYContainers(satPressuresArray, invSatOilBMuArray);
426 
427  updateSaturationPressure_(regionIdx);
428  }
429  }
430 
434  unsigned numRegions() const
435  { return inverseOilBMuTable_.size(); }
436 
440  template <class Evaluation>
441  Evaluation internalEnergy(unsigned,
442  const Evaluation&,
443  const Evaluation&,
444  const Evaluation&) const
445  {
446  throw std::runtime_error("Requested the enthalpy of oil but the thermal option is not enabled");
447  }
448 
452  template <class Evaluation>
453  Evaluation viscosity(unsigned regionIdx,
454  const Evaluation& /*temperature*/,
455  const Evaluation& pressure,
456  const Evaluation& Rs) const
457  {
458  // ATTENTION: Rs is the first axis!
459  const Evaluation& invBo = inverseOilBTable_[regionIdx].eval(Rs, pressure, /*extrapolate=*/true);
460  const Evaluation& invMuoBo = inverseOilBMuTable_[regionIdx].eval(Rs, pressure, /*extrapolate=*/true);
461 
462  return invBo/invMuoBo;
463  }
464 
468  template <class Evaluation>
469  Evaluation saturatedViscosity(unsigned regionIdx,
470  const Evaluation& /*temperature*/,
471  const Evaluation& pressure) const
472  {
473  // ATTENTION: Rs is the first axis!
474  const Evaluation& invBo = inverseSaturatedOilBTable_[regionIdx].eval(pressure, /*extrapolate=*/true);
475  const Evaluation& invMuoBo = inverseSaturatedOilBMuTable_[regionIdx].eval(pressure, /*extrapolate=*/true);
476 
477  return invBo/invMuoBo;
478  }
479 
483  template <class Evaluation>
484  Evaluation inverseFormationVolumeFactor(unsigned regionIdx,
485  const Evaluation& /*temperature*/,
486  const Evaluation& pressure,
487  const Evaluation& Rs) const
488  {
489  // ATTENTION: Rs is represented by the _first_ axis!
490  return inverseOilBTable_[regionIdx].eval(Rs, pressure, /*extrapolate=*/true);
491  }
492 
496  template <class Evaluation>
497  Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx,
498  const Evaluation& /*temperature*/,
499  const Evaluation& pressure) const
500  {
501  // ATTENTION: Rs is represented by the _first_ axis!
502  return inverseSaturatedOilBTable_[regionIdx].eval(pressure, /*extrapolate=*/true);
503  }
504 
508  template <class Evaluation>
509  Evaluation saturatedGasDissolutionFactor(unsigned regionIdx,
510  const Evaluation& /*temperature*/,
511  const Evaluation& pressure) const
512  { return saturatedGasDissolutionFactorTable_[regionIdx].eval(pressure, /*extrapolate=*/true); }
513 
521  template <class Evaluation>
522  Evaluation saturatedGasDissolutionFactor(unsigned regionIdx,
523  const Evaluation& /*temperature*/,
524  const Evaluation& pressure,
525  const Evaluation& oilSaturation,
526  Evaluation maxOilSaturation) const
527  {
528  Evaluation tmp =
529  saturatedGasDissolutionFactorTable_[regionIdx].eval(pressure, /*extrapolate=*/true);
530 
531  // apply the vaporization parameters for the gas phase (cf. the Eclipse VAPPARS
532  // keyword)
533  maxOilSaturation = min(maxOilSaturation, Scalar(1.0));
534  if (vapPar2_ > 0.0 && maxOilSaturation > 0.01 && oilSaturation < maxOilSaturation) {
535  constexpr const Scalar eps = 0.001;
536  const Evaluation& So = max(oilSaturation, eps);
537  tmp *= max(1e-3, pow(So/maxOilSaturation, vapPar2_));
538  }
539 
540  return tmp;
541  }
542 
549  template <class Evaluation>
550  Evaluation saturationPressure(unsigned regionIdx,
551  const Evaluation&,
552  const Evaluation& Rs) const
553  {
554  typedef MathToolbox<Evaluation> Toolbox;
555 
556  const auto& RsTable = saturatedGasDissolutionFactorTable_[regionIdx];
557  constexpr const Scalar eps = std::numeric_limits<typename Toolbox::Scalar>::epsilon()*1e6;
558 
559  // use the saturation pressure function to get a pretty good initial value
560  Evaluation pSat = saturationPressure_[regionIdx].eval(Rs, /*extrapolate=*/true);
561 
562  // Newton method to do the remaining work. If the initial
563  // value is good, this should only take two to three
564  // iterations...
565  bool onProbation = false;
566  for (int i = 0; i < 20; ++i) {
567  const Evaluation& f = RsTable.eval(pSat, /*extrapolate=*/true) - Rs;
568  const Evaluation& fPrime = RsTable.evalDerivative(pSat, /*extrapolate=*/true);
569 
570  // If the derivative is "zero" Newton will not converge,
571  // so simply return our initial guess.
572  if (std::abs(scalarValue(fPrime)) < 1.0e-30) {
573  return pSat;
574  }
575 
576  const Evaluation& delta = f/fPrime;
577 
578  pSat -= delta;
579 
580  if (pSat < 0.0) {
581  // if the pressure is lower than 0 Pascals, we set it back to 0. if this
582  // happens twice, we give up and just return 0 Pa...
583  if (onProbation)
584  return 0.0;
585 
586  onProbation = true;
587  pSat = 0.0;
588  }
589 
590  if (std::abs(scalarValue(delta)) < std::abs(scalarValue(pSat))*eps)
591  return pSat;
592  }
593 
594  std::stringstream errlog;
595  errlog << "Finding saturation pressure did not converge:"
596  << " pSat = " << pSat
597  << ", Rs = " << Rs;
598 #if HAVE_OPM_COMMON
599  OpmLog::debug("Live oil saturation pressure", errlog.str());
600 #endif
601  throw NumericalIssue(errlog.str());
602  }
603 
604  template <class Evaluation>
605  Evaluation diffusionCoefficient(const Evaluation& /*temperature*/,
606  const Evaluation& /*pressure*/,
607  unsigned /*compIdx*/) const
608  {
609  throw std::runtime_error("Not implemented: The PVT model does not provide a diffusionCoefficient()");
610  }
611  const Scalar gasReferenceDensity(unsigned regionIdx) const
612  { return gasReferenceDensity_[regionIdx]; }
613 
614  const Scalar oilReferenceDensity(unsigned regionIdx) const
615  { return oilReferenceDensity_[regionIdx]; }
616 
617  const std::vector<TabulatedTwoDFunction>& inverseOilBTable() const
618  { return inverseOilBTable_; }
619 
620  const std::vector<TabulatedTwoDFunction>& oilMuTable() const
621  { return oilMuTable_; }
622 
623  const std::vector<TabulatedTwoDFunction>& inverseOilBMuTable() const
624  { return inverseOilBMuTable_; }
625 
626  const std::vector<TabulatedOneDFunction>& saturatedOilMuTable() const
627  { return saturatedOilMuTable_; }
628 
629  const std::vector<TabulatedOneDFunction>& inverseSaturatedOilBTable() const
630  { return inverseSaturatedOilBTable_; }
631 
632  const std::vector<TabulatedOneDFunction>& inverseSaturatedOilBMuTable() const
633  { return inverseSaturatedOilBMuTable_; }
634 
635  const std::vector<TabulatedOneDFunction>& saturatedGasDissolutionFactorTable() const
636  { return saturatedGasDissolutionFactorTable_; }
637 
638  const std::vector<TabulatedOneDFunction>& saturationPressure() const
639  { return saturationPressure_; }
640 
641  Scalar vapPar2() const
642  { return vapPar2_; }
643 
644  bool operator==(const LiveOilPvt<Scalar>& data) const
645  {
646  return this->gasReferenceDensity_ == data.gasReferenceDensity_ &&
647  this->oilReferenceDensity_ == data.oilReferenceDensity_ &&
648  this->inverseOilBTable() == data.inverseOilBTable() &&
649  this->oilMuTable() == data.oilMuTable() &&
650  this->inverseOilBMuTable() == data.inverseOilBMuTable() &&
651  this->saturatedOilMuTable() == data.saturatedOilMuTable() &&
652  this->inverseSaturatedOilBTable() == data.inverseSaturatedOilBTable() &&
653  this->inverseSaturatedOilBMuTable() == data.inverseSaturatedOilBMuTable() &&
654  this->saturatedGasDissolutionFactorTable() == data.saturatedGasDissolutionFactorTable() &&
655  this->vapPar2() == data.vapPar2();
656  }
657 
658 private:
659  void updateSaturationPressure_(unsigned regionIdx)
660  {
661  typedef std::pair<Scalar, Scalar> Pair;
662  const auto& gasDissolutionFac = saturatedGasDissolutionFactorTable_[regionIdx];
663 
664  // create the function representing saturation pressure depending of the mass
665  // fraction in gas
666  size_t n = gasDissolutionFac.numSamples();
667  Scalar delta = (gasDissolutionFac.xMax() - gasDissolutionFac.xMin())/Scalar(n + 1);
668 
669  SamplingPoints pSatSamplePoints;
670  Scalar Rs = 0;
671  for (size_t i=0; i <= n; ++ i) {
672  Scalar pSat = gasDissolutionFac.xMin() + Scalar(i)*delta;
673  Rs = saturatedGasDissolutionFactor(regionIdx,
674  /*temperature=*/Scalar(1e30),
675  pSat);
676 
677  Pair val(Rs, pSat);
678  pSatSamplePoints.push_back(val);
679  }
680 
681  //Prune duplicate Rs values (can occur, and will cause problems in further interpolation)
682  auto x_coord_comparator = [](const Pair& a, const Pair& b) { return a.first == b.first; };
683  auto last = std::unique(pSatSamplePoints.begin(), pSatSamplePoints.end(), x_coord_comparator);
684  if (std::distance(pSatSamplePoints.begin(), last) > 1) // only remove them if there are more than two points
685  pSatSamplePoints.erase(last, pSatSamplePoints.end());
686 
687  saturationPressure_[regionIdx].setContainerOfTuples(pSatSamplePoints);
688  }
689 
690  std::vector<Scalar> gasReferenceDensity_;
691  std::vector<Scalar> oilReferenceDensity_;
692  std::vector<TabulatedTwoDFunction> inverseOilBTable_;
693  std::vector<TabulatedTwoDFunction> oilMuTable_;
694  std::vector<TabulatedTwoDFunction> inverseOilBMuTable_;
695  std::vector<TabulatedOneDFunction> saturatedOilMuTable_;
696  std::vector<TabulatedOneDFunction> inverseSaturatedOilBTable_;
697  std::vector<TabulatedOneDFunction> inverseSaturatedOilBMuTable_;
698  std::vector<TabulatedOneDFunction> saturatedGasDissolutionFactorTable_;
699  std::vector<TabulatedOneDFunction> saturationPressure_;
700 
701  Scalar vapPar2_;
702 };
703 
704 } // namespace Opm
705 
706 #endif
A traits class which provides basic mathematical functions for arbitrary scalar floating point values...
Implements a linearly interpolated scalar function that depends on one variable.
Implements a scalar function that depends on two variables and which is sampled uniformly in the X di...
This class represents the Pressure-Volume-Temperature relations of the oil phas with dissolved gas.
Definition: LiveOilPvt.hpp:52
Evaluation internalEnergy(unsigned, const Evaluation &, const Evaluation &, const Evaluation &) const
Returns the specific enthalpy [J/kg] of oil given a set of parameters.
Definition: LiveOilPvt.hpp:441
void setSaturatedOilFormationVolumeFactor(unsigned regionIdx, const SamplingPoints &samplePoints)
Initialize the function for the oil formation volume factor.
Definition: LiveOilPvt.hpp:303
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:484
Evaluation saturatedGasDissolutionFactor(unsigned regionIdx, const Evaluation &, const Evaluation &pressure, const Evaluation &oilSaturation, Evaluation maxOilSaturation) const
Returns the gas dissolution factor [m^3/m^3] of the oil phase.
Definition: LiveOilPvt.hpp:522
Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &, const Evaluation &pressure) const
Returns the formation volume factor [-] of the fluid phase.
Definition: LiveOilPvt.hpp:497
void setSaturatedOilViscosity(unsigned regionIdx, const SamplingPoints &samplePoints)
Initialize the phase viscosity for gas saturated oil.
Definition: LiveOilPvt.hpp:357
Evaluation saturationPressure(unsigned regionIdx, const Evaluation &, 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:550
void setSaturatedOilGasDissolutionFactor(unsigned regionIdx, const SamplingPoints &samplePoints)
Initialize the function for the gas dissolution factor .
Definition: LiveOilPvt.hpp:291
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:453
void setReferenceDensities(unsigned regionIdx, Scalar rhoRefOil, Scalar rhoRefGas, Scalar)
Initialize the reference densities of all fluids for a given PVT region.
Definition: LiveOilPvt.hpp:277
void setOilViscosity(unsigned regionIdx, const TabulatedTwoDFunction &muo)
Initialize the viscosity of the oil phase.
Definition: LiveOilPvt.hpp:347
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:509
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:469
unsigned numRegions() const
Return the number of PVT regions which are considered by this PVT-object.
Definition: LiveOilPvt.hpp:434
void initEnd()
Finish initializing the oil phase PVT properties.
Definition: LiveOilPvt.hpp:385
void setInverseOilFormationVolumeFactor(unsigned regionIdx, const TabulatedTwoDFunction &invBo)
Initialize the function for the oil formation volume factor.
Definition: LiveOilPvt.hpp:339
Definition: Exceptions.hpp:46
Implements a linearly interpolated scalar function that depends on one variable.
Definition: Tabulated1DFunction.hpp:48
Implements a scalar function that depends on two variables and which is sampled uniformly in the X di...
Definition: UniformXTabulated2DFunction.hpp:54
Definition: MathToolbox.hpp:50