From fa572188fa040cb36c63e48119e78b398e4a7492 Mon Sep 17 00:00:00 2001
From: Arne Morten Kvarving <arne.morten.kvarving@sintef.no>
Date: Wed, 12 Feb 2020 15:04:21 +0100
Subject: [PATCH] changed: avoid use of boost::multi_array in VFPProdTable

---
 .../EclipseState/Schedule/VFPProdTable.hpp    |  21 +--
 .../EclipseState/Schedule/VFPProdTable.cpp    | 125 ++++++++----------
 tests/parser/TableManagerTests.cpp            |  21 ++-
 3 files changed, 76 insertions(+), 91 deletions(-)

diff --git a/opm/parser/eclipse/EclipseState/Schedule/VFPProdTable.hpp b/opm/parser/eclipse/EclipseState/Schedule/VFPProdTable.hpp
index 495f90efc..916918e82 100644
--- a/opm/parser/eclipse/EclipseState/Schedule/VFPProdTable.hpp
+++ b/opm/parser/eclipse/EclipseState/Schedule/VFPProdTable.hpp
@@ -21,7 +21,8 @@
 #define OPM_PARSER_ECLIPSE_ECLIPSESTATE_TABLES_VFPPRODTABLE_HPP_
 
 
-#include <boost/multi_array.hpp>
+#include <array>
+#include <vector>
 
 namespace Opm {
 
@@ -34,8 +35,7 @@ namespace Opm {
  */
 class VFPProdTable {
 public:
-    typedef boost::multi_array<double, 5> array_type;
-    typedef boost::array<array_type::index, 5> extents;
+    typedef std::vector<double> array_type;
 
     enum FLO_TYPE {
         FLO_OIL=1,
@@ -130,12 +130,10 @@ public:
     }
 
     /**
-     * Returns the data of the table itself. The data is ordered so that
+     * Returns the data of the table itself. For ordered access
+     * use operator()(thp_idx, wfr_idx, gfr_idx, alq_idx, flo_idx)
      *
-     * table = getTable();
-     * bhp = table[thp_idx][wfr_idx][gfr_idx][alq_idx][flo_idx];
-     *
-     * gives the bottom hole pressure value in the table for the coordinate
+     * This gives the bottom hole pressure value in the table for the coordinate
      * given by
      * flo_axis = getFloAxis();
      * thp_axis = getTHPAxis();
@@ -150,7 +148,10 @@ public:
     }
 
     bool operator==(const VFPProdTable& data) const;
-    VFPProdTable& operator=(const VFPProdTable& data);
+
+    std::array<size_t,5> shape() const;
+
+    double operator()(size_t thp_idx, size_t wfr_idx, size_t gfr_idx, size_t alq_idx, size_t flo_idx) const;
 
 private:
 
@@ -171,6 +172,8 @@ private:
 
     void check(const DeckKeyword& table, const double factor);
 
+    double& operator()(size_t thp_idx, size_t wfr_idx, size_t gfr_idx, size_t alq_idx, size_t flo_idx);
+
     static void scaleValues(std::vector<double>& values,
                             const double& scaling_factor);
 
diff --git a/src/opm/parser/eclipse/EclipseState/Schedule/VFPProdTable.cpp b/src/opm/parser/eclipse/EclipseState/Schedule/VFPProdTable.cpp
index 744c594a3..cca2e5ec8 100644
--- a/src/opm/parser/eclipse/EclipseState/Schedule/VFPProdTable.cpp
+++ b/src/opm/parser/eclipse/EclipseState/Schedule/VFPProdTable.cpp
@@ -18,6 +18,8 @@
 */
 
 #include <algorithm>
+#include <cassert>
+#include <cmath>
 #include <iostream>
 
 #include <opm/common/OpmLog/OpmLog.hpp>
@@ -159,14 +161,6 @@ VFPProdTable::VFPProdTable(int table_num,
     m_wfr_data = wfr_data;
     m_gfr_data = gfr_data;
     m_alq_data = alq_data;
-
-    extents shape;
-    shape[0] = data.shape()[0];
-    shape[1] = data.shape()[1];
-    shape[2] = data.shape()[2];
-    shape[3] = data.shape()[3];
-    shape[4] = data.shape()[4];
-    m_data.resize(shape);
     m_data = data;
 
     //check();
@@ -276,14 +270,8 @@ VFPProdTable::VFPProdTable( const DeckKeyword& table, const UnitSystem& deck_uni
     size_t ng = m_gfr_data.size();
     size_t na = m_alq_data.size();
     size_t nf = m_flo_data.size();
-    extents shape;
-    shape[0] = nt;
-    shape[1] = nw;
-    shape[2] = ng;
-    shape[3] = na;
-    shape[4] = nf;
-    m_data.resize(shape);
-    std::fill_n(m_data.data(), m_data.num_elements(), std::nan("0"));
+    m_data.resize(nt*nw*ng*na*nf);
+    std::fill_n(m_data.data(), m_data.size(), std::nan("0"));
 
     //Check that size of table matches size of axis:
     if (table.size() != nt*nw*ng*na + 6) {
@@ -315,7 +303,7 @@ VFPProdTable::VFPProdTable( const DeckKeyword& table, const UnitSystem& deck_uni
                         << t << "," << w << "," << g << "," << a << "," << f
                         << "]=" << bhp_tht[f] << " too large" << std::endl;
             }
-            m_data[t][w][g][a][f] = table_scaling_factor*bhp_tht[f];
+            (*this)(t,w,g,a,f) = table_scaling_factor*bhp_tht[f];
         }
     }
 
@@ -355,54 +343,42 @@ void VFPProdTable::check(const DeckKeyword& keyword, const double table_scaling_
     assert(std::is_sorted(m_gfr_data.begin(), m_gfr_data.end()));
     assert(std::is_sorted(m_alq_data.begin(), m_alq_data.end()));
 
+    size_t nt = m_thp_data.size();
+    size_t nw = m_wfr_data.size();
+    size_t ng = m_gfr_data.size();
+    size_t na = m_alq_data.size();
+    size_t nf = m_flo_data.size();
+
     //Check data size matches axes
-    assert(m_data.num_dimensions() == 5);
-    assert(m_data.shape()[0] == m_thp_data.size());
-    assert(m_data.shape()[1] == m_wfr_data.size());
-    assert(m_data.shape()[2] == m_gfr_data.size());
-    assert(m_data.shape()[3] == m_alq_data.size());
-    assert(m_data.shape()[4] == m_flo_data.size());
+    assert(m_data.size() == nt*nw*ng*na*nf);
 
-
-    //Check that all elements have been set
+    //Check that bhp(thp) is a monotonic increasing function.
+    //If this is not the case, we might not be able to determine
+    //the thp from the bhp easily
     typedef array_type::size_type size_type;
-    for (size_type t=0; t<m_data.shape()[0]; ++t) {
-        for (size_type w=0; w<m_data.shape()[1]; ++w) {
-            for (size_type g=0; g<m_data.shape()[2]; ++g) {
-                for (size_type a=0; a<m_data.shape()[3]; ++a) {
-                    for (size_type f=0; f<m_data.shape()[4]; ++f) {
-                        if (std::isnan(m_data[t][w][g][a][f])) {
+    std::string points;
+    for (size_type w = 0; w < nw; ++w) {
+        for (size_type g = 0; g < ng; ++g) {
+            for (size_type a = 0; a < na; ++a) {
+                for (size_type f = 0; f < nf; ++f) {
+                    double bhp_last = (*this)(0,w,g,a,f);
+                    for (size_type t = 0; t < nt; ++t) {
+                        //Check that all elements have been set
+                        if (std::isnan((*this)(t,w,g,a,f))) {
                             //TODO: Replace with proper log message
                             std::cerr << "VFPPROD element ["
                                     << t << "," << w << "," << g << "," << a << "," << f
                                     << "] not set!" << std::endl;
                             throw std::invalid_argument("Missing VFPPROD value");
                         }
-                    }
-                }
-            }
-        }
-    }
-
-
-    //Check that bhp(thp) is a monotonic increasing function.
-    //If this is not the case, we might not be able to determine
-    //the thp from the bhp easily
-    std::string points = "";
-    for (size_type w=0; w<m_data.shape()[1]; ++w) {
-        for (size_type g=0; g<m_data.shape()[2]; ++g) {
-            for (size_type a=0; a<m_data.shape()[3]; ++a) {
-                for (size_type f=0; f<m_data.shape()[4]; ++f) {
-                    double bhp_last = m_data[0][w][g][a][f];
-                    for (size_type t=0; t<m_data.shape()[0]; ++t) {
-                        if (m_data[t][w][g][a][f] < bhp_last) {
+                        if ((*this)(t,w,g,a,f) < bhp_last) {
                             points += "At point (FLOW, THP, WFR, GFR, ALQ) = "
                                     + std::to_string(f) + " " + std::to_string(t) + " "
                                     + std::to_string(w) + " " + std::to_string(g) + " "
                                     + std::to_string(a) + " at BHP = "
-                                    + std::to_string(m_data[t][w][g][a][f] / table_scaling_factor) + "\n";
+                                    + std::to_string((*this)(t,w,g,a,f) / table_scaling_factor) + "\n";
                         }
-                        bhp_last = m_data[t][w][g][a][f];
+                        bhp_last = (*this)(t,w,g,a,f);
                     }
                 }
             }
@@ -543,25 +519,34 @@ bool VFPProdTable::operator==(const VFPProdTable& data) const {
 }
 
 
-VFPProdTable& VFPProdTable::operator=(const VFPProdTable& data) {
-    m_table_num = data.m_table_num;
-    m_datum_depth = data.m_datum_depth;
-    m_flo_type = data.m_flo_type;
-    m_wfr_type = data.m_wfr_type;
-    m_gfr_type = data.m_gfr_type;
-    m_alq_type = data.m_alq_type;
-    m_flo_data = data.m_flo_data;
-    m_thp_data = data.m_thp_data;
-    m_wfr_data = data.m_wfr_data;
-    m_gfr_data = data.m_gfr_data;
-    m_alq_data = data.m_alq_data;
-    extents shape;
-    for (size_t i = 0; i < 5; ++i)
-        shape[i] = data.m_data.shape()[i];
-    m_data.resize(shape);
-    for (size_t i = 0; i < data.m_data.num_elements(); ++i)
-        *(m_data.data() + i) = *(data.m_data.data() + i);
-    return *this;
+double VFPProdTable::operator()(size_t thp_idx, size_t wfr_idx, size_t gfr_idx, size_t alq_idx, size_t flo_idx) const {
+    size_t nw = m_wfr_data.size();
+    size_t ng = m_gfr_data.size();
+    size_t na = m_alq_data.size();
+    size_t nf = m_flo_data.size();
+
+    return m_data[thp_idx*nw*ng*na*nf + wfr_idx*ng*na*nf + gfr_idx*na*nf + alq_idx*nf + flo_idx];
+}
+
+
+double& VFPProdTable::operator()(size_t thp_idx, size_t wfr_idx, size_t gfr_idx, size_t alq_idx, size_t flo_idx) {
+    size_t nw = m_wfr_data.size();
+    size_t ng = m_gfr_data.size();
+    size_t na = m_alq_data.size();
+    size_t nf = m_flo_data.size();
+
+    return m_data[thp_idx*nw*ng*na*nf + wfr_idx*ng*na*nf + gfr_idx*na*nf + alq_idx*nf + flo_idx];
+}
+
+
+std::array<size_t,5> VFPProdTable::shape() const {
+    size_t nt = m_thp_data.size();
+    size_t nw = m_wfr_data.size();
+    size_t ng = m_gfr_data.size();
+    size_t na = m_alq_data.size();
+    size_t nf = m_flo_data.size();
+
+    return {nt, nw, ng, na, nf};
 }
 
 } //Namespace opm
diff --git a/tests/parser/TableManagerTests.cpp b/tests/parser/TableManagerTests.cpp
index 099f20d48..8dd63162f 100644
--- a/tests/parser/TableManagerTests.cpp
+++ b/tests/parser/TableManagerTests.cpp
@@ -737,8 +737,7 @@ VFPPROD \n\
     //The data itself
     {
         typedef Opm::VFPProdTable::array_type::size_type size_type;
-        const Opm::VFPProdTable::array_type& data = vfpprodTable.getTable();
-        const size_type* size = data.shape();
+        const auto size = vfpprodTable.shape();
 
         BOOST_CHECK_EQUAL(size[0], 2);
         BOOST_CHECK_EQUAL(size[1], 2);
@@ -750,13 +749,13 @@ VFPPROD \n\
         double conversion_factor = 100000.0;
 
         double index = 0.5;
-        for (size_type a=0; a<size[3]; ++a) {
-            for (size_type g=0; g<size[2]; ++g) {
-                for (size_type w=0; w<size[1]; ++w) {
-                    for (size_type t=0; t<size[0]; ++t) {
-                        for (size_type f=0; f<size[4]; ++f) {
+        for (size_type a = 0; a < size[3]; ++a) {
+            for (size_type g = 0;  g < size[2]; ++g) {
+                for (size_type w = 0; w < size[1]; ++w) {
+                    for (size_type t = 0; t < size[0]; ++t) {
+                        for (size_type f = 0; f < size[4]; ++f) {
                             index += 1.0;
-                            BOOST_CHECK_EQUAL(data[t][w][g][a][f], index*conversion_factor);
+                            BOOST_CHECK_EQUAL(const_cast<const VFPProdTable&>(vfpprodTable)(t,w,g,a,f), index*conversion_factor);
                         }
                     }
                 }
@@ -855,15 +854,13 @@ VFPPROD \n\
 
     //The data itself
     {
-        typedef Opm::VFPProdTable::array_type::size_type size_type;
-        const Opm::VFPProdTable::array_type& data = vfpprodTable.getTable();
-        const size_type* size = data.shape();
+        const auto size = vfpprodTable.shape();
 
         //Table given as BHP => barsa. Convert to pascal
         double conversion_factor = 100000.0;
 
         BOOST_CHECK_EQUAL(size[0]*size[1]*size[2]*size[3]*size[4], 1);
-        BOOST_CHECK_EQUAL(data[0][0][0][0][0], 1.5*conversion_factor);
+        BOOST_CHECK_EQUAL(const_cast<const VFPProdTable&>(vfpprodTable)(0,0,0,0,0), 1.5*conversion_factor);
     }
 }