OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Refactoring and updated tests

Browse Source
This commit is contained in:
babrodtk 2015-07-09 10:24:47 +02:00
parent 0467af953c
commit 79410685ca
4 changed files with 233 additions and 139 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
opm/autodiff/BlackoilModelBase_impl.hpp
View File

@ -1379,7 +1379,7 @@ namespace detail {
const double& thp = well_controls_iget_target(wc, current);
const double& alq = well_controls_iget_alq(wc, current);
bhp_targets (w) = vfp_properties_->prod_bhp(vfp, aqua[w], liquid[w], vapour[w], thp, alq);
bhp_targets (w) = vfp_properties_->getProd()->bhp(vfp, aqua[w], liquid[w], vapour[w], thp, alq);
rate_targets(w) = -1e100;
}
break;
@ -1787,7 +1787,7 @@ namespace detail {
double alq = well_controls_iget_alq(wc, ctrl_index);
int table_id = well_controls_iget_vfp(wc, ctrl_index);
well_state.thp()[w] = vfp_properties_->prod_thp(table_id, aqua, liquid, vapour, bhp[w], alq);
well_state.thp()[w] = vfp_properties_->getProd()->thp(table_id, aqua, liquid, vapour, bhp[w], alq);
//Assume only one THP control specified for each well
break;

78
opm/autodiff/VFPProperties.cpp
View File

@ -34,44 +34,72 @@ VFPProperties::VFPProperties() {
VFPProperties::VFPProperties(const VFPInjTable* inj_table, const VFPProdTable* prod_table) {
if (inj_table != NULL) {
m_inj_tables[inj_table->getTableNum()] = inj_table;
//FIXME: Implement VFPInjProperties
OPM_THROW(std::logic_error, "VFPInjProperties not implemented yet");
}
if (prod_table != NULL) {
m_prod_tables[prod_table->getTableNum()] = prod_table;
m_prod.reset(new VFPProdProperties(prod_table));
}
}
VFPProperties::VFPProperties(const std::map<int, VFPInjTable>& inj_tables,
const std::map<int, VFPProdTable>& prod_tables) {
init(inj_tables);
init(prod_tables);
//FIXME: Implement VFPInjProperties
OPM_THROW(std::logic_error, "VFPInjProperties not implemented yet");
m_prod.reset(new VFPProdProperties(prod_tables));
}
VFPProperties::VFPProperties(const std::map<int, VFPInjTable>& inj_tables) {
init(inj_tables);
//FIXME: Implement VFPInjProperties
OPM_THROW(std::logic_error, "VFPInjProperties not implemented yet");
}
VFPProperties::VFPProperties(const std::map<int, VFPProdTable>& prod_tables) {
init(prod_tables);
m_prod.reset(new VFPProdProperties(prod_tables));
}
void VFPProperties::init(const std::map<int, VFPInjTable>& inj_tables) {
//Populate injection table pointers.
for (const auto& table : inj_tables) {
m_inj_tables[table.first] = &table.second;
}
VFPProdProperties::VFPProdProperties() {
}
void VFPProperties::init(const std::map<int, VFPProdTable>& prod_tables) {
VFPProdProperties::VFPProdProperties(const VFPProdTable* table){
m_tables[table->getTableNum()] = table;
}
VFPProdProperties::VFPProdProperties(const std::map<int, VFPProdTable>& tables) {
init(tables);
}
void VFPProdProperties::init(const std::map<int, VFPProdTable>& prod_tables) {
//Populate production table pointers
for (const auto& table : prod_tables) {
m_prod_tables[table.first] = &table.second;
m_tables[table.first] = &table.second;
}
}
VFPProperties::ADB::V VFPProperties::prod_bhp(int table_id,
VFPProdProperties::ADB::V VFPProdProperties::bhp(int table_id,
const Wells& wells,
const ADB::V& qs,
const ADB::V& thp,
@ -86,10 +114,10 @@ VFPProperties::ADB::V VFPProperties::prod_bhp(int table_id,
const ADB::V& o = subset(qs, Span(nw, 1, BlackoilPhases::Liquid*nw));
const ADB::V& g = subset(qs, Span(nw, 1, BlackoilPhases::Vapour*nw));
return prod_bhp(table_id, w, o, g, thp, alq);
return bhp(table_id, w, o, g, thp, alq);
}
VFPProperties::ADB::V VFPProperties::prod_bhp(int table_id,
VFPProdProperties::ADB::V VFPProdProperties::bhp(int table_id,
const ADB::V& aqua,
const ADB::V& liquid,
const ADB::V& vapour,
@ -107,12 +135,12 @@ VFPProperties::ADB::V VFPProperties::prod_bhp(int table_id,
ADB::V bhp_vals;
bhp_vals.resize(nw);
for (int i=0; i<nw; ++i) {
bhp_vals[i] = prod_bhp(table_id, aqua[i], liquid[i], vapour[i], thp[i], alq[i]);
bhp_vals[i] = bhp(table_id, aqua[i], liquid[i], vapour[i], thp[i], alq[i]);
}
return bhp_vals;
}
double VFPProperties::prod_bhp(int table_id,
double VFPProdProperties::bhp(int table_id,
const double& aqua,
const double& liquid,
const double& vapour,
@ -137,7 +165,7 @@ double VFPProperties::prod_bhp(int table_id,
}
double VFPProperties::prod_thp(int table_id,
double VFPProdProperties::thp(int table_id,
const double& aqua,
const double& liquid,
const double& vapour,
@ -279,9 +307,9 @@ double VFPProperties::prod_thp(int table_id,
const VFPProdTable* VFPProperties::getProdTable(int table_id) const {
auto entry = m_prod_tables.find(table_id);
if (entry == m_prod_tables.end()) {
const VFPProdTable* VFPProdProperties::getProdTable(int table_id) const {
auto entry = m_tables.find(table_id);
if (entry == m_tables.end()) {
OPM_THROW(std::invalid_argument, "Nonexistent table " << table_id << " referenced.");
}
else {
@ -289,7 +317,7 @@ const VFPProdTable* VFPProperties::getProdTable(int table_id) const {
}
}
VFPProperties::InterpData VFPProperties::find_interp_data(const double& value, const std::vector<double>& values) {
VFPProdProperties::InterpData VFPProdProperties::find_interp_data(const double& value, const std::vector<double>& values) {
InterpData retval;
//First element greater than or equal to value
@ -330,7 +358,7 @@ VFPProperties::InterpData VFPProperties::find_interp_data(const double& value, c
#pragma GCC optimize ("unroll-loops")
#endif
double VFPProperties::interpolate(const VFPProdTable::array_type& array,
double VFPProdProperties::interpolate(const VFPProdTable::array_type& array,
const InterpData& flo_i,
const InterpData& thp_i,
const InterpData& wfr_i,
@ -408,7 +436,7 @@ double VFPProperties::interpolate(const VFPProdTable::array_type& array,
#endif
double VFPProperties::find_x(const double& x0,
double VFPProdProperties::find_x(const double& x0,
const double& x1,
const double& y0,
const double& y1,

72
opm/autodiff/VFPProperties.hpp
View File

@ -31,18 +31,15 @@
namespace Opm {
class VFPProdProperties;
class VFPInjProperties;
/**
* Class which linearly interpolates BHP as a function of rate, tubing head pressure,
* water fraction, gas fraction, and artificial lift for production VFP tables, and similarly
* the BHP as a function of the rate and tubing head pressure.
* A thin wrapper class that holds one VFPProdProperties and one
* VFPInjProperties object.
*/
class VFPProperties {
public:
typedef AutoDiffBlock<double> ADB;
/**
* Empty constructor
*/
VFPProperties();
/**
@ -76,6 +73,47 @@ public:
*/
VFPProperties(const std::map<int, VFPProdTable>& prod_tables);
const VFPInjProperties* getInj() const {
return m_inj.get();
}
const VFPProdProperties* getProd() const {
return m_prod.get();
}
private:
std::shared_ptr<VFPInjProperties> m_inj;
std::shared_ptr<VFPProdProperties> m_prod;
};
/**
* Class which linearly interpolates BHP as a function of rate, tubing head pressure,
* water fraction, gas fraction, and artificial lift for production VFP tables, and similarly
* the BHP as a function of the rate and tubing head pressure.
*/
class VFPProdProperties {
public:
typedef AutoDiffBlock<double> ADB;
/**
* Empty constructor
*/
VFPProdProperties();
/**
* Constructor
* Takes *no* ownership of data.
* @param prod_table A *single* VFPPROD table
*/
VFPProdProperties(const VFPProdTable* prod_table);
/**
* Constructor
* Takes *no* ownership of data.
* @param prod_tables A map of different VFPPROD tables.
*/
VFPProdProperties(const std::map<int, VFPProdTable>& prod_tables);
/**
* Linear interpolation of bhp as function of the input parameters.
* @param table_id Table number to use
@ -87,7 +125,7 @@ public:
* @return The bottom hole pressure, interpolated/extrapolated linearly using
* the above parameters from the values in the input table.
*/
ADB::V prod_bhp(int table_id,
ADB::V bhp(int table_id,
const Wells& wells,
const ADB::V& qs,
const ADB::V& thp,
@ -106,7 +144,7 @@ public:
* the above parameters from the values in the input table, for each entry in the
* input ADB objects.
*/
ADB::V prod_bhp(int table_id,
ADB::V bhp(int table_id,
const ADB::V& aqua,
const ADB::V& liquid,
const ADB::V& vapour,
@ -125,7 +163,7 @@ public:
* @return The bottom hole pressure, interpolated/extrapolated linearly using
* the above parameters from the values in the input table.
*/
double prod_bhp(int table_id,
double bhp(int table_id,
const double& aqua,
const double& liquid,
const double& vapour,
@ -144,7 +182,7 @@ public:
* @return The tubing hole pressure, interpolated/extrapolated linearly using
* the above parameters from the values in the input table.
*/
double prod_thp(int table_id,
double thp(int table_id,
const double& aqua,
const double& liquid,
const double& vapour,
@ -188,8 +226,8 @@ public:
//Water-oil ratio = water / oil
return aqua / liquid;
case VFPProdTable::WFR_WCT:
//Water cut = water / (water + oil + gas)
return aqua / (aqua + liquid + vapour);
//Water cut = water / (water + oil)
return aqua / (aqua + liquid);
case VFPProdTable::WFR_WGR:
//Water-gas ratio = water / gas
return aqua / vapour;
@ -225,8 +263,7 @@ public:
private:
// Map which connects the table number with the table itself
std::map<int, const VFPProdTable*> m_prod_tables;
std::map<int, const VFPInjTable*> m_inj_tables;
std::map<int, const VFPProdTable*> m_tables;
/**
* Helper struct for linear interpolation
@ -264,9 +301,8 @@ private:
/**
* Initialization routines
* Initialization routine
*/
void init(const std::map<int, VFPInjTable>& inj_tables);
void init(const std::map<int, VFPProdTable>& prod_tables);
/**

218
tests/test_vfpproperties.cpp
View File

@ -43,8 +43,8 @@
const double max_d_tol = 1.0e-12;
const double sad_tol = 1.0e-9;
const double max_d_tol = 1.0e-10;
const double sad_tol = 1.0e-8;
@ -54,7 +54,7 @@ const double sad_tol = 1.0e-9;
* values data is given at
*/
struct TrivialFixture {
typedef Opm::VFPProperties::ADB ADB;
typedef Opm::VFPProdProperties::ADB ADB;
TrivialFixture() : thp_axis{0.0, 1.0},
wfr_axis{0.0, 0.5, 1.0},
@ -108,7 +108,7 @@ struct TrivialFixture {
double u = l / static_cast<double>(nu-1);
for (int m=0; m<nv; ++m) {
double v = m / static_cast<double>(nv-1);
// table[thp_idx][wfr_idx][gfr_idx][alq_idx][flo_idx];
data[i][j][k][l][m] = x + 2*y + 3*z + 4*u + 5*v;
}
}
@ -158,10 +158,10 @@ struct TrivialFixture {
data);
//Initialize properties that use the table
properties.reset(new Opm::VFPProperties(NULL, &table));
properties.reset(new Opm::VFPProdProperties(&table));
}
std::shared_ptr<Opm::VFPProperties> properties;
std::shared_ptr<Opm::VFPProdProperties> properties;
Opm::VFPProdTable table;
private:
@ -194,10 +194,12 @@ BOOST_AUTO_TEST_CASE(GetTable)
//Create wells
const int nphases = 3;
const int nwells = 5;
const int nwells = 1;
const int nperfs = 1;
std::shared_ptr<Wells> wells(create_wells(nphases, nwells, nperfs),
destroy_wells);
const int cells[] = {5};
add_well(INJECTOR, 100, 1, NULL, cells, NULL, NULL, wells.get());
//Create interpolation points
double aqua_d = 1.5;
@ -212,22 +214,22 @@ BOOST_AUTO_TEST_CASE(GetTable)
ADB::V thp_adb = ADB::V::Constant(1, thp_d);
ADB::V alq_adb = ADB::V::Constant(1, alq_d);
ADB::V qs_adb;
ADB::V qs_adb(3);
qs_adb << aqua_adb, liquid_adb, vapour_adb;
//Check that different versions of the prod_bph function work
ADB::V a = properties->prod_bhp(1, aqua_adb, liquid_adb, vapour_adb, thp_adb, alq_adb);
double b = properties->prod_bhp(1, aqua_d, liquid_d, vapour_d, thp_d, alq_d);
ADB::V c = properties->prod_bhp(1, *wells, qs_adb, thp_adb, alq_adb);
ADB::V a = properties->bhp(1, aqua_adb, liquid_adb, vapour_adb, thp_adb, alq_adb);
double b = properties->bhp(1, aqua_d, liquid_d, vapour_d, thp_d, alq_d);
ADB::V c = properties->bhp(1, *wells, qs_adb, thp_adb, alq_adb);
//Check that results are actually equal
BOOST_CHECK_EQUAL(a[0], b);
BOOST_CHECK_EQUAL(a[0], c[0]);
//Table 2 does not exist.
BOOST_CHECK_THROW(properties->prod_bhp(2, *wells, qs_adb, thp_adb, alq_adb), std::invalid_argument);
BOOST_CHECK_THROW(properties->bhp(2, *wells, qs_adb, thp_adb, alq_adb), std::invalid_argument);
}
/**
@ -242,11 +244,11 @@ BOOST_AUTO_TEST_CASE(InterpolateZero)
//Check interpolation
double sum = 0.0;
int n=5;
for (int i=0; i<n; ++i) {
for (int i=1; i<n; ++i) {
const double x = i / static_cast<double>(n-1);
for (int j=0; j<n; ++j) {
for (int j=1; j<n; ++j) {
const double y = j / static_cast<double>(n-1);
for (int k=0; k<n; ++k) {
for (int k=1; k<n; ++k) {
const double z = k / static_cast<double>(n-1);
for (int l=0; l<n; ++l) {
const double u = l / static_cast<double>(n-1);
@ -254,7 +256,7 @@ BOOST_AUTO_TEST_CASE(InterpolateZero)
const double v = m / static_cast<double>(n-1);
//Note order of arguments!
sum += properties->prod_bhp(1, v, x, y, z, u);
sum += properties->bhp(1, v, x, y, z, u);
}
}
}
@ -277,11 +279,11 @@ BOOST_AUTO_TEST_CASE(InterpolateOne)
//Check interpolation
double sum = 0.0;
int n=5;
for (int i=0; i<n; ++i) {
for (int i=1; i<n; ++i) {
const double x = i / static_cast<double>(n-1);
for (int j=0; j<n; ++j) {
for (int j=1; j<n; ++j) {
const double y = j / static_cast<double>(n-1);
for (int k=0; k<n; ++k) {
for (int k=1; k<n; ++k) {
const double z = k / static_cast<double>(n-1);
for (int l=0; l<n; ++l) {
const double u = l / static_cast<double>(n-1);
@ -289,14 +291,14 @@ BOOST_AUTO_TEST_CASE(InterpolateOne)
const double v = m / static_cast<double>(n-1);
//Note order of arguments!
sum += properties->prod_bhp(1, v, x, y, z, u);
sum += properties->bhp(1, v, x, y, z, u);
}
}
}
}
}
double reference = n*n*n*n*n;
double reference = (n-1)*(n-1)*(n-1)*n*n;
BOOST_CHECK_EQUAL(sum, reference);
}
@ -318,19 +320,25 @@ BOOST_AUTO_TEST_CASE(InterpolatePlane)
int n=5;
for (int i=0; i<=n; ++i) {
const double x = i / static_cast<double>(n);
for (int j=0; j<=n; ++j) {
const double y = j / static_cast<double>(n);
for (int k=0; k<=n; ++k) {
const double z = k / static_cast<double>(n);
for (int j=1; j<=n; ++j) {
const double aqua = j / static_cast<double>(n);
for (int k=1; k<=n; ++k) {
const double vapour = k / static_cast<double>(n);
for (int l=0; l<=n; ++l) {
const double u = l / static_cast<double>(n);
for (int m=0; m<=n; ++m) {
const double v = m / static_cast<double>(n);
double reference = x + 2*y + 3*z + 4*u + 5*v;
for (int m=1; m<=n; ++m) {
const double liquid = m / static_cast<double>(n);
//Find values that should be in table
double v = properties->getFlo(aqua, liquid, vapour, table.getFloType());
double y = properties->getWFR(aqua, liquid, vapour, table.getWFRType());
double z = properties->getGFR(aqua, liquid, vapour, table.getGFRType());
double reference = x + 2*y + 3*z+ 4*u + 5*v;
reference_sum += reference;
//Note order of arguments!
double value = properties->prod_bhp(1, v, x, y, z, u);
//Note order of arguments! id, aqua, liquid, vapour, thp , alq
double value = properties->bhp(1, aqua, liquid, vapour, x, u);
sum += value;
double abs_diff = std::abs(value - reference);
@ -365,21 +373,28 @@ BOOST_AUTO_TEST_CASE(ExtrapolatePlane)
double sad = 0.0; // Sum absolute difference
double max_d = 0.0; // Maximum difference
int n=1;
for (int i=-5; i<=n+5; ++i) {
int o=5;
for (int i=0; i<=n+o; ++i) {
const double x = i / static_cast<double>(n);
for (int j=-5; j<=n+5; ++j) {
const double y = j / static_cast<double>(n);
for (int k=-5; k<=n+5; ++k) {
const double z = k / static_cast<double>(n);
for (int l=-5; l<=n+5; ++l) {
for (int j=1; j<=n+o; ++j) {
const double aqua = j / static_cast<double>(n);
for (int k=1; k<=n+o; ++k) {
const double vapour = k / static_cast<double>(n);
for (int l=0; l<=n+o; ++l) {
const double u = l / static_cast<double>(n);
for (int m=-5; m<=n+5; ++m) {
const double v = m / static_cast<double>(n);
double reference = x + 2*y + 3*z + 4*u + 5*v;
for (int m=1; m<=n+o; ++m) {
const double liquid = m / static_cast<double>(n);
//Find values that should be in table
double v = properties->getFlo(aqua, liquid, vapour, table.getFloType());
double y = properties->getWFR(aqua, liquid, vapour, table.getWFRType());
double z = properties->getGFR(aqua, liquid, vapour, table.getGFRType());
double reference = x + 2*y + 3*z+ 4*u + 5*v;
reference_sum += reference;
//Note order of arguments!
double value = properties->prod_bhp(1, v, x, y, z, u);
//Note order of arguments! id, aqua, liquid, vapour, thp , alq
double value = properties->bhp(1, aqua, liquid, vapour, x, u);
sum += value;
double abs_diff = std::abs(value - reference);
@ -411,11 +426,11 @@ BOOST_AUTO_TEST_CASE(ExtrapolatePlaneADB)
//Temporary variables used to represent independent wells
const int num_wells = 5;
ADB::V thp(num_wells);
ADB::V wfr(num_wells);
ADB::V gfr(num_wells);
ADB::V alq(num_wells);
ADB::V flo(num_wells);
ADB::V x_v(num_wells);
ADB::V aqua_v(num_wells);
ADB::V vapour_v(num_wells);
ADB::V u_v(num_wells);
ADB::V liquid_v(num_wells);
//Check linear extrapolation (i.e., using values of x, y, etc. outside our interpolant domain)
double sum = 0.0;
@ -423,30 +438,37 @@ BOOST_AUTO_TEST_CASE(ExtrapolatePlaneADB)
double sad = 0.0; // Sum absolute difference
double max_d = 0.0; // Maximum difference
int n=1;
for (int i=-5; i<=n+5; ++i) {
int o=5;
for (int i=0; i<=n+o; ++i) {
const double x = i / static_cast<double>(n);
for (int j=-5; j<=n+5; ++j) {
const double y = j / static_cast<double>(n);
for (int k=-5; k<=n+5; ++k) {
const double z = k / static_cast<double>(n);
for (int l=-5; l<=n+5; ++l) {
for (int j=1; j<=n+o; ++j) {
const double aqua = j / static_cast<double>(n);
for (int k=1; k<=n+o; ++k) {
const double vapour = k / static_cast<double>(n);
for (int l=0; l<=n+o; ++l) {
const double u = l / static_cast<double>(n);
for (int m=-5; m<=n+5; ++m) {
const double v = m / static_cast<double>(n);
for (int m=1; m<=n+o; ++m) {
const double liquid = m / static_cast<double>(n);
for (unsigned int w=0; w<num_wells; ++w) {
thp[w] = x*w;
wfr[w] = y*w;
gfr[w] = z*w;
alq[w] = u*w;
flo[w] = v*w;
x_v[w] = x*(w+1);
aqua_v[w] = aqua*(w+1);
vapour_v[w] = vapour*(w+1);
u_v[w] = u*(w+1);
liquid_v[w] = liquid*(w+1);
}
ADB::V bhp_val = properties->prod_bhp(1, flo, thp, wfr, gfr, alq);
ADB::V bhp_val = properties->bhp(1, aqua_v, liquid_v, vapour_v, x_v, u_v);
double value = 0.0;
double reference = 0.0;
for (int w=0; w < num_wells; ++w) {
reference = x*w + 2*y*w + 3*z*w + 4*u*w + 5*v*w;
//Find values that should be in table
double v = properties->getFlo(aqua*(w+1), liquid*(w+1), vapour*(w+1), table.getFloType());
double y = properties->getWFR(aqua*(w+1), liquid*(w+1), vapour*(w+1), table.getWFRType());
double z = properties->getGFR(aqua*(w+1), liquid*(w+1), vapour*(w+1), table.getGFRType());
reference = x*(w+1) + 2*y + 3*z + 4*u*(w+1) + 5*v;
value = bhp_val[w];
sum += value;
@ -518,7 +540,7 @@ BOOST_AUTO_TEST_CASE(InterpolateADBAndQs)
}
//Call the bhp function
ADB::V bhp = properties->prod_bhp(1, *wells, qs, thp, alq);
ADB::V bhp = properties->bhp(1, *wells, qs, thp, alq);
//Calculate reference
//First, find the three phases
@ -578,7 +600,7 @@ BOOST_AUTO_TEST_SUITE_END() // Trivial tests
struct ConversionFixture {
typedef Opm::VFPProperties::ADB ADB;
typedef Opm::VFPProdProperties::ADB ADB;
ConversionFixture() :
num_wells(5),
@ -625,19 +647,19 @@ BOOST_AUTO_TEST_CASE(getFlo)
}
{
ADB::V flo = Opm::VFPProperties::getFlo(aqua, liquid, vapour, Opm::VFPProdTable::FLO_OIL);
ADB::V flo = Opm::VFPProdProperties::getFlo(aqua, liquid, vapour, Opm::VFPProdTable::FLO_OIL);
const double* computed = &flo[0];
BOOST_CHECK_EQUAL_COLLECTIONS(ref_flo_oil.begin(), ref_flo_oil.end(), computed, computed+num_wells);
}
{
ADB::V flo = Opm::VFPProperties::getFlo(aqua, liquid, vapour, Opm::VFPProdTable::FLO_LIQ);
ADB::V flo = Opm::VFPProdProperties::getFlo(aqua, liquid, vapour, Opm::VFPProdTable::FLO_LIQ);
const double* computed = &flo[0];
BOOST_CHECK_EQUAL_COLLECTIONS(ref_flo_liq.begin(), ref_flo_liq.end(), computed, computed+num_wells);
}
{
ADB::V flo = Opm::VFPProperties::getFlo(aqua, liquid, vapour, Opm::VFPProdTable::FLO_GAS);
ADB::V flo = Opm::VFPProdProperties::getFlo(aqua, liquid, vapour, Opm::VFPProdTable::FLO_GAS);
const double* computed = &flo[0];
BOOST_CHECK_EQUAL_COLLECTIONS(ref_flo_gas.begin(), ref_flo_gas.end(), computed, computed+num_wells);
}
@ -652,24 +674,24 @@ BOOST_AUTO_TEST_CASE(getWFR)
std::vector<double> ref_wfr_wgr(num_wells);
for (int i=0; i<num_wells; ++i) {
ref_wfr_wor[i] = aqua[i] / liquid[i];
ref_wfr_wct[i] = aqua[i] / (aqua[i] + liquid[i] + vapour[i]);
ref_wfr_wct[i] = aqua[i] / (aqua[i] + liquid[i]);
ref_wfr_wgr[i] = aqua[i] / vapour[i];
}
{
ADB::V flo = Opm::VFPProperties::getWFR(aqua, liquid, vapour, Opm::VFPProdTable::WFR_WOR);
ADB::V flo = Opm::VFPProdProperties::getWFR(aqua, liquid, vapour, Opm::VFPProdTable::WFR_WOR);
const double* computed = &flo[0];
BOOST_CHECK_EQUAL_COLLECTIONS(ref_wfr_wor.begin(), ref_wfr_wor.end(), computed, computed+num_wells);
}
{
ADB::V flo = Opm::VFPProperties::getWFR(aqua, liquid, vapour, Opm::VFPProdTable::WFR_WCT);
ADB::V flo = Opm::VFPProdProperties::getWFR(aqua, liquid, vapour, Opm::VFPProdTable::WFR_WCT);
const double* computed = &flo[0];
BOOST_CHECK_EQUAL_COLLECTIONS(ref_wfr_wct.begin(), ref_wfr_wct.end(), computed, computed+num_wells);
}
{
ADB::V flo = Opm::VFPProperties::getWFR(aqua, liquid, vapour, Opm::VFPProdTable::WFR_WGR);
ADB::V flo = Opm::VFPProdProperties::getWFR(aqua, liquid, vapour, Opm::VFPProdTable::WFR_WGR);
const double* computed = &flo[0];
BOOST_CHECK_EQUAL_COLLECTIONS(ref_wfr_wgr.begin(), ref_wfr_wgr.end(), computed, computed+num_wells);
}
@ -689,19 +711,19 @@ BOOST_AUTO_TEST_CASE(getGFR)
}
{
ADB::V flo = Opm::VFPProperties::getGFR(aqua, liquid, vapour, Opm::VFPProdTable::GFR_GOR);
ADB::V flo = Opm::VFPProdProperties::getGFR(aqua, liquid, vapour, Opm::VFPProdTable::GFR_GOR);
const double* computed = &flo[0];
BOOST_CHECK_EQUAL_COLLECTIONS(ref_gfr_gor.begin(), ref_gfr_gor.end(), computed, computed+num_wells);
}
{
ADB::V flo = Opm::VFPProperties::getGFR(aqua, liquid, vapour, Opm::VFPProdTable::GFR_GLR);
ADB::V flo = Opm::VFPProdProperties::getGFR(aqua, liquid, vapour, Opm::VFPProdTable::GFR_GLR);
const double* computed = &flo[0];
BOOST_CHECK_EQUAL_COLLECTIONS(ref_gfr_glr.begin(), ref_gfr_glr.end(), computed, computed+num_wells);
}
{
ADB::V flo = Opm::VFPProperties::getGFR(aqua, liquid, vapour, Opm::VFPProdTable::GFR_OGR);
ADB::V flo = Opm::VFPProdProperties::getGFR(aqua, liquid, vapour, Opm::VFPProdTable::GFR_OGR);
const double* computed = &flo[0];
BOOST_CHECK_EQUAL_COLLECTIONS(ref_gfr_ogr.begin(), ref_gfr_ogr.end(), computed, computed+num_wells);
}
@ -751,9 +773,9 @@ BOOST_AUTO_TEST_CASE(ParseInterpolateRealisticVFPPROD)
BOOST_CHECK_EQUAL(deck->numKeywords("VFPPROD"), 1);
Opm::VFPProdTable table;
table.init(deck->getKeyword("VFPPROD", 1), units);
table.init(deck->getKeyword("VFPPROD", 0), units);
Opm::VFPProperties properties(NULL, &table);
Opm::VFPProdProperties properties(&table);
//Do some rudimentary testing
//Get the BHP as a function of rate, thp, wfr, gfr, alq
@ -763,10 +785,6 @@ BOOST_AUTO_TEST_CASE(ParseInterpolateRealisticVFPPROD)
double thp[] = {16.010000000000002, 22.438571428571429, 28.867142857142859, 35.295714285714283, 41.724285714285713, 48.152857142857144, 54.581428571428575, 61.009999999999998};
int n = sizeof(liq) / sizeof(liq[0]);
#ifdef verbose
std::cout << std::fixed << std::setprecision(17);
std::cout << "a=[..." << std::endl;
#endif
int i = 0;
double sad = 0.0; //Sum of absolute difference
double max_d = 0.0; //Maximum difference
@ -775,7 +793,6 @@ BOOST_AUTO_TEST_CASE(ParseInterpolateRealisticVFPPROD)
for (int g=0; g<n; ++g) {
//for (unsigned int a=0; a<n; ++a) { //n==1, skip this loop
for (int f=0; f<n; ++f) {
//Liq given as SM3/day => convert to SM3/second
double f_i = liq[f]*1.1574074074074073e-05;
@ -791,27 +808,40 @@ BOOST_AUTO_TEST_CASE(ParseInterpolateRealisticVFPPROD)
//ALQ unit not relevant in this case
double a_i = 0.0;
//Value given as pascal, convert to barsa for comparison with reference
double value_i = properties.prod_bhp(32, f_i, t_i, w_i, g_i, a_i) * 10.0e-6;
//Now reconstruct possible aqua, liquid,
//and vapour phase compositions that satisfy the above
double abs_diff = std::abs(value_i - reference[i]);
sad += abs_diff;
max_d = std::max(max_d, abs_diff);
//liq[f] = aqua + liquid
//wct[w] = aqua / (aqua + liquid)
//gor[g] = vapour / liquid
//
// aqua = wct[w] * liq[f]
// liquid = liq[f] - aqua
// vapour = gor[g] * liquid
double aqua = w_i * f_i;
double liquid = f_i - aqua;
double vapour = g_i * liquid;
if (aqua == 0.0 || liquid == 0.0) {
//FIXME: This skips some corner cases, but will fail in current
//implementation, since getWFR(...), getGFR(...), getFlo(...)
//might perform division by zero
}
else {
//Value given as pascal, convert to barsa for comparison with reference
double value_i = properties.bhp(32, aqua, liquid, vapour, t_i, a_i) * 10.0e-6;
double abs_diff = std::abs(value_i - reference[i]);
sad += abs_diff;
max_d = std::max(max_d, abs_diff);
}
++i;
#ifdef verbose
std::cout << ((f==0)?"":", ") << value << std::flush;
#endif
}
#ifdef verbose
std::cout << " ..." << std::endl;
#endif
//}
}
}
}
#ifdef verbose
std::cout << "];" << std::endl;
#endif
BOOST_CHECK_SMALL(max_d, max_d_tol);
BOOST_CHECK_SMALL(sad, sad_tol);