python bindings, adding support for UDA type

opm/parser/eclipse/Deck/DeckItem.hpp

@@ -137,6 +137,13 @@ namespace Opm {
         const std::vector<Dimension>& activeDimensions() const;
         const std::vector<Dimension>& defaultDimensions() const;
 
+        bool is_uda() { return  (type == get_type< UDAValue >()); };
+        bool is_double() { return  type == get_type< double >(); };
+        bool is_int() { return  type == get_type< int >() ; };
+        bool is_string() { return  type == get_type< std::string >(); };
+
+        UDAValue& get_uda() { return uval[0]; };
+
     private:
         mutable std::vector< double > dval;
         std::vector< int > ival;

opm/parser/eclipse/Deck/UDAValue.hpp

@@ -64,6 +64,9 @@ public:
 
     bool operator==(const UDAValue& other) const;
     bool operator!=(const UDAValue& other) const;
+
+    bool is_numeric() { return numeric_value; }
+
 private:
     bool numeric_value;
     double double_value;

python/cxx/deck_keyword.cpp

@@ -110,8 +110,30 @@ py::array_t<double> get_SI_array(const DeckKeyword& kw) {
     return convert::numpy_array( kw.getSIDoubleData() );
 }
 
+bool uda_item_is_numeric(DeckItem *  item)
+{
+    if( !item->is_uda() )
+        throw std::logic_error("deck item doesn't support user defined quantities");
+
+    UDAValue uda = item->get_uda();
+
+    return uda.is_numeric();
 }
 
+double get_uda_double(DeckItem *  item)
+{
+    UDAValue uda = item->get_uda();
+    return uda.get<double>();
+}
+
+std::string get_uda_str(DeckItem *  item)
+{
+    UDAValue uda = item->get_uda();
+    return uda.get<std::string>();
+}
+
+
+}
 void python::common::export_DeckKeyword(py::module& module) {
     py::class_< DeckKeyword >( module, "DeckKeyword")
         .def(py::init<const ParserKeyword& >())
@@ -186,6 +208,10 @@ void python::common::export_DeckKeyword(py::module& module) {
 
     py::class_< DeckItem >(module, "DeckItem")
         .def( "__len__", &DeckItem::data_size )
+        .def("is_uda", &DeckItem::is_uda)
+        .def("is_double", &DeckItem::is_double)
+        .def("is_int", &DeckItem::is_int)
+        .def("is_string", &DeckItem::is_string)
         .def("get_str", &DeckItem::get<std::string>)
         .def("get_int", &DeckItem::get<int>)
         .def("get_raw", &DeckItem::get<double>)
@@ -193,8 +219,11 @@ void python::common::export_DeckKeyword(py::module& module) {
         .def("get_data_list", &item_to_pylist)
         .def("get_raw_data_list", &raw_data_to_pylist)
         .def("get_SI_data_list", &SI_data_to_pylist)
-        .def("has_value", &DeckItem::hasValue)
+        .def("__has_value", &DeckItem::hasValue)
-        .def("defaulted", &DeckItem::defaultApplied)
+        .def("__defaulted", &DeckItem::defaultApplied)
+        .def("__is_numberic", &uda_item_is_numeric)
+        .def("__uda_double", &get_uda_double)
+        .def("__uda_str", &get_uda_str)
         ;
 
 

python/python/opm/_common.py

@@ -13,6 +13,7 @@ from .libopmcommon_python import action
 
 from .libopmcommon_python import Parser, ParseContext
 from .libopmcommon_python import DeckKeyword
+from .libopmcommon_python import DeckItem
 
 from .libopmcommon_python import EclipseState
 from .libopmcommon_python import FieldProperties

python/python/opm/io/deck/__init__.py

@@ -1 +1,40 @@
 from opm._common import DeckKeyword
+from opm._common import DeckItem
+
+
+# in prinsiple it should be possible to use the  has_value(int) function
+# on not only the first element. However, in 99% of the use cases it is the
+# first element which is of interesst. Hence for python bindings this is
+# hardcoded to the first element
+
+@property
+def defaulted_deckitem(self):
+    return self.__defaulted(0)
+
+
+@property
+def has_value_deckitem(self):
+    return self.__has_value(0)
+
+
+@property
+def get_item_deckitem(self):
+
+    if self.is_int():
+        return self.get_int(0)
+    elif self.is_string():
+        return self.get_str(0)
+    elif self.is_double():
+        return self.get_raw(0)
+    elif self.is_uda():
+        if self.__is_numberic():
+            return self.__uda_double()
+        else:
+            return self.__uda_str()
+    else:
+        raise ValueError("Deck Item, unknown type")
+
+
+setattr(DeckItem, "defaulted", defaulted_deckitem)
+setattr(DeckItem, "valid", has_value_deckitem)
+setattr(DeckItem, "value", get_item_deckitem)

python/tests/data/SPE1CASE1.DATA

@@ -0,0 +1,480 @@
+-- This reservoir simulation deck is made available under the Open Database
+-- License: http://opendatacommons.org/licenses/odbl/1.0/. Any rights in
+-- individual contents of the database are licensed under the Database Contents
+-- License: http://opendatacommons.org/licenses/dbcl/1.0/
+
+-- Copyright (C) 2015 Statoil
+
+-- This simulation is based on the data given in 
+-- 'Comparison of Solutions to a Three-Dimensional
+-- Black-Oil Reservoir Simulation Problem' by Aziz S. Odeh,
+-- Journal of Petroleum Technology, January 1981
+
+
+---------------------------------------------------------------------------
+------------------------ SPE1 - CASE 1 ------------------------------------
+---------------------------------------------------------------------------
+
+RUNSPEC
+-- -------------------------------------------------------------------------
+
+TITLE
+   SPE1 - CASE 1
+
+DIMENS
+   10 10 3 /
+
+-- The number of equilibration regions is inferred from the EQLDIMS
+-- keyword.
+EQLDIMS
+/
+
+-- The number of PVTW tables is inferred from the TABDIMS keyword;
+-- when no data is included in the keyword the default values are used.
+TABDIMS
+/
+
+REGDIMS
+ 2 /
+
+
+OIL
+GAS
+WATER
+DISGAS
+-- As seen from figure 4 in Odeh, GOR is increasing with time,
+-- which means that dissolved gas is present
+
+FIELD
+
+START
+   1 'JAN' 2015 /
+
+WELLDIMS
+-- Item 1: maximum number of wells in the model
+-- 	   - there are two wells in the problem; injector and producer
+-- Item 2: maximum number of grid blocks connected to any one well
+-- 	   - must be one as the wells are located at specific grid blocks
+-- Item 3: maximum number of groups in the model
+-- 	   - we are dealing with only one 'group'
+-- Item 4: maximum number of wells in any one group
+-- 	   - there must be two wells in a group as there are two wells in total
+   2 1 1 2 /
+
+UNIFOUT
+
+
+UDQDIMS
+ 16 16 0 16 16 0 0 50 0 20 /
+
+-- Dimensions for the user defined arguments facility
+-- number of keyword arguments in which UDQs replace numerical values
+-- ratained for back-compatibility
+-- total number of unique instances in which a UDQ is used in a keyword argument
+UDADIMS
+ 10   1*  10 /
+
+
+GRID
+
+-- The INIT keyword is used to request an .INIT file. The .INIT file
+-- is written before the simulation actually starts, and contains grid
+-- properties and saturation tables as inferred from the input
+-- deck. There are no other keywords which can be used to configure
+-- exactly what is written to the .INIT file.
+INIT
+
+
+-- -------------------------------------------------------------------------
+NOECHO
+
+DX 
+-- There are in total 300 cells with length 1000ft in x-direction	
+   	300*1000 /
+DY
+-- There are in total 300 cells with length 1000ft in y-direction	
+	300*1000 /
+DZ
+-- The layers are 20, 30 and 50 ft thick, in each layer there are 100 cells
+	100*20 100*30 100*50 /
+
+TOPS
+-- The depth of the top of each grid block
+	100*8325 /
+
+PORO
+-- Constant porosity of 0.3 throughout all 300 grid cells
+   	300*0.3 /
+
+PERMX
+-- The layers have perm. 500mD, 50mD and 200mD, respectively.
+	100*500 100*50 100*200 /
+
+PERMY
+-- Equal to PERMX
+	100*500 100*50 100*200 /
+
+PERMZ
+-- Cannot find perm. in z-direction in Odeh's paper
+-- For the time being, we will assume PERMZ equal to PERMX and PERMY:
+	100*500 100*50 100*200 /
+ECHO
+
+PROPS
+-- -------------------------------------------------------------------------
+
+PVTW
+-- Item 1: pressure reference (psia)
+-- Item 2: water FVF (rb per bbl or rb per stb)
+-- Item 3: water compressibility (psi^{-1})
+-- Item 4: water viscosity (cp)
+-- Item 5: water 'viscosibility' (psi^{-1})
+
+-- Using values from Norne:
+-- In METRIC units:
+-- 	277.0 1.038 4.67E-5 0.318 0.0 /
+-- In FIELD units:
+    	4017.55 1.038 3.22E-6 0.318 0.0 /
+
+ROCK
+-- Item 1: reference pressure (psia)
+-- Item 2: rock compressibility (psi^{-1})
+
+-- Using values from table 1 in Odeh:
+	14.7 3E-6 /
+
+SWOF
+-- Column 1: water saturation
+--   	     - this has been set to (almost) equally spaced values from 0.12 to 1
+-- Column 2: water relative permeability
+--   	     - generated from the Corey-type approx. formula
+--	       the coeffisient is set to 10e-5, S_{orw}=0 and S_{wi}=0.12
+-- Column 3: oil relative permeability when only oil and water are present
+--	     - we will use the same values as in column 3 in SGOF.
+-- 	       This is not really correct, but since only the first 
+--	       two values are of importance, this does not really matter
+-- Column 4: water-oil capillary pressure (psi) 
+
+0.12	0    		 	1	0
+0.18	4.64876033057851E-008	1	0
+0.24	0.000000186		0.997	0
+0.3	4.18388429752066E-007	0.98	0
+0.36	7.43801652892562E-007	0.7	0
+0.42	1.16219008264463E-006	0.35	0
+0.48	1.67355371900826E-006	0.2	0
+0.54	2.27789256198347E-006	0.09	0
+0.6	2.97520661157025E-006	0.021	0
+0.66	3.7654958677686E-006	0.01	0
+0.72	4.64876033057851E-006	0.001	0
+0.78	0.000005625		0.0001	0
+0.84	6.69421487603306E-006	0	0
+0.91	8.05914256198347E-006	0	0
+1	0.00001			0	0 /
+
+
+SGOF
+-- Column 1: gas saturation
+-- Column 2: gas relative permeability
+-- Column 3: oil relative permeability when oil, gas and connate water are present
+-- Column 4: oil-gas capillary pressure (psi)
+-- 	     - stated to be zero in Odeh's paper
+
+-- Values in column 1-3 are taken from table 3 in Odeh's paper:
+0	0	1	0
+0.001	0	1	0
+0.02	0	0.997	0
+0.05	0.005	0.980	0
+0.12	0.025	0.700	0
+0.2	0.075	0.350	0
+0.25	0.125	0.200	0
+0.3	0.190	0.090	0
+0.4	0.410	0.021	0
+0.45	0.60	0.010	0
+0.5	0.72	0.001	0
+0.6	0.87	0.0001	0
+0.7	0.94	0.000	0
+0.85	0.98	0.000	0 
+0.88	0.984	0.000	0 /
+--1.00	1.0	0.000	0 /
+-- Warning from Eclipse: first sat. value in SWOF + last sat. value in SGOF
+-- 	   		 must not be greater than 1, but Eclipse still runs
+-- Flow needs the sum to be excactly 1 so I added a row with gas sat. =  0.88
+-- The corresponding krg value was estimated by assuming linear rel. between
+-- gas sat. and krw. between gas sat. 0.85 and 1.00 (the last two values given)
+
+DENSITY
+-- Density (lb per ft³) at surface cond. of 
+-- oil, water and gas, respectively (in that order)
+
+-- Using values from Norne:
+-- In METRIC units:
+--      859.5 1033.0 0.854 /
+-- In FIELD units:
+      	53.66 64.49 0.0533 /
+
+PVDG
+-- Column 1: gas phase pressure (psia)
+-- Column 2: gas formation volume factor (rb per Mscf)
+-- 	     - in Odeh's paper the units are said to be given in rb per bbl, 
+-- 	       but this is assumed to be a mistake: FVF-values in Odeh's paper 
+--	       are given in rb per scf, not rb per bbl. This will be in 
+--	       agreement with conventions
+-- Column 3: gas viscosity (cP)
+
+-- Using values from lower right table in Odeh's table 2:
+14.700	166.666	0.008000
+264.70	12.0930	0.009600
+514.70	6.27400	0.011200
+1014.7	3.19700	0.014000
+2014.7	1.61400	0.018900
+2514.7	1.29400	0.020800
+3014.7	1.08000	0.022800
+4014.7	0.81100	0.026800
+5014.7	0.64900	0.030900
+9014.7	0.38600	0.047000 /
+
+PVTO
+-- Column 1: dissolved gas-oil ratio (Mscf per stb)
+-- Column 2: bubble point pressure (psia)
+-- Column 3: oil FVF for saturated oil (rb per stb)
+-- Column 4: oil viscosity for saturated oil (cP)
+
+-- Use values from top left table in Odeh's table 2:
+0.0010	14.7	1.0620	1.0400 /
+0.0905	264.7	1.1500	0.9750 /
+0.1800	514.7	1.2070	0.9100 /
+0.3710	1014.7	1.2950	0.8300 /
+0.6360	2014.7	1.4350	0.6950 /
+0.7750	2514.7	1.5000	0.6410 /
+0.9300	3014.7	1.5650	0.5940 /
+1.2700	4014.7	1.6950	0.5100 
+	9014.7	1.5790	0.7400 /
+1.6180	5014.7	1.8270	0.4490 
+	9014.7	1.7370	0.6310 /	
+-- It is required to enter data for undersaturated oil for the highest GOR
+-- (i.e. the last row) in the PVTO table.
+-- In order to fulfill this requirement, values for oil FVF and viscosity
+-- at 9014.7psia and GOR=1.618 for undersaturated oil have been approximated:
+-- It has been assumed that there is a linear relation between the GOR
+-- and the FVF when keeping the pressure constant at 9014.7psia.
+-- From Odeh we know that (at 9014.7psia) the FVF is 2.357 at GOR=2.984
+-- for saturated oil and that the FVF is 1.579 at GOR=1.27 for undersaturated oil,
+-- so it is possible to use the assumption described above. 
+-- An equivalent approximation for the viscosity has been used.
+/
+
+REGIONS
+
+EQLNUM
+  300*1 /
+
+
+FIPNUM
+  100*1 100* 100*2 /
+
+
+SOLUTION
+-- -------------------------------------------------------------------------
+
+EQUIL
+-- Item 1: datum depth (ft)
+-- Item 2: pressure at datum depth (psia)
+-- 	   - Odeh's table 1 says that initial reservoir pressure is 
+-- 	     4800 psi at 8400ft, which explains choice of item 1 and 2
+-- Item 3: depth of water-oil contact (ft)
+-- 	   - chosen to be directly under the reservoir
+-- Item 4: oil-water capillary pressure at the water oil contact (psi)
+-- 	   - given to be 0 in Odeh's paper
+-- Item 5: depth of gas-oil contact (ft)
+-- 	   - chosen to be directly above the reservoir
+-- Item 6: gas-oil capillary pressure at gas-oil contact (psi)
+-- 	   - given to be 0 in Odeh's paper
+-- Item 7: RSVD-table
+-- Item 8: RVVD-table
+-- Item 9: Set to 0 as this is the only value supported by OPM
+
+-- Item #: 1 2    3    4 5    6 7 8 9
+	8400 4800 8450 0 8300 0 1 0 0 /
+
+RSVD
+-- Dissolved GOR is initially constant with depth through the reservoir.
+-- The reason is that the initial reservoir pressure given is higher 
+---than the bubble point presssure of 4014.7psia, meaning that there is no 
+-- free gas initially present.
+8300 1.270
+8450 1.270 /
+
+SUMMARY
+-- -------------------------------------------------------------------------	 
+
+-- 1a) Oil rate vs time
+FOPR
+-- Field Oil Production Rate
+
+-- 1b) GOR vs time
+WGOR
+-- Well Gas-Oil Ratio
+   'PROD'
+/
+-- Using FGOR instead of WGOR:PROD results in the same graph
+FGOR
+
+-- 2a) Pressures of the cell where the injector and producer are located
+BPR
+1  1  1 /
+10 10 3 /
+/
+
+-- 2b) Gas saturation at grid points given in Odeh's paper
+BGSAT
+1  1  1 /
+1  1  2 /
+1  1  3 /
+10 1  1 /
+10 1  2 /
+10 1  3 /
+10 10 1 /
+10 10 2 /
+10 10 3 /
+/
+
+-- In order to compare Eclipse with Flow:
+WBHP
+  'INJ'
+  'PROD'
+/
+WGIR
+  'INJ'
+  'PROD'
+/
+WGIT
+  'INJ'
+  'PROD'
+/
+WGPR
+  'INJ'
+  'PROD'
+/
+WGPT
+  'INJ'
+  'PROD'
+/
+WOIR
+  'INJ'
+  'PROD'
+/
+WOIT
+  'INJ'
+  'PROD'
+/
+WOPR
+  'INJ'
+  'PROD'
+/
+WOPT
+  'INJ'
+  'PROD'
+/
+WWIR
+  'INJ'
+  'PROD'
+/
+WWIT
+  'INJ'
+  'PROD'
+/
+WWPR
+  'INJ'
+  'PROD'
+/
+WWPT
+  'INJ'
+  'PROD'
+/
+
+WUOPRL
+ 'PROD' /
+
+
+
+
+SCHEDULE
+-- -------------------------------------------------------------------------
+
+
+UDQ
+-- test 
+--oil & liquid capacities at GEFAC = 0.8995
+DEFINE WUOPRL (20000 - TIME * 2.5) /
+-- units
+UNITS  WUOPRL SM3/DAY /
+--
+/
+
+
+RPTSCHED
+	'PRES' 'SGAS' 'RS' 'WELLS' 'WELSPECS' /
+
+RPTRST
+	'BASIC=1' /
+
+
+-- If no resolution (i.e. case 1), the two following lines must be added:
+DRSDT
+ 0 /
+-- if DRSDT is set to 0, GOR cannot rise and free gas does not 
+-- dissolve in undersaturated oil -> constant bubble point pressure
+
+WELSPECS
+-- Item #: 1	 2	3	4	5	 6
+	'PROD'	'G1'	10	10	8400	'OIL' /
+	'INJ'	'G1'	1	1	8335	'GAS' /
+/
+-- Coordinates in item 3-4 are retrieved from Odeh's figure 1 and 2
+-- Note that the depth at the midpoint of the well grid blocks
+-- has been used as reference depth for bottom hole pressure in item 5
+
+COMPDAT
+-- Item #: 1	2	3	4	5	6	7	8	9
+	'PROD'	10	10	3	3	'OPEN'	1*	1*	0.5 /
+	'INJ'	1	1	1	1	'OPEN'	1*	1*	0.5 /
+/
+-- Coordinates in item 2-5 are retreived from Odeh's figure 1 and 2 
+-- Item 9 is the well bore internal diameter, 
+-- the radius is given to be 0.25ft in Odeh's paper
+
+
+WCONPROD
+-- Item #:1	2      3     4	   5  9
+	'PROD' 'OPEN' 'ORAT' WUOPRL 1* 1.5E5 2* 1000 /
+/
+-- It is stated in Odeh's paper that the maximum oil prod. rate
+-- is 20 000stb per day which explains the choice of value in item 4.
+-- The items > 4 are defaulted with the exception of item  9,
+-- the BHP lower limit, which is given to be 1000psia in Odeh's paper
+
+WCONINJE
+-- Item #:1	 2	 3	 4	5      6  7
+	'INJ'	'GAS'	'OPEN'	'RATE'	100000 1* 9014 /
+/
+-- Stated in Odeh that gas inj. rate (item 5) is 100MMscf per day
+-- BHP upper limit (item 7) should not be exceeding the highest
+-- pressure in the PVT table=9014.7psia (default is 100 000psia)
+
+TSTEP
+--Advance the simulater once a month for TEN years:
+31 28 31 30 31 30 31 31 30 31 30 31 
+31 28 31 30 31 30 31 31 30 31 30 31 
+31 28 31 30 31 30 31 31 30 31 30 31 
+31 28 31 30 31 30 31 31 30 31 30 31 
+31 28 31 30 31 30 31 31 30 31 30 31 
+31 28 31 30 31 30 31 31 30 31 30 31 
+31 28 31 30 31 30 31 31 30 31 30 31 
+31 28 31 30 31 30 31 31 30 31 30 31 
+31 28 31 30 31 30 31 31 30 31 30 31 
+31 28 31 30 31 30 31 31 30 31 30 31 /
+
+--Advance the simulator once a year for TEN years:
+--10*365 /
+
+END

python/tests/test_parse_deck.py

@@ -3,8 +3,15 @@ import json
 import opm
 import opm.io
 import os.path
+import numpy as np
+
+try:
+    from tests.utils import test_path
+except ImportError:
+    from utils import test_path
 
 from opm.io.parser import Parser, ParseContext
+from opm.io.deck import DeckKeyword
 
 class TestParse(unittest.TestCase):
 
@@ -83,6 +90,65 @@ FIPNUM
         self.assertIn( 'TESTKEY2', deck )
 
 
+    def test_parser_deckItems(self):
+
+        parser = Parser()
+
+        error_recovery = [("PARSE_RANDOM_SLASH", opm.io.action.ignore),
+                          ("PARSE_EXTRA_RECORDS", opm.io.action.ignore)]
+
+        context = ParseContext(error_recovery)
+
+        self.deck_spe1case1 = parser.parse(test_path("data/SPE1CASE1.DATA"), context)
+
+        dkw_compdate = self.deck_spe1case1["COMPDAT"]
+
+        self.assertTrue( dkw_compdate[0][0].is_string() )
+        self.assertFalse( dkw_compdate[0][1].is_string() )
+
+        self.assertTrue( dkw_compdate[0][1].is_int() )
+        self.assertFalse( dkw_compdate[0][1].is_double() )
+
+        self.assertTrue( dkw_compdate[0][8].is_double() )
+
+        self.assertTrue(dkw_compdate[0][0].value == "PROD")
+
+        conI = dkw_compdate[0][1].value
+        conJ = dkw_compdate[0][2].value
+        conK = dkw_compdate[0][3].value
+
+        self.assertEqual(dkw_compdate[0][5].value, "OPEN")
+
+        self.assertTrue((conI, conJ, conK) == (10,10,3))
+
+        self.assertFalse( dkw_compdate[0][7].valid )
+        self.assertTrue( dkw_compdate[0][7].defaulted )
+
+        self.assertEqual( dkw_compdate[0][6].value, 0)
+
+        self.assertEqual( dkw_compdate[0][8].value, 0.5)
+
+        dkw_wconprod = self.deck_spe1case1["WCONPROD"]
+
+        welln= dkw_wconprod[0][0].value
+        self.assertEqual(dkw_wconprod[0][2].value, "ORAT")
+        self.assertEqual(dkw_wconprod[0][3].value, "WUOPRL")
+        self.assertEqual(dkw_wconprod[0][5].value, 1.5e5)
+
+        dkw_permx = self.deck_spe1case1["PERMX"]
+        permx =  dkw_permx.get_raw_array()
+        self.assertEqual(len(permx), 300)
+        self.assertTrue(isinstance(permx, np.ndarray))
+        self.assertEqual(permx.dtype, "float64")
+
+        dkw_eqlnum = self.deck_spe1case1["EQLNUM"]
+        eqlnum =  dkw_eqlnum.get_int_array()
+
+        self.assertEqual(len(eqlnum), 300)
+        self.assertTrue(isinstance(eqlnum, np.ndarray))
+        self.assertEqual(eqlnum.dtype, "int32")
+
+
 
 if __name__ == "__main__":
     unittest.main()