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

Initial integration of VFPProdTables

Browse Source
This commit is contained in:
babrodtk 2015-07-07 15:49:10 +02:00
parent 3260e978da
commit d27403b427
5 changed files with 317 additions and 246 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

44
opm/autodiff/BlackoilModelBase_impl.hpp
View File

@ -32,6 +32,7 @@
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
#include <opm/autodiff/GeoProps.hpp>
#include <opm/autodiff/WellDensitySegmented.hpp>
#include <opm/autodiff/VFPProperties.hpp>
#include <opm/core/grid.h>
#include <opm/core/linalg/LinearSolverInterface.hpp>
@ -1096,6 +1097,7 @@ namespace detail {
}
bool constraintBroken(const std::vector<double>& bhp,
const std::vector<double>& thp,
const std::vector<double>& well_phase_flow_rate,
const int well,
const int num_phases,
@ -1117,6 +1119,11 @@ namespace detail {
broken = bhp[well] > target;
break;
case THP:
//FIXME: Not correct
broken = thp[well] > target;
break;
case RESERVOIR_RATE: // Intentional fall-through
case SURFACE_RATE:
broken = rateToCompare(well_phase_flow_rate,
@ -1133,6 +1140,11 @@ namespace detail {
broken = bhp[well] < target;
break;
case THP:
//FIXME: Not correct
broken = thp[well] < target;
break;
case RESERVOIR_RATE: // Intentional fall-through
case SURFACE_RATE:
// Note that the rates compared below are negative,
@ -1162,7 +1174,7 @@ namespace detail {
{
if( ! wellsActive() ) return ;
std::string modestring[3] = { "BHP", "RESERVOIR_RATE", "SURFACE_RATE" };
std::string modestring[4] = { "BHP", "THP", "RESERVOIR_RATE", "SURFACE_RATE" };
// Find, for each well, if any constraints are broken. If so,
// switch control to first broken constraint.
const int np = wells().number_of_phases;
@ -1185,7 +1197,7 @@ namespace detail {
// inequality constraint, and therefore skipped.
continue;
}
if (detail::constraintBroken(xw.bhp(), xw.wellRates(), w, np, wells().type[w], wc, ctrl_index)) {
if (detail::constraintBroken(xw.bhp(), xw.thp(), xw.wellRates(), w, np, wells().type[w], wc, ctrl_index)) {
// ctrl_index will be the index of the broken constraint after the loop.
break;
}
@ -1203,7 +1215,7 @@ namespace detail {
}
// Updating well state and primary variables.
// Target values are used as initial conditions for BHP and SURFACE_RATE
// Target values are used as initial conditions for BHP, THP, and SURFACE_RATE
const double target = well_controls_iget_target(wc, current);
const double* distr = well_controls_iget_distr(wc, current);
switch (well_controls_iget_type(wc, current)) {
@ -1211,6 +1223,10 @@ namespace detail {
xw.bhp()[w] = target;
break;
case THP:
xw.thp()[w] = target;
break;
case RESERVOIR_RATE:
// No direct change to any observable quantity at
// surface condition. In this case, use existing
@ -1330,6 +1346,14 @@ namespace detail {
const int np = wells().number_of_phases;
const int nw = wells().number_of_wells;
const ADB& aqua_adb = subset(state.qs, Span(nw, 1, BlackoilPhases::Aqua*nw));
const ADB& liquid_adb = subset(state.qs, Span(nw, 1, BlackoilPhases::Liquid*nw));
const ADB& vapour_adb = subset(state.qs, Span(nw, 1, BlackoilPhases::Vapour*nw));
const V& aqua = aqua_adb.value();
const V& liquid = liquid_adb.value();
const V& vapour = vapour_adb.value();
V bhp_targets = V::Zero(nw);
V rate_targets = V::Zero(nw);
M rate_distr(nw, np*nw);
@ -1347,7 +1371,17 @@ namespace detail {
rate_targets(w) = -1e100;
}
break;
//ARB: case THP:
case THP:
{
const int vfp = well_controls_iget_vfp(wc, current);
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);
rate_targets(w) = -1e100;
}
break;
case RESERVOIR_RATE: // Intentional fall-through
case SURFACE_RATE:
@ -1371,7 +1405,7 @@ namespace detail {
}
const ADB bhp_residual = state.bhp - bhp_targets;
const ADB rate_residual = rate_distr * state.qs - rate_targets;
//const ADB thp_residual = state.bhp - vfpprop.bhp(thp_ctrl, state.qs, alq);
//wells
// Choose bhp residual for positive bhp targets.
Selector<double> bhp_selector(bhp_targets);

6
opm/autodiff/SimulatorBase_impl.hpp
View File

@ -54,7 +54,7 @@ namespace Opm
rateConverter_(props_, std::vector<int>(AutoDiffGrid::numCells(grid_), 0)),
threshold_pressures_by_face_(threshold_pressures_by_face),
is_parallel_run_( false ),
vfpProperties_(eclipse_state->getVFPProdTables())
vfpProperties_(eclipse_state->getVFPInjTables(), eclipse_state->getVFPProdTables())
{
// Misc init.
const int num_cells = AutoDiffGrid::numCells(grid);
@ -440,14 +440,14 @@ namespace Opm
well_controls_assert_number_of_phases(ctrl, int(np));
const int ok_resv =
well_controls_add_new(RESERVOIR_RATE, target,
well_controls_add_new(RESERVOIR_RATE, target, -1e100, -1e100,
& distr[0], ctrl);
// For WCONHIST/RESV the BHP limit is set to 1 atm.
// TODO: Make it possible to modify the BHP limit using
// the WELTARG keyword
const int ok_bhp =
well_controls_add_new(BHP, unit::convert::from(1.0, unit::atm),
well_controls_add_new(BHP, unit::convert::from(1.0, unit::atm), -1e100, -1e100,
NULL, ctrl);
if (ok_resv != 0 && ok_bhp != 0) {

171
opm/autodiff/VFPProperties.cpp
View File

@ -23,7 +23,6 @@
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <algorithm>
@ -72,129 +71,83 @@ void VFPProperties::init(const std::map<int, VFPProdTable>& prod_tables) {
}
}
double VFPProperties::prod_bhp(int table, const double& flo, const double& thp, const double& wfr, const double& gfr, const double& alq) {
if (m_prod_tables.find(table) == m_prod_tables.end()) {
OPM_THROW(std::invalid_argument, "Nonexistent table " << table << " referenced.");
}
const auto* tab = m_prod_tables[table];
//First, find the values to interpolate between
auto flo_i = find_interp_data(flo, tab->getFloAxis());
auto thp_i = find_interp_data(thp, tab->getTHPAxis());
auto wfr_i = find_interp_data(wfr, tab->getWFRAxis());
auto gfr_i = find_interp_data(gfr, tab->getGFRAxis());
auto alq_i = find_interp_data(alq, tab->getALQAxis());
//Then perform the interpolation itself
return interpolate(tab->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i);
}
VFPProperties::ADB VFPProperties::prod_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();
const ADB::V& g_v = gfr.value();
const ADB::V& a_v = alq.value();
const int nw = f_v.size();
//Compute the BHP for each well independently
ADB::V bhp_vals;
bhp_vals.resize(nw);
for (int i=0; i<nw; ++i) {
bhp_vals[i] = prod_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);
}
VFPProperties::ADB VFPProperties::prod_bhp(int table, const Wells& wells, const ADB& qs, const ADB& thp, const ADB& alq) {
if (m_prod_tables.find(table) == m_prod_tables.end()) {
OPM_THROW(std::invalid_argument, "Nonexistant table " << table << " referenced.");
}
VFPProperties::ADB::V VFPProperties::prod_bhp(int table_id,
const Wells& wells,
const ADB::V& qs,
const ADB::V& thp,
const ADB::V& alq) const {
const int np = wells.number_of_phases;
const int nw = wells.number_of_wells;
//Short-hands for water / oil / gas phases
//TODO enable support for two-phase.
assert(np == 3);
const ADB& w = subset(qs, Span(nw, 1, BlackoilPhases::Aqua*nw));
const ADB& o = subset(qs, Span(nw, 1, BlackoilPhases::Liquid*nw));
const ADB& g = subset(qs, Span(nw, 1, BlackoilPhases::Vapour*nw));
const ADB::V& w = subset(qs, Span(nw, 1, BlackoilPhases::Aqua*nw));
const ADB::V& o = subset(qs, Span(nw, 1, BlackoilPhases::Liquid*nw));
const ADB::V& g = subset(qs, Span(nw, 1, BlackoilPhases::Vapour*nw));
const auto* tab = m_prod_tables[table];
ADB flo = getFlo(w, o, g, tab->getFloType());
ADB wfr = getWFR(w, o, g, tab->getWFRType());
ADB gfr = getGFR(w, o, g, tab->getGFRType());
return prod_bhp(table_id, w, o, g, thp, alq);
}
//TODO: Check ALQ type here?
VFPProperties::ADB::V VFPProperties::prod_bhp(int table_id,
const ADB::V& aqua,
const ADB::V& liquid,
const ADB::V& vapour,
const ADB::V& thp,
const ADB::V& alq) const {
const int nw = thp.size();
return prod_bhp(table, flo, thp, wfr, gfr, alq);
assert(aqua.size() == nw);
assert(liquid.size() == nw);
assert(vapour.size() == nw);
assert(thp.size() == nw);
assert(alq.size() == nw);
//Compute the BHP for each well independently
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]);
}
return bhp_vals;
}
double VFPProperties::prod_bhp(int table_id,
const double& aqua,
const double& liquid,
const double& vapour,
const double& thp,
const double& alq) const {
const VFPProdTable* table = getProdTable(table_id);
//Find interpolation variables
double flo = getFlo(aqua, liquid, vapour, table->getFloType());
double wfr = getWFR(aqua, liquid, vapour, table->getWFRType());
double gfr = getGFR(aqua, liquid, vapour, table->getGFRType());
//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());
//Then perform the interpolation itself
return interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i);
}
VFPProperties::ADB VFPProperties::getFlo(const ADB& aqua, const ADB& liquid, const ADB& vapour,
const VFPProdTable::FLO_TYPE& type) {
switch (type) {
case VFPProdTable::FLO_OIL:
//Oil = liquid phase
return liquid;
case VFPProdTable::FLO_LIQ:
//Liquid = aqua + liquid phases
return aqua + liquid;
case VFPProdTable::FLO_GAS:
//Gas = vapor phase
return vapour;
case VFPProdTable::FLO_INVALID: //Intentional fall-through
default:
OPM_THROW(std::logic_error, "Invalid FLO_TYPE: '" << type << "'");
return ADB::null();
const VFPProdTable* VFPProperties::getProdTable(int table_id) const {
auto entry = m_prod_tables.find(table_id);
if (entry == m_prod_tables.end()) {
OPM_THROW(std::invalid_argument, "Nonexistent table " << table_id << " referenced.");
}
else {
return entry->second;
}
}
VFPProperties::ADB VFPProperties::getWFR(const ADB& aqua, const ADB& liquid, const ADB& vapour,
const VFPProdTable::WFR_TYPE& type) {
switch(type) {
case VFPProdTable::WFR_WOR:
//Water-oil ratio = water / oil
return aqua / liquid;
case VFPProdTable::WFR_WCT:
//Water cut = water / (water + oil + gas)
return aqua / (aqua + liquid + vapour);
case VFPProdTable::WFR_WGR:
//Water-gas ratio = water / gas
return aqua / vapour;
case VFPProdTable::WFR_INVALID: //Intentional fall-through
default:
OPM_THROW(std::logic_error, "Invalid WFR_TYPE: '" << type << "'");
return ADB::null();
}
}
VFPProperties::ADB VFPProperties::getGFR(const ADB& aqua, const ADB& liquid, const ADB& vapour,
const VFPProdTable::GFR_TYPE& type) {
switch(type) {
case VFPProdTable::GFR_GOR:
// Gas-oil ratio = gas / oil
return vapour / liquid;
case VFPProdTable::GFR_GLR:
// Gas-liquid ratio = gas / (oil + water)
return vapour / (liquid + aqua);
case VFPProdTable::GFR_OGR:
// Oil-gas ratio = oil / gas
return liquid / vapour;
case VFPProdTable::GFR_INVALID: //Intentional fall-through
default:
OPM_THROW(std::logic_error, "Invalid GFR_TYPE: '" << type << "'");
return ADB::null();
}
}
VFPProperties::InterpData VFPProperties::find_interp_data(const double& value, const std::vector<double>& values) {
InterpData retval;

133
opm/autodiff/VFPProperties.hpp
View File

@ -21,6 +21,7 @@
#define OPM_AUTODIFF_VFPPROPERTIES_HPP_
#include <opm/core/wells.h>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/VFPProdTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/VFPInjTable.hpp>
@ -77,7 +78,7 @@ public:
/**
* Linear interpolation of bhp as function of the input parameters.
* @param table Table number to use
* @param table_id Table number to use
* @param wells Wells structure with information about wells in qs
* @param qs Flow quantities
* @param thp Tubing head pressure
@ -86,50 +87,50 @@ public:
* @return The bottom hole pressure, interpolated/extrapolated linearly using
* the above parameters from the values in the input table.
*/
ADB prod_bhp(int table,
ADB::V prod_bhp(int table_id,
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
* @param gfr Gas-oil ratio, gas-liquid ratio, or oil-gas ratio
* @param alq Artificial lift or other parameter
*
* @return The bottom hole pressure, interpolated/extrapolated linearly using
* the above parameters from the values in the input table.
*/
double prod_bhp(int table,
const double& flo,
const double& thp,
const double& wfr,
const double& gfr,
const double& alq);
const ADB::V& qs,
const ADB::V& thp,
const ADB::V& alq) const;
/**
* 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 aqua Water phase
* @param liquid Oil phase
* @param vapour Gas phase
* @param thp Tubing head pressure
* @param wfr Water-oil ratio, water cut, or water-gas ratio
* @param gfr Gas-oil ratio, gas-liquid ratio, or oil-gas ratio
* @param alq Artificial lift or other parameter
*
* @return The bottom hole pressure, interpolated/extrapolated linearly using
* the above parameters from the values in the input table, for each entry in the
* input ADB objects.
*/
ADB prod_bhp(int table,
const ADB& flo,
const ADB& thp,
const ADB& wfr,
const ADB& gfr,
const ADB& alq);
ADB::V prod_bhp(int table_id,
const ADB::V& aqua,
const ADB::V& liquid,
const ADB::V& vapour,
const ADB::V& thp,
const ADB::V& alq) const;
/**
* Linear interpolation of bhp as a function of the input parameters
* @param table_id Table number to use
* @param aqua Water phase
* @param liquid Oil phase
* @param vapour Gas phase
* @param thp Tubing head pressure
* @param alq Artificial lift or other parameter
*
* @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,
const double& aqua,
const double& liquid,
const double& vapour,
const double& thp,
const double& alq) const;
//FIXME: ARB: Implement inj_bhp to match the prod_bhp's, but for injection wells.
@ -137,22 +138,71 @@ public:
* Computes the flo parameter according to the flo_type_
* @return Production rate of oil, gas or liquid.
*/
static ADB getFlo(const ADB& aqua, const ADB& liquid, const ADB& vapour,
const VFPProdTable::FLO_TYPE& type);
template <typename T>
static T getFlo(const T& aqua, const T& liquid, const T& vapour,
const VFPProdTable::FLO_TYPE& type) {
switch (type) {
case VFPProdTable::FLO_OIL:
//Oil = liquid phase
return liquid;
case VFPProdTable::FLO_LIQ:
//Liquid = aqua + liquid phases
return aqua + liquid;
case VFPProdTable::FLO_GAS:
//Gas = vapor phase
return vapour;
case VFPProdTable::FLO_INVALID: //Intentional fall-through
default:
OPM_THROW(std::logic_error, "Invalid FLO_TYPE: '" << type << "'");
}
}
/**
* Computes the wfr parameter according to the wfr_type_
* @return Production rate of oil, gas or liquid.
*/
static ADB getWFR(const ADB& aqua, const ADB& liquid, const ADB& vapour,
const VFPProdTable::WFR_TYPE& type);
template <typename T>
static T getWFR(const T& aqua, const T& liquid, const T& vapour,
const VFPProdTable::WFR_TYPE& type) {
switch(type) {
case VFPProdTable::WFR_WOR:
//Water-oil ratio = water / oil
return aqua / liquid;
case VFPProdTable::WFR_WCT:
//Water cut = water / (water + oil + gas)
return aqua / (aqua + liquid + vapour);
case VFPProdTable::WFR_WGR:
//Water-gas ratio = water / gas
return aqua / vapour;
case VFPProdTable::WFR_INVALID: //Intentional fall-through
default:
OPM_THROW(std::logic_error, "Invalid WFR_TYPE: '" << type << "'");
}
}
/**
* Computes the gfr parameter according to the gfr_type_
* @return Production rate of oil, gas or liquid.
*/
static ADB getGFR(const ADB& aqua, const ADB& liquid, const ADB& vapour,
const VFPProdTable::GFR_TYPE& type);
template <typename T>
static T getGFR(const T& aqua, const T& liquid, const T& vapour,
const VFPProdTable::GFR_TYPE& type) {
switch(type) {
case VFPProdTable::GFR_GOR:
// Gas-oil ratio = gas / oil
return vapour / liquid;
case VFPProdTable::GFR_GLR:
// Gas-liquid ratio = gas / (oil + water)
return vapour / (liquid + aqua);
case VFPProdTable::GFR_OGR:
// Oil-gas ratio = oil / gas
return liquid / vapour;
case VFPProdTable::GFR_INVALID: //Intentional fall-through
default:
OPM_THROW(std::logic_error, "Invalid GFR_TYPE: '" << type << "'");
}
}
private:
// Map which connects the table number with the table itself
@ -188,6 +238,11 @@ private:
*/
void init(const std::map<int, VFPInjTable>& inj_tables);
void init(const std::map<int, VFPProdTable>& prod_tables);
/**
* Misc helper functions
*/
const VFPProdTable* getProdTable(int table_id) const;
};
}

209
tests/test_vfpproperties.cpp
View File

@ -117,6 +117,28 @@ struct TrivialFixture {
}
}
/**
* Fills our interpolation data with "random" values
*/
void fillDataRandom() {
unsigned long randx = 42;
for (int i=0; i<nx; ++i) {
for (int j=0; j<ny; ++j) {
for (int k=0; k<nz; ++k) {
for (int l=0; l<nu; ++l) {
for (int m=0; m<nv; ++m) {
data[i][j][k][l][m] = randx;
randx = randx*1103515245 + 12345;
}
}
}
}
}
}
void initTable() {
}
@ -140,6 +162,7 @@ struct TrivialFixture {
}
std::shared_ptr<Opm::VFPProperties> properties;
Opm::VFPProdTable table;
private:
const std::vector<double> thp_axis;
@ -154,7 +177,6 @@ private:
int nv;
Opm::VFPProdTable::extents size;
Opm::VFPProdTable::array_type data;
Opm::VFPProdTable table;
};
@ -167,13 +189,45 @@ BOOST_FIXTURE_TEST_SUITE( TrivialTests, TrivialFixture )
BOOST_AUTO_TEST_CASE(GetTable)
{
fillDataRandom();
initProperties();
//Table 1 has been initialized
properties->prod_bhp(1, 0.0, 0.0, 0.0, 0.0, 0.0);
//Create wells
const int nphases = 3;
const int nwells = 5;
const int nperfs = 1;
std::shared_ptr<Wells> wells(create_wells(nphases, nwells, nperfs),
destroy_wells);
//Create interpolation points
double aqua_d = 1.5;
double liquid_d = 2.5;
double vapour_d = 3.5;
double thp_d = 4.5;
double alq_d = 5.5;
ADB::V aqua_adb = ADB::V::Constant(1, aqua_d);
ADB::V liquid_adb = ADB::V::Constant(1, liquid_d);
ADB::V vapour_adb = ADB::V::Constant(1, vapour_d);
ADB::V thp_adb = ADB::V::Constant(1, thp_d);
ADB::V alq_adb = ADB::V::Constant(1, alq_d);
ADB::V qs_adb;
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);
//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, 0.0, 0.0, 0.0, 0.0, 0.0), std::invalid_argument);
BOOST_CHECK_THROW(properties->prod_bhp(2, *wells, qs_adb, thp_adb, alq_adb), std::invalid_argument);
}
/**
@ -357,11 +411,11 @@ BOOST_AUTO_TEST_CASE(ExtrapolatePlaneADB)
//Temporary variables used to represent independent wells
const int num_wells = 5;
ADB::V xx(num_wells);
ADB::V yy(num_wells);
ADB::V zz(num_wells);
ADB::V uu(num_wells);
ADB::V vv(num_wells);
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);
//Check linear extrapolation (i.e., using values of x, y, etc. outside our interpolant domain)
double sum = 0.0;
@ -380,26 +434,20 @@ BOOST_AUTO_TEST_CASE(ExtrapolatePlaneADB)
for (int m=-5; m<=n+5; ++m) {
const double v = m / static_cast<double>(n);
for (unsigned int w=0; w<num_wells; ++w) {
xx[w] = x*w;
yy[w] = y*w;
zz[w] = z*w;
uu[w] = u*w;
vv[w] = v*w;
thp[w] = x*w;
wfr[w] = y*w;
gfr[w] = z*w;
alq[w] = u*w;
flo[w] = v*w;
}
ADB flo = ADB::constant(vv);
ADB thp = ADB::constant(xx);
ADB wfr = ADB::constant(yy);
ADB gfr = ADB::constant(zz);
ADB alq = ADB::constant(uu);
ADB bhp_val = properties->prod_bhp(1, flo, thp, wfr, gfr, alq);
ADB::V bhp_val = properties->prod_bhp(1, flo, thp, wfr, gfr, alq);
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;
value = bhp_val.value()[w];
value = bhp_val[w];
sum += value;
reference_sum += reference;
@ -450,30 +498,27 @@ BOOST_AUTO_TEST_CASE(InterpolateADBAndQs)
}
//Create some artificial flow values for our wells between 0 and 1
ADB::V qs_vals(nphases*nwells);
ADB::V qs(nphases*nwells);
for (int j=0; j<nphases; ++j) {
for (int i=0; i<nwells; ++i) {
qs_vals[j*nwells+i] = (j*nwells+i) / static_cast<double>(nwells*nphases-1.0);
qs[j*nwells+i] = (j*nwells+i) / static_cast<double>(nwells*nphases-1.0);
}
}
ADB qs = ADB::constant(qs_vals);
//Create the THP for each well
ADB::V thp_vals(nwells);
ADB::V thp(nwells);
for (int i=0; i<nwells; ++i) {
thp_vals[i] = (i) / static_cast<double>(nwells-1.0);
thp[i] = (i) / static_cast<double>(nwells-1.0);
}
ADB thp = ADB::constant(thp_vals);
//Create the ALQ for each well
ADB::V alq_vals(nwells);
ADB::V alq(nwells);
for (int i=0; i<nwells; ++i) {
alq_vals[i] = 0.0;
alq[i] = 0.0;
}
ADB alq = ADB::constant(alq_vals);
//Call the bhp function
ADB bhp = properties->prod_bhp(1, *wells, qs, thp, alq);
ADB::V bhp = properties->prod_bhp(1, *wells, qs, thp, alq);
//Calculate reference
//First, find the three phases
@ -481,9 +526,9 @@ BOOST_AUTO_TEST_CASE(InterpolateADBAndQs)
std::vector<double> oil(nwells);
std::vector<double> gas(nwells);
for (int i=0; i<nwells; ++i) {
water[i] = qs_vals[i];
oil[i] = qs_vals[nwells+i];
gas[i] = qs_vals[2*nwells+i];
water[i] = qs[i];
oil[i] = qs[nwells+i];
gas[i] = qs[2*nwells+i];
}
//Compute reference value
@ -492,16 +537,15 @@ BOOST_AUTO_TEST_CASE(InterpolateADBAndQs)
double flo = oil[i];
double wor = water[i]/oil[i];
double gor = gas[i]/oil[i];
reference[i] = thp_vals[i] + 2*wor + 3*gor + 4*alq_vals[i] + 5*flo;
reference[i] = thp[i] + 2*wor + 3*gor + 4*alq[i] + 5*flo;
}
//Check that interpolation matches
const ADB::V& values = bhp.value();
BOOST_REQUIRE_EQUAL(values.size(), nwells);
BOOST_REQUIRE_EQUAL(bhp.size(), nwells);
double sad = 0.0;
double max_d = 0.0;
for (int i=0; i<nwells; ++i) {
double value = values[i];
double value = bhp[i];
double ref = reference[i];
double abs_diff = std::abs(value-ref);
sad += abs_diff;
@ -538,27 +582,16 @@ struct ConversionFixture {
ConversionFixture() :
num_wells(5),
d_aqua(num_wells),
d_liquid(num_wells),
d_vapour(num_wells),
aqua(ADB::null()),
liquid(ADB::null()),
vapour(ADB::null())
aqua(num_wells),
liquid(num_wells),
vapour(num_wells)
{
for (int i=0; i<num_wells; ++i) {
d_aqua[i] = 300+num_wells*15;
d_liquid[i] = 500+num_wells*15;
d_vapour[i] = 700+num_wells*15;
aqua[i] = 300+num_wells*15;
liquid[i] = 500+num_wells*15;
vapour[i] = 700+num_wells*15;
}
//Copy the double vectors above to ADBs
ADB::V v_aqua = Eigen::Map<ADB::V>(&d_aqua[0], num_wells, 1);
ADB::V v_liquid = Eigen::Map<ADB::V>(&d_liquid[0], num_wells, 1);
ADB::V v_vapour = Eigen::Map<ADB::V>(&d_vapour[0], num_wells, 1);
aqua = ADB::constant(v_aqua);
liquid = ADB::constant(v_liquid);
vapour = ADB::constant(v_vapour);
}
~ConversionFixture() {
@ -567,13 +600,9 @@ struct ConversionFixture {
int num_wells;
std::vector<double> d_aqua;
std::vector<double> d_liquid;
std::vector<double> d_vapour;
ADB aqua;
ADB liquid;
ADB vapour;
ADB::V aqua;
ADB::V liquid;
ADB::V vapour;
};
@ -590,26 +619,26 @@ BOOST_AUTO_TEST_CASE(getFlo)
std::vector<double> ref_flo_liq(num_wells);
std::vector<double> ref_flo_gas(num_wells);
for (int i=0; i<num_wells; ++i) {
ref_flo_oil[i] = d_liquid[i];
ref_flo_liq[i] = d_aqua[i] + d_liquid[i];
ref_flo_gas[i] = d_vapour[i];
ref_flo_oil[i] = liquid[i];
ref_flo_liq[i] = aqua[i] + liquid[i];
ref_flo_gas[i] = vapour[i];
}
{
ADB flo = Opm::VFPProperties::getFlo(aqua, liquid, vapour, Opm::VFPProdTable::FLO_OIL);
const double* computed = &flo.value()[0];
ADB::V flo = Opm::VFPProperties::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 flo = Opm::VFPProperties::getFlo(aqua, liquid, vapour, Opm::VFPProdTable::FLO_LIQ);
const double* computed = &flo.value()[0];
ADB::V flo = Opm::VFPProperties::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 flo = Opm::VFPProperties::getFlo(aqua, liquid, vapour, Opm::VFPProdTable::FLO_GAS);
const double* computed = &flo.value()[0];
ADB::V flo = Opm::VFPProperties::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);
}
}
@ -622,26 +651,26 @@ BOOST_AUTO_TEST_CASE(getWFR)
std::vector<double> ref_wfr_wct(num_wells);
std::vector<double> ref_wfr_wgr(num_wells);
for (int i=0; i<num_wells; ++i) {
ref_wfr_wor[i] = d_aqua[i] / d_liquid[i];
ref_wfr_wct[i] = d_aqua[i] / (d_aqua[i] + d_liquid[i] + d_vapour[i]);
ref_wfr_wgr[i] = d_aqua[i] / d_vapour[i];
ref_wfr_wor[i] = aqua[i] / liquid[i];
ref_wfr_wct[i] = aqua[i] / (aqua[i] + liquid[i] + vapour[i]);
ref_wfr_wgr[i] = aqua[i] / vapour[i];
}
{
ADB flo = Opm::VFPProperties::getWFR(aqua, liquid, vapour, Opm::VFPProdTable::WFR_WOR);
const double* computed = &flo.value()[0];
ADB::V flo = Opm::VFPProperties::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 flo = Opm::VFPProperties::getWFR(aqua, liquid, vapour, Opm::VFPProdTable::WFR_WCT);
const double* computed = &flo.value()[0];
ADB::V flo = Opm::VFPProperties::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 flo = Opm::VFPProperties::getWFR(aqua, liquid, vapour, Opm::VFPProdTable::WFR_WGR);
const double* computed = &flo.value()[0];
ADB::V flo = Opm::VFPProperties::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);
}
}
@ -654,26 +683,26 @@ BOOST_AUTO_TEST_CASE(getGFR)
std::vector<double> ref_gfr_glr(num_wells);
std::vector<double> ref_gfr_ogr(num_wells);
for (int i=0; i<num_wells; ++i) {
ref_gfr_gor[i] = d_vapour[i] / d_liquid[i];
ref_gfr_glr[i] = d_vapour[i] / (d_liquid[i] + d_aqua[i]);
ref_gfr_ogr[i] = d_liquid[i] / d_vapour[i];
ref_gfr_gor[i] = vapour[i] / liquid[i];
ref_gfr_glr[i] = vapour[i] / (liquid[i] + aqua[i]);
ref_gfr_ogr[i] = liquid[i] / vapour[i];
}
{
ADB flo = Opm::VFPProperties::getGFR(aqua, liquid, vapour, Opm::VFPProdTable::GFR_GOR);
const double* computed = &flo.value()[0];
ADB::V flo = Opm::VFPProperties::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 flo = Opm::VFPProperties::getGFR(aqua, liquid, vapour, Opm::VFPProdTable::GFR_GLR);
const double* computed = &flo.value()[0];
ADB::V flo = Opm::VFPProperties::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 flo = Opm::VFPProperties::getGFR(aqua, liquid, vapour, Opm::VFPProdTable::GFR_OGR);
const double* computed = &flo.value()[0];
ADB::V flo = Opm::VFPProperties::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);
}
}