diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index 792be4c7f..8bb210d6c 100644
--- a/CMakeLists_files.cmake
+++ b/CMakeLists_files.cmake
@@ -61,6 +61,7 @@ list (APPEND MAIN_SOURCE_FILES
opm/simulators/wells/GroupState.cpp
opm/simulators/wells/ParallelWellInfo.cpp
opm/simulators/wells/TargetCalculator.cpp
+ opm/simulators/wells/VFPHelpers.cpp
opm/simulators/wells/VFPProdProperties.cpp
opm/simulators/wells/VFPInjProperties.cpp
opm/simulators/wells/WellGroupHelpers.cpp
diff --git a/opm/simulators/wells/VFPHelpers.cpp b/opm/simulators/wells/VFPHelpers.cpp
new file mode 100644
index 000000000..4c1f092b7
--- /dev/null
+++ b/opm/simulators/wells/VFPHelpers.cpp
@@ -0,0 +1,656 @@
+/*
+ Copyright 2015 SINTEF ICT, Applied Mathematics.
+
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see .
+*/
+
+#include
+#include
+
+#include
+
+#include
+
+#include
+#include
+
+#include
+#include
+#include
+
+namespace {
+
+/**
+ * Helper function that finds x for a given value of y for a line
+ * *NOTE ORDER OF ARGUMENTS*
+ */
+double findX(const double& x0,
+ const double& x1,
+ const double& y0,
+ const double& y1,
+ const double& y)
+{
+ const double dx = x1 - x0;
+ const double dy = y1 - y0;
+
+ /**
+ * y = y0 + (dy / dx) * (x - x0)
+ * => x = x0 + (y - y0) * (dx / dy)
+ *
+ * If dy is zero, use x1 as the value.
+ */
+
+ double x = 0.0;
+
+ if (dy != 0.0) {
+ x = x0 + (y-y0) * (dx/dy);
+ }
+ else {
+ x = x1;
+ }
+
+ return x;
+}
+
+/**
+ * Returns zero if input value is negative
+ */
+template
+static T chopNegativeValues(const T& value) {
+ return Opm::max(0.0, value);
+}
+
+}
+
+namespace Opm {
+namespace detail {
+
+InterpData findInterpData(const double& value_in, const std::vector& values)
+{
+ InterpData retval;
+
+ const int nvalues = values.size();
+
+ // chopping the value to be zero, which means we do not
+ // extrapolate the table towards nagative ranges
+ const double value = value_in < 0.? 0. : value_in;
+
+ //If we only have one value in our vector, return that
+ if (values.size() == 1) {
+ retval.ind_[0] = 0;
+ retval.ind_[1] = 0;
+ retval.inv_dist_ = 0.0;
+ retval.factor_ = 0.0;
+ }
+ // Else search in the vector
+ else {
+ //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= value) {
+ retval.ind_[0] = i-1;
+ retval.ind_[1] = i;
+ break;
+ }
+ }
+ }
+
+ const double start = values[retval.ind_[0]];
+ const double end = values[retval.ind_[1]];
+
+ //Find interpolation ratio
+ 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 / (end-start);
+ retval.factor_ = (value-start) * retval.inv_dist_;
+ }
+ else {
+ retval.inv_dist_ = 0.0;
+ retval.factor_ = 0.0;
+ }
+ }
+
+ // Disallow extrapolation with higher factor than 3.0.
+ // The factor 3.0 has been chosen because it works well
+ // with certain testcases, and may not be optimal.
+ if (retval.factor_ > 3.0) {
+ retval.factor_ = 3.0;
+ }
+
+ return retval;
+}
+
+VFPEvaluation operator+(VFPEvaluation lhs, const VFPEvaluation& rhs)
+{
+ lhs.value += rhs.value;
+ lhs.dthp += rhs.dthp;
+ lhs.dwfr += rhs.dwfr;
+ lhs.dgfr += rhs.dgfr;
+ lhs.dalq += rhs.dalq;
+ lhs.dflo += rhs.dflo;
+ return lhs;
+}
+
+VFPEvaluation operator-(VFPEvaluation lhs, const VFPEvaluation& rhs)
+{
+ lhs.value -= rhs.value;
+ lhs.dthp -= rhs.dthp;
+ lhs.dwfr -= rhs.dwfr;
+ lhs.dgfr -= rhs.dgfr;
+ lhs.dalq -= rhs.dalq;
+ lhs.dflo -= rhs.dflo;
+ return lhs;
+}
+
+VFPEvaluation operator*(double lhs, const VFPEvaluation& rhs)
+{
+ VFPEvaluation retval;
+ retval.value = rhs.value * lhs;
+ retval.dthp = rhs.dthp * lhs;
+ retval.dwfr = rhs.dwfr * lhs;
+ retval.dgfr = rhs.dgfr * lhs;
+ retval.dalq = rhs.dalq * lhs;
+ retval.dflo = rhs.dflo * lhs;
+ return retval;
+}
+
+VFPEvaluation interpolate(const VFPProdTable& table,
+ const InterpData& flo_i,
+ const InterpData& thp_i,
+ const InterpData& wfr_i,
+ const InterpData& gfr_i,
+ const InterpData& alq_i)
+{
+ //Values and derivatives in a 5D hypercube
+ VFPEvaluation nn[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
+ //interpolation below.
+ //The following ladder of for loops will presumably be unrolled by a reasonable compiler.
+ 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) {
+ for (int f=0; f<=1; ++f) {
+ //Shorthands for indexing
+ 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].value = table(ti,wi,gi,ai,fi);
+ }
+ }
+ }
+ }
+ }
+
+ //Calculate derivatives
+ //Note that the derivative of the two end points of a line aligned with the
+ //"axis of the derivative" are equal
+ 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) {
+ nn[0][i][j][k][l].dthp = (nn[1][i][j][k][l].value - nn[0][i][j][k][l].value) * thp_i.inv_dist_;
+ nn[i][0][j][k][l].dwfr = (nn[i][1][j][k][l].value - nn[i][0][j][k][l].value) * wfr_i.inv_dist_;
+ nn[i][j][0][k][l].dgfr = (nn[i][j][1][k][l].value - nn[i][j][0][k][l].value) * gfr_i.inv_dist_;
+ nn[i][j][k][0][l].dalq = (nn[i][j][k][1][l].value - nn[i][j][k][0][l].value) * alq_i.inv_dist_;
+ nn[i][j][k][l][0].dflo = (nn[i][j][k][l][1].value - nn[i][j][k][l][0].value) * flo_i.inv_dist_;
+
+ nn[1][i][j][k][l].dthp = nn[0][i][j][k][l].dthp;
+ nn[i][1][j][k][l].dwfr = nn[i][0][j][k][l].dwfr;
+ nn[i][j][1][k][l].dgfr = nn[i][j][0][k][l].dgfr;
+ nn[i][j][k][1][l].dalq = nn[i][j][k][0][l].dalq;
+ nn[i][j][k][l][1].dflo = nn[i][j][k][l][0].dflo;
+ }
+ }
+ }
+ }
+
+ double t1, t2; //interpolation variables, so that t1 = (1-t) and t2 = 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.
+ t2 = flo_i.factor_;
+ t1 = (1.0-t2);
+ 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] = t1*nn[t][w][g][a][0] + t2*nn[t][w][g][a][1];
+ }
+ }
+ }
+ }
+
+ t2 = alq_i.factor_;
+ t1 = (1.0-t2);
+ 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] = t1*nn[t][w][g][0][0] + t2*nn[t][w][g][1][0];
+ }
+ }
+ }
+
+ t2 = gfr_i.factor_;
+ t1 = (1.0-t2);
+ for (int t=0; t<=1; ++t) {
+ for (int w=0; w<=1; ++w) {
+ nn[t][w][0][0][0] = t1*nn[t][w][0][0][0] + t2*nn[t][w][1][0][0];
+ }
+ }
+
+ t2 = wfr_i.factor_;
+ t1 = (1.0-t2);
+ for (int t=0; t<=1; ++t) {
+ nn[t][0][0][0][0] = t1*nn[t][0][0][0][0] + t2*nn[t][1][0][0][0];
+ }
+
+ t2 = thp_i.factor_;
+ t1 = (1.0-t2);
+ nn[0][0][0][0][0] = t1*nn[0][0][0][0][0] + t2*nn[1][0][0][0][0];
+
+ return nn[0][0][0][0][0];
+}
+
+VFPEvaluation interpolate(const VFPInjTable& table,
+ const InterpData& flo_i,
+ const InterpData& thp_i)
+{
+ //Values and derivatives in a 2D plane
+ VFPEvaluation nn[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.
+ for (int t=0; t<=1; ++t) {
+ for (int f=0; f<=1; ++f) {
+ //Shorthands for indexing
+ const int ti = thp_i.ind_[t];
+ const int fi = flo_i.ind_[f];
+
+ //Copy element
+ nn[t][f].value = table(ti,fi);
+ }
+ }
+
+ //Calculate derivatives
+ //Note that the derivative of the two end points of a line aligned with the
+ //"axis of the derivative" are equal
+ for (int i=0; i<=1; ++i) {
+ nn[0][i].dthp = (nn[1][i].value - nn[0][i].value) * thp_i.inv_dist_;
+ nn[i][0].dwfr = -1e100;
+ nn[i][0].dgfr = -1e100;
+ nn[i][0].dalq = -1e100;
+ nn[i][0].dflo = (nn[i][1].value - nn[i][0].value) * flo_i.inv_dist_;
+
+ nn[1][i].dthp = nn[0][i].dthp;
+ nn[i][1].dwfr = nn[i][0].dwfr;
+ nn[i][1].dgfr = nn[i][0].dgfr;
+ nn[i][1].dalq = nn[i][0].dalq;
+ nn[i][1].dflo = nn[i][0].dflo;
+ }
+
+ double t1, t2; //interpolation variables, so that t1 = (1-t) and t2 = t.
+
+ // Go from 2D to 1D
+ t2 = flo_i.factor_;
+ t1 = (1.0-t2);
+ nn[0][0] = t1*nn[0][0] + t2*nn[0][1];
+ nn[1][0] = t1*nn[1][0] + t2*nn[1][1];
+
+ // Go from line to point on line
+ t2 = thp_i.factor_;
+ t1 = (1.0-t2);
+ nn[0][0] = t1*nn[0][0] + t2*nn[1][0];
+
+ return nn[0][0];
+}
+
+VFPEvaluation bhp(const VFPProdTable& table,
+ const double& aqua,
+ const double& liquid,
+ const double& vapour,
+ const double& thp,
+ const double& alq)
+{
+ //Find interpolation variables
+ double flo = detail::getFlo(table, aqua, liquid, vapour);
+ double wfr = detail::getWFR(table, aqua, liquid, vapour);
+ double gfr = detail::getGFR(table, aqua, liquid, vapour);
+
+ //First, find the values to interpolate between
+ //Recall that flo is negative in Opm, so switch sign.
+ auto flo_i = detail::findInterpData(-flo, table.getFloAxis());
+ auto thp_i = detail::findInterpData( thp, table.getTHPAxis());
+ auto wfr_i = detail::findInterpData( wfr, table.getWFRAxis());
+ auto gfr_i = detail::findInterpData( gfr, table.getGFRAxis());
+ auto alq_i = detail::findInterpData( alq, table.getALQAxis());
+
+ detail::VFPEvaluation retval = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i);
+
+ return retval;
+}
+
+VFPEvaluation bhp(const VFPInjTable& table,
+ const double& aqua,
+ const double& liquid,
+ const double& vapour,
+ const double& thp)
+{
+ //Find interpolation variables
+ double flo = detail::getFlo(table, aqua, liquid, vapour);
+
+ //First, find the values to interpolate between
+ auto flo_i = detail::findInterpData(flo, table.getFloAxis());
+ auto thp_i = detail::findInterpData(thp, table.getTHPAxis());
+
+ //Then perform the interpolation itself
+ detail::VFPEvaluation retval = detail::interpolate(table, flo_i, thp_i);
+
+ return retval;
+}
+
+double findTHP(const std::vector& bhp_array,
+ const std::vector& thp_array,
+ double bhp)
+{
+ int nthp = thp_array.size();
+
+ double thp = -1e100;
+
+ //Check that our thp axis is sorted
+ assert(std::is_sorted(thp_array.begin(), thp_array.end()));
+
+ /**
+ * Our *interpolated* bhp_array will be montonic increasing for increasing
+ * THP if our input BHP values are monotonic increasing for increasing
+ * THP values. However, if we have to *extrapolate* along any of the other
+ * axes, this guarantee holds no more, and bhp_array may be "random"
+ */
+ if (std::is_sorted(bhp_array.begin(), bhp_array.end())) {
+ //Target bhp less than all values in array, extrapolate
+ if (bhp <= bhp_array[0]) {
+ //TODO: LOG extrapolation
+ const double& x0 = thp_array[0];
+ const double& x1 = thp_array[1];
+ const double& y0 = bhp_array[0];
+ const double& y1 = bhp_array[1];
+ thp = findX(x0, x1, y0, y1, bhp);
+ }
+ //Target bhp greater than all values in array, extrapolate
+ else if (bhp > bhp_array[nthp-1]) {
+ //TODO: LOG extrapolation
+ const double& x0 = thp_array[nthp-2];
+ const double& x1 = thp_array[nthp-1];
+ const double& y0 = bhp_array[nthp-2];
+ const double& y1 = bhp_array[nthp-1];
+ thp = findX(x0, x1, y0, y1, bhp);
+ }
+ //Target bhp within table ranges, interpolate
+ else {
+ //Loop over the values and find min(bhp_array(thp)) == bhp
+ //so that we maximize the rate.
+
+ //Find i so that bhp_array[i-1] <= bhp <= bhp_array[i];
+ //Assuming a small number of values in bhp_array, this should be quite
+ //efficient. Other strategies might be bisection, etc.
+ int i=0;
+ bool found = false;
+ for (; i(found); //Silence compiler warning
+
+ const double& x0 = thp_array[i ];
+ const double& x1 = thp_array[i+1];
+ const double& y0 = bhp_array[i ];
+ const double& y1 = bhp_array[i+1];
+ thp = findX(x0, x1, y0, y1, bhp);
+ }
+ }
+ //bhp_array not sorted, raw search.
+ else {
+ //Find i so that bhp_array[i-1] <= bhp <= bhp_array[i];
+ //Since the BHP values might not be sorted, first search within
+ //our interpolation values, and then try to extrapolate.
+ int i=0;
+ bool found = false;
+ for (; i bhp_array[nthp-1]) {
+ //TODO: LOG extrapolation
+ const double& x0 = thp_array[nthp-2];
+ const double& x1 = thp_array[nthp-1];
+ const double& y0 = bhp_array[nthp-2];
+ const double& y1 = bhp_array[nthp-1];
+ thp = findX(x0, x1, y0, y1, bhp);
+ }
+ else {
+ OPM_THROW(std::logic_error, "Programmer error: Unable to find THP in THP array");
+ }
+ }
+
+ return thp;
+}
+
+template
+T getFlo(const VFPProdTable& table,
+ const T& aqua,
+ const T& liquid,
+ const T& vapour)
+{
+ auto type = table.getFloType();
+ switch (type) {
+ case VFPProdTable::FLO_TYPE::FLO_OIL:
+ //Oil = liquid phase
+ return liquid;
+ case VFPProdTable::FLO_TYPE::FLO_LIQ:
+ //Liquid = aqua + liquid phases
+ return aqua + liquid;
+ case VFPProdTable::FLO_TYPE::FLO_GAS:
+ //Gas = vapor phase
+ return vapour;
+ default:
+ throw std::logic_error("Invalid FLO_TYPE");
+ }
+}
+
+template
+T getFlo(const VFPInjTable& table,
+ const T& aqua,
+ const T& liquid,
+ const T& vapour)
+{
+ auto type = table.getFloType();
+ switch (type) {
+ case VFPInjTable::FLO_TYPE::FLO_OIL:
+ //Oil = liquid phase
+ return liquid;
+ case VFPInjTable::FLO_TYPE::FLO_WAT:
+ //Liquid = aqua phase
+ return aqua;
+ case VFPInjTable::FLO_TYPE::FLO_GAS:
+ //Gas = vapor phase
+ return vapour;
+ default:
+ throw std::logic_error("Invalid FLO_TYPE");
+ }
+}
+
+static constexpr double threshold = 1e-12;
+
+template
+T getWFR(const VFPProdTable& table,
+ const T& aqua,
+ const T& liquid,
+ const T& vapour)
+{
+ auto type = table.getWFRType();
+ switch(type) {
+ case VFPProdTable::WFR_TYPE::WFR_WOR: {
+ //Water-oil ratio = water / oil
+ return chopNegativeValues(-aqua) / max(threshold, chopNegativeValues(-liquid));
+ }
+ case VFPProdTable::WFR_TYPE::WFR_WCT:
+ //Water cut = water / (water + oil)
+ return chopNegativeValues(-aqua) / max(threshold, chopNegativeValues(-aqua - liquid));
+ case VFPProdTable::WFR_TYPE::WFR_WGR:
+ //Water-gas ratio = water / gas
+ return chopNegativeValues(-aqua) / max(threshold, chopNegativeValues(-vapour));
+ default:
+ throw std::logic_error("Invalid WFR_TYPE");
+ }
+}
+
+template
+T getGFR(const VFPProdTable& table,
+ const T& aqua,
+ const T& liquid,
+ const T& vapour)
+{
+ auto type = table.getGFRType();
+ switch(type) {
+ case VFPProdTable::GFR_TYPE::GFR_GOR:
+ // Gas-oil ratio = gas / oil
+ return chopNegativeValues(-vapour) / max(threshold, chopNegativeValues(-liquid));
+ case VFPProdTable::GFR_TYPE::GFR_GLR:
+ // Gas-liquid ratio = gas / (oil + water)
+ return chopNegativeValues(-vapour) / max(threshold, chopNegativeValues(-liquid - aqua));
+ case VFPProdTable::GFR_TYPE::GFR_OGR:
+ // Oil-gas ratio = oil / gas
+ return chopNegativeValues(-liquid) / max(threshold, chopNegativeValues(-vapour));
+ default:
+ throw std::logic_error("Invalid GFR_TYPE");
+ }
+}
+
+template
+const T& getTable(const std::map> tables, int table_id)
+{
+ auto entry = tables.find(table_id);
+ if (entry == tables.end()) {
+ OPM_THROW(std::invalid_argument, "Nonexistent VFP table " << table_id << " referenced.");
+ }
+ else {
+ return entry->second.get();
+ }
+}
+
+template <>
+VFPProdTable::FLO_TYPE getType(const VFPProdTable& table)
+{
+ return table.getFloType();
+}
+
+template <>
+VFPProdTable::WFR_TYPE getType(const VFPProdTable& table)
+{
+ return table.getWFRType();
+}
+
+template <>
+VFPProdTable::GFR_TYPE getType(const VFPProdTable& table)
+{
+ return table.getGFRType();
+}
+
+/**
+ * Returns the type variable for FLO for injection tables
+ */
+template <>
+VFPInjTable::FLO_TYPE getType(const VFPInjTable& table)
+{
+ return table.getFloType();
+}
+
+template const VFPInjTable& getTable(const std::map>, int);
+template const VFPProdTable& getTable(const std::map>, int);
+
+#define INSTANCE(...) \
+ template __VA_ARGS__ getFlo(const VFPInjTable&, const __VA_ARGS__&, const __VA_ARGS__&, const __VA_ARGS__&); \
+ template __VA_ARGS__ getFlo(const VFPProdTable&, const __VA_ARGS__&, const __VA_ARGS__&, const __VA_ARGS__&); \
+ template __VA_ARGS__ getGFR(const VFPProdTable&, const __VA_ARGS__&, const __VA_ARGS__&, const __VA_ARGS__&); \
+ template __VA_ARGS__ getWFR(const VFPProdTable&, const __VA_ARGS__&, const __VA_ARGS__&, const __VA_ARGS__&);
+
+INSTANCE(double)
+INSTANCE(DenseAd::Evaluation)
+INSTANCE(DenseAd::Evaluation)
+INSTANCE(DenseAd::Evaluation)
+INSTANCE(DenseAd::Evaluation)
+INSTANCE(DenseAd::Evaluation)
+INSTANCE(DenseAd::Evaluation)
+INSTANCE(DenseAd::Evaluation)
+INSTANCE(DenseAd::Evaluation)
+INSTANCE(DenseAd::Evaluation)
+INSTANCE(DenseAd::Evaluation)
+INSTANCE(DenseAd::Evaluation)
+INSTANCE(DenseAd::Evaluation)
+INSTANCE(DenseAd::Evaluation)
+
+} // namespace detail
+} // namespace Opm
diff --git a/opm/simulators/wells/VFPHelpers.hpp b/opm/simulators/wells/VFPHelpers.hpp
index 4812f5f6e..a1c17cf90 100644
--- a/opm/simulators/wells/VFPHelpers.hpp
+++ b/opm/simulators/wells/VFPHelpers.hpp
@@ -21,56 +21,31 @@
#ifndef OPM_AUTODIFF_VFPHELPERS_HPP_
#define OPM_AUTODIFF_VFPHELPERS_HPP_
-#include
-
-#include
-#include
-#include
-#include
-#include
-#include
-#include
+#include
+#include
+#include