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

Rewrote function finding interpolation data

Browse Source
This commit is contained in:
babrodtk 2015-08-19 12:37:54 +02:00
parent 62c2373a9d
commit f3553313d2
3 changed files with 55 additions and 24 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
opm/autodiff/SimulatorBase.hpp
View File

@ -31,7 +31,6 @@
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/RateConverter.hpp>
#include <opm/autodiff/VFPProperties.hpp>
#include <opm/core/grid.h>
#include <opm/core/wells.h>

48
opm/autodiff/VFPHelpers.hpp
View File

@ -203,12 +203,15 @@ struct InterpData {
/**
* Helper function to find indices etc. for linear interpolation
* Helper function to find indices etc. for linear interpolation and extrapolation
* @param value Value to find in values
* @param values Sorted list of values to search for value in.
* @return Data required to find the interpolated value
*/
inline InterpData findInterpData(const double& value, const std::vector<double>& values) {
InterpData retval;
const double abs_value = value;//std::abs(value);
const int nvalues = values.size();
//If we only have one value in our vector, return that
if (values.size() == 1) {
@ -219,25 +222,36 @@ inline InterpData findInterpData(const double& value, const std::vector<double>&
}
// Else search in the vector
else {
//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, abs_value);
//If value is less than all values, use first interval
if (value < values.front()) {
retval.ind_[0] = 0;
retval.ind_[1] = 1;
}
//If value is greater than all values, use last interval
else if (value >= values.back()) {
retval.ind_[0] = nvalues-2;
retval.ind_[1] = nvalues-1;
}
else {
//Search internal intervals
for (int i=1; i<nvalues; ++i) {
if (values[i] >= value) {
retval.ind_[0] = i-1;
retval.ind_[1] = i;
break;
}
}
}
//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();
const double start = values[retval.ind_[0]];
const double end = values[retval.ind_[1]];
//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,
if (end > start) {
//FIXME: Possible source for floating point error here if value and floor are large,
//but very close to each other
retval.inv_dist_ = 1.0 / dist;
retval.factor_ = (abs_value-*floor_iter) * retval.inv_dist_;
retval.inv_dist_ = 1.0 / (end-start);
retval.factor_ = (value-start) * retval.inv_dist_;
}
else {
retval.inv_dist_ = 0.0;

30
tests/test_vfpproperties.cpp
View File

@ -67,25 +67,43 @@ BOOST_AUTO_TEST_SUITE( HelperTests )
BOOST_AUTO_TEST_CASE(findInterpData)
{
std::vector<double> values = {1, 5, 7, 9, 11, 15};
double interpolate = 6.87;
double extrapolate_left = -1.89;
double extrapolate_right = 32.1;
double exact = 9.0;
double interpolate = 6.0;
double extrapolate_left = -1.0;
double extrapolate_right = 19;
double first = 1;
double last = 15;
Opm::detail::InterpData eval0 = Opm::detail::findInterpData(exact, values);
Opm::detail::InterpData eval1 = Opm::detail::findInterpData(interpolate, values);
Opm::detail::InterpData eval2 = Opm::detail::findInterpData(extrapolate_left, values);
Opm::detail::InterpData eval3 = Opm::detail::findInterpData(extrapolate_right, values);
Opm::detail::InterpData eval4 = Opm::detail::findInterpData(first, values);
Opm::detail::InterpData eval5 = Opm::detail::findInterpData(last, values);
BOOST_CHECK_EQUAL(eval0.ind_[0], 2);
BOOST_CHECK_EQUAL(eval0.ind_[1], 3);
BOOST_CHECK_EQUAL(eval0.factor_, 1.0);
BOOST_CHECK_EQUAL(eval1.ind_[0], 1);
BOOST_CHECK_EQUAL(eval1.ind_[1], 2);
BOOST_CHECK_EQUAL(eval1.factor_, (interpolate-values[1]) / (values[2] - values[1]));
BOOST_CHECK_EQUAL(eval1.factor_, 0.5);
BOOST_CHECK_EQUAL(eval2.ind_[0], 0);
BOOST_CHECK_EQUAL(eval2.ind_[1], 1);
BOOST_CHECK_EQUAL(eval2.factor_, (extrapolate_left-values[0]) / (values[1] - values[0]));
BOOST_CHECK_EQUAL(eval2.factor_, -0.5);
BOOST_CHECK_EQUAL(eval3.ind_[0], 4);
BOOST_CHECK_EQUAL(eval3.ind_[1], 5);
BOOST_CHECK_EQUAL(eval3.factor_, (extrapolate_right-values[4]) / (values[5] - values[4]));
BOOST_CHECK_EQUAL(eval3.factor_, 2.0);
BOOST_CHECK_EQUAL(eval4.ind_[0], 0);
BOOST_CHECK_EQUAL(eval4.ind_[1], 1);
BOOST_CHECK_EQUAL(eval4.factor_, 0.0);
BOOST_CHECK_EQUAL(eval5.ind_[0], 4);
BOOST_CHECK_EQUAL(eval5.ind_[1], 5);
BOOST_CHECK_EQUAL(eval5.factor_, 1.0);
}
BOOST_AUTO_TEST_SUITE_END() // HelperTests