diff --git a/opm/autodiff/VFPProperties.hpp b/opm/autodiff/VFPProperties.hpp index f6c468903..7cec4de92 100644 --- a/opm/autodiff/VFPProperties.hpp +++ b/opm/autodiff/VFPProperties.hpp @@ -284,26 +284,26 @@ namespace Opm { 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(), values.end()-1, value); + auto ceil_iter = std::lower_bound(values.begin()+1, values.end()-1, value); - //Find last element smaller than or equal to range - auto floor_iter = ceil_iter; - if (*floor_iter == value) { - // floor_iter == ceil_iter == value - } - else if (floor_iter > values.begin()) { - // floor_iter <= value <= ceil_iter - --floor_iter; - } + //Find last element smaller than range + auto floor_iter = ceil_iter-1; - //Now set these in the retval struct - retval.ind_[0] = floor_iter - values.begin(); - retval.ind_[1] = ceil_iter - values.begin(); + //Find the indices + int a = floor_iter - values.begin(); + int b = ceil_iter - values.begin(); + int max_size = static_cast(values.size())-1; + + //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) { + assert(dist >= 0.0); + 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; @@ -315,14 +315,15 @@ namespace Opm { return retval; } +#pragma GCC push_options +#pragma GCC optimize ("unroll-loops") double interpolate(const InterpData& flo_i, const InterpData& thp_i, const InterpData& wfr_i, const InterpData& gfr_i, const InterpData& alq_i) { - //extents shape({{2, 2, 2, 2, 2}}); - //array_type nn(shape); double nn[2][2][2][2][2]; //Pick out nearest neighbors (nn) to our evaluation point //The following ladder of for loops will presumably be unrolled by a reasonable compiler. + //If needed, this can be manually unrolled //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 //interpolation below. @@ -332,11 +333,11 @@ namespace Opm { for (int a=0; a<=1; ++a) { for (int f=0; f<=1; ++f) { //Shorthands for indexing - int ti = thp_i.ind_[t]; - int wi = wfr_i.ind_[w]; - int gi = gfr_i.ind_[g]; - int ai = alq_i.ind_[a]; - int fi = flo_i.ind_[f]; + const int ti = thp_i.ind_[t]; + const int wi = wfr_i.ind_[w]; + const int gi = gfr_i.ind_[g]; + const int ai = alq_i.ind_[a]; + const int fi = flo_i.ind_[f]; //Copy element nn[t][w][g][a][f] = data_[ti][wi][gi][ai][fi]; @@ -385,6 +386,7 @@ namespace Opm { tf = thp_i.factor_; return (1.0-tf)*nn[0][0][0][0][0] + tf*nn[1][0][0][0][0]; } +#pragma GCC pop_options //unroll loops //"Header" variables int table_num_; diff --git a/tests/test_vfpproperties.cpp b/tests/test_vfpproperties.cpp index 257715a67..3ebfe6ac1 100644 --- a/tests/test_vfpproperties.cpp +++ b/tests/test_vfpproperties.cpp @@ -288,6 +288,101 @@ BOOST_AUTO_TEST_CASE(InterpolatePlane) +/** + * Test that we can generate some dummy one-data, and + * interpolate something the analytic solution + */ +BOOST_AUTO_TEST_CASE(ExtrapolatePlane) +{ + //All of our axes go from 0 to 1 + const std::vector& thp_axis{0.0, 1.0}; + const std::vector& wfr_axis{0.0, 0.5, 1.0}; + const std::vector& gfr_axis{0.0, 0.25, 0.5, 0.75, 1}; + const std::vector& alq_axis{0.0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1}; + const std::vector& flo_axis{0.0, 0.0625, 0.125, 0.1875, 0.25, 0.3125, 0.375, 0.4375, + 0.5, 0.5625, 0.625, 0.6875, 0.75, 0.8125, 0.875, 0.9375, 1}; + + int nx = static_cast(thp_axis.size()); + int ny = static_cast(wfr_axis.size()); + int nz = static_cast(gfr_axis.size()); + int nu = static_cast(alq_axis.size()); + int nv = static_cast(flo_axis.size()); + + //Allocate data and fill with trivial data (zeros) + Opm::VFPProperties::extents size = {{ nx, ny, nz, nu, nv }}; + Opm::VFPProperties::array_type data(size); + for (int i=0; i(nx-1); + for (int j=0; j(ny-1); + for (int k=0; k(nz-1); + for (int l=0; l(nu-1); + for (int m=0; m(nv-1); + + data[i][j][k][l][m] = x + 2*y + 3*z + 4*u + 5*v; + } + } + } + } + } + + //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); + + //Check linear extrapolation (i.e., using values of x, y, etc. outside our interpolant domain) + double sum = 0.0; + double reference_sum = 0.0; + double sad = 0.0; // Sum absolute difference + double max_d = 0.0; // Maximum difference + int n=1; + for (int i=-5; i<=n+5; ++i) { + const double x = i / static_cast(n); + for (int j=-5; j<=n+5; ++j) { + const double y = j / static_cast(n); + for (int k=-5; k<=n+5; ++k) { + const double z = k / static_cast(n); + for (int l=-5; l<=n+5; ++l) { + const double u = l / static_cast(n); + for (int m=-5; m<=n+5; ++m) { + const double v = m / static_cast(n); + double reference = x + 2*y + 3*z + 4*u + 5*v; + reference_sum += reference; + + //Note order of arguments! + double value = table.bhp(v, x, y, z, u); + sum += value; + + double abs_diff = std::abs(value - reference); + + sad += std::abs(abs_diff); + max_d = std::max(max_d, abs_diff); + } + } + } + } + } + + BOOST_CHECK_CLOSE(sum, reference_sum, 0.0001); + BOOST_CHECK_SMALL(max_d, 0.0005); + BOOST_CHECK_SMALL(sad, 0.1); +} + + + /** * Tests that we can actually parse some input data, and interpolate within that space