OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Make all codes use the same linear interpolation routines.

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen 2013-03-22 15:28:16 +01:00
parent dcf88207e3
commit a8097317a5
6 changed files with 44 additions and 54 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
opm/core/io/eclipse/EclipseGridParserHelpers.hpp
View File

@ -410,7 +410,7 @@ namespace
} else if (defaultInterpolation) {
// Interpolate
for (int i=0; i<int(indx.size()); ++i) {
table[k][indx[i]] = linearInterpolationExtrap(xv, yv, x[i]);
table[k][indx[i]] = linearInterpolation(xv, yv, x[i]);
}
} else {
// Interpolate

28
opm/core/props/pvt/SinglePvtLiveGas.cpp
View File

@ -37,8 +37,8 @@
namespace Opm
{
using Opm::linearInterpolationExtrap;
using Opm::linearInterpolDerivative;
using Opm::linearInterpolation;
using Opm::linearInterpolationDerivative;
//------------------------------------------------------------------------
// Member functions
@ -184,7 +184,7 @@ namespace Opm
// To handle no-gas case.
return 0.0;
}
double satR = linearInterpolationExtrap(saturated_gas_table_[0],
double satR = linearInterpolation(saturated_gas_table_[0],
saturated_gas_table_[3], press);
double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
if (satR < maxR ) {
@ -205,13 +205,13 @@ namespace Opm
dRdpval = 0.0;
return;
}
double satR = linearInterpolationExtrap(saturated_gas_table_[0],
double satR = linearInterpolation(saturated_gas_table_[0],
saturated_gas_table_[3], press);
double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
if (satR < maxR ) {
// Saturated case
Rval = satR;
dRdpval = linearInterpolDerivative(saturated_gas_table_[0],
dRdpval = linearInterpolationDerivative(saturated_gas_table_[0],
saturated_gas_table_[3],
press);
} else {
@ -227,13 +227,13 @@ namespace Opm
const bool deriv) const
{
int section;
double Rval = linearInterpolationExtrap(saturated_gas_table_[0],
double Rval = linearInterpolation(saturated_gas_table_[0],
saturated_gas_table_[3], press,
section);
double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
if (deriv) {
if (Rval < maxR ) { // Saturated case
return linearInterpolDerivative(saturated_gas_table_[0],
return linearInterpolationDerivative(saturated_gas_table_[0],
saturated_gas_table_[item],
press);
} else { // Undersaturated case
@ -246,11 +246,11 @@ namespace Opm
return val;
}
double val1 =
linearInterpolationExtrap(undersat_gas_tables_[is][0],
linearInterpolation(undersat_gas_tables_[is][0],
undersat_gas_tables_[is][item],
maxR);
double val2 =
linearInterpolationExtrap(undersat_gas_tables_[is+1][0],
linearInterpolation(undersat_gas_tables_[is+1][0],
undersat_gas_tables_[is+1][item],
maxR);
double val = (val2 - val1)/
@ -259,14 +259,14 @@ namespace Opm
}
} else {
if (Rval < maxR ) { // Saturated case
return linearInterpolationExtrap(saturated_gas_table_[0],
return linearInterpolation(saturated_gas_table_[0],
saturated_gas_table_[item],
press);
} else { // Undersaturated case
int is = section;
// Extrapolate from first table section
if (is == 0 && press < saturated_gas_table_[0][0]) {
return linearInterpolationExtrap(undersat_gas_tables_[0][0],
return linearInterpolation(undersat_gas_tables_[0][0],
undersat_gas_tables_[0][item],
maxR);
}
@ -274,7 +274,7 @@ namespace Opm
// Extrapolate from last table section
int ltp = saturated_gas_table_[0].size() - 1;
if (is+1 == ltp && press > saturated_gas_table_[0][ltp]) {
return linearInterpolationExtrap(undersat_gas_tables_[ltp][0],
return linearInterpolation(undersat_gas_tables_[ltp][0],
undersat_gas_tables_[ltp][item],
maxR);
}
@ -290,11 +290,11 @@ namespace Opm
return val;
}
double val1 =
linearInterpolationExtrap(undersat_gas_tables_[is][0],
linearInterpolation(undersat_gas_tables_[is][0],
undersat_gas_tables_[is][item],
maxR);
double val2 =
linearInterpolationExtrap(undersat_gas_tables_[is+1][0],
linearInterpolation(undersat_gas_tables_[is+1][0],
undersat_gas_tables_[is+1][item],
maxR);
double val = val1 + w*(val2 - val1);

24
opm/core/props/pvt/SinglePvtLiveOil.cpp
View File

@ -26,8 +26,8 @@
namespace Opm
{
using Opm::linearInterpolationExtrap;
using Opm::linearInterpolDerivative;
using Opm::linearInterpolation;
using Opm::linearInterpolationDerivative;
using Opm::tableIndex;
@ -254,7 +254,7 @@ namespace Opm
if (surfvol[phase_pos_[Vapour]] == 0.0) {
return 0.0;
}
double Rval = linearInterpolationExtrap(saturated_oil_table_[0],
double Rval = linearInterpolation(saturated_oil_table_[0],
saturated_oil_table_[3], press);
double maxR = surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
if (Rval < maxR ) { // Saturated case
@ -272,12 +272,12 @@ namespace Opm
dRdpval = 0.0;
return;
}
Rval = linearInterpolationExtrap(saturated_oil_table_[0],
Rval = linearInterpolation(saturated_oil_table_[0],
saturated_oil_table_[3], press);
double maxR = surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
if (Rval < maxR ) {
// Saturated case
dRdpval = linearInterpolDerivative(saturated_oil_table_[0],
dRdpval = linearInterpolationDerivative(saturated_oil_table_[0],
saturated_oil_table_[3],
press);
} else {
@ -294,13 +294,13 @@ namespace Opm
const bool deriv) const
{
int section;
double Rval = linearInterpolationExtrap(saturated_oil_table_[0],
double Rval = linearInterpolation(saturated_oil_table_[0],
saturated_oil_table_[3],
press, section);
double maxR = (surfvol[phase_pos_[Liquid]] == 0.0) ? 0.0 : surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
if (deriv) {
if (Rval < maxR ) { // Saturated case
return linearInterpolDerivative(saturated_oil_table_[0],
return linearInterpolationDerivative(saturated_oil_table_[0],
saturated_oil_table_[item],
press);
} else { // Undersaturated case
@ -310,11 +310,11 @@ namespace Opm
ASSERT(undersat_oil_tables_[is][0].size() >= 2);
ASSERT(undersat_oil_tables_[is+1][0].size() >= 2);
double val1 =
linearInterpolDerivative(undersat_oil_tables_[is][0],
linearInterpolationDerivative(undersat_oil_tables_[is][0],
undersat_oil_tables_[is][item],
press);
double val2 =
linearInterpolDerivative(undersat_oil_tables_[is+1][0],
linearInterpolationDerivative(undersat_oil_tables_[is+1][0],
undersat_oil_tables_[is+1][item],
press);
double val = val1 + w*(val2 - val1);
@ -322,7 +322,7 @@ namespace Opm
}
} else {
if (Rval < maxR ) { // Saturated case
return linearInterpolationExtrap(saturated_oil_table_[0],
return linearInterpolation(saturated_oil_table_[0],
saturated_oil_table_[item],
press);
} else { // Undersaturated case
@ -333,11 +333,11 @@ namespace Opm
ASSERT(undersat_oil_tables_[is][0].size() >= 2);
ASSERT(undersat_oil_tables_[is+1][0].size() >= 2);
double val1 =
linearInterpolationExtrap(undersat_oil_tables_[is][0],
linearInterpolation(undersat_oil_tables_[is][0],
undersat_oil_tables_[is][item],
press);
double val2 =
linearInterpolationExtrap(undersat_oil_tables_[is+1][0],
linearInterpolation(undersat_oil_tables_[is+1][0],
undersat_oil_tables_[is+1][item],
press);
double val = val1 + w*(val2 - val1);

6
opm/core/props/rock/RockCompressibility.cpp
View File

@ -70,7 +70,7 @@ namespace Opm
return (1.0 + cpnorm + 0.5*cpnorm*cpnorm);
} else {
// return Opm::linearInterpolation(p_, poromult_, pressure);
return Opm::linearInterpolationExtrap(p_, poromult_, pressure);
return Opm::linearInterpolation(p_, poromult_, pressure);
}
}
@ -81,8 +81,8 @@ namespace Opm
} else {
//const double poromult = Opm::linearInterpolation(p_, poromult_, pressure);
//const double dporomultdp = Opm::linearInterpolationDerivative(p_, poromult_, pressure);
const double poromult = Opm::linearInterpolationExtrap(p_, poromult_, pressure);
const double dporomultdp = Opm::linearInterpolationDerivativeExtrap(p_, poromult_, pressure);
const double poromult = Opm::linearInterpolation(p_, poromult_, pressure);
const double dporomultdp = Opm::linearInterpolationDerivative(p_, poromult_, pressure);
return dporomultdp/poromult;
}

31
opm/core/utility/linInt.hpp
View File

@ -66,31 +66,18 @@ namespace Opm
return jl;
}
inline double linearInterpolation(const std::vector<double>& xv,
const std::vector<double>& yv, double x)
{
// Returns end point if x is outside xv
std::vector<double>::const_iterator lb = lower_bound(xv.begin(), xv.end(), x);
int ix2 = lb - xv.begin();
if (ix2 == 0) {
return yv[0];
} else if (ix2 == int(xv.size())) {
return yv[ix2-1];
}
int ix1 = ix2 - 1;
return (yv[ix2] - yv[ix1])/(xv[ix2] - xv[ix1])*(x - xv[ix1]) + yv[ix1];
}
inline double linearInterpolDerivative(const std::vector<double>& xv,
const std::vector<double>& yv, double x)
inline double linearInterpolationDerivative(const std::vector<double>& xv,
const std::vector<double>& yv, double x)
{
// Extrapolates if x is outside xv
int ix1 = tableIndex(xv, x);
int ix2 = ix1 + 1;
return (yv[ix2] - yv[ix1])/(xv[ix2] - xv[ix1]);
}
inline double linearInterpolationExtrap(const std::vector<double>& xv,
const std::vector<double>& yv, double x)
inline double linearInterpolation(const std::vector<double>& xv,
const std::vector<double>& yv, double x)
{
// Extrapolates if x is outside xv
int ix1 = tableIndex(xv, x);
@ -98,9 +85,9 @@ namespace Opm
return (yv[ix2] - yv[ix1])/(xv[ix2] - xv[ix1])*(x - xv[ix1]) + yv[ix1];
}
inline double linearInterpolationExtrap(const std::vector<double>& xv,
const std::vector<double>& yv,
double x, int& ix1)
inline double linearInterpolation(const std::vector<double>& xv,
const std::vector<double>& yv,
double x, int& ix1)
{
// Extrapolates if x is outside xv
ix1 = tableIndex(xv, x);
@ -108,6 +95,8 @@ namespace Opm
return (yv[ix2] - yv[ix1])/(xv[ix2] - xv[ix1])*(x - xv[ix1]) + yv[ix1];
}
} // namespace Opm

7
opm/core/utility/linearInterpolation.hpp
View File

@ -35,7 +35,8 @@
#ifndef OPM_LINEARINTERPOLATION_HEADER
#define OPM_LINEARINTERPOLATION_HEADER
#include <opm/core/utility/linInt.hpp>
#if 0
#include <vector>
#include <algorithm>
@ -103,7 +104,7 @@ namespace Opm
}
template <typename T>
T linearInterpolationExtrap(const std::vector<double>& xv,
T linearInterpolation(const std::vector<double>& xv,
const std::vector<T>& yv,
double x)
{
@ -123,5 +124,5 @@ namespace Opm
} // namespace Opm
#endif
#endif // OPM_LINEARINTERPOLATION_HEADER