OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge remote-tracking branch 'upstream/master' into opm-parser-integrate

Browse Source
This commit is contained in:
Joakim Hove 2014-01-29 11:02:21 +01:00
commit bfbc62655d
26 changed files with 694 additions and 185 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
CMakeLists_files.cmake
View File

@ -177,7 +177,8 @@ list (APPEND TEST_SOURCE_FILES
list (APPEND TEST_DATA_FILES
tests/extratestdata.xml
tests/testdata.xml
tests/testFluid.DATA
tests/liveoil.DATA
tests/wetgas.DATA
tests/testBlackoilState1.DATA
tests/testBlackoilState2.DATA
tests/wells_manager_data.data

6
cmake/Modules/FindSuperLU.cmake
View File

@ -54,7 +54,7 @@ endif()
# print message if there was still no blas found!
if(NOT BLAS_FOUND AND NOT SUPERLU_BLAS_LIBRARY)
message("BLAS not found but required for SuperLU")
message(STATUS "BLAS not found but required for SuperLU")
return()
endif()
list(APPEND CMAKE_REQUIRED_LIBRARIES "${SUPERLU_BLAS_LIBRARY}")
@ -69,7 +69,7 @@ if (NOT SUPERLU_INCLUDE_DIR)
)
endif()
if(NOT SUPERLU_INCLUDE_DIR)
message("Directory with the SuperLU include files not found")
message(STATUS "Directory with the SuperLU include files not found")
return()
endif()
list(APPEND CMAKE_REQUIRED_INCLUDES "${SUPERLU_INCLUDE_DIR}")
@ -84,7 +84,7 @@ if (NOT SUPERLU_LIBRARY)
)
endif()
if(NOT SUPERLU_LIBRARY)
message("Directory with the SuperLU library not found")
message(STATUS "Directory with the SuperLU library not found")
return()
endif()
list(APPEND CMAKE_REQUIRED_LIBRARIES "${SUPERLU_LIBRARY}")

1
cmake/Modules/Finddune-common.cmake
View File

@ -50,6 +50,7 @@ int main (void) {
HAVE_ARRAY;
HAVE_BOOST_MAKE_SHARED_HPP;
HAVE_BOOST_SHARED_PTR_HPP;
HAVE_DUNE_BOOST;
HAVE_GMP;
HAVE_MAKE_SHARED;
HAVE_MPI;

1
cmake/Modules/Finddune-istl.cmake
View File

@ -49,6 +49,7 @@ int main (void) {
HAVE_MPI;
HAVE_PARMETIS;
HAVE_SUPERLU;
HAVE_UMFPACK;
SUPERLU_MIN_VERSION_4_3;
SUPERLU_POST_2005_VERSION
")

26
cmake/Modules/OpmAliases.cmake Normal file
View File

@ -0,0 +1,26 @@
# - Alias probed variables for compatibility with DUNE buildsystem
#
# DUNE build system sets some variables which have different names
# in the CMake modules we are using; this module set those variable
# so they can be exported to config.h visible to DUNE headers
function (set_aliases)
# hardcoded list of "dune-var opm-var" pairs, where the components
# are separated by space
set (aliases
"HAVE_UMFPACK HAVE_SUITESPARSE_UMFPACK_H"
"HAVE_DUNE_BOOST HAVE_BOOST"
)
foreach (alias IN LISTS aliases)
# convert entry "X Y" into a list "X;Y", then pick apart
string (REGEX REPLACE "\ +" ";" tuple "${alias}")
list (GET tuple 0 var)
list (GET tuple 1 name)
# write this alias to cache
set (${var} ${${name}} PARENT_SCOPE)
endforeach (alias)
endfunction (set_aliases)
# always call this when the module is imported
set_aliases ()

16
cmake/Modules/OpmFind.cmake
View File

@ -137,16 +137,22 @@ macro (find_and_append_package_to prefix name)
# using config mode is better than using module (aka. find) mode
# because then the package has already done all its probes and
# stored them in the config file for us
if (NOT (${name}_FOUND OR ${NAME}_FOUND))
if (NOT DEFINED ${name}_FOUND AND NOT DEFINED ${NAME}_FOUND)
if (${name}_DIR)
message (STATUS "Finding package ${name} using config mode")
find_package (${name} ${ARGN} NO_MODULE PATHS ${${name}_DIR} NO_DEFAULT_PATH)
else (${name}_DIR)
else ()
message (STATUS "Finding package ${name} using module mode")
find_package (${name} ${ARGN})
endif (${name}_DIR)
endif (NOT (${name}_FOUND OR ${NAME}_FOUND))
endif (CMAKE_DISABLE_FIND_PACKAGE_${name})
endif ()
endif ()
if (NOT DEFINED ${name}_FOUND)
set (${name}_FOUND "${${NAME}_FOUND}")
endif ()
if (NOT DEFINED ${NAME}_FOUND)
set (${NAME}_FOUND "${${name}_FOUND}")
endif ()
endif ()
# the variable "NAME" may be replaced during find_package (as this is
# now a macro, and not a function anymore), so we must reinitialize

3
cmake/Modules/OpmLibMain.cmake
View File

@ -89,6 +89,9 @@ endif (COMMAND prereqs_hook)
include (OpmFind)
find_and_append_package_list_to (${project} ${${project}_DEPS})
# set aliases to probed variables
include (OpmAliases)
# remove the dependency on the testing framework from the main library;
# it is not possible to query for Boost twice with different components.
list (REMOVE_ITEM "${project}_LIBRARIES" "${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}")

6
cmake/Modules/OpmPackage.cmake
View File

@ -344,11 +344,11 @@ macro (find_opm_package module deps header lib defs prog conf)
if ((NOT (${module}_INCLUDE_DIR ${_and_lib_var} AND HAVE_${MODULE}))
AND (_${module}_required OR NOT _${module}_quiet))
if (DEFINED ${module}_DIR)
message ("${module}_DIR = ${${module}_DIR}")
message (STATUS "${module}_DIR = ${${module}_DIR}")
elseif (DEFINED ${module}_ROOT)
message ("${module}_ROOT = ${${module}_ROOT}")
message (STATUS "${module}_ROOT = ${${module}_ROOT}")
elseif (DEFINED ${MODULE}_ROOT)
message ("${MODULE}_ROOT = ${${MODULE}_ROOT}")
message (STATUS "${MODULE}_ROOT = ${${MODULE}_ROOT}")
endif (DEFINED ${module}_DIR)
endif ((NOT (${module}_INCLUDE_DIR ${_and_lib_var} AND HAVE_${MODULE}))
AND (_${module}_required OR NOT _${module}_quiet))

6
cmake/OPM-CMake.md
View File

@ -423,6 +423,12 @@ Wrapper script which emulates an autotools front-end, making the build
system usable with dunecontrol. There is one in the project root directory
which just forwards everything to the main script in `cmake/Scripts`.
<tr>
<td> OpmAliases
<td>
Copy variables which are probed by our find modules to the names which
are expected by DUNE.
<tr>
<td> OpmCompile
<td>

3
opm/core/io/eclipse/EclipseGridParser.cpp
View File

@ -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;

20
opm/core/props/pvt/SinglePvtConstCompr.hpp
View File

@ -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,

19
opm/core/props/pvt/SinglePvtDead.cpp
View File

@ -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,

18
opm/core/props/pvt/SinglePvtDead.hpp
View File

@ -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,

19
opm/core/props/pvt/SinglePvtDeadSpline.cpp
View File

@ -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,

18
opm/core/props/pvt/SinglePvtDeadSpline.hpp
View File

@ -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,

18
opm/core/props/pvt/SinglePvtInterface.hpp
View File

@ -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.

156
opm/core/props/pvt/SinglePvtLiveGas.cpp
View File

@ -61,7 +61,7 @@ namespace Opm
for (int i=0; i<sz; ++i) {
saturated_gas_table_[0][i] = pvtg[region_number][i][0]; // p
saturated_gas_table_[1][i] = pvtg[region_number][i][2]; // Bg
saturated_gas_table_[1][i] = 1.0/pvtg[region_number][i][2]; // 1/Bg
saturated_gas_table_[2][i] = pvtg[region_number][i][3]; // mu_g
saturated_gas_table_[3][i] = pvtg[region_number][i][1]; // Rv
}
@ -75,7 +75,7 @@ namespace Opm
undersat_gas_tables_[i][2].resize(tsize);
for (int j=0, k=0; j<tsize; ++j) {
undersat_gas_tables_[i][0][j] = pvtg[region_number][i][++k]; // Rv
undersat_gas_tables_[i][1][j] = pvtg[region_number][i][++k]; // Bg
undersat_gas_tables_[i][1][j] = 1.0/pvtg[region_number][i][++k]; // 1/Bg
undersat_gas_tables_[i][2][j] = pvtg[region_number][i][++k]; // mu_g
}
}
@ -118,7 +118,13 @@ namespace Opm
double* output_dmudp,
double* output_dmudr) const
{
OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented");
for (int i = 0; i < n; ++i) {
const PhasePresence& cnd = cond[i];
output_mu[i] = miscible_gas(p[i], r[i], cnd,2, 0);
output_dmudp[i] = miscible_gas(p[i], r[i], cnd, 2, 1);
output_dmudr[i] = miscible_gas(p[i], r[i], cnd, 2, 2);
}
}
@ -171,16 +177,39 @@ namespace Opm
double* output_dbdr) const
{
OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented");
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
const PhasePresence& cnd = cond[i];
output_b[i] = miscible_gas(p[i], r[i], cnd, 1, 0);
output_dbdp[i] = miscible_gas(p[i], r[i], cnd, 1, 1);
output_dbdr[i] = miscible_gas(p[i], r[i], cnd, 1, 2);
}
}
/// Gas resolution and its derivatives at bublepoint as a function of p.
void SinglePvtLiveGas::rbub(const int n,
void SinglePvtLiveGas::rvSat(const int n,
const double* p,
double* output_rbub,
double* output_drbubdp) const
double* output_rvSat,
double* output_drvSatdp) const
{
OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented");
for (int i = 0; i < n; ++i) {
output_rvSat[i] = linearInterpolation(saturated_gas_table_[0],
saturated_gas_table_[3],p[i]);
output_drvSatdp[i] = linearInterpolationDerivative(saturated_gas_table_[0],
saturated_gas_table_[3],p[i]);
}
}
void SinglePvtLiveGas::rsSat(const int n,
const double* /*p*/,
double* output_rsSat,
double* output_drsSatdp) const
{
std::fill(output_rsSat, output_rsSat + n, 0.0);
std::fill(output_drsSatdp, output_drsSatdp + n, 0.0);
}
/// Solution factor as a function of p and z.
@ -219,7 +248,7 @@ namespace Opm
// To handle no-gas case.
return 1.0;
}
return miscible_gas(press, surfvol, 1, false);
return 1.0/miscible_gas(press, surfvol, 1, false);
}
void SinglePvtLiveGas::evalBDeriv(const double press, const double* surfvol,
@ -231,8 +260,8 @@ namespace Opm
dBdpval = 0.0;
return;
}
Bval = miscible_gas(press, surfvol, 1, false);
dBdpval = miscible_gas(press, surfvol, 1, true);
Bval = evalB(press, surfvol);
dBdpval = -Bval*Bval*miscible_gas(press, surfvol, 1, true);
}
double SinglePvtLiveGas::evalR(const double press, const double* surfvol) const
@ -360,5 +389,110 @@ namespace Opm
}
}
double SinglePvtLiveGas::miscible_gas(const double press,
const double r,
const PhasePresence& cond,
const int item,
const int deriv) const
{
const bool isSat = cond.hasFreeOil();
// Derivative w.r.t p
if (deriv == 1) {
if (isSat) { // Saturated case
return linearInterpolationDerivative(saturated_gas_table_[0],
saturated_gas_table_[item],
press);
} else { // Undersaturated case
int is = tableIndex(saturated_gas_table_[0], press);
if (undersat_gas_tables_[is][0].size() < 2) {
double val = (saturated_gas_table_[item][is+1]
- saturated_gas_table_[item][is]) /
(saturated_gas_table_[0][is+1] -
saturated_gas_table_[0][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 = (val2 - val1)/
(saturated_gas_table_[0][is+1] - saturated_gas_table_[0][is]);
return val;
}
} else if (deriv == 2){
if (isSat) {
return 0;
} else {
int is = tableIndex(saturated_gas_table_[0], press);
double w = (press - saturated_gas_table_[0][is]) /
(saturated_gas_table_[0][is+1] - saturated_gas_table_[0][is]);
assert(undersat_gas_tables_[is][0].size() >= 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

23
opm/core/props/pvt/SinglePvtLiveGas.hpp
View File

@ -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<std::vector<double> > saturated_gas_table_;
std::vector<std::vector<std::vector<double> > > undersat_gas_tables_;

19
opm/core/props/pvt/SinglePvtLiveOil.cpp
View File

@ -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,

18
opm/core/props/pvt/SinglePvtLiveOil.hpp
View File

@ -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,

1
opm/core/simulator/BlackoilState.cpp
View File

@ -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
}

3
opm/core/simulator/BlackoilState.hpp
View File

@ -47,13 +47,16 @@ namespace Opm
std::vector<double>& surfacevol () { return surfvol_; }
std::vector<double>& gasoilratio () { return gor_ ; }
std::vector<double>& rv () {return rv_ ; }
const std::vector<double>& surfacevol () const { return surfvol_; }
const std::vector<double>& gasoilratio () const { return gor_ ; }
const std::vector<double>& rv () const {return rv_ ; }
private:
std::vector<double> surfvol_;
std::vector<double> gor_ ;
std::vector<double> rv_ ;
};
} // namespace Opm

42
opm/core/simulator/initState_impl.hpp
View File

@ -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 <class Props, class State>
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<double>& rs = state.gasoilratio();
const std::vector<double>& 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<double> allA_a(nc*np*np);
std::vector<double> allA_l(nc*np*np);
std::vector<double> allA_v(nc*np*np);
@ -641,6 +641,17 @@ namespace Opm
allcells[c] = c;
}
std::vector<double> capPressures(nc*np);
props.capPress(nc,&state.saturation()[0],&allcells[0],&capPressures[0],NULL);
std::vector<double> Pw(nc);
std::vector<double> 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<double>& 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.");
}
}

36
tests/liveoil.DATA Normal file
View File

@ -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 /

321
tests/test_blackoilfluid.cpp
View File

@ -28,24 +28,18 @@
using namespace Opm;
using namespace std;
BOOST_AUTO_TEST_CASE(test_blackoilfluid)
{
std::vector<std::shared_ptr<SinglePvtInterface> > 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<std::shared_ptr<SinglePvtInterface> > props_;
PhaseUsage phase_usage_ = phaseUsageFromDeck(deck);
enum PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2 };
int samples = 0;
std::vector<std::shared_ptr<SinglePvtInterface> > 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<double> p, std::vector<double> r,std::vector<double> z,
std::vector<std::shared_ptr<SinglePvtInterface> > props_, std::vector<Opm::PhasePresence> condition)
{
std::vector<double> mu(n);
std::vector<double> dmudp(n);
std::vector<double> dmudr(n);
std::vector<double> 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<double> p, std::vector<double> r,std::vector<double> z,
std::vector<std::shared_ptr<SinglePvtInterface> > props_, std::vector<Opm::PhasePresence> condition)
{
// test b
std::vector<double> b(n);
std::vector<double> B(n);
std::vector<double> invB(n);
std::vector<double> dinvBdp(n);
std::vector<double> dBdp(n);
std::vector<double> dbdr(n);
std::vector<double> 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<double> p, std::vector<std::shared_ptr<SinglePvtInterface> > props_){
// test bublepoint pressure
std::vector<double> rs(n);
std::vector<double> 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<double> p, std::vector<std::shared_ptr<SinglePvtInterface> > props_){
// test rv saturated
std::vector<double> rv(n);
std::vector<double> 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<std::shared_ptr<SinglePvtInterface> > 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<double> p(n);
std::vector<double> r(n);
std::vector<Opm::PhasePresence> condition(n);
std::vector<PhasePresence> condition(n);
std::vector<double> z(n * np);
std::vector<double> mu(n);
std::vector<double> dmudp(n);
std::vector<double> dmudr(n);
std::vector<double> 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;
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<std::shared_ptr<SinglePvtInterface> > 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<double> p(n);
std::vector<double> r(n);
std::vector<PhasePresence> condition(n);
std::vector<double> 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) {
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);
z[0+i*np] = 0; z[1+i*np] = r[i];
z[2+i*np] = 1;
}
// test b
std::vector<double> b(n);
std::vector<double> B(n);
std::vector<double> invB(n);
std::vector<double> dinvBdp(n);
std::vector<double> dBdp(n);
std::vector<double> dbdr(n);
std::vector<double> 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<double> rbub(n);
std::vector<double> 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_);
}

71
tests/wetgas.DATA Normal file
View File

@ -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 /
/