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

Updated VFPProperties to support a vector VFPProdTable's

Browse Source
This commit is contained in:
André R. Brodtkorb 2015-06-26 15:27:58 +02:00 committed by babrodtk
parent fe7b5f2f6f
commit 066e54bbfc
3 changed files with 114 additions and 43 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

66
opm/autodiff/VFPProperties.cpp
View File

@ -30,26 +30,43 @@
namespace Opm {
VFPProperties::VFPProperties(VFPProdTable* table) : table_(table) {
VFPProperties::VFPProperties() {
}
void VFPProperties::init(const VFPProdTable* table) {
m_tables[table->getTableNum()] = table;
}
void VFPProperties::init(const std::vector<VFPProdTable>& tables) {
//Loop through all tables, and add to our map
for (unsigned int i=0; i<tables.size(); ++i) {
int table_number = tables[i].getTableNum();
if (m_tables.find(table_number) != m_tables.end()) {
OPM_THROW(std::invalid_argument, "Duplicate table numbers found for VFPPROD");
}
else {
m_tables[table_number] = &tables[i];
}
}
}
double VFPProperties::bhp(const double& flo, const double& thp, const double& wfr, const double& gfr, const double& alq) {
double VFPProperties::bhp(int table, const double& flo, const double& thp, const double& wfr, const double& gfr, const double& alq) {
if (m_tables.find(table) == m_tables.end()) {
OPM_THROW(std::invalid_argument, "Nonexistant table " << table << " referenced.");
}
//First, find the values to interpolate between
auto flo_i = find_interp_data(flo, table_->getFloAxis());
auto thp_i = find_interp_data(thp, table_->getTHPAxis());
auto wfr_i = find_interp_data(wfr, table_->getWFRAxis());
auto gfr_i = find_interp_data(gfr, table_->getGFRAxis());
auto alq_i = find_interp_data(alq, table_->getALQAxis());
auto flo_i = find_interp_data(flo, m_tables[table]->getFloAxis());
auto thp_i = find_interp_data(thp, m_tables[table]->getTHPAxis());
auto wfr_i = find_interp_data(wfr, m_tables[table]->getWFRAxis());
auto gfr_i = find_interp_data(gfr, m_tables[table]->getGFRAxis());
auto alq_i = find_interp_data(alq, m_tables[table]->getALQAxis());
//Then perform the interpolation itself
return interpolate(flo_i, thp_i, wfr_i, gfr_i, alq_i);
return interpolate(m_tables[table]->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i);
}
VFPProperties::ADB VFPProperties::bhp(const ADB& flo, const ADB& thp, const ADB& wfr, const ADB& gfr, const ADB& alq) {
VFPProperties::ADB VFPProperties::bhp(int table, const ADB& flo, const ADB& thp, const ADB& wfr, const ADB& gfr, const ADB& alq) {
const ADB::V& f_v = flo.value();
const ADB::V& t_v = thp.value();
const ADB::V& w_v = wfr.value();
@ -62,7 +79,7 @@ VFPProperties::ADB VFPProperties::bhp(const ADB& flo, const ADB& thp, const ADB&
ADB::V bhp_vals;
bhp_vals.resize(nw);
for (int i=0; i<nw; ++i) {
bhp_vals[i] = bhp(f_v[i], t_v[i], w_v[i], g_v[i], a_v[i]);
bhp_vals[i] = bhp(table, f_v[i], t_v[i], w_v[i], g_v[i], a_v[i]);
}
//Create an ADB constant value.
return ADB::constant(bhp_vals);
@ -70,7 +87,11 @@ VFPProperties::ADB VFPProperties::bhp(const ADB& flo, const ADB& thp, const ADB&
VFPProperties::ADB VFPProperties::bhp(const Wells& wells, const ADB& qs, const ADB& thp, const ADB& alq) {
VFPProperties::ADB VFPProperties::bhp(int table, const Wells& wells, const ADB& qs, const ADB& thp, const ADB& alq) {
if (m_tables.find(table) == m_tables.end()) {
OPM_THROW(std::invalid_argument, "Nonexistant table " << table << " referenced.");
}
const int np = wells.number_of_phases;
const int nw = wells.number_of_wells;
@ -81,13 +102,13 @@ VFPProperties::ADB VFPProperties::bhp(const Wells& wells, const ADB& qs, const A
const ADB& o = subset(qs, Span(nw, 1, BlackoilPhases::Liquid*nw));
const ADB& g = subset(qs, Span(nw, 1, BlackoilPhases::Vapour*nw));
ADB flo = getFlo(w, o, g, table_->getFloType());
ADB wfr = getWFR(w, o, g, table_->getWFRType());
ADB gfr = getGFR(w, o, g, table_->getGFRType());
ADB flo = getFlo(w, o, g, m_tables[table]->getFloType());
ADB wfr = getWFR(w, o, g, m_tables[table]->getWFRType());
ADB gfr = getGFR(w, o, g, m_tables[table]->getGFRType());
//TODO: Check ALQ type here?
return bhp(flo, thp, wfr, gfr, alq);
return bhp(table, flo, thp, wfr, gfr, alq);
}
@ -191,8 +212,12 @@ VFPProperties::InterpData VFPProperties::find_interp_data(const double& value, c
#pragma GCC optimize ("unroll-loops")
#endif
double VFPProperties::interpolate(const InterpData& flo_i, const InterpData& thp_i,
const InterpData& wfr_i, const InterpData& gfr_i, const InterpData& alq_i) {
double VFPProperties::interpolate(const VFPProdTable::array_type& array,
const InterpData& flo_i,
const InterpData& thp_i,
const InterpData& wfr_i,
const InterpData& gfr_i,
const InterpData& alq_i) {
double nn[2][2][2][2][2];
//Pick out nearest neighbors (nn) to our evaluation point
@ -200,7 +225,6 @@ double VFPProperties::interpolate(const InterpData& flo_i, const InterpData& thp
//we copy to (nn) will fit better in cache than the full original table for the
//interpolation below.
//The following ladder of for loops will presumably be unrolled by a reasonable compiler.
const VFPProdTable::array_type& data = table_->getTable();
for (int t=0; t<=1; ++t) {
for (int w=0; w<=1; ++w) {
for (int g=0; g<=1; ++g) {
@ -214,7 +238,7 @@ double VFPProperties::interpolate(const InterpData& flo_i, const InterpData& thp
const int fi = flo_i.ind_[f];
//Copy element
nn[t][w][g][a][f] = data[ti][wi][gi][ai][fi];
nn[t][w][g][a][f] = array[ti][wi][gi][ai][fi];
}
}
}

52
opm/autodiff/VFPProperties.hpp
View File

@ -25,6 +25,8 @@
#include <opm/parser/eclipse/EclipseState/Tables/VFPProdTable.hpp>
#include <boost/multi_array.hpp>
#include <map>
namespace Opm {
class VFPProdTable;
@ -38,12 +40,25 @@ public:
typedef AutoDiffBlock<double> ADB;
/**
* Constructor, takes *no* ownership of data.
* Constructor
*/
VFPProperties(VFPProdTable* table);
VFPProperties();
/**
* Initialization routine for a single production table, takes *no* ownership of data.
* @param table A single VFPPROD table
*/
void init(const VFPProdTable* table);
/**
* Initialization routine, takes *no* ownership of data.
* @param tables A vector of different VFPPROD tables. Must have unique table numbers
*/
void init(const std::vector<VFPProdTable>& tables);
/**
* Linear interpolation of bhp as function of the input parameters.
* @param table Table number to use
* @param wells Wells structure with information about wells in qs
* @param qs Flow quantities
* @param thp Tubing head pressure
@ -52,10 +67,15 @@ public:
* @return The bottom hole pressure, interpolated/extrapolated linearly using
* the above parameters from the values in the input table.
*/
ADB bhp(const Wells& wells, const ADB& qs, const ADB& thp, const ADB& alq);
ADB bhp(int table,
const Wells& wells,
const ADB& qs,
const ADB& thp,
const ADB& alq);
/**
* Linear interpolation of bhp as a function of the input parameters
* @param table Table number to use
* @param flo Production rate of oil, gas or liquid
* @param thp Tubing head pressure
* @param wfr Water-oil ratio, water cut, or water-gas ratio
@ -65,10 +85,16 @@ public:
* @return The bottom hole pressure, interpolated/extrapolated linearly using
* the above parameters from the values in the input table.
*/
double bhp(const double& flo, const double& thp, const double& wfr, const double& gfr, const double& alq);
double bhp(int table,
const double& flo,
const double& thp,
const double& wfr,
const double& gfr,
const double& alq);
/**
* Linear interpolation of bhp as a function of the input parameters given as ADBs
* @param table Table number to use
* @param flo Production rate of oil, gas or liquid
* @param thp Tubing head pressure
* @param wfr Water-oil ratio, water cut, or water-gas ratio
@ -79,7 +105,12 @@ public:
* the above parameters from the values in the input table, for each entry in the
* input ADB objects.
*/
ADB bhp(const ADB& flo, const ADB& thp, const ADB& wfr, const ADB& gfr, const ADB& alq);
ADB bhp(int table,
const ADB& flo,
const ADB& thp,
const ADB& wfr,
const ADB& gfr,
const ADB& alq);
/**
* Computes the flo parameter according to the flo_type_
@ -103,7 +134,8 @@ public:
const VFPProdTable::GFR_TYPE& type);
private:
VFPProdTable* table_;
// Map which connects the table number with the table itself
std::map<int, const VFPProdTable*> m_tables;
/**
* Helper struct for linear interpolation
@ -122,8 +154,12 @@ private:
/**
* Helper function which interpolates data using the indices etc. given in the inputs.
*/
double interpolate(const InterpData& flo_i, const InterpData& thp_i,
const InterpData& wfr_i, const InterpData& gfr_i, const InterpData& alq_i);
static double interpolate(const VFPProdTable::array_type& array,
const InterpData& flo_i,
const InterpData& thp_i,
const InterpData& wfr_i,
const InterpData& gfr_i,
const InterpData& alq_i);
};

39
tests/test_vfpproperties.cpp
View File

@ -136,7 +136,8 @@ struct TrivialFixture {
data);
//Initialize properties that use the table
properties.reset(new Opm::VFPProperties(&table));
properties.reset(new Opm::VFPProperties());
properties->init(&table);
}
std::shared_ptr<Opm::VFPProperties> properties;
@ -164,6 +165,18 @@ private:
//Set F to be our test suite fixture for our "trivial" tests
BOOST_FIXTURE_TEST_SUITE( TrivialTests, TrivialFixture )
BOOST_AUTO_TEST_CASE(GetTable)
{
initProperties();
//Table 1 has been initialized
properties->bhp(1, 0.0, 0.0, 0.0, 0.0, 0.0);
//Table 2 does not exist.
BOOST_CHECK_THROW(properties->bhp(2, 0.0, 0.0, 0.0, 0.0, 0.0), std::invalid_argument);
}
/**
* Test that we can generate some dummy zero-data,
* interpolate using doubles as input, and compare against the analytic solution
@ -188,7 +201,7 @@ BOOST_AUTO_TEST_CASE(InterpolateZero)
const double v = m / static_cast<double>(n-1);
//Note order of arguments!
sum += properties->bhp(v, x, y, z, u);
sum += properties->bhp(1, v, x, y, z, u);
}
}
}
@ -223,7 +236,7 @@ BOOST_AUTO_TEST_CASE(InterpolateOne)
const double v = m / static_cast<double>(n-1);
//Note order of arguments!
sum += properties->bhp(v, x, y, z, u);
sum += properties->bhp(1, v, x, y, z, u);
}
}
}
@ -264,7 +277,7 @@ BOOST_AUTO_TEST_CASE(InterpolatePlane)
reference_sum += reference;
//Note order of arguments!
double value = properties->bhp(v, x, y, z, u);
double value = properties->bhp(1, v, x, y, z, u);
sum += value;
double abs_diff = std::abs(value - reference);
@ -313,7 +326,7 @@ BOOST_AUTO_TEST_CASE(ExtrapolatePlane)
reference_sum += reference;
//Note order of arguments!
double value = properties->bhp(v, x, y, z, u);
double value = properties->bhp(1, v, x, y, z, u);
sum += value;
double abs_diff = std::abs(value - reference);
@ -381,7 +394,7 @@ BOOST_AUTO_TEST_CASE(ExtrapolatePlaneADB)
ADB gfr = ADB::constant(zz);
ADB alq = ADB::constant(uu);
ADB bhp_val = properties->bhp(flo, thp, wfr, gfr, alq);
ADB bhp_val = properties->bhp(1, flo, thp, wfr, gfr, alq);
double value = 0.0;
double reference = 0.0;
@ -461,7 +474,7 @@ BOOST_AUTO_TEST_CASE(InterpolateADBAndQs)
ADB alq = ADB::constant(alq_vals);
//Call the bhp function
ADB bhp = properties->bhp(*wells, qs, thp, alq);
ADB bhp = properties->bhp(1, *wells, qs, thp, alq);
//Calculate reference
//First, find the three phases
@ -701,16 +714,15 @@ BOOST_AUTO_TEST_CASE(ParseInterpolateRealisticVFPPROD)
std::shared_ptr<Opm::UnitSystem> units(Opm::UnitSystem::newMETRIC());
Opm::ParserPtr parser(new Opm::Parser());
boost::filesystem::path file("tests/VFPPROD1");
boost::filesystem::path file("tests/VFPPROD2");
deck = parser->parseFile(file.string());
Opm::checkDeck(deck);
BOOST_REQUIRE(deck->hasKeyword("VFPPROD"));
BOOST_CHECK_EQUAL(deck->numKeywords("VFPPROD"), 2);
BOOST_CHECK_EQUAL(deck->numKeywords("VFPPROD"), 1);
std::vector<Opm::VFPProdTable> tables;
std::vector<Opm::VFPProperties> properties;
int num_tables = deck->numKeywords("VFPPROD");
for (int i=0; i<num_tables; ++i) {
@ -719,9 +731,8 @@ BOOST_AUTO_TEST_CASE(ParseInterpolateRealisticVFPPROD)
tables[i].init(keyword, units);
}
for (int i=0; i<num_tables; ++i) {
properties.push_back(Opm::VFPProperties(&(tables[i])));
}
Opm::VFPProperties properties;
properties.init(tables);
//Do some rudimentary testing
//Get the BHP as a function of rate, thp, wfr, gfr, alq
@ -760,7 +771,7 @@ BOOST_AUTO_TEST_CASE(ParseInterpolateRealisticVFPPROD)
double a_i = 0.0;
//Value given as pascal, convert to barsa for comparison with reference
double value_i = properties[0].bhp(f_i, t_i, w_i, g_i, a_i) * 10.0e-6;
double value_i = properties.bhp(32, f_i, t_i, w_i, g_i, a_i) * 10.0e-6;
double abs_diff = std::abs(value_i - reference[i]);
sad += abs_diff;