| OpmCompile
|
diff --git a/opm/core/io/eclipse/EclipseGridParser.cpp b/opm/core/io/eclipse/EclipseGridParser.cpp
index c33d6fd7..f0106cc3 100644
--- a/opm/core/io/eclipse/EclipseGridParser.cpp
+++ b/opm/core/io/eclipse/EclipseGridParser.cpp
@@ -97,6 +97,7 @@ namespace EclipseKeywords
string("SHEARMOD"), string("POISSONMOD"), string("PWAVEMOD"),
string("MULTPV"), string("PRESSURE"), string("SGAS"),
string("SWAT"), string("SOIL"), string("RS"),
+ string("RV"),
string("DXV"), string("DYV"), string("DZV"),
string("DEPTHZ"), string("TOPS"), string("MAPAXES"),
string("SWCR"), string("SWL"), string("SWU"),
@@ -570,6 +571,8 @@ void EclipseGridParser::convertToSI()
unit = units_.pressure;
} else if (key == "RS") {
unit = units_.gasvol_s / units_.liqvol_s;
+ } else if (key == "RV") {
+ unit = units_.liqvol_s / units_.gasvol_s;
} else if (key == "MAPAXES") {
OPM_MESSAGE("Not applying units to MAPAXES yet!");
unit = 1.0;
diff --git a/opm/core/props/pvt/SinglePvtConstCompr.hpp b/opm/core/props/pvt/SinglePvtConstCompr.hpp
index db73fc88..e97a7242 100644
--- a/opm/core/props/pvt/SinglePvtConstCompr.hpp
+++ b/opm/core/props/pvt/SinglePvtConstCompr.hpp
@@ -224,14 +224,24 @@ namespace Opm
}
- virtual void rbub(const int n,
+ virtual void rsSat(const int n,
const double* /*p*/,
- double* output_rbub,
- double* output_drbubdp) const
+ double* output_rsSat,
+ double* output_drsSatdp) const
{
- std::fill(output_rbub, output_rbub + n, 0.0);
- std::fill(output_drbubdp, output_drbubdp + n, 0.0);
+ std::fill(output_rsSat, output_rsSat + n, 0.0);
+ std::fill(output_drsSatdp, output_drsSatdp + n, 0.0);
+ }
+
+ virtual void rvSat(const int n,
+ const double* /*p*/,
+ double* output_rvSat,
+ double* output_drvSatdp) const
+
+ {
+ std::fill(output_rvSat, output_rvSat + n, 0.0);
+ std::fill(output_drvSatdp, output_drvSatdp + n, 0.0);
}
virtual void R(const int n,
diff --git a/opm/core/props/pvt/SinglePvtDead.cpp b/opm/core/props/pvt/SinglePvtDead.cpp
index 664c38f2..2db23dbe 100644
--- a/opm/core/props/pvt/SinglePvtDead.cpp
+++ b/opm/core/props/pvt/SinglePvtDead.cpp
@@ -176,13 +176,22 @@ namespace Opm
}
- void SinglePvtDead::rbub(const int n,
+ void SinglePvtDead::rsSat(const int n,
const double* /*p*/,
- double* output_rbub,
- double* output_drbubdp) const
+ double* output_rsSat,
+ double* output_drsSatdp) const
{
- std::fill(output_rbub, output_rbub + n, 0.0);
- std::fill(output_drbubdp, output_drbubdp + n, 0.0);
+ std::fill(output_rsSat, output_rsSat + n, 0.0);
+ std::fill(output_drsSatdp, output_drsSatdp + n, 0.0);
+ }
+
+ void SinglePvtDead::rvSat(const int n,
+ const double* /*p*/,
+ double* output_rvSat,
+ double* output_drvSatdp) const
+ {
+ std::fill(output_rvSat, output_rvSat + n, 0.0);
+ std::fill(output_drvSatdp, output_drvSatdp + n, 0.0);
}
void SinglePvtDead::R(const int n,
diff --git a/opm/core/props/pvt/SinglePvtDead.hpp b/opm/core/props/pvt/SinglePvtDead.hpp
index 8b36c5fc..2c40a951 100644
--- a/opm/core/props/pvt/SinglePvtDead.hpp
+++ b/opm/core/props/pvt/SinglePvtDead.hpp
@@ -49,7 +49,7 @@ namespace Opm
double* output_mu) const;
/// Viscosity and its derivatives as a function of p and r.
- /// The fluid is considered saturated if r >= rbub(p).
+ /// The fluid is considered saturated if r >= rsSat(p).
virtual void mu(const int n,
const double* p,
const double* r,
@@ -81,7 +81,7 @@ namespace Opm
double* output_dBdp) const;
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
- /// The fluid is considered saturated if r >= rbub(p).
+ /// The fluid is considered saturated if r >= rsSat(p).
virtual void b(const int n,
const double* p,
const double* r,
@@ -100,11 +100,17 @@ namespace Opm
double* output_dbdr) const;
- /// Gas resolution and its derivatives at bublepoint as a function of p.
- virtual void rbub(const int n,
+ /// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p.
+ virtual void rsSat(const int n,
const double* p,
- double* output_rbub,
- double* output_drbubdp) const;
+ double* output_rsSat,
+ double* output_drsSatdp) const;
+
+ /// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
+ virtual void rvSat(const int n,
+ const double* p,
+ double* output_rvSat,
+ double* output_drvSatdp) const;
/// Solution factor as a function of p and z.
virtual void R(const int n,
diff --git a/opm/core/props/pvt/SinglePvtDeadSpline.cpp b/opm/core/props/pvt/SinglePvtDeadSpline.cpp
index 5db3521d..a0873e0b 100644
--- a/opm/core/props/pvt/SinglePvtDeadSpline.cpp
+++ b/opm/core/props/pvt/SinglePvtDeadSpline.cpp
@@ -178,13 +178,22 @@ namespace Opm
}
- void SinglePvtDeadSpline::rbub(const int n,
+ void SinglePvtDeadSpline::rsSat(const int n,
const double* /*p*/,
- double* output_rbub,
- double* output_drbubdp) const
+ double* output_rsSat,
+ double* output_drsSatdp) const
{
- std::fill(output_rbub, output_rbub + n, 0.0);
- std::fill(output_drbubdp, output_drbubdp + n, 0.0);
+ std::fill(output_rsSat, output_rsSat + n, 0.0);
+ std::fill(output_drsSatdp, output_drsSatdp + n, 0.0);
+ }
+
+ void SinglePvtDeadSpline::rvSat(const int n,
+ const double* /*p*/,
+ double* output_rvSat,
+ double* output_drvSatdp) const
+ {
+ std::fill(output_rvSat, output_rvSat + n, 0.0);
+ std::fill(output_drvSatdp, output_drvSatdp + n, 0.0);
}
void SinglePvtDeadSpline::R(const int n,
diff --git a/opm/core/props/pvt/SinglePvtDeadSpline.hpp b/opm/core/props/pvt/SinglePvtDeadSpline.hpp
index b6af388d..337c8a80 100644
--- a/opm/core/props/pvt/SinglePvtDeadSpline.hpp
+++ b/opm/core/props/pvt/SinglePvtDeadSpline.hpp
@@ -50,7 +50,7 @@ namespace Opm
double* output_mu) const;
/// Viscosity and its derivatives as a function of p and r.
- /// The fluid is considered saturated if r >= rbub(p).
+ /// The fluid is considered saturated if r >= rsSat(p).
virtual void mu(const int n,
const double* p,
const double* r,
@@ -82,7 +82,7 @@ namespace Opm
double* output_dBdp) const;
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
- /// The fluid is considered saturated if r >= rbub(p).
+ /// The fluid is considered saturated if r >= rsSat(p).
virtual void b(const int n,
const double* p,
const double* r,
@@ -100,11 +100,17 @@ namespace Opm
double* output_dbdp,
double* output_dbdr) const;
- /// Gas resolution and its derivatives at bublepoint as a function of p.
- virtual void rbub(const int n,
+ /// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p.
+ virtual void rsSat(const int n,
const double* p,
- double* output_rbub,
- double* output_drbubdp) const;
+ double* output_rsSat,
+ double* output_drsSatdp) const;
+
+ /// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
+ virtual void rvSat(const int n,
+ const double* p,
+ double* output_rvSat,
+ double* output_drvSatdp) const;
/// Solution factor as a function of p and z.
virtual void R(const int n,
diff --git a/opm/core/props/pvt/SinglePvtInterface.hpp b/opm/core/props/pvt/SinglePvtInterface.hpp
index 1deddb5e..83d3febb 100644
--- a/opm/core/props/pvt/SinglePvtInterface.hpp
+++ b/opm/core/props/pvt/SinglePvtInterface.hpp
@@ -56,7 +56,7 @@ namespace Opm
double* output_mu) const = 0;
/// Viscosity as a function of p and r.
- /// The fluid is considered saturated if r >= rbub(p).
+ /// The fluid is considered saturated if r >= rsSat(p).
virtual void mu(const int n,
const double* p,
const double* r,
@@ -88,7 +88,7 @@ namespace Opm
double* output_dBdp) const = 0;
/// The inverse of the volume factor b = 1 / B as a function of p and r.
- /// The fluid is considered saturated if r >= rbub(p).
+ /// The fluid is considered saturated if r >= rsSat(p).
virtual void b(const int n,
const double* p,
const double* r,
@@ -106,11 +106,17 @@ namespace Opm
double* output_dbdp,
double* output_dpdr) const = 0;
- /// Gas resolution at bublepoint as a function of pressure
- virtual void rbub(const int n,
+ /// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p.
+ virtual void rsSat(const int n,
const double* p,
- double* output_rbub,
- double* output_drbubdp) const = 0;
+ double* output_rsSat,
+ double* output_drsSatdp) const = 0;
+
+ /// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
+ virtual void rvSat(const int n,
+ const double* p,
+ double* output_rvSat,
+ double* output_drvSatdp) const = 0;
/// Solution factor as a function of p and z.
diff --git a/opm/core/props/pvt/SinglePvtLiveGas.cpp b/opm/core/props/pvt/SinglePvtLiveGas.cpp
index 2852237d..639ee3ac 100644
--- a/opm/core/props/pvt/SinglePvtLiveGas.cpp
+++ b/opm/core/props/pvt/SinglePvtLiveGas.cpp
@@ -61,7 +61,7 @@ namespace Opm
for (int i=0; i= 2);
+ assert(undersat_gas_tables_[is+1][0].size() >= 2);
+ double val1 =
+ linearInterpolationDerivative(undersat_gas_tables_[is][0],
+ undersat_gas_tables_[is][item],
+ r);
+ double val2 =
+ linearInterpolationDerivative(undersat_gas_tables_[is+1][0],
+ undersat_gas_tables_[is+1][item],
+ r);
+
+ double val = val1 + w * (val2 - val1);
+ return val;
+
+ }
+ } else {
+ if (isSat) { // Saturated case
+ return linearInterpolation(saturated_gas_table_[0],
+ saturated_gas_table_[item],
+ press);
+ } else { // Undersaturated case
+ int is = tableIndex(saturated_gas_table_[0], press);
+ // Extrapolate from first table section
+ if (is == 0 && press < saturated_gas_table_[0][0]) {
+ return linearInterpolation(undersat_gas_tables_[0][0],
+ undersat_gas_tables_[0][item],
+ r);
+ }
+
+ // Extrapolate from last table section
+ //int ltp = saturated_gas_table_[0].size() - 1;
+ //if (is+1 == ltp && press > saturated_gas_table_[0][ltp]) {
+ // return linearInterpolation(undersat_gas_tables_[ltp][0],
+ // undersat_gas_tables_[ltp][item],
+ // r);
+ //}
+
+ // Interpolate between table sections
+ double w = (press - saturated_gas_table_[0][is]) /
+ (saturated_gas_table_[0][is+1] -
+ saturated_gas_table_[0][is]);
+ if (undersat_gas_tables_[is][0].size() < 2) {
+ double val = saturated_gas_table_[item][is] +
+ w*(saturated_gas_table_[item][is+1] -
+ saturated_gas_table_[item][is]);
+ return val;
+ }
+ double val1 =
+ linearInterpolation(undersat_gas_tables_[is][0],
+ undersat_gas_tables_[is][item],
+ r);
+ double val2 =
+ linearInterpolation(undersat_gas_tables_[is+1][0],
+ undersat_gas_tables_[is+1][item],
+ r);
+ double val = val1 + w*(val2 - val1);
+ return val;
+ }
+ }
+ }
+
+
+
} // namespace Opm
diff --git a/opm/core/props/pvt/SinglePvtLiveGas.hpp b/opm/core/props/pvt/SinglePvtLiveGas.hpp
index faa22254..de595393 100644
--- a/opm/core/props/pvt/SinglePvtLiveGas.hpp
+++ b/opm/core/props/pvt/SinglePvtLiveGas.hpp
@@ -47,7 +47,7 @@ namespace Opm
double* output_mu) const;
/// Viscosity and its derivatives as a function of p and r.
- /// The fluid is considered saturated if r >= rbub(p).
+ /// The fluid is considered saturated if r >= rsSat(p).
virtual void mu(const int n,
const double* p,
const double* r,
@@ -79,7 +79,7 @@ namespace Opm
double* output_dBdp) const;
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
- /// The fluid is considered saturated if r >= rbub(p).
+ /// The fluid is considered saturated if r >= rsSat(p).
virtual void b(const int n,
const double* p,
const double* r,
@@ -99,11 +99,17 @@ namespace Opm
- /// Gas resolution and its derivatives at bublepoint as a function of p.
- virtual void rbub(const int n,
+ /// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p.
+ virtual void rsSat(const int n,
const double* p,
- double* output_rbub,
- double* output_drbubdp) const;
+ double* output_rsSat,
+ double* output_drsSatdp) const;
+
+ /// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
+ virtual void rvSat(const int n,
+ const double* p,
+ double* output_rvSat,
+ double* output_drvSatdp) const;
/// Solution factor as a function of p and z.
virtual void R(const int n,
@@ -129,6 +135,11 @@ namespace Opm
const double* surfvol,
const int item,
const bool deriv = false) const;
+ double miscible_gas(const double press,
+ const double r,
+ const PhasePresence& cond,
+ const int item,
+ const int deriv = 0) const;
// PVT properties of wet gas (with vaporised oil)
std::vector > saturated_gas_table_;
std::vector > > undersat_gas_tables_;
diff --git a/opm/core/props/pvt/SinglePvtLiveOil.cpp b/opm/core/props/pvt/SinglePvtLiveOil.cpp
index def278e1..d18252ee 100644
--- a/opm/core/props/pvt/SinglePvtLiveOil.cpp
+++ b/opm/core/props/pvt/SinglePvtLiveOil.cpp
@@ -225,21 +225,30 @@ namespace Opm
}
}
- void SinglePvtLiveOil::rbub(const int n,
+ void SinglePvtLiveOil::rsSat(const int n,
const double* p,
- double* output_rbub,
- double* output_drbubdp) const
+ double* output_rsSat,
+ double* output_drsSatdp) const
{
for (int i = 0; i < n; ++i) {
- output_rbub[i] = linearInterpolation(saturated_oil_table_[0],
+ output_rsSat[i] = linearInterpolation(saturated_oil_table_[0],
saturated_oil_table_[3],p[i]);
- output_drbubdp[i] = linearInterpolationDerivative(saturated_oil_table_[0],
+ output_drsSatdp[i] = linearInterpolationDerivative(saturated_oil_table_[0],
saturated_oil_table_[3],p[i]);
}
}
+ void SinglePvtLiveOil::rvSat(const int n,
+ const double* /*p*/,
+ double* output_rvSat,
+ double* output_drvSatdp) const
+ {
+ std::fill(output_rvSat, output_rvSat + n, 0.0);
+ std::fill(output_drvSatdp, output_drvSatdp + n, 0.0);
+ }
+
/// Solution factor as a function of p and z.
void SinglePvtLiveOil::R(const int n,
const double* p,
diff --git a/opm/core/props/pvt/SinglePvtLiveOil.hpp b/opm/core/props/pvt/SinglePvtLiveOil.hpp
index 4bfa9a76..d0e5fb70 100644
--- a/opm/core/props/pvt/SinglePvtLiveOil.hpp
+++ b/opm/core/props/pvt/SinglePvtLiveOil.hpp
@@ -48,7 +48,7 @@ namespace Opm
double* output_mu) const;
/// Viscosity and its derivatives as a function of p and r.
- /// The fluid is considered saturated if r >= rbub(p).
+ /// The fluid is considered saturated if r >= rsSat(p).
virtual void mu(const int n,
const double* p,
const double* r,
@@ -80,7 +80,7 @@ namespace Opm
double* output_dBdp) const;
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
- /// The fluid is considered saturated if r >= rbub(p).
+ /// The fluid is considered saturated if r >= rsSat(p).
virtual void b(const int n,
const double* p,
const double* r,
@@ -98,11 +98,17 @@ namespace Opm
double* output_dbdp,
double* output_dbdr) const;
- /// Gas resolution and its derivatives at bublepoint as a function of p.
- virtual void rbub(const int n,
+ /// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p.
+ virtual void rsSat(const int n,
const double* p,
- double* output_rbub,
- double* output_drbubdp) const;
+ double* output_rsSat,
+ double* output_drsSatdp) const;
+
+ /// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
+ virtual void rvSat(const int n,
+ const double* p,
+ double* output_rvSat,
+ double* output_drvSatdp) const;
/// Solution factor as a function of p and z.
virtual void R(const int n,
diff --git a/opm/core/simulator/BlackoilState.cpp b/opm/core/simulator/BlackoilState.cpp
index c8c3afa8..0b864a17 100644
--- a/opm/core/simulator/BlackoilState.cpp
+++ b/opm/core/simulator/BlackoilState.cpp
@@ -7,6 +7,7 @@ void
BlackoilState::init(const UnstructuredGrid& g, int num_phases) {
SimulatorState::init(g, num_phases);
gor_.resize(g.number_of_cells, 0.) ;
+ rv_.resize(g.number_of_cells,0.);
// surfvol_ intentionally empty, left to initBlackoilSurfvol
}
diff --git a/opm/core/simulator/BlackoilState.hpp b/opm/core/simulator/BlackoilState.hpp
index 5cf8156e..c6fb7779 100644
--- a/opm/core/simulator/BlackoilState.hpp
+++ b/opm/core/simulator/BlackoilState.hpp
@@ -47,13 +47,16 @@ namespace Opm
std::vector& surfacevol () { return surfvol_; }
std::vector& gasoilratio () { return gor_ ; }
+ std::vector& rv () {return rv_ ; }
const std::vector& surfacevol () const { return surfvol_; }
const std::vector& gasoilratio () const { return gor_ ; }
+ const std::vector& rv () const {return rv_ ; }
private:
std::vector surfvol_;
std::vector gor_ ;
+ std::vector rv_ ;
};
} // namespace Opm
diff --git a/opm/core/simulator/initState_impl.hpp b/opm/core/simulator/initState_impl.hpp
index 13d7728b..12dbe7b5 100644
--- a/opm/core/simulator/initState_impl.hpp
+++ b/opm/core/simulator/initState_impl.hpp
@@ -609,26 +609,26 @@ namespace Opm
}
}
/// Initialize surface volume from pressure and saturation by z = As.
- /// Here the gas/oil ratio is used to compute an intial z for the
- /// computation of A.
+ /// Here the solution gas/oil ratio or vapor oil/gas ratio is used to
+ /// compute an intial z for the computation of A.
template
- void initBlackoilSurfvolUsingRS(const UnstructuredGrid& grid,
+ void initBlackoilSurfvolUsingRSorRV(const UnstructuredGrid& grid,
const Props& props,
State& state)
{
if (props.numPhases() != 3) {
OPM_THROW(std::runtime_error, "initBlackoilSurfvol() is only supported in three-phase simulations.");
}
-
const std::vector& rs = state.gasoilratio();
+ const std::vector& rv = state.rv();
//make input for computation of the A matrix
state.surfacevol() = state.saturation();
- //const PhaseUsage pu = phaseUsageFromDeck(deck);
const PhaseUsage pu = props.phaseUsage();
const int np = props.numPhases();
const int nc = grid.number_of_cells;
+
std::vector allA_a(nc*np*np);
std::vector allA_l(nc*np*np);
std::vector allA_v(nc*np*np);
@@ -641,6 +641,17 @@ namespace Opm
allcells[c] = c;
}
+ std::vector capPressures(nc*np);
+ props.capPress(nc,&state.saturation()[0],&allcells[0],&capPressures[0],NULL);
+
+ std::vector Pw(nc);
+ std::vector Pg(nc);
+
+ for (int c = 0; c < nc; ++c){
+ Pw[c] = state.pressure()[c] + capPressures[c*np + BlackoilPhases::Aqua];
+ Pg[c] = state.pressure()[c] + capPressures[c*np + BlackoilPhases::Vapour];
+ }
+
double z_tmp;
@@ -687,7 +698,7 @@ namespace Opm
if(state.saturation()[np*c + p] > 0)
z_tmp = 1e10;
else
- z_tmp = rs[c];
+ z_tmp = rv[c];
}
else if(p == BlackoilPhases::Vapour)
z_tmp = 1;
@@ -717,7 +728,9 @@ namespace Opm
z[2] += A_v[2 + np*col]*s[col];
}
+ double ztmp = z[2];
z[2] += z[1]*rs[c];
+ z[1] += ztmp*rv[c];
}
}
@@ -739,10 +752,21 @@ namespace Opm
int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
state.gasoilratio()[c] = rs_deck[c_deck];
}
- initBlackoilSurfvolUsingRS(grid, props, state);
+ initBlackoilSurfvolUsingRSorRV(grid, props, state);
computeSaturation(props,state);
- } else {
- OPM_THROW(std::runtime_error, "Temporarily, we require the RS field.");
+ } else if (deck.hasField("RV")){
+ const std::vector& rv_deck = deck.getFloatingPointValue("RV");
+ const int num_cells = grid.number_of_cells;
+ for (int c = 0; c < num_cells; ++c) {
+ int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
+ state.rv()[c] = rv_deck[c_deck];
+ }
+ initBlackoilSurfvolUsingRSorRV(grid, props, state);
+ computeSaturation(props,state);
+ }
+
+ else {
+ OPM_THROW(std::runtime_error, "Temporarily, we require the RS or the RV field.");
}
}
diff --git a/tests/liveoil.DATA b/tests/liveoil.DATA
new file mode 100644
index 00000000..c6a04b40
--- /dev/null
+++ b/tests/liveoil.DATA
@@ -0,0 +1,36 @@
+OIL
+GAS
+WATER
+
+FIELD
+
+PVTO
+-- Rs Pbub Bo Vo
+ .0 14.7 1.0000 1.20 /
+ .165 400. 1.0120 1.17 /
+ .335 800. 1.0255 1.14 /
+ .500 1200. 1.0380 1.11 /
+ .665 1600. 1.0510 1.08 /
+ .828 2000. 1.0630 1.06 /
+ .985 2400. 1.0750 1.03 /
+ 1.130 2800. 1.0870 1.00 /
+ 1.270 3200. 1.0985 .98 /
+ 1.390 3600. 1.1100 .95 /
+ 1.500 4000. 1.1200 .94
+ 5000. 1.1189 .94 /
+ /
+
+PVDG
+-- Pg Bg Vg
+ 14.7 178.08 .0125
+ 400. 5.4777 .0130
+ 800. 2.7392 .0135
+ 1200. 1.8198 .0140
+ 1600. 1.3648 .0145
+ 2000. 1.0957 .0150
+ 2400. 0.9099 .0155
+ 2800. 0.7799 .0160
+ 3200. 0.6871 .0165
+ 3600. 0.6035 .0170
+ 4000. 0.5432 .0175 /
+
diff --git a/tests/test_blackoilfluid.cpp b/tests/test_blackoilfluid.cpp
index 9816284a..a1491f13 100644
--- a/tests/test_blackoilfluid.cpp
+++ b/tests/test_blackoilfluid.cpp
@@ -28,24 +28,18 @@
using namespace Opm;
using namespace std;
-BOOST_AUTO_TEST_CASE(test_blackoilfluid)
-{
+
+std::vector > getProps(const EclipseGridParser deck,PhaseUsage phase_usage_){
- // read eclipse deck
- const string filename = "testFluid.DATA";
- cout << "Reading deck: " << filename << endl;
- const EclipseGridParser deck (filename);
-
- // setup pvt interface
- std::vector > props_;
- PhaseUsage phase_usage_ = phaseUsageFromDeck(deck);
enum PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2 };
int samples = 0;
+ std::vector > props_;
// Set the properties.
props_.resize(phase_usage_.num_phases);
+
// Water PVT
if (phase_usage_.phase_used[Aqua]) {
if (deck.hasField("PVTW")) {
@@ -55,6 +49,7 @@ BOOST_AUTO_TEST_CASE(test_blackoilfluid)
props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(0.5*Opm::prefix::centi*Opm::unit::Poise));
}
}
+
// Oil PVT
if (phase_usage_.phase_used[Liquid]) {
if (deck.hasField("PVDO")) {
@@ -87,6 +82,146 @@ BOOST_AUTO_TEST_CASE(test_blackoilfluid)
}
}
+ return props_;
+}
+
+void testmu(const double reltol, int n, int np, std::vector p, std::vector r,std::vector z,
+ std::vector > props_, std::vector condition)
+{
+ std::vector mu(n);
+ std::vector dmudp(n);
+ std::vector dmudr(n);
+ std::vector mu_new(n);
+ double dmudp_diff;
+ double dmudr_diff;
+ double dmudp_diff_u;
+ double dmudr_diff_u;
+
+ // test mu
+ for (int phase = 0; phase < np; ++phase) {
+ props_[phase]->mu(n, &p[0], &r[0], &condition[0], &mu_new[0], &dmudp[0], &dmudr[0]);
+ props_[phase]->mu(n, &p[0], &z[0], &mu[0]);
+ dmudp_diff = (mu_new[1]-mu_new[0])/(p[1]-p[0]);
+ dmudr_diff = (mu_new[2]-mu_new[0])/(r[2]-r[0]);
+ dmudp_diff_u = (mu_new[4]-mu_new[3])/(p[4]-p[3]);
+ dmudr_diff_u = (mu_new[5]-mu_new[3])/(r[5]-r[3]);
+
+ for (int i = 0; i < n; ++i){
+ BOOST_CHECK_CLOSE(mu_new[i],mu[i],reltol);
+ }
+
+ // saturated case
+ BOOST_CHECK_CLOSE(dmudp_diff,dmudp[0],reltol);
+ BOOST_CHECK_CLOSE(dmudr_diff,dmudr[0],reltol);
+
+ // unsaturated case
+ BOOST_CHECK_CLOSE(dmudp_diff_u,dmudp[3],reltol);
+ BOOST_CHECK_CLOSE(dmudr_diff_u,dmudr[3],reltol);
+
+ }
+}
+
+void testb(const double reltol, int n, int np, std::vector p, std::vector r,std::vector z,
+ std::vector > props_, std::vector condition)
+{
+ // test b
+ std::vector b(n);
+ std::vector B(n);
+ std::vector invB(n);
+ std::vector dinvBdp(n);
+ std::vector dBdp(n);
+ std::vector dbdr(n);
+ std::vector dbdp(n);
+ double dbdp_diff;
+ double dbdr_diff;
+ double dbdp_diff_u;
+ double dbdr_diff_u;
+
+ for (int phase = 0; phase < np; ++phase) {
+ props_[phase]->b(n, &p[0], &r[0], &condition[0], &b[0], &dbdp[0], &dbdr[0]);
+ props_[phase]->dBdp(n, &p[0], &z[0], &B[0], &dBdp[0]);
+ dbdp_diff = (b[1]-b[0])/(p[1]-p[0]);
+ dbdr_diff = (b[2]-b[0])/(r[2]-r[0]);
+ dbdp_diff_u = (b[4]-b[3])/(p[4]-p[3]);
+ dbdr_diff_u = (b[5]-b[3])/(r[5]-r[3]);
+ for (int i = 0; i < n; ++i){
+ invB[i] = 1/B[i];
+ dinvBdp[i] = -1/pow(B[i],2) * dBdp[i];
+ }
+
+ for (int i = 0; i < n; ++i){
+ BOOST_CHECK_CLOSE(invB[i],b[i] , reltol);
+ BOOST_CHECK_CLOSE(dinvBdp[i],dbdp[i] , reltol);
+
+ }
+ // saturated case
+ BOOST_CHECK_CLOSE(dbdp_diff,dbdp[0], reltol);
+ BOOST_CHECK_CLOSE(dbdr_diff,dbdr[0], reltol);
+
+ // unsaturated case
+ BOOST_CHECK_CLOSE(dbdp_diff_u,dbdp[3], reltol);
+ BOOST_CHECK_CLOSE(dbdr_diff_u,dbdr[3], reltol);
+ }
+}
+
+void testrsSat(double reltol, int n, int np, std::vector p, std::vector > props_){
+ // test bublepoint pressure
+ std::vector rs(n);
+ std::vector drsdp(n);
+ double drsdp_diff;
+ double drsdp_diff_u;
+
+ for (int phase = 0; phase < np; ++phase) {
+ props_[phase] ->rsSat(n, &p[0], &rs[0], &drsdp[0]);
+
+ drsdp_diff = (rs[1]-rs[0])/(p[1]-p[0]);
+ drsdp_diff_u = (rs[4]-rs[3])/(p[4]-p[3]);
+
+ // saturated case
+ BOOST_CHECK_CLOSE(drsdp_diff,drsdp[0], reltol);
+
+ // unsaturad case
+ BOOST_CHECK_CLOSE(drsdp_diff_u,drsdp[3], reltol);
+
+ }
+}
+
+void testrvSat(double reltol, int n, int np, std::vector p, std::vector > props_){
+ // test rv saturated
+ std::vector rv(n);
+ std::vector drvdp(n);
+ double drvdp_diff;
+ double drvdp_diff_u;
+
+ for (int phase = 0; phase < np; ++phase) {
+ props_[phase] ->rvSat(n, &p[0], &rv[0], &drvdp[0]);
+
+ drvdp_diff = (rv[1]-rv[0])/(p[1]-p[0]);
+ drvdp_diff_u = (rv[4]-rv[3])/(p[4]-p[3]);
+
+ // saturated case
+ BOOST_CHECK_CLOSE(drvdp_diff,drvdp[0], reltol);
+
+ // unsaturad case
+ BOOST_CHECK_CLOSE(drvdp_diff_u,drvdp[3], reltol);
+
+ }
+}
+
+BOOST_AUTO_TEST_CASE(test_liveoil)
+{
+
+
+ // read eclipse deck
+ const string filename = "liveoil.DATA";
+ cout << "Reading deck: " << filename << endl;
+ const EclipseGridParser deck (filename);
+
+ // setup pvt interface
+ PhaseUsage phase_usage_ = phaseUsageFromDeck(deck);
+ std::vector > props_ = getProps(deck,phase_usage_);
+
+
// setup a test case. We will check 6 [p,r] pairs and compare them to both the [p,z] interface and a finite difference
// approximation of the derivatives.
const int n = 6;
@@ -97,25 +232,16 @@ BOOST_AUTO_TEST_CASE(test_blackoilfluid)
std::vector p(n);
std::vector r(n);
- std::vector condition(n);
+ std::vector condition(n);
std::vector z(n * np);
- std::vector mu(n);
- std::vector dmudp(n);
- std::vector dmudr(n);
- std::vector mu_new(n);
- double dmudp_diff;
- double dmudr_diff;
- double dmudp_diff_u;
- double dmudr_diff_u;
-
// Used for forward difference calculations
- const double h_p = 100000;
+ const double h_p = 1e4;
const double h_r = 1;
// saturated
- p[0] = 10000000;
+ p[0] = 1e7;
p[1] = p[0] + h_p;
p[2] = p[0];
@@ -144,90 +270,85 @@ BOOST_AUTO_TEST_CASE(test_blackoilfluid)
}
- // test mu
- for (int phase = 1; phase < 2; ++phase) {
- props_[phase]->mu(n, &p[0], &r[0], &condition[0], &mu_new[0], &dmudp[0], &dmudr[0]);
- props_[phase]->mu(n, &p[0], &z[0], &mu[0]);
- dmudp_diff = (mu_new[1]-mu_new[0])/h_p;
- dmudr_diff = (mu_new[2]-mu_new[0])/h_r;
- dmudp_diff_u = (mu_new[4]-mu_new[3])/h_p;
- dmudr_diff_u = (mu_new[5]-mu_new[3])/h_r;
-
- for (int i = 0; i < n; ++i){
- BOOST_CHECK_CLOSE(mu_new[i],mu[i],reltol);
- }
- // saturated case
- BOOST_CHECK_CLOSE(dmudp_diff,dmudp[0],reltol);
- BOOST_CHECK_CLOSE(dmudr_diff,dmudr[0],reltol);
-
- // unsaturated case
- BOOST_CHECK_CLOSE(dmudp_diff_u,dmudp[3] , reltol);
- BOOST_CHECK_CLOSE(dmudr_diff_u,dmudr[3] , reltol);
-
- }
-
- // test b
- std::vector b(n);
- std::vector B(n);
- std::vector invB(n);
- std::vector dinvBdp(n);
- std::vector dBdp(n);
- std::vector dbdr(n);
- std::vector dbdp(n);
- double dbdp_diff;
- double dbdr_diff;
- double dbdp_diff_u;
- double dbdr_diff_u;
-
- for (int phase = 1; phase < 2; ++phase) {
- props_[phase]->b(n, &p[0], &r[0], &condition[0], &b[0], &dbdp[0], &dbdr[0]);
- //props_[phase]->B(n, p, z, B);
- props_[phase]->dBdp(n, &p[0], &z[0], &B[0], &dBdp[0]);
- dbdp_diff = (b[1]-b[0])/h_p;
- dbdr_diff = (b[2]-b[0])/h_r;
- dbdp_diff_u = (b[4]-b[3])/h_p;
- dbdr_diff_u = (b[5]-b[3])/h_r;
- for (int i = 0; i < n; ++i){
- invB[i] = 1/B[i];
- dinvBdp[i] = -1/pow(B[i],2) * dBdp[i];
- }
-
- for (int i = 0; i < n; ++i){
- BOOST_CHECK_CLOSE(invB[i],b[i] , reltol);
- BOOST_CHECK_CLOSE(dinvBdp[i],dbdp[i] , reltol);
-
- }
- // saturated case
- BOOST_CHECK_CLOSE(dbdp_diff,dbdp[0], reltol);
- BOOST_CHECK_CLOSE(dbdr_diff,dbdr[0], reltol);
-
- // unsaturated case
- BOOST_CHECK_CLOSE(dbdp_diff_u,dbdp[3], reltol);
- BOOST_CHECK_CLOSE(dbdr_diff_u,dbdr[3], reltol);
- }
-
- // test bublepoint pressure
- std::vector rbub(n);
- std::vector drbubdp(n);
- double drbubdp_diff;
- double drbubdp_diff_u;
-
- for (int phase = 1; phase < 2; ++phase) {
- props_[phase] ->rbub(n, &p[0], &rbub[0], &drbubdp[0]);
-
- drbubdp_diff = (rbub[1]-rbub[0])/h_p;
- drbubdp_diff_u = (rbub[4]-rbub[3])/h_p;
-
- // saturated case
- BOOST_CHECK_CLOSE(drbubdp_diff,drbubdp[0], reltol);
-
- // unsaturad case
- BOOST_CHECK_CLOSE(drbubdp_diff_u,drbubdp[3], reltol);
-
- }
+ testmu(reltol, n, np, p, r,z, props_, condition);
+ testb(reltol,n,np,p,r,z,props_,condition);
+ testrsSat(reltol,n,np,p,props_);
+ testrvSat(reltol,n,np,p,props_);
}
+
+BOOST_AUTO_TEST_CASE(test_wetgas)
+{
+
+
+ // read eclipse deck
+ const string filename = "wetgas.DATA";
+ cout << "Reading deck: " << filename << endl;
+ const EclipseGridParser deck (filename);
+
+ // setup pvt interface
+ PhaseUsage phase_usage_ = phaseUsageFromDeck(deck);
+ std::vector > props_ = getProps(deck,phase_usage_);
+
+
+ // setup a test case. We will check 6 [p,r] pairs and compare them to both the [p,z] interface and a finite difference
+ // approximation of the derivatives.
+ const int n = 6;
+ const int np = phase_usage_.num_phases;
+
+ // the tolerance for acceptable difference in values
+ const double reltol = 1e-9;
+
+ std::vector p(n);
+ std::vector r(n);
+ std::vector condition(n);
+ std::vector z(n * np);
+
+
+ // Used for forward difference calculations
+ const double h_p = 1e4;
+ const double h_r = 1e-7;
+
+ // saturated
+ p[0] = 1e7;
+ p[1] = p[0] + h_p;
+ p[2] = p[0];
+
+ r[0] = 5e-5;
+ r[1] = 5e-5;
+ r[2] = 5e-5 + h_r;
+
+ condition[0].setFreeOil();
+ condition[1].setFreeOil();
+ condition[2].setFreeOil();
+
+
+ // undersaturated
+ p[3] = p[0];
+ p[4] = p[1];
+ p[5] = p[2];
+
+ r[3] = 1e-5;
+ r[4] = 1e-5;
+ r[5] = 1e-5 +h_r;
+
+ // Corresponing z factors, used to compare with the [p,z] interface
+ for (int i = 0; i < n; ++i) {
+ z[0+i*np] = 0; z[1+i*np] = r[i];
+ z[2+i*np] = 1;
+
+ }
+
+ testmu(reltol, n, np, p, r,z, props_, condition);
+
+ testb(reltol,n,np,p,r,z,props_,condition);
+
+ testrsSat(reltol,n,np,p,props_);
+
+ testrvSat(reltol,n,np,p,props_);
+
+}
diff --git a/tests/wetgas.DATA b/tests/wetgas.DATA
new file mode 100644
index 00000000..4fad316e
--- /dev/null
+++ b/tests/wetgas.DATA
@@ -0,0 +1,71 @@
+WATER
+OIL
+GAS
+
+FIELD
+
+-- PVT PROPERTIES OF DRY GAS (NO VAPOURISED OIL)
+-- FROM SPE3 Blackoil Kleppe
+--
+-- 'Pressure' 'Oil FVF' 'Oil Visc'
+PVDO
+ 1214.7000 1.0632 0.3668
+ 1814.7000 1.0518 0.4241
+ 2414.7000 1.0418 0.5018
+ 3014.7000 1.0332 0.6068
+ 3214.7000 1.0308 0.6461
+ 3364.7000 1.0291 0.6753
+ 3414.7000 1.0285 0.6852
+ 3443.8831 1.0282 0.6912
+/
+
+-- Wet Gas PVT Properties (Vapourised Oil)
+-- Column Properties are:
+-- 'Gas Pressure' 'Gas OGR' 'Gas FVF' 'Gas Visc'
+-- Units: psia stb /Mscf rb /Mscf cp
+PVTG
+ 1214.7000 0.0013130 2.2799 0.0149
+ 0 2.2815 0.01488/
+ 1814.7000 0.00353 1.4401 0.01791
+ 0.001313 1.4429 0.01782
+ 0 1.4445 0.01735 /
+ 2414.7000 0.01102 1.0438 0.02328
+ 0.00353 1.0495 0.02267
+ 0.001313 1.0512 0.0225
+ 0 1.0522 0.02240 /
+ 3014.7000 0.0331 0.8456 0.0318
+ 0.01102 0.8489 0.02924
+ 0.00353 0.8500 0.02844
+ 0.001313 0.8503 0.02820
+ 0 0.8505 0.02807 /
+ 3214.7000 0.0454 0.8082 0.03539
+ 0.0331 0.8080 0.03371
+ 0.01102 0.8075 0.03113
+ 0.00353 0.8073 0.03029
+ 0.001313 0.8073 0.03004
+ 0 0.8073 0.02989 /
+ 3364.7000 0.05670 0.7875 0.0384
+ 0.04540 0.7860 0.03667
+ 0.03310 0.7843 0.03515
+ 0.01102 0.7814 0.03429
+ 0.00353 0.7804 0.03162
+ 0.001313 0.7801 0.03136
+ 0 0.7799 0.03121 /
+ 3416.7575 0.0612 0.7816 0.03955
+ 0.0567 0.7809 0.0386
+ 0.0454 0.7789 0.03717
+ 0.0331 0.7768 0.03564
+ 0.01102 0.7731 0.03296
+ 0.00353 0.7718 0.03207
+ 0.001313 0.7714 0.03181
+ 0 0.7712 0.03166 /
+ 3449.3322 0.0642 0.7783 0.0403
+ 0.0612 0.7777 0.0395
+ 0.0567 0.7769 0.03892
+ 0.0454 0.7747 0.03748
+ 0.0331 0.7723 0.03594
+ 0.01102 0.7681 0.03325
+ 0.00353 0.7666 0.03236
+ 0.001313 0.7662 0.0321
+ 0 0.7660 0.03194 /
+/
|