My Project
WetHumidGasPvt.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_HUMID_GAS_PVT_HPP
28 #define OPM_WET_HUMID_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 
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>
52 {
53  using SamplingPoints = std::vector<std::pair<Scalar, Scalar>>;
54 
55 public:
58 
60  {
61  vapPar1_ = 0.0;
62  }
63 
64  WetHumidGasPvt(const std::vector<Scalar>& gasReferenceDensity,
65  const std::vector<Scalar>& oilReferenceDensity,
66  const std::vector<Scalar>& waterReferenceDensity,
67  const std::vector<TabulatedTwoDFunction>& inverseGasBRvwSat,
68  const std::vector<TabulatedTwoDFunction>& inverseGasBRvSat,
69  const std::vector<TabulatedOneDFunction>& inverseSaturatedGasB,
70  const std::vector<TabulatedTwoDFunction>& gasMuRvwSat,
71  const std::vector<TabulatedTwoDFunction>& gasMuRvSat,
72  const std::vector<TabulatedTwoDFunction>& inverseGasBMuRvwSat,
73  const std::vector<TabulatedTwoDFunction>& inverseGasBMuRvSat,
74  const std::vector<TabulatedOneDFunction>& inverseSaturatedGasBMu,
75  const std::vector<TabulatedOneDFunction>& saturatedWaterVaporizationFactorTable,
76  const std::vector<TabulatedOneDFunction>& saturatedOilVaporizationFactorTable,
77  const std::vector<TabulatedOneDFunction>& saturationPressure,
78  Scalar vapPar1)
79  : gasReferenceDensity_(gasReferenceDensity)
80  , oilReferenceDensity_(oilReferenceDensity)
81  , waterReferenceDensity_(waterReferenceDensity)
82  , inverseGasBRvwSat_(inverseGasBRvwSat) // inverse of Bg evaluated at saturated water-gas ratio (Rvw) values; pvtg
83  , inverseGasBRvSat_(inverseGasBRvSat) // inverse of Bg evaluated at saturated oil-gas ratio (Rv) values; pvtgw
84  , inverseSaturatedGasB_(inverseSaturatedGasB) // evaluated at saturated water-gas ratio (Rvw) and oil-gas ratio (Rv) values; pvtgw
85  , gasMuRvwSat_(gasMuRvwSat) // Mug evaluated at saturated water-gas ratio (Rvw) values; pvtg
86  , gasMuRvSat_(gasMuRvSat) // Mug evaluated at saturated oil-gas ratio (Rv) values; pvtgw
87  , inverseGasBMuRvwSat_(inverseGasBMuRvwSat) // Bg^-1*Mug evaluated at saturated water-gas ratio (Rvw) values; pvtg
88  , inverseGasBMuRvSat_(inverseGasBMuRvSat) // Bg^-1*Mug evaluated at saturated oil-gas ratio (Rv) values; pvtgw
89  , inverseSaturatedGasBMu_(inverseSaturatedGasBMu) //pvtgw
90  , saturatedWaterVaporizationFactorTable_(saturatedWaterVaporizationFactorTable) //pvtgw
91  , saturatedOilVaporizationFactorTable_(saturatedOilVaporizationFactorTable) //pvtg
92  , saturationPressure_(saturationPressure)
93  , vapPar1_(vapPar1)
94  {
95  }
96 
97 
98 #if HAVE_ECL_INPUT
104  void initFromState(const EclipseState& eclState, const Schedule& schedule)
105  {
106  const auto& pvtgwTables = eclState.getTableManager().getPvtgwTables();
107  const auto& pvtgTables = eclState.getTableManager().getPvtgTables();
108  const auto& densityTable = eclState.getTableManager().getDensityTable();
109 
110  assert(pvtgwTables.size() == densityTable.size());
111  assert(pvtgTables.size() == densityTable.size());
112 
113  size_t numRegions = pvtgwTables.size();
114  setNumRegions(numRegions);
115 
116  for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
117  Scalar rhoRefO = densityTable[regionIdx].oil;
118  Scalar rhoRefG = densityTable[regionIdx].gas;
119  Scalar rhoRefW = densityTable[regionIdx].water;
120 
121  setReferenceDensities(regionIdx, rhoRefO, rhoRefG, rhoRefW);
122  }
123  // Table PVTGW
124  for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
125  const auto& pvtgwTable = pvtgwTables[regionIdx];
126 
127  const auto& saturatedTable = pvtgwTable.getSaturatedTable();
128  assert(saturatedTable.numRows() > 1);
129 
130  // PVTGW table contains values at saturated Rv
131  auto& gasMuRvSat = gasMuRvSat_[regionIdx];
132  auto& invGasBRvSat = inverseGasBRvSat_[regionIdx];
133  auto& invSatGasB = inverseSaturatedGasB_[regionIdx];
134  auto& invSatGasBMu = inverseSaturatedGasBMu_[regionIdx];
135  auto& waterVaporizationFac = saturatedWaterVaporizationFactorTable_[regionIdx];
136 
137  waterVaporizationFac.setXYArrays(saturatedTable.numRows(),
138  saturatedTable.getColumn("PG"),
139  saturatedTable.getColumn("RW"));
140 
141  std::vector<Scalar> invSatGasBArray;
142  std::vector<Scalar> invSatGasBMuArray;
143 
144  // extract the table for the gas viscosity and formation volume factors
145  for (unsigned outerIdx = 0; outerIdx < saturatedTable.numRows(); ++ outerIdx) {
146  Scalar pg = saturatedTable.get("PG" , outerIdx);
147  Scalar B = saturatedTable.get("BG" , outerIdx);
148  Scalar mu = saturatedTable.get("MUG" , outerIdx);
149 
150  invGasBRvSat.appendXPos(pg);
151  gasMuRvSat.appendXPos(pg);
152 
153  invSatGasBArray.push_back(1.0/B);
154  invSatGasBMuArray.push_back(1.0/(mu*B));
155 
156  assert(invGasBRvSat.numX() == outerIdx + 1);
157  assert(gasMuRvSat.numX() == outerIdx + 1);
158 
159  const auto& underSaturatedTable = pvtgwTable.getUnderSaturatedTable(outerIdx);
160  size_t numRows = underSaturatedTable.numRows();
161  for (size_t innerIdx = 0; innerIdx < numRows; ++ innerIdx) {
162  Scalar Rw = underSaturatedTable.get("RW" , innerIdx);
163  Scalar Bg = underSaturatedTable.get("BG" , innerIdx);
164  Scalar mug = underSaturatedTable.get("MUG" , innerIdx);
165 
166  invGasBRvSat.appendSamplePoint(outerIdx, Rw, 1.0/Bg);
167  gasMuRvSat.appendSamplePoint(outerIdx, Rw, mug);
168  }
169  }
170 
171  {
172  std::vector<double> tmpPressure = saturatedTable.getColumn("PG").vectorCopy( );
173 
174  invSatGasB.setXYContainers(tmpPressure, invSatGasBArray);
175  invSatGasBMu.setXYContainers(tmpPressure, invSatGasBMuArray);
176  }
177 
178  // make sure to have at least two sample points per gas pressure value
179  for (unsigned xIdx = 0; xIdx < invGasBRvSat.numX(); ++xIdx) {
180  // a single sample point is definitely needed
181  assert(invGasBRvSat.numY(xIdx) > 0);
182 
183  // everything is fine if the current table has two or more sampling points
184  // for a given mole fraction
185  if (invGasBRvSat.numY(xIdx) > 1)
186  continue;
187 
188  // find the master table which will be used as a template to extend the
189  // current line. We define master table as the first table which has values
190  // for undersaturated gas...
191  size_t masterTableIdx = xIdx + 1;
192  for (; masterTableIdx < saturatedTable.numRows(); ++masterTableIdx)
193  {
194  if (pvtgwTable.getUnderSaturatedTable(masterTableIdx).numRows() > 1)
195  break;
196  }
197 
198  if (masterTableIdx >= saturatedTable.numRows())
199  throw std::runtime_error("PVTGW tables are invalid: The last table must exhibit at least one "
200  "entry for undersaturated gas!");
201 
202 
203  // extend the current table using the master table.
204  extendPvtgwTable_(regionIdx,
205  xIdx,
206  pvtgwTable.getUnderSaturatedTable(xIdx),
207  pvtgwTable.getUnderSaturatedTable(masterTableIdx));
208  }
209  }
210 
211  // Table PVTG
212  for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
213  const auto& pvtgTable = pvtgTables[regionIdx];
214 
215  const auto& saturatedTable = pvtgTable.getSaturatedTable();
216  assert(saturatedTable.numRows() > 1);
217  // PVTG table contains values at saturated Rvw
218  auto& gasMuRvwSat = gasMuRvwSat_[regionIdx];
219  auto& invGasBRvwSat = inverseGasBRvwSat_[regionIdx];
220  auto& invSatGasB = inverseSaturatedGasB_[regionIdx];
221  auto& invSatGasBMu = inverseSaturatedGasBMu_[regionIdx];
222  auto& oilVaporizationFac = saturatedOilVaporizationFactorTable_[regionIdx];
223 
224  oilVaporizationFac.setXYArrays(saturatedTable.numRows(),
225  saturatedTable.getColumn("PG"),
226  saturatedTable.getColumn("RV"));
227 
228  std::vector<Scalar> invSatGasBArray;
229  std::vector<Scalar> invSatGasBMuArray;
230 
232  for (unsigned outerIdx = 0; outerIdx < saturatedTable.numRows(); ++ outerIdx) {
233  Scalar pg = saturatedTable.get("PG" , outerIdx);
234  Scalar B = saturatedTable.get("BG" , outerIdx);
235  Scalar mu = saturatedTable.get("MUG" , outerIdx);
236 
237  invGasBRvwSat.appendXPos(pg);
238  gasMuRvwSat.appendXPos(pg);
239 
240  invSatGasBArray.push_back(1.0/B);
241  invSatGasBMuArray.push_back(1.0/(mu*B));
242 
243  assert(invGasBRvwSat.numX() == outerIdx + 1);
244  assert(gasMuRvwSat.numX() == outerIdx + 1);
245 
246  const auto& underSaturatedTable = pvtgTable.getUnderSaturatedTable(outerIdx);
247  size_t numRows = underSaturatedTable.numRows();
248  for (size_t innerIdx = 0; innerIdx < numRows; ++ innerIdx) {
249  Scalar Rv = underSaturatedTable.get("RV" , innerIdx);
250  Scalar Bg = underSaturatedTable.get("BG" , innerIdx);
251  Scalar mug = underSaturatedTable.get("MUG" , innerIdx);
252 
253  invGasBRvwSat.appendSamplePoint(outerIdx, Rv, 1.0/Bg);
254  gasMuRvwSat.appendSamplePoint(outerIdx, Rv, mug);
255  }
256  }
257 
258  {
259  std::vector<double> tmpPressure = saturatedTable.getColumn("PG").vectorCopy( );
260 
261  invSatGasB.setXYContainers(tmpPressure, invSatGasBArray);
262  invSatGasBMu.setXYContainers(tmpPressure, invSatGasBMuArray);
263  }
264 
265  // make sure to have at least two sample points per gas pressure value
266  for (unsigned xIdx = 0; xIdx < invGasBRvwSat.numX(); ++xIdx) {
267  // a single sample point is definitely needed
268  assert(invGasBRvwSat.numY(xIdx) > 0);
269 
270  // everything is fine if the current table has two or more sampling points
271  // for a given mole fraction
272  if (invGasBRvwSat.numY(xIdx) > 1)
273  continue;
274 
275  // find the master table which will be used as a template to extend the
276  // current line. We define master table as the first table which has values
277  // for undersaturated gas...
278  size_t masterTableIdx = xIdx + 1;
279  for (; masterTableIdx < saturatedTable.numRows(); ++masterTableIdx)
280  {
281  if (pvtgTable.getUnderSaturatedTable(masterTableIdx).numRows() > 1)
282  break;
283  }
284 
285  if (masterTableIdx >= saturatedTable.numRows())
286  throw std::runtime_error("PVTG tables are invalid: The last table must exhibit at least one "
287  "entry for undersaturated gas!");
288 
289 
290  // extend the current table using the master table.
291  extendPvtgTable_(regionIdx,
292  xIdx,
293  pvtgTable.getUnderSaturatedTable(xIdx),
294  pvtgTable.getUnderSaturatedTable(masterTableIdx));
295  }
296  } //end PVTGW
297  vapPar1_ = 0.0;
298  const auto& oilVap = schedule[0].oilvap();
299  if (oilVap.getType() == OilVaporizationProperties::OilVaporization::VAPPARS) {
300  vapPar1_ = oilVap.vap1();
301  }
302 
303  initEnd();
304  }
305 
306 private:
307  void extendPvtgwTable_(unsigned regionIdx,
308  unsigned xIdx,
309  const SimpleTable& curTable,
310  const SimpleTable& masterTable)
311  {
312  std::vector<double> RvArray = curTable.getColumn("RW").vectorCopy();
313  std::vector<double> gasBArray = curTable.getColumn("BG").vectorCopy();
314  std::vector<double> gasMuArray = curTable.getColumn("MUG").vectorCopy();
315 
316  auto& invGasBRvSat = inverseGasBRvSat_[regionIdx];
317  auto& gasMuRvSat = gasMuRvSat_[regionIdx];
318 
319  for (size_t newRowIdx = 1; newRowIdx < masterTable.numRows(); ++ newRowIdx) {
320  const auto& RVColumn = masterTable.getColumn("RW");
321  const auto& BGColumn = masterTable.getColumn("BG");
322  const auto& viscosityColumn = masterTable.getColumn("MUG");
323 
324  // compute the gas pressure for the new entry
325  Scalar diffRv = RVColumn[newRowIdx] - RVColumn[newRowIdx - 1];
326  Scalar newRv = RvArray.back() + diffRv;
327 
328  // calculate the compressibility of the master table
329  Scalar B1 = BGColumn[newRowIdx];
330  Scalar B2 = BGColumn[newRowIdx - 1];
331  Scalar x = (B1 - B2)/( (B1 + B2)/2.0 );
332 
333  // calculate the gas formation volume factor which exhibits the same
334  // "compressibility" for the new value of Rv
335  Scalar newBg = gasBArray.back()*(1.0 + x/2.0)/(1.0 - x/2.0);
336 
337  // calculate the "viscosibility" of the master table
338  Scalar mu1 = viscosityColumn[newRowIdx];
339  Scalar mu2 = viscosityColumn[newRowIdx - 1];
340  Scalar xMu = (mu1 - mu2)/( (mu1 + mu2)/2.0 );
341 
342  // calculate the gas formation volume factor which exhibits the same
343  // compressibility for the new pressure
344  Scalar newMug = gasMuArray.back()*(1.0 + xMu/2)/(1.0 - xMu/2.0);
345 
346  // append the new values to the arrays which we use to compute the additional
347  // values ...
348  RvArray.push_back(newRv);
349  gasBArray.push_back(newBg);
350  gasMuArray.push_back(newMug);
351 
352  // ... and register them with the internal table objects
353  invGasBRvSat.appendSamplePoint(xIdx, newRv, 1.0/newBg);
354  gasMuRvSat.appendSamplePoint(xIdx, newRv, newMug);
355  }
356  }
357  void extendPvtgTable_(unsigned regionIdx,
358  unsigned xIdx,
359  const SimpleTable& curTable,
360  const SimpleTable& masterTable)
361  {
362  std::vector<double> RvArray = curTable.getColumn("RV").vectorCopy();
363  std::vector<double> gasBArray = curTable.getColumn("BG").vectorCopy();
364  std::vector<double> gasMuArray = curTable.getColumn("MUG").vectorCopy();
365 
366  auto& invGasBRvwSat= inverseGasBRvwSat_[regionIdx];
367  auto& gasMuRvwSat = gasMuRvwSat_[regionIdx];
368 
369  for (size_t newRowIdx = 1; newRowIdx < masterTable.numRows(); ++ newRowIdx) {
370  const auto& RVColumn = masterTable.getColumn("RV");
371  const auto& BGColumn = masterTable.getColumn("BG");
372  const auto& viscosityColumn = masterTable.getColumn("MUG");
373 
374  // compute the gas pressure for the new entry
375  Scalar diffRv = RVColumn[newRowIdx] - RVColumn[newRowIdx - 1];
376  Scalar newRv = RvArray.back() + diffRv;
377 
378  // calculate the compressibility of the master table
379  Scalar B1 = BGColumn[newRowIdx];
380  Scalar B2 = BGColumn[newRowIdx - 1];
381  Scalar x = (B1 - B2)/( (B1 + B2)/2.0 );
382 
383  // calculate the gas formation volume factor which exhibits the same
384  // "compressibility" for the new value of Rv
385  Scalar newBg = gasBArray.back()*(1.0 + x/2.0)/(1.0 - x/2.0);
386 
387  // calculate the "viscosibility" of the master table
388  Scalar mu1 = viscosityColumn[newRowIdx];
389  Scalar mu2 = viscosityColumn[newRowIdx - 1];
390  Scalar xMu = (mu1 - mu2)/( (mu1 + mu2)/2.0 );
391 
392  // calculate the gas formation volume factor which exhibits the same
393  // compressibility for the new pressure
394  Scalar newMug = gasMuArray.back()*(1.0 + xMu/2)/(1.0 - xMu/2.0);
395 
396  // append the new values to the arrays which we use to compute the additional
397  // values ...
398  RvArray.push_back(newRv);
399  gasBArray.push_back(newBg);
400  gasMuArray.push_back(newMug);
401 
402  // ... and register them with the internal table objects
403  invGasBRvwSat.appendSamplePoint(xIdx, newRv, 1.0/newBg);
404  gasMuRvwSat.appendSamplePoint(xIdx, newRv, newMug);
405  }
406  }
407 
408 public:
409 #endif // HAVE_ECL_INPUT
410 
411  void setNumRegions(size_t numRegions)
412  {
413  waterReferenceDensity_.resize(numRegions);
414  oilReferenceDensity_.resize(numRegions);
415  gasReferenceDensity_.resize(numRegions);
416  inverseGasBRvwSat_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
417  inverseGasBRvSat_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
418  inverseGasBMuRvwSat_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
419  inverseGasBMuRvSat_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
420  inverseSaturatedGasB_.resize(numRegions);
421  inverseSaturatedGasBMu_.resize(numRegions);
422  gasMuRvwSat_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
423  gasMuRvSat_.resize(numRegions, TabulatedTwoDFunction{TabulatedTwoDFunction::InterpolationPolicy::RightExtreme});
424  saturatedWaterVaporizationFactorTable_.resize(numRegions);
425  saturatedOilVaporizationFactorTable_.resize(numRegions);
426  saturationPressure_.resize(numRegions);
427  }
428 
432  void setReferenceDensities(unsigned regionIdx,
433  Scalar rhoRefOil,
434  Scalar rhoRefGas,
435  Scalar rhoRefWater)
436  {
437  waterReferenceDensity_[regionIdx] = rhoRefWater;
438  oilReferenceDensity_[regionIdx] = rhoRefOil;
439  gasReferenceDensity_[regionIdx] = rhoRefGas;
440  }
441 
447  void setSaturatedGasWaterVaporizationFactor(unsigned regionIdx, const SamplingPoints& samplePoints)
448  { saturatedWaterVaporizationFactorTable_[regionIdx].setContainerOfTuples(samplePoints); }
449 
455  void setSaturatedGasOilVaporizationFactor(unsigned regionIdx, const SamplingPoints& samplePoints)
456  { saturatedOilVaporizationFactorTable_[regionIdx].setContainerOfTuples(samplePoints); }
457 
458 
462  void initEnd()
463  {
464 
465  //PVTGW
466  // calculate the final 2D functions which are used for interpolation.
467  size_t numRegions = gasMuRvSat_.size();
468  for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
469  // calculate the table which stores the inverse of the product of the gas
470  // formation volume factor and the gas viscosity
471  const auto& gasMuRvSat = gasMuRvSat_[regionIdx];
472  const auto& invGasBRvSat = inverseGasBRvSat_[regionIdx];
473  assert(gasMuRvSat.numX() == invGasBRvSat.numX());
474 
475  auto& invGasBMuRvSat = inverseGasBMuRvSat_[regionIdx];
476  auto& invSatGasB = inverseSaturatedGasB_[regionIdx];
477  auto& invSatGasBMu = inverseSaturatedGasBMu_[regionIdx];
478 
479  std::vector<Scalar> satPressuresArray;
480  std::vector<Scalar> invSatGasBArray;
481  std::vector<Scalar> invSatGasBMuArray;
482  for (size_t pIdx = 0; pIdx < gasMuRvSat.numX(); ++pIdx) {
483  invGasBMuRvSat.appendXPos(gasMuRvSat.xAt(pIdx));
484 
485  assert(gasMuRvSat.numY(pIdx) == invGasBRvSat.numY(pIdx));
486 
487  size_t numRw = gasMuRvSat.numY(pIdx);
488  for (size_t RwIdx = 0; RwIdx < numRw; ++RwIdx)
489  invGasBMuRvSat.appendSamplePoint(pIdx,
490  gasMuRvSat.yAt(pIdx, RwIdx),
491  invGasBRvSat.valueAt(pIdx, RwIdx)
492  / gasMuRvSat.valueAt(pIdx, RwIdx));
493 
494  // the sampling points in UniformXTabulated2DFunction are always sorted
495  // in ascending order. Thus, the value for saturated gas is the last one
496  // (i.e., the one with the largest Rw value)
497  satPressuresArray.push_back(gasMuRvSat.xAt(pIdx));
498  invSatGasBArray.push_back(invGasBRvSat.valueAt(pIdx, numRw - 1));
499  invSatGasBMuArray.push_back(invGasBMuRvSat.valueAt(pIdx, numRw - 1));
500  }
501 
502  invSatGasB.setXYContainers(satPressuresArray, invSatGasBArray);
503  invSatGasBMu.setXYContainers(satPressuresArray, invSatGasBMuArray);
504  }
505 
506  //PVTG
507  // calculate the final 2D functions which are used for interpolation.
508  //size_t numRegions = gasMuRvwSat_.size();
509  for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
510  // calculate the table which stores the inverse of the product of the gas
511  // formation volume factor and the gas viscosity
512  const auto& gasMuRvwSat = gasMuRvwSat_[regionIdx];
513  const auto& invGasBRvwSat = inverseGasBRvwSat_[regionIdx];
514  assert(gasMuRvwSat.numX() == invGasBRvwSat.numX());
515 
516  auto& invGasBMuRvwSat = inverseGasBMuRvwSat_[regionIdx];
517  auto& invSatGasB = inverseSaturatedGasB_[regionIdx];
518  auto& invSatGasBMu = inverseSaturatedGasBMu_[regionIdx];
519 
520  std::vector<Scalar> satPressuresArray;
521  std::vector<Scalar> invSatGasBArray;
522  std::vector<Scalar> invSatGasBMuArray;
523  for (size_t pIdx = 0; pIdx < gasMuRvwSat.numX(); ++pIdx) {
524  invGasBMuRvwSat.appendXPos(gasMuRvwSat.xAt(pIdx));
525 
526  assert(gasMuRvwSat.numY(pIdx) == invGasBRvwSat.numY(pIdx));
527 
528  size_t numRw = gasMuRvwSat.numY(pIdx);
529  for (size_t RwIdx = 0; RwIdx < numRw; ++RwIdx)
530  invGasBMuRvwSat.appendSamplePoint(pIdx,
531  gasMuRvwSat.yAt(pIdx, RwIdx),
532  invGasBRvwSat.valueAt(pIdx, RwIdx)
533  / gasMuRvwSat.valueAt(pIdx, RwIdx));
534 
535  // the sampling points in UniformXTabulated2DFunction are always sorted
536  // in ascending order. Thus, the value for saturated gas is the last one
537  // (i.e., the one with the largest Rw value)
538  satPressuresArray.push_back(gasMuRvwSat.xAt(pIdx));
539  invSatGasBArray.push_back(invGasBRvwSat.valueAt(pIdx, numRw - 1));
540  invSatGasBMuArray.push_back(invGasBMuRvwSat.valueAt(pIdx, numRw - 1));
541  }
542 
543  invSatGasB.setXYContainers(satPressuresArray, invSatGasBArray);
544  invSatGasBMu.setXYContainers(satPressuresArray, invSatGasBMuArray);
545 
546  updateSaturationPressure_(regionIdx);
547  }
548  }
549 
553  unsigned numRegions() const
554  { return gasReferenceDensity_.size(); }
555 
559  template <class Evaluation>
560  Evaluation internalEnergy(unsigned,
561  const Evaluation&,
562  const Evaluation&,
563  const Evaluation&) const
564  {
565  throw std::runtime_error("Requested the enthalpy of gas but the thermal option is not enabled");
566  }
567 
571  template <class Evaluation>
572  Evaluation viscosity(unsigned regionIdx,
573  const Evaluation& /*temperature*/,
574  const Evaluation& pressure,
575  const Evaluation& Rv,
576  const Evaluation& Rvw) const
577  {
578  const Evaluation& temperature = 1E30;
579 
580  if (Rv >= (1.0 - 1e-10)*saturatedOilVaporizationFactor(regionIdx, temperature, pressure)) {
581  const Evaluation& invBg = inverseGasBRvSat_[regionIdx].eval(pressure, Rvw, /*extrapolate=*/true);
582  const Evaluation& invMugBg = inverseGasBMuRvSat_[regionIdx].eval(pressure, Rvw, /*extrapolate=*/true);
583  return invBg/invMugBg;
584  }
585  else {
586  // for Rv undersaturated viscosity is evaluated at saturated Rvw values
587  const Evaluation& invBg = inverseGasBRvwSat_[regionIdx].eval(pressure, Rv, /*extrapolate=*/true);
588  const Evaluation& invMugBg = inverseGasBMuRvwSat_[regionIdx].eval(pressure, Rv, /*extrapolate=*/true);
589  return invBg/invMugBg;
590  }
591  }
592 
596  template <class Evaluation>
597  Evaluation saturatedViscosity(unsigned regionIdx,
598  const Evaluation& /*temperature*/,
599  const Evaluation& pressure) const
600  {
601  const Evaluation& invBg = inverseSaturatedGasB_[regionIdx].eval(pressure, /*extrapolate=*/true);
602  const Evaluation& invMugBg = inverseSaturatedGasBMu_[regionIdx].eval(pressure, /*extrapolate=*/true);
603 
604  return invBg/invMugBg;
605  }
606 
610  // template <class Evaluation>
611  // Evaluation inverseFormationVolumeFactor(unsigned regionIdx,
612  // const Evaluation& /*temperature*/,
613  // const Evaluation& pressure,
614  // const Evaluation& Rw) const
615  // { return inverseGasB_[regionIdx].eval(pressure, Rw, /*extrapolate=*/true); }
616 
617  template <class Evaluation>
618  Evaluation inverseFormationVolumeFactor(unsigned regionIdx,
619  const Evaluation& /*temperature*/,
620  const Evaluation& pressure,
621  const Evaluation& Rv,
622  const Evaluation& Rvw) const
623  {
624  const Evaluation& temperature = 1E30;
625 
626  if (Rv >= (1.0 - 1e-10)*saturatedOilVaporizationFactor(regionIdx, temperature, pressure)) {
627  return inverseGasBRvSat_[regionIdx].eval(pressure, Rvw, /*extrapolate=*/true);
628  }
629  else {
630  // for Rv undersaturated Bg^-1 is evaluated at saturated Rvw values
631  return inverseGasBRvwSat_[regionIdx].eval(pressure, Rv, /*extrapolate=*/true);
632  }
633 
634  }
635 
636 
640  template <class Evaluation>
641  Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx,
642  const Evaluation& /*temperature*/,
643  const Evaluation& pressure) const
644  { return inverseSaturatedGasB_[regionIdx].eval(pressure, /*extrapolate=*/true); }
645 
649  template <class Evaluation>
650  Evaluation saturatedWaterVaporizationFactor(unsigned regionIdx,
651  const Evaluation& /*temperature*/,
652  const Evaluation& pressure) const
653  {
654  return saturatedWaterVaporizationFactorTable_[regionIdx].eval(pressure, /*extrapolate=*/true);
655  }
656 
657  template <class Evaluation>
658  Evaluation saturatedOilVaporizationFactor(unsigned regionIdx,
659  const Evaluation& /*temperature*/,
660  const Evaluation& pressure) const
661  {
662  return saturatedOilVaporizationFactorTable_[regionIdx].eval(pressure, /*extrapolate=*/true);
663  }
664 
672  template <class Evaluation>
673  Evaluation saturatedOilVaporizationFactor(unsigned regionIdx,
674  const Evaluation& /*temperature*/,
675  const Evaluation& pressure,
676  const Evaluation& oilSaturation,
677  Evaluation maxOilSaturation) const
678  {
679  Evaluation tmp =
680  saturatedOilVaporizationFactorTable_[regionIdx].eval(pressure, /*extrapolate=*/true);
681 
682  // apply the vaporization parameters for the gas phase (cf. the Eclipse VAPPARS
683  // keyword)
684  maxOilSaturation = min(maxOilSaturation, Scalar(1.0));
685  if (vapPar1_ > 0.0 && maxOilSaturation > 0.01 && oilSaturation < maxOilSaturation) {
686  constexpr const Scalar eps = 0.001;
687  const Evaluation& So = max(oilSaturation, eps);
688  tmp *= max(1e-3, pow(So/maxOilSaturation, vapPar1_));
689  }
690 
691  return tmp;
692  }
693 
701  //PJPE assume dependence on Rv
702  template <class Evaluation>
703  Evaluation saturationPressure(unsigned regionIdx,
704  const Evaluation&,
705  const Evaluation& Rw) const
706  {
707  using Toolbox = MathToolbox<Evaluation>;
708 
709  const auto& RwTable = saturatedWaterVaporizationFactorTable_[regionIdx];
710  constexpr const Scalar eps = std::numeric_limits<typename Toolbox::Scalar>::epsilon()*1e6;
711 
712  // use the tabulated saturation pressure function to get a pretty good initial value
713  Evaluation pSat = saturationPressure_[regionIdx].eval(Rw, /*extrapolate=*/true);
714 
715  // Newton method to do the remaining work. If the initial
716  // value is good, this should only take two to three
717  // iterations...
718  bool onProbation = false;
719  for (unsigned i = 0; i < 20; ++i) {
720  const Evaluation& f = RwTable.eval(pSat, /*extrapolate=*/true) - Rw;
721  const Evaluation& fPrime = RwTable.evalDerivative(pSat, /*extrapolate=*/true);
722 
723  // If the derivative is "zero" Newton will not converge,
724  // so simply return our initial guess.
725  if (std::abs(scalarValue(fPrime)) < 1.0e-30) {
726  return pSat;
727  }
728 
729  const Evaluation& delta = f/fPrime;
730 
731  pSat -= delta;
732 
733  if (pSat < 0.0) {
734  // if the pressure is lower than 0 Pascals, we set it back to 0. if this
735  // happens twice, we give up and just return 0 Pa...
736  if (onProbation)
737  return 0.0;
738 
739  onProbation = true;
740  pSat = 0.0;
741  }
742 
743  if (std::abs(scalarValue(delta)) < std::abs(scalarValue(pSat))*eps)
744  return pSat;
745  }
746 
747  std::stringstream errlog;
748  errlog << "Finding saturation pressure did not converge:"
749  << " pSat = " << pSat
750  << ", Rw = " << Rw;
751 #if HAVE_OPM_COMMON
752  OpmLog::debug("Wet gas saturation pressure", errlog.str());
753 #endif
754  throw NumericalIssue(errlog.str());
755  }
756 
757  template <class Evaluation>
758  Evaluation diffusionCoefficient(const Evaluation& /*temperature*/,
759  const Evaluation& /*pressure*/,
760  unsigned /*compIdx*/) const
761  {
762  throw std::runtime_error("Not implemented: The PVT model does not provide a diffusionCoefficient()");
763  }
764 
765  const Scalar gasReferenceDensity(unsigned regionIdx) const
766  { return gasReferenceDensity_[regionIdx]; }
767 
768  const Scalar oilReferenceDensity(unsigned regionIdx) const
769  { return oilReferenceDensity_[regionIdx]; }
770 
771  const Scalar waterReferenceDensity(unsigned regionIdx) const
772  { return waterReferenceDensity_[regionIdx]; }
773 
774  const std::vector<TabulatedTwoDFunction>& inverseGasB() const {
775  return inverseGasBRvSat_;
776  }
777 
778  const std::vector<TabulatedOneDFunction>& inverseSaturatedGasB() const {
779  return inverseSaturatedGasB_;
780  }
781 
782  const std::vector<TabulatedTwoDFunction>& gasMu() const {
783  return gasMuRvSat_;
784  }
785 
786  const std::vector<TabulatedTwoDFunction>& inverseGasBMu() const {
787  return inverseGasBMuRvSat_;
788  }
789 
790  const std::vector<TabulatedOneDFunction>& inverseSaturatedGasBMu() const {
791  return inverseSaturatedGasBMu_;
792  }
793 
794  const std::vector<TabulatedOneDFunction>& saturatedWaterVaporizationFactorTable() const {
795  return saturatedWaterVaporizationFactorTable_;
796  }
797 
798  const std::vector<TabulatedOneDFunction>& saturatedOilVaporizationFactorTable() const {
799  return saturatedOilVaporizationFactorTable_;
800  }
801 
802 
803  const std::vector<TabulatedOneDFunction>& saturationPressure() const {
804  return saturationPressure_;
805  }
806 
807  Scalar vapPar1() const {
808  return vapPar1_;
809  }
810 
811  bool operator==(const WetHumidGasPvt<Scalar>& data) const
812  {
813  return this->gasReferenceDensity_ == data.gasReferenceDensity_ &&
814  this->oilReferenceDensity_ == data.oilReferenceDensity_ &&
815  this->waterReferenceDensity_ == data.waterReferenceDensity_ &&
816  this->inverseGasB() == data.inverseGasB() &&
817  this->inverseSaturatedGasB() == data.inverseSaturatedGasB() &&
818  this->gasMu() == data.gasMu() &&
819  this->inverseGasBMu() == data.inverseGasBMu() &&
820  this->inverseSaturatedGasBMu() == data.inverseSaturatedGasBMu() &&
821  this->saturatedWaterVaporizationFactorTable() == data.saturatedWaterVaporizationFactorTable() &&
822  this->saturationPressure() == data.saturationPressure() &&
823  this->vapPar1() == data.vapPar1();
824  }
825 
826 private:
827  void updateSaturationPressure_(unsigned regionIdx)
828  {
829  const auto& oilVaporizationFac = saturatedOilVaporizationFactorTable_[regionIdx];
830 
831  // create the taublated function representing saturation pressure depending of
832  // Rv
833  size_t n = oilVaporizationFac.numSamples();
834  Scalar delta = (oilVaporizationFac.xMax() - oilVaporizationFac.xMin())/Scalar(n + 1);
835 
836  SamplingPoints pSatSamplePoints;
837  Scalar Rv = 0;
838  for (size_t i = 0; i <= n; ++ i) {
839  Scalar pSat = oilVaporizationFac.xMin() + Scalar(i)*delta;
840  Rv = saturatedOilVaporizationFactor(regionIdx, /*temperature=*/Scalar(1e30), pSat);
841 
842  pSatSamplePoints.emplace_back(Rv, pSat);
843  }
844 
845  //Prune duplicate Rv values (can occur, and will cause problems in further interpolation)
846  auto x_coord_comparator = [](const auto& a, const auto& b) { return a.first == b.first; };
847  auto last = std::unique(pSatSamplePoints.begin(), pSatSamplePoints.end(), x_coord_comparator);
848  pSatSamplePoints.erase(last, pSatSamplePoints.end());
849 
850  saturationPressure_[regionIdx].setContainerOfTuples(pSatSamplePoints);
851  }
852 
853  std::vector<Scalar> gasReferenceDensity_;
854  std::vector<Scalar> waterReferenceDensity_;
855  std::vector<Scalar> oilReferenceDensity_;
856  std::vector<TabulatedTwoDFunction> inverseGasBRvwSat_;
857  std::vector<TabulatedTwoDFunction> inverseGasBRvSat_;
858  std::vector<TabulatedOneDFunction> inverseSaturatedGasB_;
859  std::vector<TabulatedTwoDFunction> gasMuRvwSat_;
860  std::vector<TabulatedTwoDFunction> gasMuRvSat_;
861  std::vector<TabulatedTwoDFunction> inverseGasBMuRvwSat_;
862  std::vector<TabulatedTwoDFunction> inverseGasBMuRvSat_;
863  std::vector<TabulatedOneDFunction> inverseSaturatedGasBMu_;
864  std::vector<TabulatedOneDFunction> saturatedWaterVaporizationFactorTable_;
865  std::vector<TabulatedOneDFunction> saturatedOilVaporizationFactorTable_;
866  std::vector<TabulatedOneDFunction> saturationPressure_;
867 
868  Scalar vapPar1_;
869 };
870 
871 } // namespace Opm
872 
873 #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
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 phase with vaporized oil a...
Definition: WetHumidGasPvt.hpp:52
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: WetHumidGasPvt.hpp:703
unsigned numRegions() const
Return the number of PVT regions which are considered by this PVT-object.
Definition: WetHumidGasPvt.hpp:553
void setSaturatedGasOilVaporizationFactor(unsigned regionIdx, const SamplingPoints &samplePoints)
Initialize the function for the oil vaporization factor .
Definition: WetHumidGasPvt.hpp:455
Evaluation viscosity(unsigned regionIdx, const Evaluation &, const Evaluation &pressure, const Evaluation &Rv, const Evaluation &Rvw) const
Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
Definition: WetHumidGasPvt.hpp:572
Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &, const Evaluation &pressure) const
Returns the formation volume factor [-] of water saturated gas at a given pressure.
Definition: WetHumidGasPvt.hpp:641
Evaluation saturatedWaterVaporizationFactor(unsigned regionIdx, const Evaluation &, const Evaluation &pressure) const
Returns the water vaporization factor [m^3/m^3] of the water phase.
Definition: WetHumidGasPvt.hpp:650
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: WetHumidGasPvt.hpp:597
Evaluation internalEnergy(unsigned, const Evaluation &, const Evaluation &, const Evaluation &) const
Returns the specific enthalpy [J/kg] of gas given a set of parameters.
Definition: WetHumidGasPvt.hpp:560
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: WetHumidGasPvt.hpp:673
void setReferenceDensities(unsigned regionIdx, Scalar rhoRefOil, Scalar rhoRefGas, Scalar rhoRefWater)
Initialize the reference densities of all fluids for a given PVT region.
Definition: WetHumidGasPvt.hpp:432
Evaluation inverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &, const Evaluation &pressure, const Evaluation &Rv, const Evaluation &Rvw) const
Returns the formation volume factor [-] of the fluid phase.
Definition: WetHumidGasPvt.hpp:618
void initEnd()
Finish initializing the gas phase PVT properties.
Definition: WetHumidGasPvt.hpp:462
void setSaturatedGasWaterVaporizationFactor(unsigned regionIdx, const SamplingPoints &samplePoints)
Initialize the function for the water vaporization factor .
Definition: WetHumidGasPvt.hpp:447
Definition: MathToolbox.hpp:50