Minor fixes for PR

opm/autodiff/BlackoilModelBase_impl.hpp

@@ -1972,7 +1972,6 @@ namespace detail {
                             vapour = wr[w*np + pu.phase_pos[ Gas ] ];
                         }
 
-                        auto wc = wells().ctrls[w];
                         double alq = well_controls_iget_alq(wc, ctrl_index);
                         int table_id = well_controls_iget_vfp(wc, ctrl_index);
 

opm/autodiff/VFPHelpers.hpp

@@ -208,7 +208,7 @@ struct InterpData {
 inline InterpData findInterpData(const double& value, const std::vector<double>& values) {
     InterpData retval;
 
-    const double abs_value = std::abs(value);
+    const double abs_value = value;//std::abs(value);
 
     //If we only have one value in our vector, return that
     if (values.size() == 1) {
@@ -515,7 +515,8 @@ inline VFPEvaluation bhp(const VFPProdTable* table,
     double gfr = detail::getGFR(aqua, liquid, vapour, table->getGFRType());
 
     //First, find the values to interpolate between
-    auto flo_i = detail::findInterpData(flo, table->getFloAxis());
+    //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());
@@ -786,10 +787,10 @@ ADB gather_vars(const std::vector<const TABLE*>& well_tables,
         const ADB& values = map.find(key)->second;
 
         //Get indices to all elements that should use this ADB
-        const std::vector<int>& elems = value;
+        const std::vector<int>& current = value;
 
         //Add these elements to retval
-        retval = retval + superset(subset(values, elems), elems, values.size());
+        retval = retval + superset(subset(values, current), current, values.size());
     }
 
     return retval;

opm/autodiff/VFPInjProperties.cpp

@@ -115,7 +115,7 @@ VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector<int>& table_id,
     ADB::V dflo = ADB::V::Zero(nw);
 
     //Get the table for each well
-    std::vector<const VFPInjTable*> well_tables(nw, NULL);
+    std::vector<const VFPInjTable*> well_tables(nw, nullptr);
     for (int i=0; i<nw; ++i) {
         if (table_id[i] > 0) {
             well_tables[i] = detail::getTable(m_tables, table_id[i]);
@@ -128,7 +128,7 @@ VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector<int>& table_id,
     //Compute the BHP for each well independently
     for (int i=0; i<nw; ++i) {
         const VFPInjTable* table = well_tables[i];
-        if (table != NULL) {
+        if (table != nullptr) {
             //First, find the values to interpolate between
             auto flo_i = detail::findInterpData(flo.value()[i], table->getFloAxis());
             auto thp_i = detail::findInterpData(thp.value()[i], table->getTHPAxis());

opm/autodiff/VFPInjProperties.hpp

@@ -46,14 +46,14 @@ public:
      * Takes *no* ownership of data.
      * @param inj_table A *single* VFPINJ table
      */
-    VFPInjProperties(const VFPInjTable* inj_table);
+    explicit VFPInjProperties(const VFPInjTable* inj_table);
 
     /**
      * Constructor
      * Takes *no* ownership of data.
      * @param inj_tables A map of different VFPINJ tables.
      */
-    VFPInjProperties(const std::map<int, VFPInjTable>& inj_tables);
+    explicit VFPInjProperties(const std::map<int, VFPInjTable>& inj_tables);
 
     /**
      * Linear interpolation of bhp as function of the input parameters.
@@ -135,7 +135,7 @@ public:
     /**
      * Returns true if no vfp tables are in the current map
      */
-    inline const bool empty() const {
+    bool empty() const {
         return m_tables.empty();
     }
 

opm/autodiff/VFPProdProperties.cpp

@@ -105,7 +105,7 @@ VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector<int>& table_id,
     ADB::V dflo = ADB::V::Zero(nw);
 
     //Get the table for each well
-    std::vector<const VFPProdTable*> well_tables(nw, NULL);
+    std::vector<const VFPProdTable*> well_tables(nw, nullptr);
     for (int i=0; i<nw; ++i) {
         if (table_id[i] >= 0) {
             well_tables[i] = detail::getTable(m_tables, table_id[i]);
@@ -120,9 +120,10 @@ VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector<int>& table_id,
     //Compute the BHP for each well independently
     for (int i=0; i<nw; ++i) {
         const VFPProdTable* table = well_tables[i];
-        if (table != NULL) {
+        if (table != nullptr) {
             //First, find the values to interpolate between
-            auto flo_i = detail::findInterpData(flo.value()[i], table->getFloAxis());
+            //Value of FLO is negative in OPM for producers, but positive in VFP table
+            auto flo_i = detail::findInterpData(-flo.value()[i], table->getFloAxis());
             auto thp_i = detail::findInterpData( thp.value()[i], table->getTHPAxis());
             auto wfr_i = detail::findInterpData( wfr.value()[i], table->getWFRAxis());
             auto gfr_i = detail::findInterpData( gfr.value()[i], table->getGFRAxis());
@@ -130,38 +131,6 @@ VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector<int>& table_id,
 
             detail::VFPEvaluation bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i);
 
-            /*
-            static const int N=40;
-            std::cout << "bhp=" << bhp_val.value << ";" << std::endl;
-            std::cout << "flo=" << flo.value()[i] << ";" << std::endl;
-            std::cout << "thp=" << thp.value()[i] << ";" << std::endl;
-            std::cout << "wfr=" << wfr.value()[i] << ";" << std::endl;
-            std::cout << "gfr=" << gfr.value()[i] << ";" << std::endl;
-            std::cout << "alq=" << alq.value()[i] << ";" << std::endl;
-            std::cout << "bhp_vfp=[" << std::endl;
-            for (int j=0; j<N; ++j) {
-                const double start = table->getFloAxis().front();
-                const double end = table->getFloAxis().back();
-                const double dist = end - start;
-                double flo_d = start + (j/static_cast<double>(N-1)) * dist;
-                auto flo_i = detail::findInterpData(flo_d, table->getFloAxis());
-                detail::VFPEvaluation bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i);
-                std::cout << bhp_val.value << ",";
-            }
-            std::cout << "];" << std::endl;
-
-            std::cout << "flo_vfp=[" << std::endl;
-            for (int j=0; j<N; ++j) {
-                const double start = table->getFloAxis().front();
-                const double end = table->getFloAxis().back();
-                const double dist = end - start;
-                double flo_d = start + (j/static_cast<double>(N-1)) * dist;
-                std::cout << flo_d << ",";
-            }
-            std::cout << "];" << std::endl;
-            */
-
-
             value[i] = bhp_val.value;
             dthp[i] = bhp_val.dthp;
             dwfr[i] = bhp_val.dwfr;
@@ -247,8 +216,9 @@ double VFPProdProperties::thp(int table_id,
      * Find the function bhp_array(thp) by creating a 1D view of the data
      * by interpolating for every value of thp. This might be somewhat
      * expensive, but let us assome that nthp is small
+     * Recall that flo is negative in Opm, so switch the sign
      */
-    auto flo_i = detail::findInterpData(flo, table->getFloAxis());
+    auto flo_i = detail::findInterpData(-flo, table->getFloAxis());
     auto wfr_i = detail::findInterpData( wfr, table->getWFRAxis());
     auto gfr_i = detail::findInterpData( gfr, table->getGFRAxis());
     auto alq_i = detail::findInterpData( alq, table->getALQAxis());

opm/autodiff/VFPProdProperties.hpp

@@ -50,14 +50,14 @@ public:
      * Takes *no* ownership of data.
      * @param prod_table A *single* VFPPROD table
      */
-    VFPProdProperties(const VFPProdTable* prod_table);
+    explicit VFPProdProperties(const VFPProdTable* prod_table);
 
     /**
      * Constructor
      * Takes *no* ownership of data.
      * @param prod_tables A map of different VFPPROD tables.
      */
-    VFPProdProperties(const std::map<int, VFPProdTable>& prod_tables);
+    explicit VFPProdProperties(const std::map<int, VFPProdTable>& prod_tables);
 
     /**
      * Linear interpolation of bhp as function of the input parameters.
@@ -146,7 +146,7 @@ public:
     /**
      * Returns true if no vfp tables are in the current map
      */
-    inline const bool empty() const {
+    bool empty() const {
         return m_tables.empty();
     }
 

opm/autodiff/VFPProperties.hpp

@@ -44,7 +44,7 @@ public:
      * @param inj_table  A *single* VFPINJ table or NULL (no table)
      * @param prod_table A *single* VFPPROD table or NULL (no table)
      */
-    VFPProperties(const VFPInjTable* inj_table, const VFPProdTable* prod_table);
+    explicit VFPProperties(const VFPInjTable* inj_table, const VFPProdTable* prod_table);
 
     /**
      * Constructor
@@ -58,14 +58,14 @@ public:
     /**
      * Returns the VFP properties for injection wells
      */
-    inline const VFPInjProperties* getInj() const {
+    const VFPInjProperties* getInj() const {
         return m_inj.get();
     }
 
     /**
      * Returns the VFP properties for production wells
      */
-    inline const VFPProdProperties* getProd() const {
+    const VFPProdProperties* getProd() const {
         return m_prod.get();
     }
 

opm/autodiff/WellDensitySegmented.cpp

@@ -147,7 +147,6 @@ Opm::WellDensitySegmented::computeConnectionPressureDelta(const Wells& wells,
                                                           const std::vector<double>& z_perf,
                                                           const std::vector<double>& dens_perf,
                                                           const double gravity) {
-    const int np = wells.number_of_phases;
     const int nw = wells.number_of_wells;
     const int nperf = wells.well_connpos[nw];
 
@@ -177,7 +176,6 @@ Opm::WellDensitySegmented::computeConnectionPressureDelta(const Wells& wells,
             const double z_above = perf == wells.well_connpos[w] ? wells.depth_ref[w] : z_perf[perf - 1];
             const double dz = z_perf[perf] - z_above;
             dp_perf[perf] = dz * dens_perf[perf] * gravity;
-            //dens[wells.well_connpos[w]]
         }
     }
 

tests/test_vfpproperties.cpp

@@ -31,9 +31,12 @@
 #include <map>
 #include <sstream>
 #include <limits>
+#include <vector>
 
+#include <opm/core/utility/platform_dependent/disable_warnings.h>
 #include <boost/test/unit_test.hpp>
 #include <boost/filesystem.hpp>
+#include <opm/core/utility/platform_dependent/reenable_warnings.h>
 
 #include <opm/core/wells.h>
 #include <opm/parser/eclipse/Parser/Parser.hpp>
@@ -58,6 +61,43 @@ const double sad_tol = 1.0e-8;
 
 
 
+
+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;
+
+    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);
+
+    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(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(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_AUTO_TEST_SUITE_END() // HelperTests
+
+
+
+
+
+
+
+
+
 struct ConversionFixture {
     typedef Opm::VFPProdProperties::ADB ADB;
 
@@ -391,9 +431,9 @@ BOOST_AUTO_TEST_CASE(GetTable)
     add_well(INJECTOR, 100, 1, NULL, cells, NULL, NULL, wells.get());
 
     //Create interpolation points
-    double aqua_d   = 0.15;
+    double aqua_d   = -0.15;
-    double liquid_d = 0.25;
+    double liquid_d = -0.25;
-    double vapour_d = 0.35;
+    double vapour_d = -0.35;
     double thp_d    =  0.45;
     double alq_d    =  0.55;
 
@@ -523,13 +563,13 @@ BOOST_AUTO_TEST_CASE(InterpolatePlane)
     for (int i=0; i<=n; ++i) {
         const double thp = i / static_cast<double>(n);
         for (int j=1; j<=n; ++j) {
-            const double aqua = j / static_cast<double>(n);
+            const double aqua = -j / static_cast<double>(n);
             for (int k=1; k<=n; ++k) {
-                const double vapour = k / static_cast<double>(n);
+                const double vapour = -k / static_cast<double>(n);
                 for (int l=0; l<=n; ++l) {
                     const double alq = l / static_cast<double>(n);
                     for (int m=1; m<=n; ++m) {
-                        const double liquid = m / static_cast<double>(n);
+                        const double liquid = -m / static_cast<double>(n);
 
                         //Find values that should be in table
                         double flo = Opm::detail::getFlo(aqua, liquid, vapour, table.getFloType());
@@ -537,7 +577,7 @@ BOOST_AUTO_TEST_CASE(InterpolatePlane)
                         double gfr = Opm::detail::getGFR(aqua, liquid, vapour, table.getGFRType());
 
                         //Calculate reference
-                        double reference = thp + 2*wfr + 3*gfr+ 4*alq + 5*flo;
+                        double reference = thp + 2*wfr + 3*gfr+ 4*alq - 5*flo;
 
                         //Calculate actual
                         //Note order of arguments: id, aqua, liquid, vapour, thp, alq
@@ -545,7 +585,7 @@ BOOST_AUTO_TEST_CASE(InterpolatePlane)
 
 
                         double abs_diff = std::abs(actual - reference);
-                        double max_d = std::max(max_d, abs_diff);
+                        max_d = std::max(max_d, abs_diff);
                         sad = sad + abs_diff;
                     }
                 }
@@ -578,20 +618,20 @@ BOOST_AUTO_TEST_CASE(ExtrapolatePlane)
     for (int i=0; i<=n+o; ++i) {
         const double x = i / static_cast<double>(n);
         for (int j=1; j<=n+o; ++j) {
-            const double aqua = j / static_cast<double>(n);
+            const double aqua = -j / static_cast<double>(n);
             for (int k=1; k<=n+o; ++k) {
-                const double vapour = k / static_cast<double>(n);
+                const double vapour = -k / static_cast<double>(n);
                 for (int l=0; l<=n+o; ++l) {
                     const double u = l / static_cast<double>(n);
                     for (int m=1; m<=n+o; ++m) {
-                        const double liquid = m / static_cast<double>(n);
+                        const double liquid = -m / static_cast<double>(n);
 
                         //Find values that should be in table
                         double v = Opm::detail::getFlo(aqua, liquid, vapour, table.getFloType());
                         double y = Opm::detail::getWFR(aqua, liquid, vapour, table.getWFRType());
                         double z = Opm::detail::getGFR(aqua, liquid, vapour, table.getGFRType());
 
-                        double reference = x + 2*y + 3*z+ 4*u + 5*v;
+                        double reference = x + 2*y + 3*z+ 4*u - 5*v;
                         reference_sum += reference;
 
                         //Note order of arguments! id, aqua, liquid, vapour, thp , alq
@@ -635,13 +675,13 @@ BOOST_AUTO_TEST_CASE(ExtrapolatePlaneADB)
     for (int i=0; i<=n+o; ++i) {
         const double x = i / static_cast<double>(n);
         for (int j=1; j<=n+o; ++j) {
-            const double aqua = j / static_cast<double>(n);
+            const double aqua = -j / static_cast<double>(n);
             for (int k=1; k<=n+o; ++k) {
-                const double vapour = k / static_cast<double>(n);
+                const double vapour = -k / static_cast<double>(n);
                 for (int l=0; l<=n+o; ++l) {
                     const double u = l / static_cast<double>(n);
                     for (int m=1; m<=n+o; ++m) {
-                        const double liquid = m / static_cast<double>(n);
+                        const double liquid = -m / static_cast<double>(n);
 
                         //Temporary variables used to represent independent wells
                         const int num_wells = 5;
@@ -678,7 +718,7 @@ BOOST_AUTO_TEST_CASE(ExtrapolatePlaneADB)
                             double y = Opm::detail::getWFR(aqua*(w+1), liquid*(w+1), vapour*(w+1), table.getWFRType());
                             double z = Opm::detail::getGFR(aqua*(w+1), liquid*(w+1), vapour*(w+1), table.getGFRType());
 
-                            reference = x*(w+1) + 2*y + 3*z + 4*u*(w+1) + 5*v;
+                            reference = x*(w+1) + 2*y + 3*z + 4*u*(w+1) - 5*v;
                             value = bhp_val[w];
 
                             sum += value;
@@ -733,7 +773,7 @@ BOOST_AUTO_TEST_CASE(InterpolateADBAndQs)
     ADB::V qs_v(nphases*nwells);
     for (int j=0; j<nphases; ++j) {
         for (int i=0; i<nwells; ++i) {
-            qs_v[j*nwells+i] = (j*nwells+i) / static_cast<double>(nwells*nphases-1.0);
+            qs_v[j*nwells+i] = -(j*nwells+i) / static_cast<double>(nwells*nphases-1.0);
         }
     }
     ADB qs = ADB::constant(qs_v);
@@ -778,7 +818,7 @@ BOOST_AUTO_TEST_CASE(InterpolateADBAndQs)
         double flo = oil[i];
         double wor = water[i]/oil[i];
         double gor = gas[i]/oil[i];
-        reference[i] = thp_v[i] + 2*wor + 3*gor + 4*alq_v[i] + 5*flo;
+        reference[i] = thp_v[i] + 2*wor + 3*gor + 4*alq_v[i] - 5*flo;
     }
 
     //Check that interpolation matches
@@ -817,13 +857,13 @@ BOOST_AUTO_TEST_CASE(PartialDerivatives)
     for (int i=0; i<=n; ++i) {
         const double thp = i / static_cast<double>(n);
         for (int j=1; j<=n; ++j) {
-            const double aqua = j / static_cast<double>(n);
+            const double aqua = -j / static_cast<double>(n);
             for (int k=1; k<=n; ++k) {
-                const double vapour = k / static_cast<double>(n);
+                const double vapour = -k / static_cast<double>(n);
                 for (int l=0; l<=n; ++l) {
                     const double alq = l / static_cast<double>(n);
                     for (int m=1; m<=n; ++m) {
-                        const double liquid = m / static_cast<double>(n);
+                        const double liquid = -m / static_cast<double>(n);
 
                         //Find values that should be in table
                         double flo = Opm::detail::getFlo(aqua, liquid, vapour, table.getFloType());
@@ -832,7 +872,7 @@ BOOST_AUTO_TEST_CASE(PartialDerivatives)
 
                         //Calculate reference
                         VFPEvaluation reference;
-                        reference.value = thp + 2*wfr + 3*gfr+ 4*alq + 5*flo;
+                        reference.value = thp + 2*wfr + 3*gfr+ 4*alq - 5*flo;
                         reference.dthp = 1;
                         reference.dwfr = 2;
                         reference.dgfr = 3;
@@ -1040,7 +1080,7 @@ BOOST_AUTO_TEST_CASE(ParseInterpolateRealisticVFPPROD)
                 //for (unsigned int a=0; a<n; ++a) { //n==1, skip this loop
                     for (int f=0; f<n; ++f) {
                         //Liq given as SM3/day => convert to SM3/second
-                        double f_i = liq[f]*1.1574074074074073e-05;
+                        double f_i = -liq[f]*1.1574074074074073e-05;
 
                         //THP given as BARSA => convert to Pascal
                         double t_i = thp[t]*100000.0;