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
44namespace Opm {
45
50template <class Scalar>
52{
53 using SamplingPoints = std::vector<std::pair<Scalar, Scalar>>;
54
55public:
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
306private:
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
408public:
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
826private:
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 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