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 use the newly implemented VFPProdTable class in opm-parser (with units)

Browse Source
This commit is contained in:
André R. Brodtkorb 2015-06-25 16:30:46 +02:00 committed by babrodtk
parent 1d7f601fab
commit 179a210ad5
3 changed files with 165 additions and 638 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

356
opm/autodiff/VFPProperties.cpp
View File

@ -19,317 +19,31 @@
#include "config.h"
#include <opm/autodiff/VFPProperties.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <algorithm>
namespace Opm {
/**
* Helper function that checks if an item exists in a record, and has a
* non-zero size
*/
bool itemValid(DeckRecordConstPtr& record, const char* name) {
if (record->size() == 0) {
return false;
}
else {
DeckItemPtr item;
//TODO: Should we instead here allow the exception to propagate?
try {
item = record->getItem(name);
}
catch (...) {
return false;
}
if (item->size() > 0) {
return true;
}
else {
return false;
}
}
VFPProperties::VFPProperties(VFPProdTable* table) : table_(table) {
}
VFPProperties::VFPProperties(DeckKeywordConstPtr table) {
auto iter = table->begin();
auto header = (*iter++);
assert(itemValid(header, "TABLE"));
table_num_ = header->getItem("TABLE")->getInt(0);
assert(itemValid(header, "DATUM_DEPTH"));
datum_depth_ = header->getItem("DATUM_DEPTH")->getRawDouble(0);
//Rate type
assert(itemValid(header, "RATE_TYPE"));
std::string flo_string = header->getItem("RATE_TYPE")->getString(0);
if (flo_string == "OIL") {
flo_type_ = FLO_OIL;
}
else if (flo_string == "LIQ") {
flo_type_ = FLO_LIQ;
}
else if (flo_string == "GAS") {
flo_type_ = FLO_GAS;
}
else {
flo_type_ = FLO_INVALID;
OPM_THROW(std::runtime_error, "Invalid RATE_TYPE string: '" << flo_string << "'");
}
//Water fraction
assert(itemValid(header, "WFR"));
std::string wfr_string = header->getItem("WFR")->getString(0);
if (wfr_string == "WOR") {
wfr_type_ = WFR_WOR;
}
else if (wfr_string == "WCT") {
wfr_type_ = WFR_WCT;
}
else if (wfr_string == "WGR") {
wfr_type_ = WFR_WGR;
}
else {
wfr_type_ = WFR_INVALID;
OPM_THROW(std::runtime_error, "Invalid WFR string: '" << wfr_string << "'");
}
//Gas fraction
assert(itemValid(header, "GFR"));
std::string gfr_string = header->getItem("GFR")->getString(0);
if (gfr_string == "GOR") {
gfr_type_ = GFR_GOR;
}
else if (gfr_string == "GLR") {
gfr_type_ = GFR_GLR;
}
else if (gfr_string == "OGR") {
gfr_type_ = GFR_OGR;
}
else {
gfr_type_ = GFR_INVALID;
OPM_THROW(std::runtime_error, "Invalid GFR string: '" << gfr_string << "'");
}
//Definition of THP values, must be THP
if (itemValid(header, "PRESSURE_DEF")) {
std::string quantity_string = header->getItem("PRESSURE_DEF")->getString(0);
assert(quantity_string == "THP");
}
//Artificial lift
if (itemValid(header, "ALQ_DEF")) {
std::string alq_string = header->getItem("ALQ_DEF")->getString(0);
if (alq_string == "GRAT") {
alq_type_ = ALQ_GRAT;
}
else if (alq_string == "IGLR") {
alq_type_ = ALQ_IGLR;
}
else if (alq_string == "TGLR") {
alq_type_ = ALQ_TGLR;
}
else if (alq_string == "PUMP") {
alq_type_ = ALQ_PUMP;
}
else if (alq_string == "COMP") {
alq_type_ = ALQ_COMP;
}
else if (alq_string == "BEAN") {
alq_type_ = ALQ_BEAN;
}
else if (alq_string == "UNDEF") {
alq_type_ = ALQ_UNDEF;
}
else {
alq_type_ = ALQ_INVALID;
OPM_THROW(std::runtime_error, "Invalid ALQ_DEF string: '" << alq_string << "'");
}
}
else {
alq_type_ = ALQ_UNDEF;
}
//Units used for this table
if (itemValid(header, "UNITS")) {
//TODO: Should check that table unit matches rest of deck.
std::string unit_string = header->getItem("UNITS")->getString(0);
if (unit_string == "METRIC") {
}
else if (unit_string == "FIELD") {
}
else if (unit_string == "LAB") {
}
else if (unit_string == "PVT-M") {
}
else {
OPM_THROW(std::runtime_error, "Invalid UNITS string: '" << unit_string << "'");
}
}
else {
//Do nothing, table implicitly same unit as rest of deck
}
//Quantity in the body of the table
if (itemValid(header, "BODY_DEF")) {
std::string body_string = header->getItem("BODY_DEF")->getString(0);
if (body_string == "TEMP") {
OPM_THROW(std::logic_error, "Invalid BODY_DEF string: TEMP not supported");
}
else if (body_string == "BHP") {
}
else {
OPM_THROW(std::runtime_error, "Invalid BODY_DEF string: '" << body_string << "'");
}
}
else {
//Default to BHP
}
//Get actual rate / flow values
flo_data_ = (*iter++)->getItem("FLOW_VALUES")->getSIDoubleData();
//Get actual tubing head pressure values
thp_data_ = (*iter++)->getItem("THP_VALUES")->getSIDoubleData();
//Get actual water fraction values
wfr_data_ = (*iter++)->getItem("WFR_VALUES")->getRawDoubleData(); //FIXME: unit
//Get actual gas fraction values
gfr_data_ = (*iter++)->getItem("GFR_VALUES")->getRawDoubleData(); //FIXME: unit
//Get actual gas fraction values
alq_data_ = (*iter++)->getItem("ALQ_VALUES")->getRawDoubleData(); //FIXME: unit
//Finally, read the actual table itself.
size_t nt = thp_data_.size();
size_t nw = wfr_data_.size();
size_t ng = gfr_data_.size();
size_t na = alq_data_.size();
size_t nf = flo_data_.size();
extents shape;
shape[0] = nt;
shape[1] = nw;
shape[2] = ng;
shape[3] = na;
shape[4] = nf;
data_.resize(shape);
for (; iter!=table->end(); ++iter) {
//Get indices (subtract 1 to get 0-based index)
int t = (*iter)->getItem("THP_INDEX")->getInt(0) - 1;
int w = (*iter)->getItem("WFR_INDEX")->getInt(0) - 1;
int g = (*iter)->getItem("GFR_INDEX")->getInt(0) - 1;
int a = (*iter)->getItem("ALQ_INDEX")->getInt(0) - 1;
//Rest of values (bottom hole pressure or tubing head temperature) have index of flo value
const std::vector<double>& bhp_tht = (*iter)->getItem("VALUES")->getRawDoubleData(); //FIXME: unit
std::copy(bhp_tht.begin(), bhp_tht.end(), &data_[t][w][g][a][0]);
}
check();
}
VFPProperties::VFPProperties(int table_num,
double datum_depth,
FLO_TYPE flo_type,
WFR_TYPE wfr_type,
GFR_TYPE gfr_type,
ALQ_TYPE alq_type,
const std::vector<double>& flo_data,
const std::vector<double>& thp_data,
const std::vector<double>& wfr_data,
const std::vector<double>& gfr_data,
const std::vector<double>& alq_data,
array_type data
) :
table_num_(table_num),
datum_depth_(datum_depth),
flo_type_(flo_type),
wfr_type_(wfr_type),
gfr_type_(gfr_type),
alq_type_(alq_type),
flo_data_(flo_data),
thp_data_(thp_data),
wfr_data_(wfr_data),
gfr_data_(gfr_data),
alq_data_(alq_data),
data_(data) {
check();
}
void VFPProperties::check() {
//Table number
assert(table_num_ > 0);
//Misc types
assert(flo_type_ >= FLO_OIL && flo_type_ < FLO_INVALID);
assert(wfr_type_ >= WFR_WOR && wfr_type_ < WFR_INVALID);
assert(gfr_type_ >= GFR_GOR && gfr_type_ < GFR_INVALID);
assert(alq_type_ >= ALQ_GRAT && alq_type_ < ALQ_INVALID);
//Data axis size
assert(flo_data_.size() > 0);
assert(thp_data_.size() > 0);
assert(wfr_data_.size() > 0);
assert(gfr_data_.size() > 0);
assert(alq_data_.size() > 0);
//Data axis sorted?
assert(is_sorted(flo_data_.begin(), flo_data_.end()));
assert(is_sorted(thp_data_.begin(), thp_data_.end()));
assert(is_sorted(wfr_data_.begin(), wfr_data_.end()));
assert(is_sorted(gfr_data_.begin(), gfr_data_.end()));
assert(is_sorted(alq_data_.begin(), alq_data_.end()));
//Check data size matches axes
assert(data_.num_dimensions() == 5);
assert(data_.shape()[0] == thp_data_.size());
assert(data_.shape()[1] == wfr_data_.size());
assert(data_.shape()[2] == gfr_data_.size());
assert(data_.shape()[3] == alq_data_.size());
assert(data_.shape()[4] == flo_data_.size());
//Finally, check that all data is within reasonable ranges, defined to be up-to 1.0e10...
typedef array_type::size_type size_type;
for (size_type t=0; t<data_.shape()[0]; ++t) {
for (size_type w=0; w<data_.shape()[1]; ++w) {
for (size_type g=0; g<data_.shape()[2]; ++g) {
for (size_type a=0; a<data_.shape()[3]; ++a) {
for (size_type f=0; f<data_.shape()[4]; ++f) {
if (data_[t][w][g][a][f] > 1.0e10) {
//TODO: Replace with proper log message
std::cerr << "Too large value encountered in VFPPROD in ["
<< t << "," << w << "," << g << "," << a << "," << f << "]="
<< data_[t][w][g][a][f] << std::endl;
}
}
}
}
}
}
}
double VFPProperties::bhp(const double& flo, const double& thp, const double& wfr, const double& gfr, const double& alq) {
//First, find the values to interpolate between
auto flo_i = find_interp_data(flo, flo_data_);
auto thp_i = find_interp_data(thp, thp_data_);
auto wfr_i = find_interp_data(wfr, wfr_data_);
auto gfr_i = find_interp_data(gfr, gfr_data_);
auto alq_i = find_interp_data(alq, alq_data_);
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(flo_i, thp_i, wfr_i, gfr_i, alq_i);
@ -367,9 +81,9 @@ 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);
ADB wfr = getWFR(w, o, g);
ADB gfr = getGFR(w, o, g);
ADB flo = getFlo(w, o, g, table_->getFloType());
ADB wfr = getWFR(w, o, g, table_->getWFRType());
ADB gfr = getGFR(w, o, g, table_->getGFRType());
//TODO: Check ALQ type here?
@ -377,57 +91,60 @@ VFPProperties::ADB VFPProperties::bhp(const Wells& wells, const ADB& qs, const A
}
VFPProperties::ADB VFPProperties::getFlo(const ADB& aqua, const ADB& liquid, const ADB& vapour) {
switch (flo_type_) {
case FLO_OIL:
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 FLO_LIQ:
case VFPProdTable::FLO_LIQ:
//Liquid = aqua + liquid phases
return aqua + liquid;
case FLO_GAS:
case VFPProdTable::FLO_GAS:
//Gas = vapor phase
return vapour;
case FLO_INVALID: //Intentional fall-through
case VFPProdTable::FLO_INVALID: //Intentional fall-through
default:
OPM_THROW(std::logic_error, "Invalid FLO_TYPE: '" << flo_type_ << "'");
OPM_THROW(std::logic_error, "Invalid FLO_TYPE: '" << type << "'");
return ADB::null();
}
}
VFPProperties::ADB VFPProperties::getWFR(const ADB& aqua, const ADB& liquid, const ADB& vapour) {
switch(wfr_type_) {
case WFR_WOR:
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 WFR_WCT:
case VFPProdTable::WFR_WCT:
//Water cut = water / (water + oil + gas)
return aqua / (aqua + liquid + vapour);
case WFR_WGR:
case VFPProdTable::WFR_WGR:
//Water-gas ratio = water / gas
return aqua / vapour;
case WFR_INVALID: //Intentional fall-through
case VFPProdTable::WFR_INVALID: //Intentional fall-through
default:
OPM_THROW(std::logic_error, "Invalid WFR_TYPE: '" << wfr_type_ << "'");
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) {
switch(gfr_type_) {
case GFR_GOR:
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 GFR_GLR:
case VFPProdTable::GFR_GLR:
// Gas-liquid ratio = gas / (oil + water)
return vapour / (liquid + aqua);
case GFR_OGR:
case VFPProdTable::GFR_OGR:
// Oil-gas ratio = oil / gas
return liquid / vapour;
case GFR_INVALID: //Intentional fall-through
case VFPProdTable::GFR_INVALID: //Intentional fall-through
default:
OPM_THROW(std::logic_error, "Invalid GFR_TYPE: '" << flo_type_ << "'");
OPM_THROW(std::logic_error, "Invalid GFR_TYPE: '" << type << "'");
return ADB::null();
}
}
@ -483,6 +200,7 @@ 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) {
@ -496,7 +214,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] = data[ti][wi][gi][ai][fi];
}
}
}

110
opm/autodiff/VFPProperties.hpp
View File

@ -21,94 +21,26 @@
#define OPM_AUTODIFF_VFPPROPERTIES_HPP_
#include <opm/core/wells.h>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/VFPProdTable.hpp>
#include <boost/multi_array.hpp>
namespace Opm {
class VFPProdTable;
/**
* Class which linearly interpolates BHP as a function of rate type, tubing head pressure,
* water fraction, gas fraction, and artificial lift.
*/
class VFPProperties {
public:
typedef boost::multi_array<double, 5> array_type;
typedef boost::array<array_type::index, 5> extents;
typedef AutoDiffBlock<double> ADB;
///Rate type
enum FLO_TYPE {
FLO_OIL=1, //< Oil rate
FLO_LIQ, //< Liquid rate
FLO_GAS, //< Gas rate
//FLO_WG
//FLO_TM
FLO_INVALID
};
///Water fraction variable
enum WFR_TYPE {
WFR_WOR=11, //< Water-oil ratio
WFR_WCT, //< Water cut
WFR_WGR, //< Water-gas ratio
WFR_INVALID
};
///Gas fraction variable
enum GFR_TYPE {
GFR_GOR=21, //< Gas-oil ratio
GFR_GLR, //< Gas-liquid ratio
GFR_OGR, //< Oil-gas ratio
GFR_INVALID
};
///Artificial lift quantity
enum ALQ_TYPE {
ALQ_GRAT=31, //< Lift as injection rate
ALQ_IGLR, //< Injection gas-liquid ratio
ALQ_TGLR, //< Total gas-liquid ratio
ALQ_PUMP, //< Pump rating
ALQ_COMP, //< Compressor power
ALQ_BEAN, //< Choke diameter
ALQ_UNDEF, //< Undefined
ALQ_INVALID
};
/**
* Constructor
* @param table_num VFP table number
* @param datum_depth Reference depth for BHP
* @param flo_type Specifies what flo_data represents
* @param wfr_type Specifies what wfr_data represents
* @param gfr_type Specifies what gfr_data represents
* @param alq_type Specifies what alq_data represents
* @param flo_data Axis for flo_type
* @param thp_data Axis for thp_type
* @param wfr_data Axis for wfr_type
* @param gfr_data Axis for gfr_type
* @param alq_data Axis for alq_type
* @param data BHP to be interpolated. Given as a 5D array so that
* BHP = data[thp][wfr][gfr][alq][flo] for the indices thp, wfr, etc.
* Constructor, takes *no* ownership of data.
*/
VFPProperties(int table_num,
double datum_depth,
FLO_TYPE flo_type,
WFR_TYPE wfr_type,
GFR_TYPE gfr_type,
ALQ_TYPE alq_type,
const std::vector<double>& flo_data,
const std::vector<double>& thp_data,
const std::vector<double>& wfr_data,
const std::vector<double>& gfr_data,
const std::vector<double>& alq_data,
array_type data);
/**
* Constructor which parses a deck keyword and retrieves the relevant parts for a
* VFP table.
*/
VFPProperties(DeckKeywordConstPtr table);
VFPProperties(VFPProdTable* table);
/**
* Linear interpolation of bhp as function of the input parameters.
@ -153,26 +85,25 @@ public:
* Computes the flo parameter according to the flo_type_
* @return Production rate of oil, gas or liquid.
*/
ADB getFlo(const ADB& aqua, const ADB& liquid, const ADB& vapour);
static ADB getFlo(const ADB& aqua, const ADB& liquid, const ADB& vapour,
const VFPProdTable::FLO_TYPE& type);
/**
* Computes the wfr parameter according to the wfr_type_
* @return Production rate of oil, gas or liquid.
*/
ADB getWFR(const ADB& aqua, const ADB& liquid, const ADB& vapour);
static ADB getWFR(const ADB& aqua, const ADB& liquid, const ADB& vapour,
const VFPProdTable::WFR_TYPE& type);
/**
* Computes the gfr parameter according to the gfr_type_
* @return Production rate of oil, gas or liquid.
*/
ADB getGFR(const ADB& aqua, const ADB& liquid, const ADB& vapour);
static ADB getGFR(const ADB& aqua, const ADB& liquid, const ADB& vapour,
const VFPProdTable::GFR_TYPE& type);
private:
/**
* Debug function that runs a series of asserts to check for sanity of inputs.
* Called after constructor to check that everything looks ok.
*/
void check();
VFPProdTable* table_;
/**
* Helper struct for linear interpolation
@ -194,23 +125,6 @@ private:
double interpolate(const InterpData& flo_i, const InterpData& thp_i,
const InterpData& wfr_i, const InterpData& gfr_i, const InterpData& alq_i);
//"Header" variables
int table_num_;
double datum_depth_;
FLO_TYPE flo_type_;
WFR_TYPE wfr_type_;
GFR_TYPE gfr_type_;
ALQ_TYPE alq_type_;
//The actual table axes
std::vector<double> flo_data_;
std::vector<double> thp_data_;
std::vector<double> wfr_data_;
std::vector<double> gfr_data_;
std::vector<double> alq_data_;
//The data itself
array_type data_;
};
}

337
tests/test_vfpproperties.cpp
View File

@ -38,6 +38,7 @@
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/EclipseState/checkDeck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/VFPProdTable.hpp>
#include <opm/autodiff/VFPProperties.hpp>
@ -78,30 +79,13 @@ struct TrivialFixture {
/**
* Fills our interpolation data with zeros
*/
void fillZeros() {
void fillData(double value) {
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] = 0.0;
}
}
}
}
}
}
/**
* Fills our interpolation data with ones
*/
void fillOnes() {
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] = 1.0;
data[i][j][k][l][m] = value;
}
}
}
@ -112,7 +96,7 @@ struct TrivialFixture {
/**
* Fills our interpolation data with an ND plane
*/
void fillPlane() {
void fillDataPlane() {
for (int i=0; i<nx; ++i) {
double x = i / static_cast<double>(nx-1);
for (int j=0; j<ny; ++j) {
@ -132,6 +116,31 @@ struct TrivialFixture {
}
}
void initTable() {
}
void initProperties() {
//Initialize table
table.init(1,
1000.0,
Opm::VFPProdTable::FLO_OIL,
Opm::VFPProdTable::WFR_WOR,
Opm::VFPProdTable::GFR_GOR,
Opm::VFPProdTable::ALQ_UNDEF,
flo_axis,
thp_axis,
wfr_axis,
gfr_axis,
alq_axis,
data);
//Initialize properties that use the table
properties.reset(new Opm::VFPProperties(&table));
}
std::shared_ptr<Opm::VFPProperties> properties;
private:
const std::vector<double> thp_axis;
const std::vector<double> wfr_axis;
const std::vector<double> gfr_axis;
@ -142,8 +151,9 @@ struct TrivialFixture {
int nz;
int nu;
int nv;
Opm::VFPProperties::extents size;
Opm::VFPProperties::array_type data;
Opm::VFPProdTable::extents size;
Opm::VFPProdTable::array_type data;
Opm::VFPProdTable table;
};
@ -159,21 +169,8 @@ BOOST_FIXTURE_TEST_SUITE( TrivialTests, TrivialFixture )
*/
BOOST_AUTO_TEST_CASE(InterpolateZero)
{
fillZeros();
//Create table
Opm::VFPProperties table(1,
1000.0,
Opm::VFPProperties::FLO_OIL,
Opm::VFPProperties::WFR_WOR,
Opm::VFPProperties::GFR_GOR,
Opm::VFPProperties::ALQ_UNDEF,
flo_axis,
thp_axis,
wfr_axis,
gfr_axis,
alq_axis,
data);
fillData(0.0);
initProperties();
//Check interpolation
double sum = 0.0;
@ -190,7 +187,7 @@ BOOST_AUTO_TEST_CASE(InterpolateZero)
const double v = m / static_cast<double>(n-1);
//Note order of arguments!
sum += table.bhp(v, x, y, z, u);
sum += properties->bhp(v, x, y, z, u);
}
}
}
@ -207,21 +204,8 @@ BOOST_AUTO_TEST_CASE(InterpolateZero)
*/
BOOST_AUTO_TEST_CASE(InterpolateOne)
{
fillOnes();
//Create table
Opm::VFPProperties table(1,
1000.0,
Opm::VFPProperties::FLO_OIL,
Opm::VFPProperties::WFR_WOR,
Opm::VFPProperties::GFR_GOR,
Opm::VFPProperties::ALQ_UNDEF,
flo_axis,
thp_axis,
wfr_axis,
gfr_axis,
alq_axis,
data);
fillData(1.0);
initProperties();
//Check interpolation
double sum = 0.0;
@ -238,7 +222,7 @@ BOOST_AUTO_TEST_CASE(InterpolateOne)
const double v = m / static_cast<double>(n-1);
//Note order of arguments!
sum += table.bhp(v, x, y, z, u);
sum += properties->bhp(v, x, y, z, u);
}
}
}
@ -256,21 +240,8 @@ BOOST_AUTO_TEST_CASE(InterpolateOne)
*/
BOOST_AUTO_TEST_CASE(InterpolatePlane)
{
fillPlane();
//Create table
Opm::VFPProperties table(1,
1000.0,
Opm::VFPProperties::FLO_OIL,
Opm::VFPProperties::WFR_WOR,
Opm::VFPProperties::GFR_GOR,
Opm::VFPProperties::ALQ_UNDEF,
flo_axis,
thp_axis,
wfr_axis,
gfr_axis,
alq_axis,
data);
fillDataPlane();
initProperties();
//Check interpolation
double sum = 0.0;
@ -292,7 +263,7 @@ BOOST_AUTO_TEST_CASE(InterpolatePlane)
reference_sum += reference;
//Note order of arguments!
double value = table.bhp(v, x, y, z, u);
double value = properties->bhp(v, x, y, z, u);
sum += value;
double abs_diff = std::abs(value - reference);
@ -318,21 +289,8 @@ BOOST_AUTO_TEST_CASE(InterpolatePlane)
*/
BOOST_AUTO_TEST_CASE(ExtrapolatePlane)
{
fillPlane();
//Create table
Opm::VFPProperties table(1,
1000.0,
Opm::VFPProperties::FLO_OIL,
Opm::VFPProperties::WFR_WOR,
Opm::VFPProperties::GFR_GOR,
Opm::VFPProperties::ALQ_UNDEF,
flo_axis,
thp_axis,
wfr_axis,
gfr_axis,
alq_axis,
data);
fillDataPlane();
initProperties();
//Check linear extrapolation (i.e., using values of x, y, etc. outside our interpolant domain)
double sum = 0.0;
@ -354,7 +312,7 @@ BOOST_AUTO_TEST_CASE(ExtrapolatePlane)
reference_sum += reference;
//Note order of arguments!
double value = table.bhp(v, x, y, z, u);
double value = properties->bhp(v, x, y, z, u);
sum += value;
double abs_diff = std::abs(value - reference);
@ -381,21 +339,8 @@ BOOST_AUTO_TEST_CASE(ExtrapolatePlane)
*/
BOOST_AUTO_TEST_CASE(ExtrapolatePlaneADB)
{
fillPlane();
//Create table
Opm::VFPProperties table(1,
1000.0,
Opm::VFPProperties::FLO_OIL,
Opm::VFPProperties::WFR_WOR,
Opm::VFPProperties::GFR_GOR,
Opm::VFPProperties::ALQ_UNDEF,
flo_axis,
thp_axis,
wfr_axis,
gfr_axis,
alq_axis,
data);
fillDataPlane();
initProperties();
//Temporary variables used to represent independent wells
const int num_wells = 5;
@ -435,7 +380,7 @@ BOOST_AUTO_TEST_CASE(ExtrapolatePlaneADB)
ADB gfr = ADB::constant(zz);
ADB alq = ADB::constant(uu);
ADB bhp_val = table.bhp(flo, thp, wfr, gfr, alq);
ADB bhp_val = properties->bhp(flo, thp, wfr, gfr, alq);
double value = 0.0;
double reference = 0.0;
@ -470,21 +415,8 @@ BOOST_AUTO_TEST_CASE(ExtrapolatePlaneADB)
*/
BOOST_AUTO_TEST_CASE(Wells_Qs_thp_and_alq_interpolation)
{
fillPlane();
//Create table
Opm::VFPProperties table(1,
1000.0,
Opm::VFPProperties::FLO_OIL,
Opm::VFPProperties::WFR_WOR,
Opm::VFPProperties::GFR_GOR,
Opm::VFPProperties::ALQ_UNDEF,
flo_axis,
thp_axis,
wfr_axis,
gfr_axis,
alq_axis,
data);
fillDataPlane();
initProperties();
//Create wells
const int nphases = 3;
@ -528,7 +460,7 @@ BOOST_AUTO_TEST_CASE(Wells_Qs_thp_and_alq_interpolation)
ADB alq = ADB::constant(alq_vals);
//Call the bhp function
ADB bhp = table.bhp(*wells, qs, thp, alq);
ADB bhp = properties->bhp(*wells, qs, thp, alq);
//Calculate reference
//First, find the three phases
@ -588,46 +520,43 @@ BOOST_AUTO_TEST_SUITE_END() // Trivial tests
struct ConversionFixture : public TrivialFixture {
ConversionFixture() : TrivialFixture(),
struct ConversionFixture {
typedef Opm::VFPProperties::ADB ADB;
ConversionFixture() :
num_wells(5),
aqua_v(num_wells),
liquid_v(num_wells),
vapour_v(num_wells),
d_aqua(num_wells),
d_liquid(num_wells),
d_vapour(num_wells),
aqua(ADB::null()),
liquid(ADB::null()),
vapour(ADB::null())
{
fillPlane();
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;
}
//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() {
}
void fillPhases(double aqua_in, double liquid_in, double vapour_in) {
for (int i=0; i<num_wells; ++i) {
aqua_v[i] = aqua_in;
liquid_v[i] = liquid_in;
vapour_v[i] = vapour_in;
}
aqua = ADB::constant(aqua_v);
liquid = ADB::constant(liquid_v);
vapour = ADB::constant(vapour_v);
}
void checkEqual(const ADB& lhs, double rhs) {
const ADB::V& values = lhs.value();
for (int i=0; i<values.size(); ++i) {
BOOST_CHECK_EQUAL(static_cast<double>(values[i]), rhs);
}
}
int num_wells;
ADB::V aqua_v;
ADB::V liquid_v;
ADB::V vapour_v;
std::vector<double> d_aqua;
std::vector<double> d_liquid;
std::vector<double> d_vapour;
ADB aqua;
ADB liquid;
@ -643,57 +572,49 @@ BOOST_FIXTURE_TEST_SUITE( ConversionTests, ConversionFixture )
BOOST_AUTO_TEST_CASE(getFlo)
{
double water = 3.0;
double oil = 5.0;
double gas = 7.0;
//Compute reference solutions
std::vector<double> ref_flo_oil;
std::vector<double> ref_flo_liq;
std::vector<double> ref_flo_gas;
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];
}
Opm::VFPProperties::FLO_TYPE types[] = {
Opm::VFPProperties::FLO_OIL,
Opm::VFPProperties::FLO_LIQ,
Opm::VFPProperties::FLO_GAS
};
double reference[] = {
oil,
water+oil,
gas
};
const int num_types = sizeof(types) / sizeof(types[0]);
for (int i=0; i<num_types; ++i) {
//Create table
Opm::VFPProperties vfp_table(1,
1000.0,
types[i],
Opm::VFPProperties::WFR_WOR,
Opm::VFPProperties::GFR_GOR,
Opm::VFPProperties::ALQ_UNDEF,
flo_axis,
thp_axis,
wfr_axis,
gfr_axis,
alq_axis,
data);
fillPhases(water, oil, gas);
ADB flo = vfp_table.getFlo(aqua, liquid, vapour);
{
ADB flo = Opm::VFPProperties::getFlo(aqua, liquid, vapour, Opm::VFPProdTable::FLO_OIL);
const double* computed = &flo.value()[0];
BOOST_CHECK_EQUAL_COLLECTIONS(ref_flo_oil.begin(), ref_flo_oil.end(), computed, computed+num_wells);
}
BOOST_TEST_MESSAGE("Using FLO_TYPE=" << types[i]);
checkEqual(flo, reference[i]);
{
ADB flo = Opm::VFPProperties::getFlo(aqua, liquid, vapour, Opm::VFPProdTable::FLO_LIQ);
const double* computed = &flo.value()[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];
BOOST_CHECK_EQUAL_COLLECTIONS(ref_flo_gas.begin(), ref_flo_gas.end(), computed, computed+num_wells);
}
}
/*
BOOST_AUTO_TEST_CASE(getWFR)
{
double water = 3.0;
double oil = 5.0;
double gas = 7.0;
Opm::VFPProperties::WFR_TYPE types[] = {
Opm::VFPProperties::WFR_WOR,
Opm::VFPProperties::WFR_WCT,
Opm::VFPProperties::WFR_WGR
Opm::VFPProdTable::WFR_TYPE types[] = {
Opm::VFPProdTable::WFR_WOR,
Opm::VFPProdTable::WFR_WCT,
Opm::VFPProdTable::WFR_WGR
};
double reference[] = {
water / oil,
@ -703,22 +624,8 @@ BOOST_AUTO_TEST_CASE(getWFR)
const int num_types = sizeof(types) / sizeof(types[0]);
for (int i=0; i<num_types; ++i) {
//Create table
Opm::VFPProperties vfp_table(1,
1000.0,
Opm::VFPProperties::FLO_OIL,
types[i],
Opm::VFPProperties::GFR_GOR,
Opm::VFPProperties::ALQ_UNDEF,
flo_axis,
thp_axis,
wfr_axis,
gfr_axis,
alq_axis,
data);
fillPhases(water, oil, gas);
ADB flo = vfp_table.getWFR(aqua, liquid, vapour);
ADB flo = Opm::VFPProperties::getWFR(aqua, liquid, vapour, types[i]);
BOOST_TEST_MESSAGE("Using WFR_TYPE=" << types[i]);
checkEqual(flo, reference[i]);
@ -728,9 +635,6 @@ BOOST_AUTO_TEST_CASE(getWFR)
BOOST_AUTO_TEST_CASE(getGFR)
{
double water = 3.0;
double oil = 5.0;
double gas = 7.0;
Opm::VFPProperties::GFR_TYPE types[] = {
Opm::VFPProperties::GFR_GOR,
@ -745,22 +649,8 @@ BOOST_AUTO_TEST_CASE(getGFR)
const int num_types = sizeof(types) / sizeof(types[0]);
for (int i=0; i<num_types; ++i) {
//Create table
Opm::VFPProperties vfp_table(1,
1000.0,
Opm::VFPProperties::FLO_OIL,
Opm::VFPProperties::WFR_WOR,
types[i],
Opm::VFPProperties::ALQ_UNDEF,
flo_axis,
thp_axis,
wfr_axis,
gfr_axis,
alq_axis,
data);
fillPhases(water, oil, gas);
ADB flo = vfp_table.getGFR(aqua, liquid, vapour);
ADB flo = Opm::VFPProperties::getGFR(aqua, liquid, vapour);
BOOST_TEST_MESSAGE("Using GFR_TYPE=" << types[i]);
checkEqual(flo, reference[i]);
@ -768,6 +658,7 @@ BOOST_AUTO_TEST_CASE(getGFR)
}
*/
BOOST_AUTO_TEST_SUITE_END() // unit tests
@ -789,7 +680,6 @@ BOOST_AUTO_TEST_SUITE_END() // unit tests
BOOST_AUTO_TEST_SUITE(IntegrationTests)
extern const double reference[];
@ -800,6 +690,7 @@ extern const double reference[];
BOOST_AUTO_TEST_CASE(ParseVFPProdAndInterpolate)
{
Opm::DeckConstPtr deck;
std::shared_ptr<Opm::UnitSystem> units(Opm::UnitSystem::newMETRIC());
Opm::ParserPtr parser(new Opm::Parser());
boost::filesystem::path file("tests/VFPPROD1");
@ -808,13 +699,17 @@ BOOST_AUTO_TEST_CASE(ParseVFPProdAndInterpolate)
Opm::checkDeck(deck);
BOOST_REQUIRE(deck->hasKeyword("VFPPROD"));
BOOST_CHECK_EQUAL(deck->numKeywords("VFPPROD"), 2);
std::vector<Opm::VFPProperties> tables;
std::vector<Opm::VFPProdTable> tables;
std::vector<Opm::VFPProperties> properties;
int num_tables = deck->numKeywords("VFPPROD");
for (int i=0; i<num_tables; ++i) {
Opm::DeckKeywordConstPtr table = deck->getKeyword("VFPPROD" , i);
tables.push_back(Opm::VFPProperties(table));
Opm::DeckKeywordConstPtr keyword = deck->getKeyword("VFPPROD" , i);
tables.push_back(Opm::VFPProdTable());
tables[i].init(keyword, units);
properties.push_back(Opm::VFPProperties(&tables[i]));
}
//Do some rudimentary testing
@ -855,7 +750,7 @@ BOOST_AUTO_TEST_CASE(ParseVFPProdAndInterpolate)
//Value given as BARSA
//FIXME: should convert to Pascal when proper conversion in VFPProperties.cpp is in place
double value_i = tables[0].bhp(f_i, t_i, w_i, g_i, a_i);
double value_i = properties[0].bhp(f_i, t_i, w_i, g_i, a_i);
double abs_diff = std::abs(value_i - reference[i]);
sad += abs_diff;