Merge pull request #1164 from stefoss23/deckkeyword_vector_unit

DeckKeyword w/ vectors: takes numpy vectors.

opm/parser/eclipse/Deck/DeckKeyword.hpp

@@ -44,7 +44,7 @@ namespace Opm {
         DeckKeyword(const ParserKeyword& parserKeyword, const Location& location, const std::string& keywordName);
         DeckKeyword(const ParserKeyword& parserKeyword, const std::vector<std::vector<DeckValue>>& record_list, UnitSystem& system_active, UnitSystem& system_default);
         DeckKeyword(const ParserKeyword& parserKeyword, const std::vector<int>& data);
-        DeckKeyword(const ParserKeyword& parserKeyword, const std::vector<double>& data);
+        DeckKeyword(const ParserKeyword& parserKeyword, const std::vector<double>& data, UnitSystem& system_active, UnitSystem& system_default);
 
         const std::string& name() const;
         void setFixedSize();

python/cxx/converters.hpp

@@ -3,6 +3,7 @@
 
 #include <sstream>
 #include <pybind11/pybind11.h>
+#include <pybind11/numpy.h>
 
 namespace py = pybind11;
 
@@ -20,4 +21,31 @@ std::string str( const T& t ) {
     return stream.str();
 }
 
+namespace convert {
+
+template <class T>
+std::vector<T> vector(py::array_t<T>& input) {
+    T * input_ptr    = (T *) input.request().ptr; 
+    std::vector<T> output(input.size());
+  
+    for (int i = 0; i < input.size(); i++)
+        output[i] = input_ptr[i];
+
+    return output;
+}
+
+
+template <class T>
+py::array_t<T> numpy_array(const std::vector<T>& input) {
+    auto output =  py::array_t<T>(input.size());
+    T * py_array_ptr = (T*)output.request().ptr;
+
+    for (size_t i = 0; i < input.size(); i++)
+        py_array_ptr[i] = input[i];
+
+    return output;
+}
+
+}
+
 #endif //SUNBEAM_CONVERTERS_HPP

python/cxx/deck_keyword.cpp

@@ -97,6 +97,19 @@ void push_string_as_deck_value(std::vector<DeckValue>& record, const std::string
 
 }
 
+
+py::array_t<int> get_int_array(const DeckKeyword& kw) {
+    return convert::numpy_array( kw.getIntData() );
+}
+
+py::array_t<double> get_raw_array(const DeckKeyword& kw) {
+    return convert::numpy_array( kw.getRawDoubleData() );
+}
+
+py::array_t<double> get_SI_array(const DeckKeyword& kw) {
+    return convert::numpy_array( kw.getSIDoubleData() );
+}
+
 }
 
 void python::common::export_DeckKeyword(py::module& module) {
@@ -148,7 +161,19 @@ void python::common::export_DeckKeyword(py::module& module) {
         .def( "__getitem__", getRecord, ref_internal)
         .def( "__len__", &DeckKeyword::size )
         .def_property_readonly("name", &DeckKeyword::name )
-        ;
+
+    .def(py::init([](const ParserKeyword& parser_keyword, py::array_t<int> py_data) {               
+            return DeckKeyword(parser_keyword, convert::vector(py_data));
+        } ) )
+
+    .def(py::init([](const ParserKeyword& parser_keyword, py::array_t<double> py_data, UnitSystem& active_system, UnitSystem& default_system) {               
+            return DeckKeyword(parser_keyword, convert::vector(py_data), active_system, default_system);
+        } ) )
+        
+    .def("get_int_array", &get_int_array)
+    .def("get_raw_array", &get_raw_array)
+    .def("get_SI_array", &get_SI_array)
+         ;
 
 
     py::class_< DeckRecord >( module, "DeckRecord")

python/tests/test_parser.py

@@ -2,6 +2,8 @@ import unittest
 import os.path
 import sys
 
+import numpy as np
+
 from opm.io.parser import Parser
 from opm.io.parser import ParseContext
 from opm.io.deck import DeckKeyword
@@ -53,7 +55,7 @@ FIPNUM
         deck = parser.parse_string(string)
         deck = parser.parse_string(string, context)
 
-    def test_create_deck_kw(self):
+    def test_deck_kw_records(self):
         parser = Parser()
         deck = parser.parse_string(self.REGIONDATA)
         active_unit_system = deck.active_unit_system()
@@ -128,6 +130,26 @@ FIPNUM
             raise DeckKeyword(parser["AQANTRC"], [["1*2.2", "ABC", 8]], active_unit_system, default_unit_system)
 
 
+    def test_deck_kw_vector(self):
+        parser = Parser()
+        deck = parser.parse_string(self.REGIONDATA)
+        active_unit_system = deck.active_unit_system()
+        default_unit_system = deck.default_unit_system()
+        self.assertEqual(active_unit_system.name, "Field")
+
+        int_array = np.array([0, 1, 2, 3])
+        hbnum_kw = DeckKeyword( parser["HBNUM"], int_array)
+        assert( np.array_equal(hbnum_kw.get_int_array(), int_array) )
+
+        raw_array = np.array([1.1, 2.2, 3.3])
+        zcorn_kw = DeckKeyword( parser["ZCORN"], raw_array, active_unit_system, default_unit_system)
+        assert( np.array_equal(zcorn_kw.get_raw_array(), raw_array) )
+        si_array = zcorn_kw.get_SI_array()
+        self.assertAlmostEqual( si_array[0], 1.1 * unit_foot )
+        self.assertAlmostEqual( si_array[2], 3.3 * unit_foot )
+        
+    
+
 if __name__ == "__main__":
     unittest.main()
 

src/opm/parser/eclipse/Deck/DeckKeyword.cpp

@@ -134,23 +134,12 @@ namespace Opm {
         const ParserItem& parser_item = parser_record.get(0);
 
         setDataKeyword();
-        DeckItem item;
-        if (parser_item.dataType() == type_tag::fdouble) {
-            // The unit system treatment here is a total hack to compile
-            UnitSystem unit_system(UnitSystem::UnitType::UNIT_TYPE_METRIC);
-            const auto& active_dim = unit_system.getDimension("1");
-            const auto& default_dim = unit_system.getDimension("1");
-            item = DeckItem(parser_item.name(), double(), { active_dim }, { default_dim });
-            for (double val : data)
-                item.push_back(val);
-        }
-        else if (parser_item.dataType() == type_tag::integer) {
-            item = DeckItem(parser_item.name(), int());
-            for (int val : data)
-                item.push_back(val);
-        }
-        else
+        if (parser_item.dataType() != type_tag::integer)
             throw std::invalid_argument("Input to DeckKeyword '" + name() + "': cannot be std::vector<int>.");
+      
+        DeckItem item(parser_item.name(), int() );
+        for (int val : data)
+            item.push_back(val);
 
         DeckRecord deck_record;
         deck_record.addItem( std::move(item) );
@@ -158,7 +147,7 @@ namespace Opm {
     }
 
 
-    DeckKeyword::DeckKeyword(const ParserKeyword& parserKeyword, const std::vector<double>& data) :
+    DeckKeyword::DeckKeyword(const ParserKeyword& parserKeyword, const std::vector<double>& data, UnitSystem& system_active, UnitSystem& system_default) :
         DeckKeyword(parserKeyword)
     {
         if (!parserKeyword.isDataKeyword())
@@ -171,11 +160,14 @@ namespace Opm {
         if (parser_item.dataType() != type_tag::fdouble)
             throw std::invalid_argument("Input to DeckKeyword '" + name() + "': cannot be std::vector<double>.");
 
-        // The unit system treatment here is a total hack to compile
-        UnitSystem unit_system(UnitSystem::UnitType::UNIT_TYPE_METRIC);
-        const auto& active_dim = unit_system.getDimension("1");
-        const auto& default_dim = unit_system.getDimension("1");
-        DeckItem item(parser_item.name(), double(), { active_dim }, { default_dim });
+        auto& dim = parser_item.dimensions();
+        std::vector<Dimension> active_dimensions;
+        std::vector<Dimension> default_dimensions;
+        if (dim.size() > 0) {
+             active_dimensions.push_back( system_active.getNewDimension(dim[0]) );
+             default_dimensions.push_back( system_default.getNewDimension(dim[0]) );
+        }
+        DeckItem item(parser_item.name(), double(), active_dimensions, default_dimensions);
         for (double val : data)
             item.push_back(val);
 

tests/parser/DeckValueTests.cpp

@@ -23,6 +23,8 @@
 
 #include <boost/test/unit_test.hpp>
 
+#include <opm/common/utility/numeric/cmp.hpp>
+
 #include <opm/parser/eclipse/Units/UnitSystem.hpp>
 
 #include <opm/parser/eclipse/Parser/Parser.hpp>
@@ -142,6 +144,11 @@ BOOST_AUTO_TEST_CASE(DeckKeywordConstructor) {
 BOOST_AUTO_TEST_CASE(DeckKeywordVectorInt) {
     
    Parser parser;
+   Deck deck;
+
+   UnitSystem& unit_default = deck.getDefaultUnitSystem();
+   UnitSystem  unit_active(UnitSystem::UnitType::UNIT_TYPE_LAB);
+
    const ParserKeyword& hbnum = parser.getKeyword("HBNUM");
    const ParserKeyword& box = parser.getKeyword("BOX");
 
@@ -153,27 +160,33 @@ BOOST_AUTO_TEST_CASE(DeckKeywordVectorInt) {
    BOOST_CHECK( hbnum_kw.getIntData() == data );   
 
    std::vector<double> data_double = {1.1, 2.2};
-   BOOST_CHECK_THROW(DeckKeyword(hbnum, data_double), std::invalid_argument);
+   BOOST_CHECK_THROW(DeckKeyword(hbnum, data_double, unit_active, unit_default), std::invalid_argument);
 }
 
 
 BOOST_AUTO_TEST_CASE(DeckKeywordVectorDouble) {
     
    Parser parser;
-   const ParserKeyword& poro = parser.getKeyword("PORO");
+   Deck deck;
+
+   UnitSystem& unit_default = deck.getDefaultUnitSystem();
+   UnitSystem  unit_active(UnitSystem::UnitType::UNIT_TYPE_LAB);
+
+   const ParserKeyword& zcorn = parser.getKeyword("ZCORN");  //vector of dim length
    const ParserKeyword& box = parser.getKeyword("BOX");
 
    std::vector<double> data = {1.1, 2.2, 3.3};
-   BOOST_CHECK_THROW(DeckKeyword(box, data), std::invalid_argument);
-   DeckKeyword poro_kw(poro, data);
-   BOOST_CHECK(poro_kw.isDataKeyword());
-   BOOST_CHECK_EQUAL(poro_kw.getDataSize(), 3);
-   BOOST_CHECK( poro_kw.getRawDoubleData() == data );
+   
+   BOOST_CHECK_THROW(DeckKeyword(box, data, unit_active, unit_default), std::invalid_argument);
+   DeckKeyword zcorn_kw(zcorn, data, unit_active, unit_default);
+   BOOST_CHECK(zcorn_kw.isDataKeyword());
+   BOOST_CHECK_EQUAL(zcorn_kw.getDataSize(), 3);
+   BOOST_CHECK( zcorn_kw.getRawDoubleData() == data );
+   std::vector<double> SI_data = zcorn_kw.getSIDoubleData();
+   BOOST_CHECK( cmp::scalar_equal<double>(SI_data[0], 0.011) );
+   BOOST_CHECK( cmp::scalar_equal<double>(SI_data[1], 0.022) );
+   BOOST_CHECK( cmp::scalar_equal<double>(SI_data[2], 0.033) );
 
-   std::vector<int> data_int = {1, 2};
-   poro_kw = DeckKeyword(poro, data_int);
-   std::vector<double> raw_int = {1.0, 2.0};
-   BOOST_CHECK( poro_kw.getRawDoubleData() == raw_int );
 }