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

Added initial calculation of derivatives

Browse Source
This commit is contained in:
babrodtk 2015-08-04 11:24:23 +02:00
parent 0d36d81e51
commit f424a26651
2 changed files with 128 additions and 51 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

150
opm/autodiff/VFPProperties.cpp
View File

@ -320,32 +320,36 @@ const VFPProdTable* VFPProdProperties::getProdTable(int table_id) const {
VFPProdProperties::InterpData VFPProdProperties::find_interp_data(const double& value, const std::vector<double>& values) {
InterpData retval;
//First element greater than or equal to value
//Start with the second element, so that floor_iter does not go out of range
//Don't access out-of-range, therefore values.end()-1
auto ceil_iter = std::lower_bound(values.begin()+1, values.end()-1, value);
//Find last element smaller than range
auto floor_iter = ceil_iter-1;
//Find the indices
const int a = floor_iter - values.begin();
const int b = ceil_iter - values.begin();
const int max_size = std::max(static_cast<int>(values.size()) - 1, 0);
//Clamp indices to range of vector
retval.ind_[0] = a;
retval.ind_[1] = std::min(b, max_size);
//Find interpolation ratio
double dist = (*ceil_iter - *floor_iter);
if (dist > 0.0) {
//Possible source for floating point error here if value and floor are large,
//but very close to each other
retval.factor_ = (value-*floor_iter) / dist;
//If we only have one value in our vector, return that
if (values.size() == 1) {
retval.ind_[0] = 0;
retval.ind_[1] = 0;
retval.factor_ = 0.0;
}
// Else search in the vector
else {
retval.factor_ = 1.0;
//First element greater than or equal to value
//Start with the second element, so that floor_iter does not go out of range
//Don't access out-of-range, therefore values.end()-1
auto ceil_iter = std::lower_bound(values.begin()+1, values.end()-1, value);
//Find last element smaller than range
auto floor_iter = ceil_iter-1;
//Find the indices
retval.ind_[0] = floor_iter - values.begin();
retval.ind_[1] = ceil_iter - values.begin();
//Find interpolation ratio
double dist = (*ceil_iter - *floor_iter);
if (std::abs(dist) > 0.0) {
//Possible source for floating point error here if value and floor are large,
//but very close to each other
retval.factor_ = (value-*floor_iter) / dist;
}
else {
retval.factor_ = 0.0;
}
}
return retval;
@ -363,8 +367,16 @@ double VFPProdProperties::interpolate(const VFPProdTable::array_type& array,
const InterpData& wfr_i,
const InterpData& gfr_i,
const InterpData& alq_i) {
//Values in our hypercube
double nn[2][2][2][2][2];
//Derivatives in our hypercube
double dthp[2][2][2][2][2];
double dwfr[2][2][2][2][2];
double dgfr[2][2][2][2][2];
double dalq[2][2][2][2][2];
double dflo[2][2][2][2][2];
//Pick out nearest neighbors (nn) to our evaluation point
//This is not really required, but performance-wise it may pay off, since the 32-elements
//we copy to (nn) will fit better in cache than the full original table for the
@ -390,44 +402,98 @@ double VFPProdProperties::interpolate(const VFPProdTable::array_type& array,
}
}
//Remove dimensions one by one
// Example: going from 3D to 2D to 1D, we start by interpolating along
// the z axis first, leaving a 2D problem. Then interpolating along the y
// axis, leaving a 1D, problem, etc.
double tf = flo_i.factor_;
for (int t=0; t<=1; ++t) {
for (int w=0; w<=1; ++w) {
for (int g=0; g<=1; ++g) {
for (int a=0; a<=1; ++a) {
nn[t][w][g][a][0] = (1.0-tf)*nn[t][w][g][a][0] + tf*nn[t][w][g][a][1];
//Calculate derivatives
for (int i=0; i<=1; ++i) {
for (int j=0; j<=1; ++j) {
for (int k=0; k<=1; ++k) {
for (int l=0; l<=1; ++l) {
dthp[0][i][j][k][l] = nn[1][i][j][k][l] - nn[0][i][j][k][l];
dwfr[i][0][j][k][l] = nn[i][1][j][k][l] - nn[i][0][j][k][l];
dgfr[i][j][0][k][l] = nn[i][j][1][k][l] - nn[i][j][0][k][l];
dalq[i][j][k][0][l] = nn[i][j][k][1][l] - nn[i][j][k][0][l];
dflo[i][j][k][l][0] = nn[i][j][k][l][1] - nn[i][j][k][l][0];
//For simplicity of the rest of the interpolation code,
//we copy the derivatives in the full hypercube
dthp[1][i][j][k][l] = dthp[0][i][j][k][l];
dwfr[i][1][j][k][l] = dwfr[i][0][j][k][l];
dgfr[i][j][1][k][l] = dgfr[i][j][0][k][l];
dalq[i][j][k][1][l] = dalq[i][j][k][0][l];
dflo[i][j][k][l][1] = dflo[i][j][k][l][0];
}
}
}
}
tf = alq_i.factor_;
double a, b; //interpolation variables, so that a = (1-t) and b = t.
//Remove dimensions one by one
// Example: going from 3D to 2D to 1D, we start by interpolating along
// the z axis first, leaving a 2D problem. Then interpolating along the y
// axis, leaving a 1D, problem, etc.
b = flo_i.factor_;
a = (1.0-b);
for (int t=0; t<=1; ++t) {
for (int w=0; w<=1; ++w) {
for (int g=0; g<=1; ++g) {
nn[t][w][g][0][0] = (1.0-tf)*nn[t][w][g][0][0] + tf*nn[t][w][g][1][0];
for (int a=0; a<=1; ++a) {
nn[t][w][g][a][0] = a*nn[t][w][g][a][0] + b*nn[t][w][g][a][1];
dthp[t][w][g][a][0] = a*dthp[t][w][g][a][0] + b*dthp[t][w][g][a][1];
dwfr[t][w][g][a][0] = a*dwfr[t][w][g][a][0] + b*dwfr[t][w][g][a][1];
dgfr[t][w][g][a][0] = a*dgfr[t][w][g][a][0] + b*dgfr[t][w][g][a][1];
dalq[t][w][g][a][0] = a*dalq[t][w][g][a][0] + b*dalq[t][w][g][a][1];
dflo[t][w][g][a][0] = a*dflo[t][w][g][a][0] + b*dflo[t][w][g][a][1];
}
}
}
}
tf = gfr_i.factor_;
b = alq_i.factor_;
a = (1.0-b);
for (int t=0; t<=1; ++t) {
for (int w=0; w<=1; ++w) {
nn[t][w][0][0][0] = (1.0-tf)*nn[t][w][0][0][0] + tf*nn[t][w][1][0][0];
for (int g=0; g<=1; ++g) {
nn[t][w][g][0][0] = a*nn[t][w][g][0][0] + b*nn[t][w][g][1][0];
dthp[t][w][g][0][0] = a*dthp[t][w][g][0][0] + b*dthp[t][w][g][1][0];
dwfr[t][w][g][0][0] = a*dwfr[t][w][g][0][0] + b*dwfr[t][w][g][1][0];
dgfr[t][w][g][0][0] = a*dgfr[t][w][g][0][0] + b*dgfr[t][w][g][1][0];
dalq[t][w][g][0][0] = a*dalq[t][w][g][0][0] + b*dalq[t][w][g][1][0];
dflo[t][w][g][0][0] = a*dflo[t][w][g][0][0] + b*dflo[t][w][g][1][0];
}
}
}
tf = wfr_i.factor_;
b = gfr_i.factor_;
a = (1.0-b);
for (int t=0; t<=1; ++t) {
nn[t][0][0][0][0] = (1.0-tf)*nn[t][0][0][0][0] + tf*nn[t][1][0][0][0];
for (int w=0; w<=1; ++w) {
nn[t][w][0][0][0] = a*nn[t][w][0][0][0] + b*nn[t][w][1][0][0];
dthp[t][w][0][0][0] = a*dthp[t][w][0][0][0] + b*dthp[t][w][1][0][0];
dwfr[t][w][0][0][0] = a*dwfr[t][w][0][0][0] + b*dwfr[t][w][1][0][0];
dgfr[t][w][0][0][0] = a*dgfr[t][w][0][0][0] + b*dgfr[t][w][1][0][0];
dalq[t][w][0][0][0] = a*dalq[t][w][0][0][0] + b*dalq[t][w][1][0][0];
dflo[t][w][0][0][0] = a*dflo[t][w][0][0][0] + b*dflo[t][w][1][0][0];
}
}
tf = thp_i.factor_;
return (1.0-tf)*nn[0][0][0][0][0] + tf*nn[1][0][0][0][0];
b = wfr_i.factor_;
a = (1.0-b);
for (int t=0; t<=1; ++t) {
nn[t][0][0][0][0] = a*nn[t][0][0][0][0] + b*nn[t][1][0][0][0];
dthp[t][0][0][0][0] = a*dthp[t][0][0][0][0] + b*dthp[t][1][0][0][0];
dwfr[t][0][0][0][0] = a*dwfr[t][0][0][0][0] + b*dwfr[t][1][0][0][0];
dgfr[t][0][0][0][0] = a*dgfr[t][0][0][0][0] + b*dgfr[t][1][0][0][0];
dalq[t][0][0][0][0] = a*dalq[t][0][0][0][0] + b*dalq[t][1][0][0][0];
dflo[t][0][0][0][0] = a*dflo[t][0][0][0][0] + b*dflo[t][1][0][0][0];
}
b = thp_i.factor_;
a = (1.0-b);
return a*nn[0][0][0][0][0] + b*nn[1][0][0][0][0];
}
#ifdef __GNUC__

29
tests/test_vfpproperties.cpp
View File

@ -798,8 +798,10 @@ VFPPROD \n\
Opm::VFPProdProperties properties(&table);
const int n = 5; //Number of points to check per axis
double sad = 0.0; //Sum of absolute difference
double max_d = 0.0; //Maximum difference
double bhp_sad = 0.0; //Sum of absolute difference
double bhp_max_d = 0.0; //Maximum difference
double thp_sad = 0.0;
double thp_max_d = 0.0;
for (int w=0; w<n; ++w) { //water
for (int o=0; o<n; ++o) { //oil
for (int g=0; g<n; ++g) { //gas
@ -811,20 +813,29 @@ VFPPROD \n\
double thp = t * 456.78;
double alq = a * 42.24;
double bhp = properties.bhp(42, aqua, liquid, vapour, thp, alq);
double ref_bhp = thp;
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;
double diff = std::abs(bhp - ref_bhp);
sad += diff;
max_d = std::max(diff, max_d);
double bhp_diff = std::abs(bhp_interp - bhp_ref);
bhp_sad += bhp_diff;
bhp_max_d = std::max(bhp_diff, bhp_max_d);
double thp_diff = std::abs(thp_interp - thp_ref);
thp_sad += thp_diff;
thp_max_d = std::max(thp_diff, thp_max_d);
}
}
}
}
}
BOOST_CHECK_SMALL(max_d, max_d_tol);
BOOST_CHECK_SMALL(sad, sad_tol);
BOOST_CHECK_SMALL(bhp_max_d, max_d_tol);
BOOST_CHECK_SMALL(bhp_sad, sad_tol);
BOOST_CHECK_SMALL(thp_max_d, max_d_tol);
BOOST_CHECK_SMALL(thp_sad, sad_tol);
}
/**