My Project
WetGasPvt.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_WET_GAS_PVT_HPP
28 #define OPM_WET_GAS_PVT_HPP
29 
33 
34 #if HAVE_ECL_INPUT
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>
38 #endif
39 
40 #if HAVE_OPM_COMMON
41 #include <opm/common/OpmLog/OpmLog.hpp>
42 #endif
43 
44 namespace Opm {
45 
50 template <class Scalar>
51 class WetGasPvt
52 {
53  using SamplingPoints = std::vector<std::pair<Scalar, Scalar>>;
54 
55 public:
58 
59  WetGasPvt()
60  {
61  vapPar1_ = 0.0;
62  }
63 
64  WetGasPvt(const std::vector<Scalar>& gasReferenceDensity,
65  const std::vector<Scalar>& oilReferenceDensity,
66  const std::vector<TabulatedTwoDFunction>& inverseGasB,
67  const std::vector<TabulatedOneDFunction>& inverseSaturatedGasB,
68  const std::vector<TabulatedTwoDFunction>& gasMu,
69  const std::vector<TabulatedTwoDFunction>& inverseGasBMu,
70  const std::vector<TabulatedOneDFunction>& inverseSaturatedGasBMu,
71  const std::vector<TabulatedOneDFunction>& saturatedOilVaporizationFactorTable,
72  const std::vector<TabulatedOneDFunction>& saturationPressure,
73  Scalar vapPar1)
74  : gasReferenceDensity_(gasReferenceDensity)
75  , oilReferenceDensity_(oilReferenceDensity)
76  , inverseGasB_(inverseGasB)
77  , inverseSaturatedGasB_(inverseSaturatedGasB)
78  , gasMu_(gasMu)
79  , inverseGasBMu_(inverseGasBMu)
80  , inverseSaturatedGasBMu_(inverseSaturatedGasBMu)
81  , saturatedOilVaporizationFactorTable_(saturatedOilVaporizationFactorTable)
82  , saturationPressure_(saturationPressure)
83  , vapPar1_(vapPar1)
84  {
85  }
86 
87 
88 #if HAVE_ECL_INPUT
94  void initFromState(const EclipseState& eclState, const Schedule& schedule)
95  {
96  const auto& pvtgTables = eclState.getTableManager().getPvtgTables();
97  const auto& densityTable = eclState.getTableManager().getDensityTable();
98 
99  assert(pvtgTables.size() == densityTable.size());
100 
101  size_t numRegions = pvtgTables.size();
102  setNumRegions(numRegions);
103 
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;
108 
109  setReferenceDensities(regionIdx, rhoRefO, rhoRefG, rhoRefW);
110  }
111 
112  for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
113  const auto& pvtgTable = pvtgTables[regionIdx];
114 
115  const auto& saturatedTable = pvtgTable.getSaturatedTable();
116  assert(saturatedTable.numRows() > 1);
117 
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];
123 
124  oilVaporizationFac.setXYArrays(saturatedTable.numRows(),
125  saturatedTable.getColumn("PG"),
126  saturatedTable.getColumn("RV"));
127 
128  std::vector<Scalar> invSatGasBArray;
129  std::vector<Scalar> invSatGasBMuArray;
130 
131  // extract the table for the gas dissolution and the oil formation volume factors
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);
136 
137  invGasB.appendXPos(pg);
138  gasMu.appendXPos(pg);
139 
140  invSatGasBArray.push_back(1.0/B);
141  invSatGasBMuArray.push_back(1.0/(mu*B));
142 
143  assert(invGasB.numX() == outerIdx + 1);
144  assert(gasMu.numX() == outerIdx + 1);
145 
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);
152 
153  invGasB.appendSamplePoint(outerIdx, Rv, 1.0/Bg);
154  gasMu.appendSamplePoint(outerIdx, Rv, mug);
155  }
156  }
157 
158  {
159  std::vector<double> tmpPressure = saturatedTable.getColumn("PG").vectorCopy( );
160 
161  invSatGasB.setXYContainers(tmpPressure, invSatGasBArray);
162  invSatGasBMu.setXYContainers(tmpPressure, invSatGasBMuArray);
163  }
164 
165  // make sure to have at least two sample points per gas pressure value
166  for (unsigned xIdx = 0; xIdx < invGasB.numX(); ++xIdx) {
167  // a single sample point is definitely needed
168  assert(invGasB.numY(xIdx) > 0);
169 
170  // everything is fine if the current table has two or more sampling points
171  // for a given mole fraction
172  if (invGasB.numY(xIdx) > 1)
173  continue;
174 
175  // find the master table which will be used as a template to extend the
176  // current line. We define master table as the first table which has values
177  // for undersaturated gas...
178  size_t masterTableIdx = xIdx + 1;
179  for (; masterTableIdx < saturatedTable.numRows(); ++masterTableIdx)
180  {
181  if (pvtgTable.getUnderSaturatedTable(masterTableIdx).numRows() > 1)
182  break;
183  }
184 
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!");
188 
189 
190  // extend the current table using the master table.
191  extendPvtgTable_(regionIdx,
192  xIdx,
193  pvtgTable.getUnderSaturatedTable(xIdx),
194  pvtgTable.getUnderSaturatedTable(masterTableIdx));
195  }
196  }
197 
198  vapPar1_ = 0.0;
199  const auto& oilVap = schedule[0].oilvap();
200  if (oilVap.getType() == OilVaporizationProperties::OilVaporization::VAPPARS) {
201  vapPar1_ = oilVap.vap1();
202  }
203 
204  initEnd();
205  }
206 
207 private:
208  void extendPvtgTable_(unsigned regionIdx,
209  unsigned xIdx,
210  const SimpleTable& curTable,
211  const SimpleTable& masterTable)
212  {
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();
216 
217  auto& invGasB = inverseGasB_[regionIdx];
218  auto& gasMu = gasMu_[regionIdx];
219 
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");
224 
225  // compute the gas pressure for the new entry
226  Scalar diffRv = RVColumn[newRowIdx] - RVColumn[newRowIdx - 1];
227  Scalar newRv = RvArray.back() + diffRv;
228 
229  // calculate the compressibility of the master table
230  Scalar B1 = BGColumn[newRowIdx];
231  Scalar B2 = BGColumn[newRowIdx - 1];
232  Scalar x = (B1 - B2)/( (B1 + B2)/2.0 );
233 
234  // calculate the gas formation volume factor which exhibits the same
235  // "compressibility" for the new value of Rv
236  Scalar newBg = gasBArray.back()*(1.0 + x/2.0)/(1.0 - x/2.0);
237 
238  // calculate the "viscosibility" of the master table
239  Scalar mu1 = viscosityColumn[newRowIdx];
240  Scalar mu2 = viscosityColumn[newRowIdx - 1];
241  Scalar xMu = (mu1 - mu2)/( (mu1 + mu2)/2.0 );
242 
243  // calculate the gas formation volume factor which exhibits the same
244  // compressibility for the new pressure
245  Scalar newMug = gasMuArray.back()*(1.0 + xMu/2)/(1.0 - xMu/2.0);
246 
247  // append the new values to the arrays which we use to compute the additional
248  // values ...
249  RvArray.push_back(newRv);
250  gasBArray.push_back(newBg);
251  gasMuArray.push_back(newMug);
252 
253  // ... and register them with the internal table objects
254  invGasB.appendSamplePoint(xIdx, newRv, 1.0/newBg);
255  gasMu.appendSamplePoint(xIdx, newRv, newMug);
256  }
257  }
258 
259 public:
260 #endif // HAVE_ECL_INPUT
261 
262  void setNumRegions(size_t numRegions)
263  {
264  oilReferenceDensity_.resize(numRegions);
265  gasReferenceDensity_.resize(numRegions);
266  inverseGasB_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
267  inverseGasBMu_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
268  inverseSaturatedGasB_.resize(numRegions);
269  inverseSaturatedGasBMu_.resize(numRegions);
270  gasMu_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
271  saturatedOilVaporizationFactorTable_.resize(numRegions);
272  saturationPressure_.resize(numRegions);
273  }
274 
278  void setReferenceDensities(unsigned regionIdx,
279  Scalar rhoRefOil,
280  Scalar rhoRefGas,
281  Scalar /*rhoRefWater*/)
282  {
283  oilReferenceDensity_[regionIdx] = rhoRefOil;
284  gasReferenceDensity_[regionIdx] = rhoRefGas;
285  }
286 
292  void setSaturatedGasOilVaporizationFactor(unsigned regionIdx, const SamplingPoints& samplePoints)
293  { saturatedOilVaporizationFactorTable_[regionIdx].setContainerOfTuples(samplePoints); }
294 
304  void setSaturatedGasFormationVolumeFactor(unsigned regionIdx, const SamplingPoints& samplePoints)
305  {
306  auto& invGasB = inverseGasB_[regionIdx];
307 
308  const auto& RvTable = saturatedOilVaporizationFactorTable_[regionIdx];
309 
310  constexpr const Scalar T = 273.15 + 15.56; // [K]
311 
312  constexpr const Scalar RvMin = 0.0;
313  Scalar RvMax = RvTable.eval(saturatedOilVaporizationFactorTable_[regionIdx].xMax(), /*extrapolate=*/true);
314 
315  Scalar poMin = samplePoints.front().first;
316  Scalar poMax = samplePoints.back().first;
317 
318  constexpr const size_t nRv = 20;
319  size_t nP = samplePoints.size()*2;
320 
321  Scalar rhooRef = oilReferenceDensity_[regionIdx];
322 
323  TabulatedOneDFunction gasFormationVolumeFactor;
324  gasFormationVolumeFactor.setContainerOfTuples(samplePoints);
325 
326  updateSaturationPressure_(regionIdx);
327 
328  // calculate a table of estimated densities depending on pressure and gas mass
329  // fraction. note that this assumes oil of constant compressibility. (having said
330  // that, if only the saturated gas densities are available, there's not much
331  // choice.)
332  for (size_t RvIdx = 0; RvIdx < nRv; ++RvIdx) {
333  Scalar Rv = RvMin + (RvMax - RvMin)*RvIdx/nRv;
334 
335  invGasB.appendXPos(Rv);
336 
337  for (size_t pIdx = 0; pIdx < nP; ++pIdx) {
338  Scalar pg = poMin + (poMax - poMin)*pIdx/nP;
339 
340  Scalar poSat = saturationPressure(regionIdx, T, Rv);
341  Scalar BgSat = gasFormationVolumeFactor.eval(poSat, /*extrapolate=*/true);
342  Scalar drhoo_dp = (1.1200 - 1.1189)/((5000 - 4000)*6894.76);
343  Scalar rhoo = rhooRef/BgSat*(1 + drhoo_dp*(pg - poSat));
344 
345  Scalar Bg = rhooRef/rhoo;
346 
347  invGasB.appendSamplePoint(RvIdx, pg, 1.0/Bg);
348  }
349  }
350  }
351 
364  void setInverseGasFormationVolumeFactor(unsigned regionIdx, const TabulatedTwoDFunction& invBg)
365  { inverseGasB_[regionIdx] = invBg; }
366 
372  void setGasViscosity(unsigned regionIdx, const TabulatedTwoDFunction& mug)
373  { gasMu_[regionIdx] = mug; }
374 
382  void setSaturatedGasViscosity(unsigned regionIdx, const SamplingPoints& samplePoints )
383  {
384  auto& oilVaporizationFac = saturatedOilVaporizationFactorTable_[regionIdx];
385 
386  constexpr const Scalar RvMin = 0.0;
387  Scalar RvMax = oilVaporizationFac.eval(saturatedOilVaporizationFactorTable_[regionIdx].xMax(), /*extrapolate=*/true);
388 
389  Scalar poMin = samplePoints.front().first;
390  Scalar poMax = samplePoints.back().first;
391 
392  constexpr const size_t nRv = 20;
393  size_t nP = samplePoints.size()*2;
394 
395  TabulatedOneDFunction mugTable;
396  mugTable.setContainerOfTuples(samplePoints);
397 
398  // calculate a table of estimated densities depending on pressure and gas mass
399  // fraction
400  for (size_t RvIdx = 0; RvIdx < nRv; ++RvIdx) {
401  Scalar Rv = RvMin + (RvMax - RvMin)*RvIdx/nRv;
402 
403  gasMu_[regionIdx].appendXPos(Rv);
404 
405  for (size_t pIdx = 0; pIdx < nP; ++pIdx) {
406  Scalar pg = poMin + (poMax - poMin)*pIdx/nP;
407  Scalar mug = mugTable.eval(pg, /*extrapolate=*/true);
408 
409  gasMu_[regionIdx].appendSamplePoint(RvIdx, pg, mug);
410  }
411  }
412  }
413 
417  void initEnd()
418  {
419  // calculate the final 2D functions which are used for interpolation.
420  size_t numRegions = gasMu_.size();
421  for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
422  // calculate the table which stores the inverse of the product of the gas
423  // formation volume factor and the gas viscosity
424  const auto& gasMu = gasMu_[regionIdx];
425  const auto& invGasB = inverseGasB_[regionIdx];
426  assert(gasMu.numX() == invGasB.numX());
427 
428  auto& invGasBMu = inverseGasBMu_[regionIdx];
429  auto& invSatGasB = inverseSaturatedGasB_[regionIdx];
430  auto& invSatGasBMu = inverseSaturatedGasBMu_[regionIdx];
431 
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));
437 
438  assert(gasMu.numY(pIdx) == invGasB.numY(pIdx));
439 
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));
446 
447  // the sampling points in UniformXTabulated2DFunction are always sorted
448  // in ascending order. Thus, the value for saturated gas is the last one
449  // (i.e., the one with the largest Rv value)
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));
453  }
454 
455  invSatGasB.setXYContainers(satPressuresArray, invSatGasBArray);
456  invSatGasBMu.setXYContainers(satPressuresArray, invSatGasBMuArray);
457 
458  updateSaturationPressure_(regionIdx);
459  }
460  }
461 
465  unsigned numRegions() const
466  { return gasReferenceDensity_.size(); }
467 
471  template <class Evaluation>
472  Evaluation internalEnergy(unsigned,
473  const Evaluation&,
474  const Evaluation&,
475  const Evaluation&) const
476  {
477  throw std::runtime_error("Requested the enthalpy of gas but the thermal option is not enabled");
478  }
479 
483  template <class Evaluation>
484  Evaluation viscosity(unsigned regionIdx,
485  const Evaluation& /*temperature*/,
486  const Evaluation& pressure,
487  const Evaluation& Rv,
488  const Evaluation& /*Rvw*/) const
489  {
490  const Evaluation& invBg = inverseGasB_[regionIdx].eval(pressure, Rv, /*extrapolate=*/true);
491  const Evaluation& invMugBg = inverseGasBMu_[regionIdx].eval(pressure, Rv, /*extrapolate=*/true);
492 
493  return invBg/invMugBg;
494  }
495 
499  template <class Evaluation>
500  Evaluation saturatedViscosity(unsigned regionIdx,
501  const Evaluation& /*temperature*/,
502  const Evaluation& pressure) const
503  {
504  const Evaluation& invBg = inverseSaturatedGasB_[regionIdx].eval(pressure, /*extrapolate=*/true);
505  const Evaluation& invMugBg = inverseSaturatedGasBMu_[regionIdx].eval(pressure, /*extrapolate=*/true);
506 
507  return invBg/invMugBg;
508  }
509 
513  template <class Evaluation>
514  Evaluation inverseFormationVolumeFactor(unsigned regionIdx,
515  const Evaluation& /*temperature*/,
516  const Evaluation& pressure,
517  const Evaluation& Rv,
518  const Evaluation& /*Rvw*/) const
519  { return inverseGasB_[regionIdx].eval(pressure, Rv, /*extrapolate=*/true); }
520 
524  template <class Evaluation>
525  Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx,
526  const Evaluation& /*temperature*/,
527  const Evaluation& pressure) const
528  { return inverseSaturatedGasB_[regionIdx].eval(pressure, /*extrapolate=*/true); }
529 
533  template <class Evaluation>
534  Evaluation saturatedWaterVaporizationFactor(unsigned /*regionIdx*/,
535  const Evaluation& /*temperature*/,
536  const Evaluation& /*pressure*/) const
537  { return 0.0; /* this is non-humid gas! */ }
538 
542  template <class Evaluation>
543  Evaluation saturatedOilVaporizationFactor(unsigned regionIdx,
544  const Evaluation& /*temperature*/,
545  const Evaluation& pressure) const
546  {
547  return saturatedOilVaporizationFactorTable_[regionIdx].eval(pressure, /*extrapolate=*/true);
548  }
549 
557  template <class Evaluation>
558  Evaluation saturatedOilVaporizationFactor(unsigned regionIdx,
559  const Evaluation& /*temperature*/,
560  const Evaluation& pressure,
561  const Evaluation& oilSaturation,
562  Evaluation maxOilSaturation) const
563  {
564  Evaluation tmp =
565  saturatedOilVaporizationFactorTable_[regionIdx].eval(pressure, /*extrapolate=*/true);
566 
567  // apply the vaporization parameters for the gas phase (cf. the Eclipse VAPPARS
568  // keyword)
569  maxOilSaturation = min(maxOilSaturation, Scalar(1.0));
570  if (vapPar1_ > 0.0 && maxOilSaturation > 0.01 && oilSaturation < maxOilSaturation) {
571  constexpr const Scalar eps = 0.001;
572  const Evaluation& So = max(oilSaturation, eps);
573  tmp *= max(1e-3, pow(So/maxOilSaturation, vapPar1_));
574  }
575 
576  return tmp;
577  }
578 
589  template <class Evaluation>
590  Evaluation saturationPressure(unsigned regionIdx,
591  const Evaluation&,
592  const Evaluation& Rv) const
593  {
594  typedef MathToolbox<Evaluation> Toolbox;
595 
596  const auto& RvTable = saturatedOilVaporizationFactorTable_[regionIdx];
597  constexpr const Scalar eps = std::numeric_limits<typename Toolbox::Scalar>::epsilon()*1e6;
598 
599  // use the tabulated saturation pressure function to get a pretty good initial value
600  Evaluation pSat = saturationPressure_[regionIdx].eval(Rv, /*extrapolate=*/true);
601 
602  // Newton method to do the remaining work. If the initial
603  // value is good, this should only take two to three
604  // iterations...
605  bool onProbation = false;
606  for (unsigned i = 0; i < 20; ++i) {
607  const Evaluation& f = RvTable.eval(pSat, /*extrapolate=*/true) - Rv;
608  const Evaluation& fPrime = RvTable.evalDerivative(pSat, /*extrapolate=*/true);
609 
610  // If the derivative is "zero" Newton will not converge,
611  // so simply return our initial guess.
612  if (std::abs(scalarValue(fPrime)) < 1.0e-30) {
613  return pSat;
614  }
615 
616  const Evaluation& delta = f/fPrime;
617 
618  pSat -= delta;
619 
620  if (pSat < 0.0) {
621  // if the pressure is lower than 0 Pascals, we set it back to 0. if this
622  // happens twice, we give up and just return 0 Pa...
623  if (onProbation)
624  return 0.0;
625 
626  onProbation = true;
627  pSat = 0.0;
628  }
629 
630  if (std::abs(scalarValue(delta)) < std::abs(scalarValue(pSat))*eps)
631  return pSat;
632  }
633 
634  std::stringstream errlog;
635  errlog << "Finding saturation pressure did not converge:"
636  << " pSat = " << pSat
637  << ", Rv = " << Rv;
638 #if HAVE_OPM_COMMON
639  OpmLog::debug("Wet gas saturation pressure", errlog.str());
640 #endif
641  throw NumericalIssue(errlog.str());
642  }
643 
644  template <class Evaluation>
645  Evaluation diffusionCoefficient(const Evaluation& /*temperature*/,
646  const Evaluation& /*pressure*/,
647  unsigned /*compIdx*/) const
648  {
649  throw std::runtime_error("Not implemented: The PVT model does not provide a diffusionCoefficient()");
650  }
651 
652  const Scalar gasReferenceDensity(unsigned regionIdx) const
653  { return gasReferenceDensity_[regionIdx]; }
654 
655  const Scalar oilReferenceDensity(unsigned regionIdx) const
656  { return oilReferenceDensity_[regionIdx]; }
657 
658  const std::vector<TabulatedTwoDFunction>& inverseGasB() const {
659  return inverseGasB_;
660  }
661 
662  const std::vector<TabulatedOneDFunction>& inverseSaturatedGasB() const {
663  return inverseSaturatedGasB_;
664  }
665 
666  const std::vector<TabulatedTwoDFunction>& gasMu() const {
667  return gasMu_;
668  }
669 
670  const std::vector<TabulatedTwoDFunction>& inverseGasBMu() const {
671  return inverseGasBMu_;
672  }
673 
674  const std::vector<TabulatedOneDFunction>& inverseSaturatedGasBMu() const {
675  return inverseSaturatedGasBMu_;
676  }
677 
678  const std::vector<TabulatedOneDFunction>& saturatedOilVaporizationFactorTable() const {
679  return saturatedOilVaporizationFactorTable_;
680  }
681 
682  const std::vector<TabulatedOneDFunction>& saturationPressure() const {
683  return saturationPressure_;
684  }
685 
686  Scalar vapPar1() const {
687  return vapPar1_;
688  }
689 
690  bool operator==(const WetGasPvt<Scalar>& data) const
691  {
692  return this->gasReferenceDensity_ == data.gasReferenceDensity_ &&
693  this->oilReferenceDensity_ == data.oilReferenceDensity_ &&
694  this->inverseGasB() == data.inverseGasB() &&
695  this->inverseSaturatedGasB() == data.inverseSaturatedGasB() &&
696  this->gasMu() == data.gasMu() &&
697  this->inverseGasBMu() == data.inverseGasBMu() &&
698  this->inverseSaturatedGasBMu() == data.inverseSaturatedGasBMu() &&
699  this->saturatedOilVaporizationFactorTable() == data.saturatedOilVaporizationFactorTable() &&
700  this->saturationPressure() == data.saturationPressure() &&
701  this->vapPar1() == data.vapPar1();
702  }
703 
704 private:
705  void updateSaturationPressure_(unsigned regionIdx)
706  {
707  const auto& oilVaporizationFac = saturatedOilVaporizationFactorTable_[regionIdx];
708 
709  // create the taublated function representing saturation pressure depending of
710  // Rv
711  size_t n = oilVaporizationFac.numSamples();
712  Scalar delta = (oilVaporizationFac.xMax() - oilVaporizationFac.xMin())/Scalar(n + 1);
713 
714  SamplingPoints pSatSamplePoints;
715  Scalar Rv = 0;
716  for (size_t i = 0; i <= n; ++ i) {
717  Scalar pSat = oilVaporizationFac.xMin() + Scalar(i)*delta;
718  Rv = saturatedOilVaporizationFactor(regionIdx, /*temperature=*/Scalar(1e30), pSat);
719 
720  pSatSamplePoints.emplace_back(Rv, pSat);
721  }
722 
723  //Prune duplicate Rv values (can occur, and will cause problems in further interpolation)
724  auto x_coord_comparator = [](const auto& a, const auto& b) { return a.first == b.first; };
725  auto last = std::unique(pSatSamplePoints.begin(), pSatSamplePoints.end(), x_coord_comparator);
726  if (std::distance(pSatSamplePoints.begin(), last) > 1) // only remove them if there are more than two points
727  pSatSamplePoints.erase(last, pSatSamplePoints.end());
728 
729  saturationPressure_[regionIdx].setContainerOfTuples(pSatSamplePoints);
730  }
731 
732  std::vector<Scalar> gasReferenceDensity_;
733  std::vector<Scalar> oilReferenceDensity_;
734  std::vector<TabulatedTwoDFunction> inverseGasB_;
735  std::vector<TabulatedOneDFunction> inverseSaturatedGasB_;
736  std::vector<TabulatedTwoDFunction> gasMu_;
737  std::vector<TabulatedTwoDFunction> inverseGasBMu_;
738  std::vector<TabulatedOneDFunction> inverseSaturatedGasBMu_;
739  std::vector<TabulatedOneDFunction> saturatedOilVaporizationFactorTable_;
740  std::vector<TabulatedOneDFunction> saturationPressure_;
741 
742  Scalar vapPar1_;
743 };
744 
745 } // namespace Opm
746 
747 #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...
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
Implements a scalar function that depends on two variables and which is sampled uniformly in the X di...
Definition: UniformXTabulated2DFunction.hpp:54
This class represents the Pressure-Volume-Temperature relations of the gas phas with vaporized oil.
Definition: WetGasPvt.hpp:52
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:543
void initEnd()
Finish initializing the gas phase PVT properties.
Definition: WetGasPvt.hpp:417
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:590
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:558
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
Evaluation saturatedWaterVaporizationFactor(unsigned, const Evaluation &, const Evaluation &) const
Returns the water vaporization factor [m^3/m^3] of the gasphase.
Definition: WetGasPvt.hpp:534
unsigned numRegions() const
Return the number of PVT regions which are considered by this PVT-object.
Definition: WetGasPvt.hpp:465
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
void setSaturatedGasOilVaporizationFactor(unsigned regionIdx, const SamplingPoints &samplePoints)
Initialize the function for the oil vaporization factor .
Definition: WetGasPvt.hpp:292
void setInverseGasFormationVolumeFactor(unsigned regionIdx, const TabulatedTwoDFunction &invBg)
Initialize the function for the gas formation volume factor.
Definition: WetGasPvt.hpp:364
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: MathToolbox.hpp:50