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

Browse Source
This commit is contained in:
babrodtk 2015-08-11 12:21:06 +02:00
parent 08dd631a8d
commit 2994d1d932
5 changed files with 237 additions and 131 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
opm/autodiff/BlackoilModelBase_impl.hpp
View File

@ -1246,7 +1246,7 @@ namespace detail {
const double& alq = well_controls_iget_alq(wc, current);
//Set *BHP* target by calculating bhp from THP
xw.bhp()[w] = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq).value;
xw.bhp()[w] = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq);
break;
}

113
opm/autodiff/VFPHelpers.hpp
View File

@ -73,6 +73,7 @@ inline ADB zeroIfNan(const ADB& values) {
/**
* Computes the flo parameter according to the flo_type_
* for production tables
* @return Production rate of oil, gas or liquid.
*/
template <typename T>
@ -96,6 +97,32 @@ static T getFlo(const T& aqua, const T& liquid, const T& vapour,
/**
* Computes the flo parameter according to the flo_type_
* for injection tables
* @return Production rate of oil, gas or liquid.
*/
template <typename T>
static T getFlo(const T& aqua, const T& liquid, const T& vapour,
const VFPInjTable::FLO_TYPE& type) {
switch (type) {
case VFPInjTable::FLO_OIL:
//Oil = liquid phase
return liquid;
case VFPInjTable::FLO_WAT:
//Liquid = aqua phase
return aqua;
case VFPInjTable::FLO_GAS:
//Gas = vapor phase
return vapour;
case VFPInjTable::FLO_INVALID: //Intentional fall-through
default:
OPM_THROW(std::logic_error, "Invalid FLO_TYPE: '" << type << "'");
}
}
@ -222,6 +249,60 @@ inline InterpData findInterpData(const double& value, const std::vector<double>&
/**
* An "ADB-like" structure with a single value and a set of derivatives
*/
struct adb_like {
adb_like() : value(0.0), dthp(0.0), dwfr(0.0), dgfr(0.0), dalq(0.0), dflo(0.0) {};
double value;
double dthp;
double dwfr;
double dgfr;
double dalq;
double dflo;
};
inline adb_like operator+(
adb_like lhs,
const adb_like& rhs) {
lhs.value += rhs.value;
lhs.dthp += rhs.dthp;
lhs.dwfr += rhs.dwfr;
lhs.dgfr += rhs.dgfr;
lhs.dalq += rhs.dalq;
lhs.dflo += rhs.dflo;
return lhs;
}
inline adb_like operator-(
adb_like lhs,
const adb_like& rhs) {
lhs.value -= rhs.value;
lhs.dthp -= rhs.dthp;
lhs.dwfr -= rhs.dwfr;
lhs.dgfr -= rhs.dgfr;
lhs.dalq -= rhs.dalq;
lhs.dflo -= rhs.dflo;
return lhs;
}
inline adb_like operator*(
double lhs,
const adb_like& rhs) {
adb_like retval;
retval.value = rhs.value * lhs;
retval.dthp = rhs.dthp * lhs;
retval.dwfr = rhs.dwfr * lhs;
retval.dgfr = rhs.dgfr * lhs;
retval.dalq = rhs.dalq * lhs;
retval.dflo = rhs.dflo * lhs;
return retval;
}
/**
* Helper function which interpolates data using the indices etc. given in the inputs.
@ -230,7 +311,7 @@ inline InterpData findInterpData(const double& value, const std::vector<double>&
#pragma GCC push_options
#pragma GCC optimize ("unroll-loops")
#endif
inline VFPProdProperties::adb_like interpolate(
inline adb_like interpolate(
const VFPProdTable::array_type& array,
const InterpData& flo_i,
const InterpData& thp_i,
@ -239,7 +320,7 @@ inline VFPProdProperties::adb_like interpolate(
const InterpData& alq_i) {
//Values and derivatives in a 5D hypercube
VFPProdProperties::adb_like nn[2][2][2][2][2];
adb_like nn[2][2][2][2][2];
//Pick out nearest neighbors (nn) to our evaluation point
@ -347,6 +428,34 @@ inline VFPProdProperties::adb_like interpolate(
inline adb_like bhp(const VFPProdTable* table,
const double& aqua,
const double& liquid,
const double& vapour,
const double& thp,
const double& alq) {
//Find interpolation variables
double flo = detail::getFlo(aqua, liquid, vapour, table->getFloType());
double wfr = detail::getWFR(aqua, liquid, vapour, table->getWFRType());
double gfr = detail::getGFR(aqua, liquid, vapour, table->getGFRType());
//First, find the values to interpolate between
auto flo_i = detail::findInterpData(flo, table->getFloAxis());
auto thp_i = detail::findInterpData(thp, table->getTHPAxis());
auto wfr_i = detail::findInterpData(wfr, table->getWFRAxis());
auto gfr_i = detail::findInterpData(gfr, table->getGFRAxis());
auto alq_i = detail::findInterpData(alq, table->getALQAxis());
//Then perform the interpolation itself
detail::adb_like retval = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i);
return retval;
}
} // namespace detail

24
opm/autodiff/VFPProdProperties.cpp
View File

@ -353,7 +353,7 @@ VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector<int>& table_id,
auto gfr_i = detail::findInterpData(gfr.value()[i], table->getGFRAxis());
auto alq_i = detail::findInterpData(alq.value()[i], table->getALQAxis());
adb_like bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i);
detail::adb_like bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i);
value[i] = bhp_val.value;
dthp[i] = bhp_val.dthp;
@ -405,8 +405,7 @@ VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector<int>& table_id,
VFPProdProperties::adb_like VFPProdProperties::bhp(int table_id,
double VFPProdProperties::bhp(int table_id,
const double& aqua,
const double& liquid,
const double& vapour,
@ -414,27 +413,12 @@ VFPProdProperties::adb_like VFPProdProperties::bhp(int table_id,
const double& alq) const {
const VFPProdTable* table = getProdTable(table_id);
//Find interpolation variables
double flo = detail::getFlo(aqua, liquid, vapour, table->getFloType());
double wfr = detail::getWFR(aqua, liquid, vapour, table->getWFRType());
double gfr = detail::getGFR(aqua, liquid, vapour, table->getGFRType());
//First, find the values to interpolate between
auto flo_i = detail::findInterpData(flo, table->getFloAxis());
auto thp_i = detail::findInterpData(thp, table->getTHPAxis());
auto wfr_i = detail::findInterpData(wfr, table->getWFRAxis());
auto gfr_i = detail::findInterpData(gfr, table->getGFRAxis());
auto alq_i = detail::findInterpData(alq, table->getALQAxis());
//Then perform the interpolation itself
adb_like retval = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i);
return retval;
detail::adb_like retval = detail::bhp(table, aqua, liquid, vapour, thp, alq);
return retval.value;
}
double VFPProdProperties::thp(int table_id,
const double& aqua,
const double& liquid,

52
opm/autodiff/VFPProdProperties.hpp
View File

@ -42,19 +42,6 @@ class VFPProdProperties {
public:
typedef AutoDiffBlock<double> ADB;
/**
* An "ADB-like" structure with a single value and a set of derivatives
*/
struct adb_like {
adb_like() : value(0.0), dthp(0.0), dwfr(0.0), dgfr(0.0), dalq(0.0), dflo(0.0) {};
double value;
double dthp;
double dwfr;
double dgfr;
double dalq;
double dflo;
};
/**
* Empty constructor
*/
@ -126,7 +113,7 @@ public:
* @return The bottom hole pressure, interpolated/extrapolated linearly using
* the above parameters from the values in the input table.
*/
adb_like bhp(int table_id,
double bhp(int table_id,
const double& aqua,
const double& liquid,
const double& vapour,
@ -173,43 +160,6 @@ private:
inline VFPProdProperties::adb_like operator+(
VFPProdProperties::adb_like lhs,
const VFPProdProperties::adb_like& rhs) {
lhs.value += rhs.value;
lhs.dthp += rhs.dthp;
lhs.dwfr += rhs.dwfr;
lhs.dgfr += rhs.dgfr;
lhs.dalq += rhs.dalq;
lhs.dflo += rhs.dflo;
return lhs;
}
inline VFPProdProperties::adb_like operator-(
VFPProdProperties::adb_like lhs,
const VFPProdProperties::adb_like& rhs) {
lhs.value -= rhs.value;
lhs.dthp -= rhs.dthp;
lhs.dwfr -= rhs.dwfr;
lhs.dgfr -= rhs.dgfr;
lhs.dalq -= rhs.dalq;
lhs.dflo -= rhs.dflo;
return lhs;
}
inline VFPProdProperties::adb_like operator*(
double lhs,
const VFPProdProperties::adb_like& rhs) {
VFPProdProperties::adb_like retval;
retval.value = rhs.value * lhs;
retval.dthp = rhs.dthp * lhs;
retval.dwfr = rhs.dwfr * lhs;
retval.dgfr = rhs.dgfr * lhs;
retval.dalq = rhs.dalq * lhs;
retval.dflo = rhs.dflo * lhs;
return retval;
}
} //namespace

177
tests/test_vfpproperties.cpp
View File

@ -30,6 +30,7 @@
#include <memory>
#include <map>
#include <sstream>
#include <limits>
#include <boost/test/unit_test.hpp>
#include <boost/filesystem.hpp>
@ -228,7 +229,7 @@ BOOST_AUTO_TEST_SUITE_END() // unit tests
*/
struct TrivialFixture {
typedef Opm::VFPProdProperties::ADB ADB;
typedef Opm::VFPProdProperties::adb_like adb_like;
typedef Opm::detail::adb_like adb_like;
TrivialFixture() : table_ids(1, 1),
thp_axis{0.0, 1.0},
@ -299,13 +300,14 @@ struct TrivialFixture {
*/
inline void fillDataRandom() {
unsigned long randx = 42;
static double max_val = static_cast<double>(std::numeric_limits<unsigned long>::max());
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;
data[i][j][k][l][m] = randx / max_val;
randx = (randx*1103515245 + 12345);
}
}
}
@ -389,11 +391,11 @@ BOOST_AUTO_TEST_CASE(GetTable)
add_well(INJECTOR, 100, 1, NULL, cells, NULL, NULL, wells.get());
//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;
double aqua_d = 0.15;
double liquid_d = 0.25;
double vapour_d = 0.35;
double thp_d = 0.45;
double alq_d = 0.55;
ADB aqua_adb = createConstantScalarADB(aqua_d);
ADB liquid_adb = createConstantScalarADB(liquid_d);
@ -405,18 +407,24 @@ BOOST_AUTO_TEST_CASE(GetTable)
qs_adb_v << aqua_adb.value(), liquid_adb.value(), vapour_adb.value();
ADB qs_adb = ADB::constant(qs_adb_v);
//Check that our reference has not changed
Opm::detail::adb_like ref= Opm::detail::bhp(&table, aqua_d, liquid_d, vapour_d, thp_d, alq_d);
BOOST_CHECK_CLOSE(ref.value, 1.0923565702101556, max_d_tol);
BOOST_CHECK_CLOSE(ref.dthp, 0.13174065498177251, max_d_tol);
BOOST_CHECK_CLOSE(ref.dwfr, -1.2298177745501071, max_d_tol);
BOOST_CHECK_CLOSE(ref.dgfr, 0.82988935779290274, max_d_tol);
BOOST_CHECK_CLOSE(ref.dalq, 1.8148520254931713, max_d_tol);
BOOST_CHECK_CLOSE(ref.dflo, 9.0944843574181924, max_d_tol);
//Check that different versions of the prod_bph function work
ADB a = properties->bhp(table_ids, aqua_adb, liquid_adb, vapour_adb, thp_adb, alq_adb);
adb_like b = properties->bhp(table_ids[0], aqua_d, liquid_d, vapour_d, thp_d, alq_d);
double b =properties->bhp(table_ids[0], aqua_d, liquid_d, vapour_d, thp_d, alq_d);
ADB c = properties->bhp(table_ids, *wells, qs_adb, thp_adb, alq_adb);
//Check that results are actually equal
double d = a.value()[0];
double e = b.value;
double f = c.value()[0];
BOOST_CHECK_EQUAL(d, e);
BOOST_CHECK_EQUAL(d, f);
//Check that results are actually equal reference
BOOST_CHECK_EQUAL(a.value()[0], ref.value);
BOOST_CHECK_EQUAL(b, ref.value);
BOOST_CHECK_EQUAL(c.value()[0], ref.value);
//Table 2 does not exist.
std::vector<int> table_ids_wrong(1, 2);
@ -447,7 +455,7 @@ BOOST_AUTO_TEST_CASE(InterpolateZero)
const double v = m / static_cast<double>(n-1);
//Note order of arguments!
sum += properties->bhp(1, v, x, y, z, u).value;
sum += properties->bhp(1, v, x, y, z, u);
}
}
}
@ -482,7 +490,7 @@ BOOST_AUTO_TEST_CASE(InterpolateOne)
const double v = m / static_cast<double>(n-1);
//Note order of arguments!
const double value = properties->bhp(1, v, x, y, z, u).value;
const double value = properties->bhp(1, v, x, y, z, u);
sum += value;
}
@ -508,8 +516,8 @@ BOOST_AUTO_TEST_CASE(InterpolatePlane)
initProperties();
//Temps used to store reference and actual variables
adb_like sad;
adb_like max_d;
double sad = 0.0;
double max_d = 0.0;
//Check interpolation
for (int i=0; i<=n; ++i) {
@ -529,33 +537,15 @@ BOOST_AUTO_TEST_CASE(InterpolatePlane)
double gfr = Opm::detail::getGFR(aqua, liquid, vapour, table.getGFRType());
//Calculate reference
adb_like reference;
reference.value = thp + 2*wfr + 3*gfr+ 4*alq + 5*flo;
reference.dthp = 1;
reference.dwfr = 2;
reference.dgfr = 3;
reference.dalq = 4;
reference.dflo = 5;
double reference = thp + 2*wfr + 3*gfr+ 4*alq + 5*flo;
//Calculate actual
//Note order of arguments: id, aqua, liquid, vapour, thp, alq
adb_like actual = properties->bhp(1, aqua, liquid, vapour, thp, alq);
double actual = properties->bhp(1, aqua, liquid, vapour, thp, alq);
adb_like abs_diff = actual - reference;
abs_diff.value = std::abs(abs_diff.value);
abs_diff.dthp = std::abs(abs_diff.dthp);
abs_diff.dwfr = std::abs(abs_diff.dwfr);
abs_diff.dgfr = std::abs(abs_diff.dgfr);
abs_diff.dalq = std::abs(abs_diff.dalq);
abs_diff.dflo = std::abs(abs_diff.dflo);
max_d.value = std::max(max_d.value, abs_diff.value);
max_d.dthp = std::max(max_d.dthp, abs_diff.dthp);
max_d.dwfr = std::max(max_d.dwfr, abs_diff.dwfr);
max_d.dgfr = std::max(max_d.dgfr, abs_diff.dgfr);
max_d.dalq = std::max(max_d.dalq, abs_diff.dalq);
max_d.dflo = std::max(max_d.dflo, abs_diff.dflo);
double abs_diff = std::abs(actual - reference);
double max_d = std::max(max_d, abs_diff);
sad = sad + abs_diff;
}
}
@ -563,19 +553,8 @@ BOOST_AUTO_TEST_CASE(InterpolatePlane)
}
}
BOOST_CHECK_SMALL(max_d.value, max_d_tol);
BOOST_CHECK_SMALL(max_d.dthp, max_d_tol);
BOOST_CHECK_SMALL(max_d.dwfr, max_d_tol);
BOOST_CHECK_SMALL(max_d.dgfr, max_d_tol);
BOOST_CHECK_SMALL(max_d.dalq, max_d_tol);
BOOST_CHECK_SMALL(max_d.dflo, max_d_tol);
BOOST_CHECK_SMALL(sad.value, sad_tol);
BOOST_CHECK_SMALL(sad.dthp, sad_tol);
BOOST_CHECK_SMALL(sad.dwfr, sad_tol);
BOOST_CHECK_SMALL(sad.dgfr, sad_tol);
BOOST_CHECK_SMALL(sad.dalq, sad_tol);
BOOST_CHECK_SMALL(sad.dflo, sad_tol);
BOOST_CHECK_SMALL(max_d, max_d_tol);
BOOST_CHECK_SMALL(sad, sad_tol);
}
@ -616,7 +595,7 @@ BOOST_AUTO_TEST_CASE(ExtrapolatePlane)
reference_sum += reference;
//Note order of arguments! id, aqua, liquid, vapour, thp , alq
double value = properties->bhp(1, aqua, liquid, vapour, x, u).value;
double value = properties->bhp(1, aqua, liquid, vapour, x, u);
sum += value;
double abs_diff = std::abs(value - reference);
@ -820,6 +799,90 @@ BOOST_AUTO_TEST_CASE(InterpolateADBAndQs)
/**
* Test that the partial derivatives are reasonable
*/
BOOST_AUTO_TEST_CASE(PartialDerivatives)
{
const int n=5;
fillDataPlane();
initProperties();
//Temps used to store reference and actual variables
adb_like sad;
adb_like max_d;
//Check interpolation
for (int i=0; i<=n; ++i) {
const double thp = i / 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 alq = l / static_cast<double>(n);
for (int m=1; m<=n; ++m) {
const double liquid = m / static_cast<double>(n);
//Find values that should be in table
double flo = Opm::detail::getFlo(aqua, liquid, vapour, table.getFloType());
double wfr = Opm::detail::getWFR(aqua, liquid, vapour, table.getWFRType());
double gfr = Opm::detail::getGFR(aqua, liquid, vapour, table.getGFRType());
//Calculate reference
adb_like reference;
reference.value = thp + 2*wfr + 3*gfr+ 4*alq + 5*flo;
reference.dthp = 1;
reference.dwfr = 2;
reference.dgfr = 3;
reference.dalq = 4;
reference.dflo = 5;
//Calculate actual
//Note order of arguments: id, aqua, liquid, vapour, thp, alq
adb_like actual = Opm::detail::bhp(&table, aqua, liquid, vapour, thp, alq);
adb_like abs_diff = actual - reference;
abs_diff.value = std::abs(abs_diff.value);
abs_diff.dthp = std::abs(abs_diff.dthp);
abs_diff.dwfr = std::abs(abs_diff.dwfr);
abs_diff.dgfr = std::abs(abs_diff.dgfr);
abs_diff.dalq = std::abs(abs_diff.dalq);
abs_diff.dflo = std::abs(abs_diff.dflo);
max_d.value = std::max(max_d.value, abs_diff.value);
max_d.dthp = std::max(max_d.dthp, abs_diff.dthp);
max_d.dwfr = std::max(max_d.dwfr, abs_diff.dwfr);
max_d.dgfr = std::max(max_d.dgfr, abs_diff.dgfr);
max_d.dalq = std::max(max_d.dalq, abs_diff.dalq);
max_d.dflo = std::max(max_d.dflo, abs_diff.dflo);
sad = sad + abs_diff;
}
}
}
}
}
BOOST_CHECK_SMALL(max_d.value, max_d_tol);
BOOST_CHECK_SMALL(max_d.dthp, max_d_tol);
BOOST_CHECK_SMALL(max_d.dwfr, max_d_tol);
BOOST_CHECK_SMALL(max_d.dgfr, max_d_tol);
BOOST_CHECK_SMALL(max_d.dalq, max_d_tol);
BOOST_CHECK_SMALL(max_d.dflo, max_d_tol);
BOOST_CHECK_SMALL(sad.value, sad_tol);
BOOST_CHECK_SMALL(sad.dthp, sad_tol);
BOOST_CHECK_SMALL(sad.dwfr, sad_tol);
BOOST_CHECK_SMALL(sad.dgfr, sad_tol);
BOOST_CHECK_SMALL(sad.dalq, sad_tol);
BOOST_CHECK_SMALL(sad.dflo, sad_tol);
}
BOOST_AUTO_TEST_SUITE_END() // Trivial tests
@ -912,7 +975,7 @@ VFPPROD \n\
double thp = t * 456.78;
double alq = a * 42.24;
double bhp_interp = properties.bhp(42, aqua, liquid, vapour, thp, alq).value;
double bhp_interp = properties.bhp(42, aqua, liquid, vapour, thp, alq);
double bhp_ref = thp;
double thp_interp = properties.thp(42, aqua, liquid, vapour, bhp_ref, alq);
double thp_ref = thp;
@ -1012,7 +1075,7 @@ BOOST_AUTO_TEST_CASE(ParseInterpolateRealisticVFPPROD)
}
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).value * 10.0e-6;
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;