Added linear extrapolation

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
+            //Find last element smaller than range
-            auto floor_iter = ceil_iter;
+            auto floor_iter = ceil_iter-1;
-            if (*floor_iter == value) {
-                // floor_iter == ceil_iter == value
-            }
-            else if (floor_iter > values.begin()) {
-                // floor_iter <= value <= ceil_iter
-                --floor_iter;
-            }
 
-            //Now set these in the retval struct
+            //Find the indices
-            retval.ind_[0] = floor_iter - values.begin();
+            int a = floor_iter - values.begin();
-            retval.ind_[1] = ceil_iter - values.begin();
+            int b = ceil_iter - values.begin();
+            int max_size = static_cast<int>(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];
+                                const int ti = thp_i.ind_[t];
-                                int wi = wfr_i.ind_[w];
+                                const int wi = wfr_i.ind_[w];
-                                int gi = gfr_i.ind_[g];
+                                const int gi = gfr_i.ind_[g];
-                                int ai = alq_i.ind_[a];
+                                const int ai = alq_i.ind_[a];
-                                int fi = flo_i.ind_[f];
+                                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_;

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<double>& thp_axis{0.0, 1.0};
+    const std::vector<double>& wfr_axis{0.0, 0.5, 1.0};
+    const std::vector<double>& gfr_axis{0.0, 0.25, 0.5, 0.75, 1};
+    const std::vector<double>& alq_axis{0.0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1};
+    const std::vector<double>& 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<int>(thp_axis.size());
+    int ny = static_cast<int>(wfr_axis.size());
+    int nz = static_cast<int>(gfr_axis.size());
+    int nu = static_cast<int>(alq_axis.size());
+    int nv = static_cast<int>(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; ++i) {
+        double x = i / static_cast<double>(nx-1);
+        for (int j=0; j<ny; ++j) {
+            double y = j / static_cast<double>(ny-1);
+            for (int k=0; k<nz; ++k) {
+                double z = k / static_cast<double>(nz-1);
+                for (int l=0; l<nu; ++l) {
+                    double u = l / static_cast<double>(nu-1);
+                    for (int m=0; m<nv; ++m) {
+                        double v = m / static_cast<double>(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<double>(n);
+        for (int j=-5; j<=n+5; ++j) {
+            const double y = j / static_cast<double>(n);
+            for (int k=-5; k<=n+5; ++k) {
+                const double z = k / static_cast<double>(n);
+                for (int l=-5; l<=n+5; ++l) {
+                    const double u = l / static_cast<double>(n);
+                    for (int m=-5; m<=n+5; ++m) {
+                        const double v = m / static_cast<double>(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