Merge pull request #269 from totto82/addPBDVandPDDV

Add support for PBDV and PDDV in initstateequil.hh
This commit is contained in:
Atgeirr Flø Rasmussen 2018-01-09 11:40:34 +01:00 committed by GitHub
commit 2d9285dbab
4 changed files with 414 additions and 18 deletions

View File

@ -241,6 +241,138 @@ private:
}
};
/**
* Type that implements "dissolved gas-oil ratio"
* tabulated as a function of depth policy. Data
* typically from keyword 'PBVD'.
*/
template <class FluidSystem>
class PBVD : public RsFunction
{
public:
/**
* Constructor.
*
* \param[in] pvtRegionIdx The pvt region index
* \param[in] depth Depth nodes.
* \param[in] pbub Bubble-point pressure at @c depth.
*/
PBVD(const int pvtRegionIdx,
const std::vector<double>& depth,
const std::vector<double>& pbub)
: pvtRegionIdx_(pvtRegionIdx)
, pbubVsDepth_(depth, pbub)
{}
/**
* Function call.
*
* \param[in] depth Depth at which to calculate RS
* value.
*
* \param[in] Pressure in the cell
*
* \param[in] temp Temperature at which to calculate RS
* value.
*
* \return Dissolved gas-oil ratio (RS) at depth @c
* depth and pressure @c press.
*/
double operator()(const double depth,
const double cellPress,
const double temp,
const double satGas = 0.0) const
{
double press = cellPress;
if (satGas <= 0.0) {
if (pbubVsDepth_.xMin() > depth)
press = pbubVsDepth_.valueAt(0);
else if (pbubVsDepth_.xMax() < depth)
press = pbubVsDepth_.valueAt(pbubVsDepth_.numSamples() - 1);
else
press = pbubVsDepth_.eval(depth, /*extrapolate=*/false);
}
return satRs(std::min(press, cellPress), temp);
}
private:
typedef Opm::Tabulated1DFunction<double> PbubVsDepthFunc;
const int pvtRegionIdx_;
PbubVsDepthFunc pbubVsDepth_;
double satRs(const double press, const double temp) const
{
return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press);
}
};
/**
* Type that implements "vaporized oil-gas ratio"
* tabulated as a function of depth policy. Data
* taken from keyword 'PDVD'.
*/
template <class FluidSystem>
class PDVD : public RsFunction
{
public:
/**
* Constructor.
*
* \param[in] pvtRegionIdx The pvt region index
* \param[in] depth Depth nodes.
* \param[in] pbub Dew-point pressure at @c depth.
*/
PDVD(const int pvtRegionIdx,
const std::vector<double>& depth,
const std::vector<double>& pdew)
: pvtRegionIdx_(pvtRegionIdx)
, pdewVsDepth_(depth, pdew)
{}
/**
* Function call.
*
* \param[in] depth Depth at which to calculate RV
* value.
*
* \param[in] cellPress Pressure in the cell
*
* \param[in] temp Temperature at which to calculate RV
* value.
*
* \return Vaporized oil-gas ratio (RV) at depth @c
* depth and pressure @c press.
*/
double operator()(const double depth,
const double cellPress,
const double temp,
const double satOil = 0.0) const
{
double press = cellPress;
if (satOil <= 0.0) {
if (pdewVsDepth_.xMin() > depth)
press = pdewVsDepth_.valueAt(0);
else if (pdewVsDepth_.xMax() < depth)
press = pdewVsDepth_.valueAt(pdewVsDepth_.numSamples() - 1);
else
press = pdewVsDepth_.eval(depth, /*extrapolate=*/false);
}
return satRv(std::min(press, cellPress), temp);
}
private:
typedef Opm::Tabulated1DFunction<double> PdewVsDepthFunc;
const int pvtRegionIdx_;
PdewVsDepthFunc pdewVsDepth_;
double satRv(const double press, const double temp) const
{
return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);
}
};
/**
* Type that implements "vaporized oil-gas ratio"

View File

@ -45,6 +45,8 @@
#include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/RsvdTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/RvvdTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/PbvdTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/PdvdTable.hpp>
#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/common/data/SimulationDataContainer.hpp>
@ -974,7 +976,6 @@ public:
// Create Rs functions.
rsFunc_.reserve(rec.size());
if (FluidSystem::enableDissolvedGas()) {
const Opm::TableContainer& rsvdTables = tables.getRsvdTables();
for (size_t i = 0; i < rec.size(); ++i) {
if (eqlmap.cells(i).empty()) {
rsFunc_.push_back(std::shared_ptr<Miscibility::RsVD<FluidSystem>>());
@ -982,14 +983,26 @@ public:
}
const int pvtIdx = regionPvtIdx_[i];
if (!rec[i].liveOilInitConstantRs()) {
if (rsvdTables.size() <= 0) {
OPM_THROW(std::runtime_error, "Cannot initialise: RSVD table not available.");
const Opm::TableContainer& rsvdTables = tables.getRsvdTables();
const Opm::TableContainer& pbvdTables = tables.getPbvdTables();
if (rsvdTables.size() > 0) {
const Opm::RsvdTable& rsvdTable = rsvdTables.getTable<Opm::RsvdTable>(i);
std::vector<double> depthColumn = rsvdTable.getColumn("DEPTH").vectorCopy();
std::vector<double> rsColumn = rsvdTable.getColumn("RS").vectorCopy();
rsFunc_.push_back(std::make_shared<Miscibility::RsVD<FluidSystem>>(pvtIdx,
depthColumn, rsColumn));
} else if (pbvdTables.size() > 0) {
const Opm::PbvdTable& pbvdTable = pbvdTables.getTable<Opm::PbvdTable>(i);
std::vector<double> depthColumn = pbvdTable.getColumn("DEPTH").vectorCopy();
std::vector<double> pbubColumn = pbvdTable.getColumn("PBUB").vectorCopy();
rsFunc_.push_back(std::make_shared<Miscibility::PBVD<FluidSystem>>(pvtIdx,
depthColumn, pbubColumn));
} else {
OPM_THROW(std::runtime_error, "Cannot initialise: RSVD or PBVD table not available.");
}
const Opm::RsvdTable& rsvdTable = rsvdTables.getTable<Opm::RsvdTable>(i);
std::vector<double> depthColumn = rsvdTable.getColumn("DEPTH").vectorCopy();
std::vector<double> rsColumn = rsvdTable.getColumn("RS").vectorCopy();
rsFunc_.push_back(std::make_shared<Miscibility::RsVD<FluidSystem>>(pvtIdx,
depthColumn, rsColumn));
}
else {
if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) {
@ -1012,7 +1025,6 @@ public:
rvFunc_.reserve(rec.size());
if (FluidSystem::enableVaporizedOil()) {
const Opm::TableContainer& rvvdTables = tables.getRvvdTables();
for (size_t i = 0; i < rec.size(); ++i) {
if (eqlmap.cells(i).empty()) {
rvFunc_.push_back(std::shared_ptr<Miscibility::RvVD<FluidSystem>>());
@ -1020,16 +1032,24 @@ public:
}
const int pvtIdx = regionPvtIdx_[i];
if (!rec[i].wetGasInitConstantRv()) {
if (rvvdTables.size() <= 0) {
OPM_THROW(std::runtime_error, "Cannot initialise: RVVD table not available.");
const Opm::TableContainer& rvvdTables = tables.getRvvdTables();
const Opm::TableContainer& pdvdTables = tables.getPdvdTables();
if (rvvdTables.size() > 0) {
const Opm::RvvdTable& rvvdTable = rvvdTables.getTable<Opm::RvvdTable>(i);
std::vector<double> depthColumn = rvvdTable.getColumn("DEPTH").vectorCopy();
std::vector<double> rvColumn = rvvdTable.getColumn("RV").vectorCopy();
rvFunc_.push_back(std::make_shared<Miscibility::RvVD<FluidSystem>>(pvtIdx,
depthColumn, rvColumn));
} else if (pdvdTables.size() > 0) {
const Opm::PdvdTable& pdvdTable = pdvdTables.getTable<Opm::PdvdTable>(i);
std::vector<double> depthColumn = pdvdTable.getColumn("DEPTH").vectorCopy();
std::vector<double> pdewColumn = pdvdTable.getColumn("PDEW").vectorCopy();
rvFunc_.push_back(std::make_shared<Miscibility::PDVD<FluidSystem>>(pvtIdx,
depthColumn, pdewColumn));
} else {
OPM_THROW(std::runtime_error, "Cannot initialise: RVVD or PDCD table not available.");
}
const Opm::RvvdTable& rvvdTable = rvvdTables.getTable<Opm::RvvdTable>(i);
std::vector<double> depthColumn = rvvdTable.getColumn("DEPTH").vectorCopy();
std::vector<double> rvColumn = rvvdTable.getColumn("RV").vectorCopy();
rvFunc_.push_back(std::make_shared<Miscibility::RvVD<FluidSystem>>(pvtIdx,
depthColumn, rvColumn));
}
else {
if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) {

View File

@ -0,0 +1,152 @@
NOECHO
RUNSPEC ======
WATER
OIL
GAS
DISGAS
VAPOIL
TABDIMS
1 1 40 20 1 20 /
DIMENS
1 1 20
/
WELLDIMS
30 10 2 30 /
START
1 'JAN' 1990 /
NSTACK
25 /
EQLDIMS
-- NTEQUL
1 /
FMTOUT
FMTIN
GRID ======
DX
20*1.0
/
DY
20*1.0
/
DZ
20*5.0
/
TOPS
0.0
/
PORO
20*0.2 /
PERMZ
20*1.0
/
PERMY
20*100.0
/
PERMX
20*100.0
/
BOX
1 1 1 1 1 1 /
PROPS ======
PVTO
-- Rs Pbub Bo Vo
0 1. 1.0000 1.20 /
20 40. 1.0120 1.17 /
40 80. 1.0255 1.14 /
60 120. 1.0380 1.11 /
80 160. 1.0510 1.08 /
100 200. 1.0630 1.06 /
120 240. 1.0750 1.03 /
140 280. 1.0870 1.00 /
160 320. 1.0985 .98 /
180 360. 1.1100 .95 /
200 400. 1.1200 .94
500. 1.1189 .94 /
/
PVTG
-- Pg Rv Bg Vg
100 0.0001 0.010 0.1
0.0 0.0104 0.1 /
200 0.0004 0.005 0.2
0.0 0.0054 0.2 /
/
SWOF
0.2 0 1 0.9
1 1 0 0.1
/
SGOF
0 0 1 0.2
0.8 1 0 0.5
/
PVTW
--RefPres Bw Comp Vw Cv
1. 1.0 4.0E-5 0.96 0.0 /
ROCK
--RefPres Comp
1. 5.0E-5 /
DENSITY
700 1000 1
/
SOLUTION ======
EQUIL
45 150 50 0.25 45 0.35 1 1 0
/
PBVD
0 1.0
50 150. /
PDVD
50. 150.
100. 100 /
RPTSOL
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=2' /
SUMMARY ======
RUNSUM
SEPARATE
SCHEDULE ======
RPTSCHED
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' /
END

View File

@ -849,6 +849,97 @@ void test_DeckWithRSVDAndRVVD()
}
}
void test_DeckWithPBVDAndPDVD()
{
typedef typename TTAG(TestEquilTypeTag) TypeTag;
auto simulator = initSimulator<TypeTag>("data/equil_pbvd_and_pdvd.DATA");
const auto& eclipseState = simulator->gridManager().eclState();
Opm::GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid());
Ewoms::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(), eclipseState, simulator->gridManager().grid(), 9.80665);
const auto& pressures = comp.press();
REQUIRE(pressures.size() == 3);
REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
const int first = 0, last = grid.number_of_cells - 1;
// The relative tolerance is too loose to be very useful,
// but the answer we are checking is the result of an ODE
// solver, and it is unclear if we should check it against
// the true answer or something else.
const double reltol = 1.0e-6;
CHECK_CLOSE(pressures[0][first], 14821552, reltol);
CHECK_CLOSE(pressures[0][last], 15479828, reltol);
CHECK_CLOSE(pressures[1][first], 14911552, reltol);
CHECK_CLOSE(pressures[1][last], 15489828, reltol);
const auto& sats = comp.saturation();
// std::cout << "Saturations:\n";
// for (const auto& sat : sats) {
// for (const double s : sat) {
// std::cout << s << ' ';
// }
// std::cout << std::endl;
// }
const std::vector<double> s_opm[3]{ // opm
{ 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2426402656423233, 0.5383705390740118, 0.7844998821510003, 0.9152832369551807, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.1817779931230221, 0.08471676304481934, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
{ 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.7573597343576767, 0.4616294609259882, 0.03372212472597758, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }
};
for (int phase = 0; phase < 3; ++phase) {
REQUIRE(sats[phase].size() == s_opm[phase].size());
for (size_t i = 0; i < s_opm[phase].size(); ++i) {
//std::cout << std::setprecision(10) << sats[phase][i] << '\n';
CHECK_CLOSE(sats[phase][i], s_opm[phase][i], reltol);
}
}
const auto& rs = comp.rs();
const std::vector<double> rs_opm { // opm
74.55776480956456,
74.6008507125663,
74.6439680789467,
74.68711693934459,
74.73029732443825,
74.77350926494491,
74.81675279162118,
74.86802321984302,
74.96677993174352,
75.09034523640406,
75, 75, 75,75,75, 75, 75, 75, 75, 75 };
const auto& rv = comp.rv();
const std::vector<double> rv_opm {
0.0002488465888573874,
0.0002491051042753978,
0.0002493638084736803,
0.0002496227016360676,
0.0002498817839466295,
0.00025,
0.00025,
0.00025,
0.00025,
0.000251180039180951,
0.0002522295187440788,
0.0002275000000000001,
0.0002125,
0.0001975,
0.0001825,
0.0001675,
0.0001525,
0.0001375,
0.0001225,
0.0001075};
for (size_t i = 0; i < rv_opm.size(); ++i) {
CHECK_CLOSE(rs[i], rs_opm[i], reltol);
CHECK_CLOSE(rv[i], rv_opm[i], reltol);
}
}
void test_DeckWithSwatinit()
{
#if 0
@ -1016,6 +1107,7 @@ int main(int argc, char** argv)
test_DeckWithLiveOil();
test_DeckWithLiveGas();
test_DeckWithRSVDAndRVVD();
test_DeckWithPBVDAndPDVD();
//test_DeckWithSwatinit();
return 0;