Merge pull request #502 from OPM/opm-parser-integrate

Integrate new parser

CMakeLists_files.cmake

@@ -165,10 +165,13 @@ list (APPEND TEST_SOURCE_FILES
 	tests/test_param.cpp
 	tests/test_blackoilfluid.cpp
 	tests/test_shadow.cpp
-        tests/test_units.cpp
+	tests/test_units.cpp
-        tests/test_blackoilstate.cpp
+	tests/test_blackoilstate.cpp
-        tests/test_wellsmanager.cpp
+	tests/test_parser.cpp
-        tests/test_wellcontrols.cpp
+	tests/test_wellsmanager.cpp
+	tests/test_wellcontrols.cpp
+	tests/test_wellsgroup.cpp
+	tests/test_wellcollection.cpp
   )
 
 # originally generated with the command:
@@ -177,10 +180,13 @@ list (APPEND TEST_DATA_FILES
 	tests/extratestdata.xml
 	tests/testdata.xml
 	tests/liveoil.DATA
-  tests/wetgas.DATA
+	tests/wetgas.DATA
-  tests/testBlackoilState1.DATA
+	tests/testBlackoilState1.DATA
-  tests/testBlackoilState2.DATA
+	tests/testBlackoilState2.DATA
-  tests/wells_manager_data.data
+	tests/wells_manager_data.data
+	tests/wells_manager_data_expanded.data
+	tests/wells_manager_data_wellSTOP.data
+	tests/wells_group.data
 	)
 
 # originally generated with the command:

cmake/Modules/Findcjson.cmake

@@ -31,13 +31,14 @@ if (CMAKE_SIZEOF_VOID_P)
   math (EXPR _BITS "8 * ${CMAKE_SIZEOF_VOID_P}")
 endif (CMAKE_SIZEOF_VOID_P)
 
+string(REGEX REPLACE "${PROJECT_SOURCE_DIR}/?(.*)" "\\1"  BUILD_DIR_SUFFIX "${PROJECT_BINARY_DIR}")
+
 find_library (CJSON_LIBRARY
   NAMES "cjson"
   HINTS "${CJSON_ROOT}"
   PATHS "${PROJECT_BINARY_DIR}/../opm-parser"
-        "${PROJECT_BINARY_DIR}/../opm-parser-build"
+        "${PROJECT_BINARY_DIR}/../opm-parser${BUILD_DIR_SUFFIX}"
-        "${PROJECT_BINARY_DIR}/../../opm-parser/build"
+        "${PROJECT_BINARY_DIR}/../../opm-parser/${BUILD_DIR_SUFFIX}"
-        "${PROJECT_BINARY_DIR}/../../opm-parser/cmake-build"
   PATH_SUFFIXES "lib" "lib${_BITS}" "lib/${CMAKE_LIBRARY_ARCHITECTURE}"
                 "opm/json"
   DOC "Path to cjson library archive/shared object files"

cmake/Modules/Findopm-parser.cmake

@@ -5,9 +5,9 @@
 #
 # If found, it sets these variables:
 #
-#	HAVE_OPM_PARSER              Defined if a test program compiled
+#       HAVE_OPM_PARSER              Defined if a test program compiled
-#	OPM_PARSER_INCLUDE_DIRS      Header file directories
+#       OPM_PARSER_INCLUDE_DIRS      Header file directories
-#	OPM_PARSER_LIBRARIES         Archives and shared objects
+#       OPM_PARSER_LIBRARIES         Archives and shared objects
 
 include (FindPackageHandleStandardArgs)
 
@@ -31,6 +31,9 @@ if ((NOT OPM_PARSER_ROOT) AND OPM_ROOT)
   set (OPM_PARSER_ROOT "${OPM_ROOT}/opm-parser")
 endif ()
 
+# Detect the build dir suffix or subdirectory
+string(REGEX REPLACE "${PROJECT_SOURCE_DIR}/?(.*)" "\\1"  BUILD_DIR_SUFFIX "${PROJECT_BINARY_DIR}")
+
 # if a root is specified, then don't search in system directories
 # or in relative directories to this one
 if (OPM_PARSER_ROOT)
@@ -43,9 +46,8 @@ else ()
     "${PROJECT_SOURCE_DIR}/../opm-parser")
   set (_opm_parser_build
     "${PROJECT_BINARY_DIR}/../opm-parser"
-    "${PROJECT_BINARY_DIR}/../opm-parser-build"
+    "${PROJECT_BINARY_DIR}/../opm-parser${BUILD_DIR_SUFFIX}"
-    "${PROJECT_BINARY_DIR}/../../opm-parser/build"
+    "${PROJECT_BINARY_DIR}/../../opm-parser/${BUILD_DIR_SUFFIX}")
-    "${PROJECT_BINARY_DIR}/../../opm-parser/cmake-build")
 endif ()
 
 # use this header as signature
@@ -61,18 +63,18 @@ find_path (OPM_PARSER_INCLUDE_DIR
 # then it is probably a build directory; read the CMake cache of
 # opm-parser to figure out where the source directory is
 if ((NOT OPM_PARSER_INCLUDE_DIR) AND
-	(OPM_PARSER_ROOT AND (EXISTS "${OPM_PARSER_ROOT}/CMakeCache.txt")))
+        (OPM_PARSER_ROOT AND (EXISTS "${OPM_PARSER_ROOT}/CMakeCache.txt")))
   set (_regex "^OPMParser_SOURCE_DIR:STATIC=\(.*\)$")
   file (STRINGS
-	"${OPM_PARSER_ROOT}/CMakeCache.txt"
+        "${OPM_PARSER_ROOT}/CMakeCache.txt"
-	_cache_entry
+        _cache_entry
-	REGEX "${_regex}")
+        REGEX "${_regex}")
   string(REGEX REPLACE "${_regex}" "\\1"
-	OPM_PARSER_INCLUDE_DIR
+        OPM_PARSER_INCLUDE_DIR
-	"${_cache_entry}")
+        "${_cache_entry}")
   if (OPM_PARSER_INCLUDE_DIR)
-	set (OPM_PARSER_INCLUDE_DIR "${OPM_PARSER_INCLUDE_DIR}"
+        set (OPM_PARSER_INCLUDE_DIR "${OPM_PARSER_INCLUDE_DIR}"
-	  CACHE PATH "Path to OPM parser header files" FORCE)
+          CACHE PATH "Path to OPM parser header files" FORCE)
   endif ()
 endif ()
 
@@ -109,14 +111,16 @@ endif ()
 # get the prerequisite Boost libraries
 if (NOT Boost_FOUND)
   find_package(Boost 1.44.0
-	COMPONENTS filesystem date_time system unit_test_framework REQUIRED ${OPM_PARSER_QUIET})
+        COMPONENTS filesystem date_time system unit_test_framework REQUIRED ${OPM_PARSER_QUIET})
 endif ()
 
 # setup list of all required libraries to link with opm-parser. notice that
 # we use the plural form to get *all* the libraries needed by cjson
 set (OPM_PARSER_INCLUDE_DIRS
   ${OPM_PARSER_INCLUDE_DIR}
+  ${CJSON_INCLUDE_DIRS}
   ${Boost_INCLUDE_DIRS})
+
 set (OPM_PARSER_LIBRARIES
   ${OPM_PARSER_LIBRARY}
   ${OPM_JSON_LIBRARY}
@@ -126,7 +130,7 @@ set (OPM_PARSER_LIBRARIES
 # see if we can compile a minimum example
 # CMake logical test doesn't handle lists (sic)
 if (NOT (OPM_PARSER_INCLUDE_DIR MATCHES "-NOTFOUND"
-	  OR OPM_PARSER_LIBRARIES MATCHES "-NOTFOUND"))
+          OR OPM_PARSER_LIBRARIES MATCHES "-NOTFOUND"))
   include (CMakePushCheckState)
   include (CheckCSourceCompiles)
   cmake_push_check_state ()

cmake/Modules/OpmFind.cmake

@@ -59,6 +59,7 @@ set (_opm_proj_exemptions
   dune-istl
   dune-grid
   dune-geometry
+  opm-parser
   )
 
 # although a DUNE module, it is delivered in the OPM suite

cmake/Modules/opm-core-prereqs.cmake

@@ -3,31 +3,33 @@
 
 # defines that must be present in config.h for our headers
 set (opm-core_CONFIG_VAR
-	HAVE_ERT
+        HAVE_ERT
-	HAVE_SUITESPARSE_UMFPACK_H
+        HAVE_SUITESPARSE_UMFPACK_H
-	)
+        )
 
 # dependencies
 set (opm-core_DEPS
-	# compile with C99 support if available
+        # compile with C99 support if available
-	"C99"
+        "C99"
-	# compile with C++0x/11 support if available
+        # compile with C++0x/11 support if available
-	"CXX11Features REQUIRED"
+        "CXX11Features REQUIRED"
-	# various runtime library enhancements
+        # various runtime library enhancements
-	"Boost 1.39.0
+        "Boost 1.44.0
-		COMPONENTS date_time filesystem system unit_test_framework REQUIRED"
+                COMPONENTS date_time filesystem system unit_test_framework REQUIRED"
-	# matrix library
+        # matrix library
-	"BLAS REQUIRED"
+        "BLAS REQUIRED"
-	"LAPACK REQUIRED"
+        "LAPACK REQUIRED"
-	# Tim Davis' SuiteSparse archive
+        # Tim Davis' SuiteSparse archive
-	"SuiteSparse COMPONENTS umfpack"
+        "SuiteSparse COMPONENTS umfpack"
-	# solver
+        # solver
-	"SuperLU"
+        "SuperLU"
-	# xml processing (for config parsing)
+        # xml processing (for config parsing)
-	"TinyXML"
+        "TinyXML"
-	# Ensembles-based Reservoir Tools (ERT)
+        #Parser library
-	"ERT"
+        "opm-parser REQUIRED"
-	# DUNE dependency
+        # Ensembles-based Reservoir Tools (ERT)
-	"dune-common"
+        "ERT"
-	"dune-istl"
+        # DUNE dependency
-	)
+        "dune-common"
+        "dune-istl"
+        )

cmake/Modules/opm-parser-prereqs.cmake

@@ -13,6 +13,6 @@ set (opm-parser_DEPS
         # compile with C++0x/11 support if available
         "CXX11Features REQUIRED"
         # various runtime library enhancements
-        "Boost 1.39.0
+        "Boost 1.44.0 COMPONENTS date_time filesystem system unit_test_framework REQUIRED"
-                COMPONENTS date_time filesystem system unit_test_framework REQUIRED"
+        "cJSON"
      )

dune.module

@@ -4,4 +4,4 @@ Description: Open Porous Media Initiative Core Library
 Version: 1.1
 Label: 2013.10
 Maintainer: atgeirr@sintef.no
-Depends: dune-common (>= 2.2) dune-istl (>= 2.2)
+Depends: dune-common (>= 2.2) dune-istl (>= 2.2) opm-parser

opm/core/grid/GridManager.cpp

@@ -23,16 +23,13 @@
 #include <opm/core/grid.h>
 #include <opm/core/grid/cart_grid.h>
 #include <opm/core/grid/cornerpoint_grid.h>
+
+#include <array>
 #include <algorithm>
 #include <numeric>
 
-
-
 namespace Opm
 {
-
-
-
     /// Construct a 3d corner-point grid from a deck.
     GridManager::GridManager(const Opm::EclipseGridParser& deck)
     {
@@ -58,6 +55,30 @@ namespace Opm
         }
     }
 
+    /// Construct a 3d corner-point grid from a deck.
+    GridManager::GridManager(Opm::DeckConstPtr newParserDeck)
+    {
+        // We accept two different ways to specify the grid.
+        //    1. Corner point format.
+        //       Requires ZCORN, COORDS, DIMENS or SPECGRID, optionally
+        //       ACTNUM, optionally MAPAXES.
+        //       For this format, we will verify that DXV, DYV, DZV,
+        //       DEPTHZ and TOPS are not present.
+        //    2. Tensor grid format.
+        //       Requires DXV, DYV, DZV, optionally DEPTHZ or TOPS.
+        //       For this format, we will verify that ZCORN, COORDS
+        //       and ACTNUM are not present.
+        //       Note that for TOPS, we only allow a uniform vector of values.
+
+        if (newParserDeck->hasKeyword("ZCORN") && newParserDeck->hasKeyword("COORD")) {
+            initFromDeckCornerpoint(newParserDeck);
+        } else if (newParserDeck->hasKeyword("DXV") && newParserDeck->hasKeyword("DYV") && newParserDeck->hasKeyword("DZV")) {
+            initFromDeckTensorgrid(newParserDeck);
+        } else {
+            OPM_THROW(std::runtime_error, "Could not initialize grid from deck. "
+                  "Need either ZCORN + COORD or DXV + DYV + DZV keywords.");
+        }
+    }
 
 
 
@@ -152,6 +173,71 @@ namespace Opm
         }
     }
 
+    // Construct corner-point grid from deck.
+    void GridManager::initFromDeckCornerpoint(Opm::DeckConstPtr newParserDeck)
+    {
+        // Extract data from deck.
+        // Collect in input struct for preprocessing.
+        struct grdecl grdecl;
+        createGrdecl(newParserDeck, grdecl);
+
+        // Process grid.
+        ug_ = create_grid_cornerpoint(&grdecl, 0.0);
+        if (!ug_) {
+            OPM_THROW(std::runtime_error, "Failed to construct grid.");
+        }
+    }
+
+    void GridManager::createGrdecl(Opm::DeckConstPtr newParserDeck, struct grdecl &grdecl)
+    {
+        // Extract data from deck.
+        const std::vector<double>& zcorn = newParserDeck->getKeyword("ZCORN")->getSIDoubleData();
+        const std::vector<double>& coord = newParserDeck->getKeyword("COORD")->getSIDoubleData();
+        const int* actnum = NULL;
+        if (newParserDeck->hasKeyword("ACTNUM")) {
+            actnum = &(newParserDeck->getKeyword("ACTNUM")->getIntData()[0]);
+        }
+
+        std::array<int, 3> dims;
+        if (newParserDeck->hasKeyword("DIMENS")) {
+            Opm::DeckKeywordConstPtr dimensKeyword = newParserDeck->getKeyword("DIMENS");
+            dims[0] = dimensKeyword->getRecord(0)->getItem(0)->getInt(0);
+            dims[1] = dimensKeyword->getRecord(0)->getItem(1)->getInt(0);
+            dims[2] = dimensKeyword->getRecord(0)->getItem(2)->getInt(0);
+        } else if (newParserDeck->hasKeyword("SPECGRID")) {
+            Opm::DeckKeywordConstPtr specgridKeyword = newParserDeck->getKeyword("SPECGRID");
+            dims[0] = specgridKeyword->getRecord(0)->getItem(0)->getInt(0);
+            dims[1] = specgridKeyword->getRecord(0)->getItem(1)->getInt(0);
+            dims[2] = specgridKeyword->getRecord(0)->getItem(2)->getInt(0);
+        } else {
+            OPM_THROW(std::runtime_error, "Deck must have either DIMENS or SPECGRID.");
+        }
+
+        // Collect in input struct for preprocessing.
+
+        grdecl.zcorn = &zcorn[0];
+        grdecl.coord = &coord[0];
+        grdecl.actnum = actnum;
+        grdecl.dims[0] = dims[0];
+        grdecl.dims[1] = dims[1];
+        grdecl.dims[2] = dims[2];
+
+        if (newParserDeck->hasKeyword("MAPAXES")) {
+            Opm::DeckKeywordConstPtr mapaxesKeyword = newParserDeck->getKeyword("MAPAXES");
+            Opm::DeckRecordConstPtr mapaxesRecord = mapaxesKeyword->getRecord(0);
+
+            // memleak alert: here we need to make sure that C code
+            // can properly take ownership of the grdecl.mapaxes
+            // object. if it is not freed, it will result in a
+            // memleak...
+            double *cWtfMapaxes = static_cast<double*>(malloc(sizeof(double)*mapaxesRecord->size()));
+            for (unsigned i = 0; i < mapaxesRecord->size(); ++i)
+                cWtfMapaxes[i] = mapaxesRecord->getItem(i)->getSIDouble(0);
+            grdecl.mapaxes = cWtfMapaxes;
+        } else
+            grdecl.mapaxes = NULL;
+
+    }
 
     namespace
     {
@@ -230,6 +316,71 @@ namespace Opm
         }
     }
 
+    void GridManager::initFromDeckTensorgrid(Opm::DeckConstPtr newParserDeck)
+    {
+        // Extract logical cartesian size.
+        std::array<int, 3> dims;
+        if (newParserDeck->hasKeyword("DIMENS")) {
+            Opm::DeckKeywordConstPtr dimensKeyword = newParserDeck->getKeyword("DIMENS");
+            dims[0] = dimensKeyword->getRecord(0)->getItem(0)->getInt(0);
+            dims[1] = dimensKeyword->getRecord(0)->getItem(1)->getInt(0);
+            dims[2] = dimensKeyword->getRecord(0)->getItem(2)->getInt(0);
+        } else if (newParserDeck->hasKeyword("SPECGRID")) {
+            Opm::DeckKeywordConstPtr specgridKeyword = newParserDeck->getKeyword("SPECGRID");
+            dims[0] = specgridKeyword->getRecord(0)->getItem(0)->getInt(0);
+            dims[1] = specgridKeyword->getRecord(0)->getItem(1)->getInt(0);
+            dims[2] = specgridKeyword->getRecord(0)->getItem(2)->getInt(0);
+        } else {
+            OPM_THROW(std::runtime_error, "Deck must have either DIMENS or SPECGRID.");
+        }
+
+        // Extract coordinates (or offsets from top, in case of z).
+        const std::vector<double>& dxv = newParserDeck->getKeyword("DXV")->getSIDoubleData();
+        const std::vector<double>& dyv = newParserDeck->getKeyword("DYV")->getSIDoubleData();
+        const std::vector<double>& dzv = newParserDeck->getKeyword("DZV")->getSIDoubleData();
+        std::vector<double> x = coordsFromDeltas(dxv);
+        std::vector<double> y = coordsFromDeltas(dyv);
+        std::vector<double> z = coordsFromDeltas(dzv);
+
+        // Check that number of cells given are consistent with DIMENS/SPECGRID.
+        if (dims[0] != int(dxv.size())) {
+            OPM_THROW(std::runtime_error, "Number of DXV data points do not match DIMENS or SPECGRID.");
+        }
+        if (dims[1] != int(dyv.size())) {
+            OPM_THROW(std::runtime_error, "Number of DYV data points do not match DIMENS or SPECGRID.");
+        }
+        if (dims[2] != int(dzv.size())) {
+            OPM_THROW(std::runtime_error, "Number of DZV data points do not match DIMENS or SPECGRID.");
+        }
+
+        // Extract top corner depths, if available.
+        const double* top_depths = 0;
+        std::vector<double> top_depths_vec;
+        if (newParserDeck->hasKeyword("DEPTHZ")) {
+            const std::vector<double>& depthz = newParserDeck->getKeyword("DEPTHZ")->getSIDoubleData();
+            if (depthz.size() != x.size()*y.size()) {
+                OPM_THROW(std::runtime_error, "Incorrect size of DEPTHZ: " << depthz.size());
+            }
+            top_depths = &depthz[0];
+        } else if (newParserDeck->hasKeyword("TOPS")) {
+            // We only support constant values for TOPS.
+            // It is not 100% clear how we best can deal with
+            // varying TOPS (stair-stepping grid, or not).
+            const std::vector<double>& tops = newParserDeck->getKeyword("TOPS")->getSIDoubleData();
+            if (std::count(tops.begin(), tops.end(), tops[0]) != int(tops.size())) {
+                OPM_THROW(std::runtime_error, "We do not support nonuniform TOPS, please use ZCORN/COORDS instead.");
+            }
+            top_depths_vec.resize(x.size()*y.size(), tops[0]);
+            top_depths = &top_depths_vec[0];
+        }
+
+        // Construct grid.
+        ug_ = create_grid_tensor3d(dxv.size(), dyv.size(), dzv.size(),
+                                   &x[0], &y[0], &z[0], top_depths);
+        if (!ug_) {
+            OPM_THROW(std::runtime_error, "Failed to construct grid.");
+        }
+    }
 
 
 } // namespace Opm

opm/core/grid/GridManager.hpp

@@ -20,10 +20,12 @@
 #ifndef OPM_GRIDMANAGER_HEADER_INCLUDED
 #define OPM_GRIDMANAGER_HEADER_INCLUDED
 
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+
 #include <string>
 
 struct UnstructuredGrid;
-
+struct grdecl;
 
 namespace Opm
 {
@@ -44,6 +46,9 @@ namespace Opm
         /// Construct a 3d corner-point grid or tensor grid from a deck.
         GridManager(const Opm::EclipseGridParser& deck);
 
+        /// Construct a 3d corner-point grid or tensor grid from a deck.
+        GridManager(Opm::DeckConstPtr newParserDeck);
+
         /// Construct a 2d cartesian grid with cells of unit size.
         GridManager(int nx, int ny);
 
@@ -72,6 +77,8 @@ namespace Opm
         /// to make it clear that we are returning a C-compatible struct.
         const UnstructuredGrid* c_grid() const;
 
+        static void createGrdecl(Opm::DeckConstPtr newParserDeck, struct grdecl &grdecl);
+
     private:
         // Disable copying and assignment.
         GridManager(const GridManager& other);
@@ -79,8 +86,10 @@ namespace Opm
 
         // Construct corner-point grid from deck.
         void initFromDeckCornerpoint(const Opm::EclipseGridParser& deck);
+        void initFromDeckCornerpoint(Opm::DeckConstPtr newParserDeck);
         // Construct tensor grid from deck.
         void initFromDeckTensorgrid(const Opm::EclipseGridParser& deck);
+        void initFromDeckTensorgrid(Opm::DeckConstPtr newParserDeck);
 
         // The managed UnstructuredGrid.
         UnstructuredGrid* ug_;

opm/core/props/BlackoilPropertiesFromDeck.cpp

@@ -23,7 +23,6 @@
 
 namespace Opm
 {
-
     BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
                                                            const UnstructuredGrid& grid,
                                                            bool init_rock)
@@ -43,6 +42,25 @@ namespace Opm
         }
     }
 
+    BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck,
+                                                           const UnstructuredGrid& grid,
+                                                           bool init_rock)
+    {
+        if (init_rock){
+           rock_.init(newParserDeck, grid);
+        }
+        pvt_.init(newParserDeck, /*numSamples=*/200);
+        SaturationPropsFromDeck<SatFuncSimpleUniform>* ptr
+            = new SaturationPropsFromDeck<SatFuncSimpleUniform>();
+        satprops_.reset(ptr);
+        ptr->init(newParserDeck, grid, /*numSamples=*/200);
+
+        if (pvt_.numPhases() != satprops_->numPhases()) {
+            OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
+                  << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ").");
+        }
+    }
+
     BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
                                                            const UnstructuredGrid& grid,
                                                            const parameter::ParameterGroup& param,
@@ -107,6 +125,70 @@ namespace Opm
         }
     }
 
+    BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck,
+                                                           const UnstructuredGrid& grid,
+                                                           const parameter::ParameterGroup& param,
+                                                           bool init_rock)
+    {
+        if(init_rock){
+            rock_.init(newParserDeck, grid);
+        }
+
+        const int pvt_samples = param.getDefault("pvt_tab_size", 200);
+        pvt_.init(newParserDeck, pvt_samples);
+
+        // Unfortunate lack of pointer smartness here...
+        const int sat_samples = param.getDefault("sat_tab_size", 200);
+        std::string threephase_model = param.getDefault<std::string>("threephase_model", "simple");
+        if (newParserDeck->hasKeyword("ENDSCALE") && threephase_model != "simple") {
+            OPM_THROW(std::runtime_error, "Sorry, end point scaling currently available for the 'simple' model only.");
+        }
+        if (sat_samples > 1) {
+            if (threephase_model == "stone2") {
+                SaturationPropsFromDeck<SatFuncStone2Uniform>* ptr
+                    = new SaturationPropsFromDeck<SatFuncStone2Uniform>();
+                satprops_.reset(ptr);
+                ptr->init(newParserDeck, grid, sat_samples);
+            } else if (threephase_model == "simple") {
+                SaturationPropsFromDeck<SatFuncSimpleUniform>* ptr
+                    = new SaturationPropsFromDeck<SatFuncSimpleUniform>();
+                satprops_.reset(ptr);
+                ptr->init(newParserDeck, grid, sat_samples);
+            } else if (threephase_model == "gwseg") {
+                SaturationPropsFromDeck<SatFuncGwsegUniform>* ptr
+                    = new SaturationPropsFromDeck<SatFuncGwsegUniform>();
+                satprops_.reset(ptr);
+                ptr->init(newParserDeck, grid, sat_samples);
+            } else {
+                OPM_THROW(std::runtime_error, "Unknown threephase_model: " << threephase_model);
+            }
+        } else {
+            if (threephase_model == "stone2") {
+                SaturationPropsFromDeck<SatFuncStone2Nonuniform>* ptr
+                    = new SaturationPropsFromDeck<SatFuncStone2Nonuniform>();
+                satprops_.reset(ptr);
+                ptr->init(newParserDeck, grid, sat_samples);
+            } else if (threephase_model == "simple") {
+                SaturationPropsFromDeck<SatFuncSimpleNonuniform>* ptr
+                    = new SaturationPropsFromDeck<SatFuncSimpleNonuniform>();
+                satprops_.reset(ptr);
+                ptr->init(newParserDeck, grid, sat_samples);
+            } else if (threephase_model == "gwseg") {
+                SaturationPropsFromDeck<SatFuncGwsegNonuniform>* ptr
+                    = new SaturationPropsFromDeck<SatFuncGwsegNonuniform>();
+                satprops_.reset(ptr);
+                ptr->init(newParserDeck, grid, sat_samples);
+            } else {
+                OPM_THROW(std::runtime_error, "Unknown threephase_model: " << threephase_model);
+            }
+        }
+
+        if (pvt_.numPhases() != satprops_->numPhases()) {
+            OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
+                  << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ").");
+        }
+    }
+
     BlackoilPropertiesFromDeck::~BlackoilPropertiesFromDeck()
     {
     }

opm/core/props/BlackoilPropertiesFromDeck.hpp

@@ -27,6 +27,9 @@
 #include <opm/core/props/satfunc/SaturationPropsFromDeck.hpp>
 #include <opm/core/io/eclipse/EclipseGridParser.hpp>
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
+
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+
 #include <memory>
 
 struct UnstructuredGrid;
@@ -44,7 +47,15 @@ namespace Opm
         /// \param[in]  grid     Grid to which property object applies, needed for the
         ///                      mapping from cell indices (typically from a processed grid)
         ///                      to logical cartesian indices consistent with the deck.
-        BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
+        BlackoilPropertiesFromDeck(const EclipseGridParser &deck,
+                                   const UnstructuredGrid& grid, bool init_rock=true );
+
+        /// Initialize from deck and grid.
+        /// \param[in]  deck     Deck input parser
+        /// \param[in]  grid     Grid to which property object applies, needed for the
+        ///                      mapping from cell indices (typically from a processed grid)
+        ///                      to logical cartesian indices consistent with the deck.
+        BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck,
                                    const UnstructuredGrid& grid, bool init_rock=true );
 
         /// Initialize from deck, grid and parameters.
@@ -63,6 +74,22 @@ namespace Opm
                                    const parameter::ParameterGroup& param,
                                    bool init_rock=true);
 
+        /// Initialize from deck, grid and parameters.
+        /// \param[in]  deck     Deck input parser
+        /// \param[in]  grid     Grid to which property object applies, needed for the
+        ///                      mapping from cell indices (typically from a processed grid)
+        ///                      to logical cartesian indices consistent with the deck.
+        /// \param[in]  param    Parameters. Accepted parameters include:
+        ///                        pvt_tab_size (200)          number of uniform sample points for dead-oil pvt tables.
+        ///                        sat_tab_size (200)          number of uniform sample points for saturation tables.
+        ///                        threephase_model("simple")  three-phase relperm model (accepts "simple" and "stone2").
+        ///                      For both size parameters, a 0 or negative value indicates that no spline fitting is to
+        ///                      be done, and the input fluid data used directly for linear interpolation.
+        BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck,
+                                   const UnstructuredGrid& grid,
+                                   const parameter::ParameterGroup& param,
+                                   bool init_rock=true);
+
         /// Destructor.
         virtual ~BlackoilPropertiesFromDeck();
 

opm/core/props/phaseUsageFromDeck.hpp

@@ -24,10 +24,49 @@
 #include <opm/core/io/eclipse/EclipseGridParser.hpp>
 #include <opm/core/props/BlackoilPhases.hpp>
 
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+
 
 namespace Opm
 {
 
+    /// Looks at presence of WATER, OIL and GAS keywords in state object
+    /// to determine active phases.
+    inline PhaseUsage phaseUsageFromDeck(Opm::EclipseStateConstPtr eclipseState)
+    {
+        PhaseUsage pu;
+        std::fill(pu.phase_used, pu.phase_used + BlackoilPhases::MaxNumPhases, 0);
+
+        // Discover phase usage.
+        if (eclipseState->hasPhase(Phase::PhaseEnum::WATER)) {
+            pu.phase_used[BlackoilPhases::Aqua] = 1;
+        }
+        if (eclipseState->hasPhase(Phase::PhaseEnum::OIL)) {
+            pu.phase_used[BlackoilPhases::Liquid] = 1;
+        }
+        if (eclipseState->hasPhase(Phase::PhaseEnum::GAS)) {
+            pu.phase_used[BlackoilPhases::Vapour] = 1;
+        }
+        pu.num_phases = 0;
+        for (int i = 0; i < BlackoilPhases::MaxNumPhases; ++i) {
+            pu.phase_pos[i] = pu.num_phases;
+            pu.num_phases += pu.phase_used[i];
+        }
+
+        // Only 2 or 3 phase systems handled.
+        if (pu.num_phases < 2 || pu.num_phases > 3) {
+            OPM_THROW(std::runtime_error, "Cannot handle cases with " << pu.num_phases << " phases.");
+        }
+
+        // We need oil systems, since we do not support the keywords needed for
+        // water-gas systems.
+        if (!pu.phase_used[BlackoilPhases::Liquid]) {
+            OPM_THROW(std::runtime_error, "Cannot handle cases with no OIL, i.e. water-gas systems.");
+        }
+
+        return pu;
+    }
 
     /// Looks at presence of WATER, OIL and GAS keywords in deck
     /// to determine active phases.
@@ -66,6 +105,43 @@ namespace Opm
         return pu;
     }
 
+    /// Looks at presence of WATER, OIL and GAS keywords in deck
+    /// to determine active phases.
+    inline PhaseUsage phaseUsageFromDeck(Opm::DeckConstPtr newParserDeck)
+    {
+        PhaseUsage pu;
+        std::fill(pu.phase_used, pu.phase_used + BlackoilPhases::MaxNumPhases, 0);
+
+        // Discover phase usage.
+        if (newParserDeck->hasKeyword("WATER")) {
+            pu.phase_used[BlackoilPhases::Aqua] = 1;
+        }
+        if (newParserDeck->hasKeyword("OIL")) {
+            pu.phase_used[BlackoilPhases::Liquid] = 1;
+        }
+        if (newParserDeck->hasKeyword("GAS")) {
+            pu.phase_used[BlackoilPhases::Vapour] = 1;
+        }
+        pu.num_phases = 0;
+        for (int i = 0; i < BlackoilPhases::MaxNumPhases; ++i) {
+            pu.phase_pos[i] = pu.num_phases;
+            pu.num_phases += pu.phase_used[i];
+        }
+
+        // Only 2 or 3 phase systems handled.
+        if (pu.num_phases < 2 || pu.num_phases > 3) {
+            OPM_THROW(std::runtime_error, "Cannot handle cases with " << pu.num_phases << " phases.");
+        }
+
+        // We need oil systems, since we do not support the keywords needed for
+        // water-gas systems.
+        if (!pu.phase_used[BlackoilPhases::Liquid]) {
+            OPM_THROW(std::runtime_error, "Cannot handle cases with no OIL, i.e. water-gas systems.");
+        }
+
+        return pu;
+    }
+
 }
 
 #endif // OPM_PHASEUSAGEFROMDECK_HEADER_INCLUDED

opm/core/props/pvt/BlackoilPvtProperties.cpp

@@ -32,6 +32,9 @@
 #include <opm/core/utility/ErrorMacros.hpp>
 #include <opm/core/utility/linearInterpolation.hpp>
 
+#include <opm/parser/eclipse/Utility/PvtwTable.hpp>
+#include <opm/parser/eclipse/Utility/PvdcoTable.hpp>
+#include <opm/parser/eclipse/Deck/Deck.hpp>
 
 namespace Opm
 {
@@ -40,7 +43,6 @@ namespace Opm
     {
     }
 
-
     void BlackoilPvtProperties::init(const EclipseGridParser& deck, const int samples)
     {
         // If we need multiple regions, this class and the SinglePvt* classes must change.
@@ -113,6 +115,83 @@ namespace Opm
         }
     }
 
+    void BlackoilPvtProperties::init(Opm::DeckConstPtr newParserDeck, int samples)
+    {
+        // If we need multiple regions, this class and the SinglePvt* classes must change.
+        region_number_ = 0;
+
+        phase_usage_ = phaseUsageFromDeck(newParserDeck);
+
+        // Surface densities. Accounting for different orders in eclipse and our code.
+        if (newParserDeck->hasKeyword("DENSITY")) {
+            Opm::DeckKeywordConstPtr densityKeyword = newParserDeck->getKeyword("DENSITY");
+            const std::vector<double>& d = densityKeyword->getRecord(region_number_)->getItem(0)->getSIDoubleData();
+            enum { ECL_oil = 0, ECL_water = 1, ECL_gas = 2 };
+            if (phase_usage_.phase_used[Aqua]) {
+                densities_[phase_usage_.phase_pos[Aqua]]   = d[ECL_water];
+            }
+            if (phase_usage_.phase_used[Vapour]) {
+                densities_[phase_usage_.phase_pos[Vapour]] = d[ECL_gas];
+            }
+            if (phase_usage_.phase_used[Liquid]) {
+                densities_[phase_usage_.phase_pos[Liquid]] = d[ECL_oil];
+            }
+        } else {
+            OPM_THROW(std::runtime_error, "Input is missing DENSITY\n");
+        }
+
+        // Set the properties.
+        props_.resize(phase_usage_.num_phases);
+        // Water PVT
+        if (phase_usage_.phase_used[Aqua]) {
+            if (newParserDeck->hasKeyword("PVTW")) {
+                Opm::PvtwTable pvtwTable(newParserDeck->getKeyword("PVTW"));
+
+                props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(pvtwTable));
+            } else {
+                // Eclipse 100 default.
+                props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(0.5*Opm::prefix::centi*Opm::unit::Poise));
+            }
+        }
+        // Oil PVT
+        if (phase_usage_.phase_used[Liquid]) {
+            if (newParserDeck->hasKeyword("PVDO")) {
+                Opm::PvdoTable pvdoTable(newParserDeck->getKeyword("PVDO"), region_number_);
+
+                props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDeadSpline(pvdoTable, samples));
+            } else if (newParserDeck->hasKeyword("PVTO")) {
+                Opm::PvtoTable pvtoTable(newParserDeck->getKeyword("PVTO"), /*tableIdx=*/0);
+
+                props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtLiveOil(pvtoTable));
+            } else if (newParserDeck->hasKeyword("PVCDO")) {
+                Opm::PvdcoTable pvcdoTable(newParserDeck->getKeyword("PVCDO"));
+
+                props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtConstCompr(pvcdoTable));
+            } else {
+                OPM_THROW(std::runtime_error, "Input is missing PVDO or PVTO\n");
+            }
+        }
+        // Gas PVT
+        if (phase_usage_.phase_used[Vapour]) {
+            if (newParserDeck->hasKeyword("PVDG")) {
+                Opm::PvdgTable pvdgTable(newParserDeck->getKeyword("PVDG"), region_number_);
+
+                props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDeadSpline(pvdgTable, samples));
+            } else if (newParserDeck->hasKeyword("PVTG")) {
+                Opm::PvtgTable pvtgTable(newParserDeck->getKeyword("PVTG"), /*tableIdx=*/0);
+
+                props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtLiveGas(pvtgTable));
+            } else {
+                OPM_THROW(std::runtime_error, "Input is missing PVDG or PVTG\n");
+            }
+        }
+
+        // Must inform pvt property objects of phase structure.
+        for (int i = 0; i < phase_usage_.num_phases; ++i) {
+            props_[i]->setPhaseConfiguration(phase_usage_.num_phases, phase_usage_.phase_pos);
+        }
+    }
+
     const double* BlackoilPvtProperties::surfaceDensities() const
     {
         return densities_;

opm/core/props/pvt/BlackoilPvtProperties.hpp

@@ -20,11 +20,12 @@
 #ifndef OPM_BLACKOILPVTPROPERTIES_HEADER_INCLUDED
 #define OPM_BLACKOILPVTPROPERTIES_HEADER_INCLUDED
 
-
-
 #include <opm/core/props/pvt/SinglePvtInterface.hpp>
 #include <opm/core/props/BlackoilPhases.hpp>
 #include <opm/core/io/eclipse/EclipseGridParser.hpp>
+
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+
 #include <string>
 #include <memory>
 
@@ -55,6 +56,11 @@ namespace Opm
         ///                 data without fitting a spline.
         void init(const EclipseGridParser& deck, const int samples);
 
+
+        /// Initialize from deck.
+        /// \param deck     An input deck from the opm-parser module.
+        void init(Opm::DeckConstPtr deck, int samples);
+
         /// \return   Object describing the active phases.
         PhaseUsage phaseUsage() const;
 

opm/core/props/pvt/SinglePvtConstCompr.hpp

@@ -20,9 +20,12 @@
 #ifndef OPM_SINGLEPVTCONSTCOMPR_HEADER_INCLUDED
 #define OPM_SINGLEPVTCONSTCOMPR_HEADER_INCLUDED
 
-
 #include <opm/core/props/pvt/SinglePvtInterface.hpp>
 #include <opm/core/utility/ErrorMacros.hpp>
+
+#include <opm/parser/eclipse/Utility/PvtwTable.hpp>
+#include <opm/parser/eclipse/Utility/PvdcoTable.hpp>
+
 #include <vector>
 #include <algorithm>
 
@@ -55,6 +58,32 @@ namespace Opm
             visc_comp_ = pvtw[region_number][4];
         }
 
+        SinglePvtConstCompr(const Opm::PvtwTable &pvtwTable)
+        {
+            if (pvtwTable.numRows() != 1)
+                OPM_THROW(std::runtime_error,
+                          "The table specified by the PVTW keyword is required"
+                          "to exhibit exactly one row.");
+            ref_press_ = pvtwTable.getPressureColumn()[0];
+            ref_B_     = pvtwTable.getFormationFactorColumn()[0];
+            comp_      = pvtwTable.getCompressibilityColumn()[0];
+            viscosity_ = pvtwTable.getViscosityColumn()[0];
+            visc_comp_ = pvtwTable.getViscosibilityColumn()[0];
+        }
+
+        SinglePvtConstCompr(const Opm::PvdcoTable &pvdcoTable)
+        {
+            if (pvdcoTable.numRows() != 1)
+                OPM_THROW(std::runtime_error,
+                          "The table specified by the PVDCO keyword is required"
+                          "to exhibit exactly one row.");
+            ref_press_ = pvdcoTable.getPressureColumn()[0];
+            ref_B_     = pvdcoTable.getFormationFactorColumn()[0];
+            comp_      = pvdcoTable.getCompressibilityColumn()[0];
+            viscosity_ = pvdcoTable.getViscosityColumn()[0];
+            visc_comp_ = pvdcoTable.getViscosibilityColumn()[0];
+        }
+
         SinglePvtConstCompr(double visc)
             : ref_press_(0.0),
               ref_B_(1.0),

opm/core/props/pvt/SinglePvtDead.cpp

@@ -62,6 +62,30 @@ namespace Opm
 //            << "\n\nvisc\n\n" << viscosity_ << std::endl;
     }
 
+    /// Constructor
+    SinglePvtDead::SinglePvtDead(const Opm::PvdoTable& pvdoTable)
+    {
+        // Copy data
+        const std::vector<double>& press = pvdoTable.getPressureColumn();
+        const std::vector<double>& b = pvdoTable.getFormationFactorColumn();
+        const std::vector<double>& visc = pvdoTable.getViscosityColumn();
+
+        const int sz = b.size();
+        std::vector<double> bInv(sz);
+        for (int i = 0; i < sz; ++i) {
+            bInv[i] = 1.0 / b[i];
+        }
+        b_ = NonuniformTableLinear<double>(press, bInv);
+        viscosity_ = NonuniformTableLinear<double>(press, visc);
+
+        // Dumping the created tables.
+//         static int count = 0;
+//         std::ofstream os((std::string("dump-") + boost::lexical_cast<std::string>(count++)).c_str());
+//         os.precision(15);
+//         os << "1/B\n\n" << one_over_B_
+//            << "\n\nvisc\n\n" << viscosity_ << std::endl;
+    }
+
     // Destructor
     SinglePvtDead::~SinglePvtDead()
     {

opm/core/props/pvt/SinglePvtDead.hpp

@@ -23,6 +23,9 @@
 
 #include <opm/core/props/pvt/SinglePvtInterface.hpp>
 #include <opm/core/utility/NonuniformTableLinear.hpp>
+
+#include <opm/parser/eclipse/Utility/PvdoTable.hpp>
+
 #include <vector>
 
 namespace Opm
@@ -40,6 +43,7 @@ namespace Opm
     public:
         typedef std::vector<std::vector<std::vector<double> > > table_t;
         SinglePvtDead(const table_t& pvd_table);
+        SinglePvtDead(const Opm::PvdoTable &pvdoTable);
         virtual ~SinglePvtDead();
 
         /// Viscosity as a function of p and z.

opm/core/props/pvt/SinglePvtDeadSpline.cpp

@@ -20,6 +20,9 @@
 #include "config.h"
 #include <opm/core/props/pvt/SinglePvtDeadSpline.hpp>
 #include <opm/core/utility/buildUniformMonotoneTable.hpp>
+
+#include <opm/parser/eclipse/Utility/SimpleTable.hpp>
+
 #include <algorithm>
 
 // Extra includes for debug dumping of tables.
@@ -29,7 +32,6 @@
 
 namespace Opm
 {
-
     //------------------------------------------------------------------------
     // Member functions
     //-------------------------------------------------------------------------
@@ -54,13 +56,42 @@ namespace Opm
         }
         buildUniformMonotoneTable(press, B_inv, samples, b_);
         buildUniformMonotoneTable(press, visc, samples, viscosity_);
+    }
 
-        // Dumping the created tables.
+    SinglePvtDeadSpline::SinglePvtDeadSpline(const Opm::PvdoTable &pvdoTable, int samples)
-//         static int count = 0;
+    {
-//         std::ofstream os((std::string("dump-") + boost::lexical_cast<std::string>(count++)).c_str());
+        int numRows = pvdoTable.numRows();
-//         os.precision(15);
+
-//         os << "1/B\n\n" << one_over_B_
+        // Copy data
-//            << "\n\nvisc\n\n" << viscosity_ << std::endl;
+        const std::vector<double> &press = pvdoTable.getPressureColumn();
+        const std::vector<double> &B = pvdoTable.getFormationFactorColumn();
+        const std::vector<double> &visc = pvdoTable.getViscosityColumn();
+
+        std::vector<double> B_inv(numRows);
+        for (int i = 0; i < numRows; ++i) {
+            B_inv[i] = 1.0 / B[i];
+        }
+
+        buildUniformMonotoneTable(press, B_inv, samples, b_);
+        buildUniformMonotoneTable(press, visc, samples, viscosity_);
+    }
+
+    SinglePvtDeadSpline::SinglePvtDeadSpline(const Opm::PvdgTable &pvdgTable, int samples)
+    {
+        int numRows = pvdgTable.numRows();
+
+        // Copy data
+        const std::vector<double> &press = pvdgTable.getPressureColumn();
+        const std::vector<double> &B = pvdgTable.getFormationFactorColumn();
+        const std::vector<double> &visc = pvdgTable.getViscosityColumn();
+
+        std::vector<double> B_inv(numRows);
+        for (int i = 0; i < numRows; ++i) {
+            B_inv[i] = 1.0 / B[i];
+        }
+
+        buildUniformMonotoneTable(press, B_inv, samples, b_);
+        buildUniformMonotoneTable(press, visc, samples, viscosity_);
     }
 
     // Destructor

opm/core/props/pvt/SinglePvtDeadSpline.hpp

@@ -20,9 +20,12 @@
 #ifndef OPM_SINGLEPVTDEADSPLINE_HEADER_INCLUDED
 #define OPM_SINGLEPVTDEADSPLINE_HEADER_INCLUDED
 
-
 #include <opm/core/props/pvt/SinglePvtInterface.hpp>
 #include <opm/core/utility/UniformTableLinear.hpp>
+
+#include <opm/parser/eclipse/Utility/PvdoTable.hpp>
+#include <opm/parser/eclipse/Utility/PvdgTable.hpp>
+
 #include <vector>
 
 namespace Opm
@@ -41,6 +44,8 @@ namespace Opm
         typedef std::vector<std::vector<std::vector<double> > > table_t;
 
         SinglePvtDeadSpline(const table_t& pvd_table, const int samples);
+        SinglePvtDeadSpline(const Opm::PvdoTable &pvdoTable, int samples);
+        SinglePvtDeadSpline(const Opm::PvdgTable &pvdgTable, int samples);
         virtual ~SinglePvtDeadSpline();
 
         /// Viscosity as a function of p and z.

opm/core/props/pvt/SinglePvtLiveGas.cpp

@@ -32,8 +32,8 @@
 #include <opm/core/props/pvt/SinglePvtLiveGas.hpp>
 #include <opm/core/utility/ErrorMacros.hpp>
 #include <opm/core/utility/linearInterpolation.hpp>
-#include <algorithm>
 
+#include <algorithm>
 
 namespace Opm
 {
@@ -41,6 +41,7 @@ namespace Opm
     using Opm::linearInterpolation;
     using Opm::linearInterpolationDerivative;
 
+
     //------------------------------------------------------------------------
     // Member functions
     //-------------------------------------------------------------------------
@@ -81,6 +82,27 @@ namespace Opm
         }
     }
 
+    SinglePvtLiveGas::SinglePvtLiveGas(const Opm::PvtgTable& pvtgTable)
+    {
+        // GAS, PVTG
+        saturated_gas_table_.resize(4);
+        saturated_gas_table_[0] = pvtgTable.getOuterTable()->getPressureColumn();
+        saturated_gas_table_[1] = pvtgTable.getOuterTable()->getGasFormationFactorColumn();
+        saturated_gas_table_[2] = pvtgTable.getOuterTable()->getGasViscosityColumn();
+        saturated_gas_table_[3] = pvtgTable.getOuterTable()->getOilSolubilityColumn();
+
+        int sz = pvtgTable.getOuterTable()->numRows();
+        undersat_gas_tables_.resize(sz);
+        for (int i=0; i<sz; ++i) {
+            const auto &undersatTable = *pvtgTable.getInnerTable(i);
+
+            undersat_gas_tables_[i].resize(3);
+            undersat_gas_tables_[i][0] = undersatTable.getOilSolubilityColumn();
+            undersat_gas_tables_[i][1] = undersatTable.getGasFormationFactorColumn();
+            undersat_gas_tables_[i][2] = pvtgTable.getOuterTable()->getGasViscosityColumn();
+        }
+    }
+
     // Destructor
      SinglePvtLiveGas::~SinglePvtLiveGas()
     {

opm/core/props/pvt/SinglePvtLiveGas.hpp

@@ -21,6 +21,9 @@
 #define OPM_SINGLEPVTLIVEGAS_HEADER_INCLUDED
 
 #include <opm/core/props/pvt/SinglePvtInterface.hpp>
+
+#include <opm/parser/eclipse/Utility/PvtgTable.hpp>
+
 #include <vector>
 
 namespace Opm
@@ -38,6 +41,7 @@ namespace Opm
         typedef std::vector<std::vector<std::vector<double> > > table_t;
 
         SinglePvtLiveGas(const table_t& pvto);
+        SinglePvtLiveGas(const Opm::PvtgTable& pvtgTable);
         virtual ~SinglePvtLiveGas();
 
         /// Viscosity as a function of p and z.

opm/core/props/pvt/SinglePvtLiveOil.cpp

@@ -21,8 +21,8 @@
 #include <opm/core/props/pvt/SinglePvtLiveOil.hpp>
 #include <opm/core/utility/ErrorMacros.hpp>
 #include <opm/core/utility/linearInterpolation.hpp>
-#include <algorithm>
 
+#include <algorithm>
 
 namespace Opm
 {
@@ -31,7 +31,6 @@ namespace Opm
     using Opm::linearInterpolationDerivative;
     using Opm::tableIndex;
 
-
     //------------------------------------------------------------------------
     // Member functions
     //-------------------------------------------------------------------------
@@ -102,7 +101,73 @@ namespace Opm
                 undersat_oil_tables_[i][2].push_back(mu);
             }
         }
+    }
 
+    SinglePvtLiveOil::SinglePvtLiveOil(const Opm::PvtoTable &pvtoTable)
+    {
+        const auto saturatedPvto = pvtoTable.getOuterTable();
+
+        // OIL, PVTO
+        saturated_oil_table_.resize(4);
+        const int sz =  saturatedPvto->numRows();
+        for (int k=0; k<4; ++k) {
+            saturated_oil_table_[k].resize(sz);
+        }
+        for (int i=0; i<sz; ++i) {
+            saturated_oil_table_[0][i] = saturatedPvto->getPressureColumn()[i]; // p
+            saturated_oil_table_[1][i] = 1.0/saturatedPvto->getOilFormationFactorColumn()[i]; // 1/Bo
+            saturated_oil_table_[2][i] = saturatedPvto->getOilViscosityColumn()[i]; // mu_o
+            saturated_oil_table_[3][i] = saturatedPvto->getGasSolubilityColumn()[i]; // Rs
+        }
+
+        undersat_oil_tables_.resize(sz);
+        for (int i=0; i<sz; ++i) {
+            const auto undersaturatedPvto = pvtoTable.getInnerTable(i);
+
+            undersat_oil_tables_[i].resize(3);
+            int tsize = undersaturatedPvto->numRows();
+            undersat_oil_tables_[i][0].resize(tsize);
+            undersat_oil_tables_[i][1].resize(tsize);
+            undersat_oil_tables_[i][2].resize(tsize);
+            for (int j=0; j<tsize; ++j) {
+                undersat_oil_tables_[i][0][j] = undersaturatedPvto->getPressureColumn()[j];  // p
+                undersat_oil_tables_[i][1][j] = 1.0/undersaturatedPvto->getOilFormationFactorColumn()[j];  // 1/Bo
+                undersat_oil_tables_[i][2][j] = undersaturatedPvto->getOilViscosityColumn()[j];  // mu_o
+            }
+        }
+
+        // Complete undersaturated tables by extrapolating from existing data
+        // as is done in Eclipse and Mrst
+        int iNext = -1;
+        for (int i=0; i<sz; ++i) {
+            // Skip records already containing undersaturated data
+            if (undersat_oil_tables_[i][0].size() > 1) {
+                continue;
+            }
+            // Look ahead for next record containing undersaturated data
+            if (iNext < i) {
+                iNext = i+1;
+                while (iNext<sz && undersat_oil_tables_[iNext][0].size() < 2) {
+                    ++iNext;
+                }
+                if (iNext == sz) OPM_THROW(std::runtime_error,"Unable to complete undersaturated table.");
+            }
+            // Add undersaturated data to current record while maintaining compressibility and viscosibility
+            typedef std::vector<std::vector<std::vector<double> > >::size_type sz_t;
+            for (sz_t j=1; j<undersat_oil_tables_[iNext][0].size(); ++j) {
+                double diffPressure = undersat_oil_tables_[iNext][0][j]-undersat_oil_tables_[iNext][0][j-1];
+                double pressure = undersat_oil_tables_[i][0].back()+diffPressure;
+                undersat_oil_tables_[i][0].push_back(pressure);
+                double compr = (1.0/undersat_oil_tables_[iNext][1][j]-1.0/undersat_oil_tables_[iNext][1][j-1])
+                        / (0.5*(1.0/undersat_oil_tables_[iNext][1][j]+1.0/undersat_oil_tables_[iNext][1][j-1]));
+                double B = (1.0/undersat_oil_tables_[i][1].back())*(1.0+0.5*compr)/(1.0-0.5*compr);
+                undersat_oil_tables_[i][1].push_back(1.0/B);
+                double visc = (undersat_oil_tables_[iNext][2][j]-undersat_oil_tables_[iNext][2][j-1])
+                        / (0.5*(undersat_oil_tables_[iNext][2][j]+undersat_oil_tables_[iNext][2][j-1]));
+                double mu = (undersat_oil_tables_[i][2].back())*(1.0+0.5*visc)/(1.0-0.5*visc);
+                undersat_oil_tables_[i][2].push_back(mu);
+            }
+        }
     }
 
     /// Destructor.

opm/core/props/pvt/SinglePvtLiveOil.hpp

@@ -20,8 +20,10 @@
 #ifndef OPM_SINGLEPVTLIVEOIL_HEADER_INCLUDED
 #define OPM_SINGLEPVTLIVEOIL_HEADER_INCLUDED
 
-
 #include <opm/core/props/pvt/SinglePvtInterface.hpp>
+
+#include <opm/parser/eclipse/Utility/PvtoTable.hpp>
+
 #include <vector>
 
 namespace Opm
@@ -39,6 +41,7 @@ namespace Opm
         typedef std::vector<std::vector<std::vector<double> > > table_t;
 
         SinglePvtLiveOil(const table_t& pvto);
+        SinglePvtLiveOil(const Opm::PvtoTable &pvtoTable);
         virtual ~SinglePvtLiveOil();
 
         /// Viscosity as a function of p and z.

opm/core/props/rock/RockCompressibility.cpp

@@ -25,6 +25,9 @@
 #include <opm/core/utility/ErrorMacros.hpp>
 #include <opm/core/utility/linearInterpolation.hpp>
 
+#include <opm/parser/eclipse/Utility/RocktabTable.hpp>
+#include <opm/parser/eclipse/Utility/RockTable.hpp>
+
 #include <iostream>
 
 namespace Opm
@@ -65,6 +68,44 @@ namespace Opm
         }
     }
 
+    RockCompressibility::RockCompressibility(Opm::DeckConstPtr newParserDeck)
+        : pref_(0.0),
+          rock_comp_(0.0)
+    {
+        if (newParserDeck->hasKeyword("ROCKTAB")) {
+            Opm::DeckKeywordConstPtr rtKeyword = newParserDeck->getKeyword("ROCKTAB");
+            if (rtKeyword->size() != 1)
+                OPM_THROW(std::runtime_error, "Can only handle a single region in ROCKTAB.");
+
+            // the number of colums of the "ROCKTAB" keyword
+            // depends on the presence of the "RKTRMDIR"
+            // keyword. Messy stuff...
+            bool isDirectional = newParserDeck->hasKeyword("RKTRMDIR");
+            if (isDirectional)
+            {
+                // well, okay. we don't support non-isotropic
+                // transmissibility multipliers yet
+                OPM_THROW(std::runtime_error, "Support for non-isotropic "
+                          "transmissibility multipliers is not implemented yet.");
+            };
+
+            Opm::RocktabTable rocktabTable(rtKeyword, isDirectional);
+
+            p_ = rocktabTable.getPressureColumn();
+            poromult_ = rocktabTable.getPoreVolumeMultiplierColumn();
+            transmult_ =  rocktabTable.getTransmissibilityMultiplierColumn();
+        } else if (newParserDeck->hasKeyword("ROCK")) {
+            Opm::RockTable rockTable(newParserDeck->getKeyword("ROCK"));
+            if (rockTable.numRows() != 1)
+                OPM_THROW(std::runtime_error, "Can only handle a single region in ROCK.");
+
+            pref_ = rockTable.getPressureColumn()[0];
+            rock_comp_ = rockTable.getCompressibilityColumn()[0];
+        } else {
+            std::cout << "**** warning: no rock compressibility data found in deck (ROCK or ROCKTAB)." << std::endl;
+        }
+    }
+
     bool RockCompressibility::isActive() const
     {
         return !p_.empty() || (rock_comp_ != 0.0);

opm/core/props/rock/RockCompressibility.hpp

@@ -20,6 +20,8 @@
 #ifndef OPM_ROCKCOMPRESSIBILITY_HEADER_INCLUDED
 #define OPM_ROCKCOMPRESSIBILITY_HEADER_INCLUDED
 
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+
 #include <vector>
 
 namespace Opm
@@ -35,6 +37,10 @@ namespace Opm
         /// Looks for the keywords ROCK and ROCKTAB.
         RockCompressibility(const EclipseGridParser& deck);
 
+        /// Construct from input deck.
+        /// Looks for the keywords ROCK and ROCKTAB.
+        RockCompressibility(Opm::DeckConstPtr newParserDeck);
+
         /// Construct from parameters.
         /// Accepts the following parameters (with defaults).
         ///    rock_compressibility_pref (100.0)   [given in bar]

opm/core/props/rock/RockFromDeck.cpp

@@ -21,6 +21,9 @@
 #include "config.h"
 #include <opm/core/props/rock/RockFromDeck.hpp>
 #include <opm/core/grid.h>
+
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+
 #include <array>
 
 namespace Opm
@@ -37,6 +40,10 @@ namespace Opm
         PermeabilityKind fillTensor(const EclipseGridParser&                 parser,
                                     std::vector<const std::vector<double>*>& tensor,
                                     std::array<int,9>&                     kmap);
+        PermeabilityKind fillTensor(Opm::DeckConstPtr newParserDeck,
+                                    std::vector<const std::vector<double>*>& tensor,
+                                    std::array<int,9>&                     kmap);
+
     } // anonymous namespace
 
 
@@ -64,6 +71,15 @@ namespace Opm
         assignPermeability(deck, grid, perm_threshold);
     }
 
+    void RockFromDeck::init(Opm::DeckConstPtr newParserDeck,
+                            const UnstructuredGrid& grid)
+    {
+        assignPorosity(newParserDeck, grid);
+        permfield_valid_.assign(grid.number_of_cells, false);
+        const double perm_threshold = 0.0; // Maybe turn into parameter?
+        assignPermeability(newParserDeck, grid, perm_threshold);
+    }
+
 
     void RockFromDeck::assignPorosity(const EclipseGridParser& parser,
                                       const UnstructuredGrid& grid)
@@ -79,6 +95,21 @@ namespace Opm
         }
     }
 
+    void RockFromDeck::assignPorosity(Opm::DeckConstPtr newParserDeck,
+                                      const UnstructuredGrid& grid)
+    {
+        porosity_.assign(grid.number_of_cells, 1.0);
+        const int* gc = grid.global_cell;
+        if (newParserDeck->hasKeyword("PORO")) {
+            const std::vector<double>& poro = newParserDeck->getKeyword("PORO")->getSIDoubleData();
+            for (int c = 0; c < int(porosity_.size()); ++c) {
+                const int deck_pos = (gc == NULL) ? c : gc[c];
+                assert(0 <= c && c < (int) porosity_.size());
+                assert(0 <= deck_pos && deck_pos < (int) poro.size());
+                porosity_[c] = poro[deck_pos];
+            }
+        }
+    }
 
 
     void RockFromDeck::assignPermeability(const EclipseGridParser& parser,
@@ -135,6 +166,59 @@ namespace Opm
         }
     }
 
+    void RockFromDeck::assignPermeability(Opm::DeckConstPtr newParserDeck,
+                                          const UnstructuredGrid& grid,
+                                          double perm_threshold)
+    {
+        const int dim              = 3;
+        const int num_global_cells = grid.cartdims[0]*grid.cartdims[1]*grid.cartdims[2];
+        const int nc = grid.number_of_cells;
+
+        assert(num_global_cells > 0);
+
+        permeability_.assign(dim * dim * nc, 0.0);
+
+        std::vector<const std::vector<double>*> tensor;
+        tensor.reserve(10);
+
+        const std::vector<double> zero(num_global_cells, 0.0);
+        tensor.push_back(&zero);
+
+        std::array<int,9> kmap;
+        PermeabilityKind pkind = fillTensor(newParserDeck, tensor, kmap);
+        if (pkind == Invalid) {
+            OPM_THROW(std::runtime_error, "Invalid permeability field.");
+        }
+
+        // Assign permeability values only if such values are
+        // given in the input deck represented by 'newParserDeck'.  In
+        // other words: Don't set any (arbitrary) default values.
+        // It is infinitely better to experience a reproducible
+        // crash than subtle errors resulting from a (poorly
+        // chosen) default value...
+        //
+        if (tensor.size() > 1) {
+            const int* gc = grid.global_cell;
+            int off = 0;
+
+            for (int c = 0; c < nc; ++c, off += dim*dim) {
+                // SharedPermTensor K(dim, dim, &permeability_[off]);
+                int       kix  = 0;
+                const int glob = (gc == NULL) ? c : gc[c];
+
+                for (int i = 0; i < dim; ++i) {
+                    for (int j = 0; j < dim; ++j, ++kix) {
+                        // K(i,j) = (*tensor[kmap[kix]])[glob];
+                        permeability_[off + kix] = (*tensor[kmap[kix]])[glob];
+                    }
+                    // K(i,i) = std::max(K(i,i), perm_threshold);
+                    permeability_[off + 3*i + i] = std::max(permeability_[off + 3*i + i], perm_threshold);
+                }
+
+                permfield_valid_[c] = std::vector<unsigned char>::value_type(1);
+            }
+        }
+    }
 
     namespace {
 
@@ -204,6 +288,72 @@ namespace Opm
             return retval;
         }
 
+        /// @brief
+        ///    Classify and verify a given permeability specification
+        ///    from a structural point of view.  In particular, we
+        ///    verify that there are no off-diagonal permeability
+        ///    components such as @f$k_{xy}@f$ unless the
+        ///    corresponding diagonal components are known as well.
+        ///
+        /// @param newParserDeck [in]
+        ///    An Eclipse data parser capable of answering which
+        ///    permeability components are present in a given input
+        ///    deck.
+        ///
+        /// @return
+        ///    An enum value with the following possible values:
+        ///        ScalarPerm     only one component was given.
+        ///        DiagonalPerm   more than one component given.
+        ///        TensorPerm     at least one cross-component given.
+        ///        None           no components given.
+        ///        Invalid        invalid set of components given.
+        PermeabilityKind classifyPermeability(Opm::DeckConstPtr newParserDeck)
+        {
+            const bool xx = newParserDeck->hasKeyword("PERMX" );
+            const bool xy = newParserDeck->hasKeyword("PERMXY");
+            const bool xz = newParserDeck->hasKeyword("PERMXZ");
+
+            const bool yx = newParserDeck->hasKeyword("PERMYX");
+            const bool yy = newParserDeck->hasKeyword("PERMY" );
+            const bool yz = newParserDeck->hasKeyword("PERMYZ");
+
+            const bool zx = newParserDeck->hasKeyword("PERMZX");
+            const bool zy = newParserDeck->hasKeyword("PERMZY");
+            const bool zz = newParserDeck->hasKeyword("PERMZ" );
+
+            int num_cross_comp = xy + xz + yx + yz + zx + zy;
+            int num_comp       = xx + yy + zz + num_cross_comp;
+            PermeabilityKind retval = None;
+            if (num_cross_comp > 0) {
+                retval = TensorPerm;
+            } else {
+                if (num_comp == 1) {
+                    retval = ScalarPerm;
+                } else if (num_comp >= 2) {
+                    retval = DiagonalPerm;
+                }
+            }
+
+            bool ok = true;
+            if (num_comp > 0) {
+                // At least one tensor component specified on input.
+                // Verify that any remaining components are OK from a
+                // structural point of view.  In particular, there
+                // must not be any cross-components (e.g., k_{xy})
+                // unless the corresponding diagonal component (e.g.,
+                // k_{xx}) is present as well...
+                //
+                ok =        xx || !(xy || xz || yx || zx) ;
+                ok = ok && (yy || !(yx || yz || xy || zy));
+                ok = ok && (zz || !(zx || zy || xz || yz));
+            }
+            if (!ok) {
+                retval = Invalid;
+            }
+
+            return retval;
+        }
+
 
         /// @brief
         ///    Copy isotropic (scalar) permeability to other diagonal
@@ -333,6 +483,106 @@ namespace Opm
             return kind;
         }
 
+        /// @brief
+        ///   Extract pointers to appropriate tensor components from
+        ///   input deck.  The permeability tensor is, generally,
+        ///   @code
+        ///        [ kxx  kxy  kxz ]
+        ///    K = [ kyx  kyy  kyz ]
+        ///        [ kzx  kzy  kzz ]
+        ///   @endcode
+        ///   We store these values in a linear array using natural
+        ///   ordering with the column index cycling the most rapidly.
+        ///   In particular we use the representation
+        ///   @code
+        ///        [  0    1    2    3    4    5    6    7    8  ]
+        ///    K = [ kxx, kxy, kxz, kyx, kyy, kyz, kzx, kzy, kzz ]
+        ///   @endcode
+        ///   Moreover, we explicitly enforce symmetric tensors by
+        ///   assigning
+        ///   @code
+        ///     3     1       6     2       7     5
+        ///    kyx = kxy,    kzx = kxz,    kzy = kyz
+        ///   @endcode
+        ///   However, we make no attempt at enforcing positive
+        ///   definite tensors.
+        ///
+        /// @param [in]  parser
+        ///    An Eclipse data parser capable of answering which
+        ///    permeability components are present in a given input
+        ///    deck as well as retrieving the numerical value of each
+        ///    permeability component in each grid cell.
+        ///
+        /// @param [out] tensor
+        /// @param [out] kmap
+        PermeabilityKind fillTensor(Opm::DeckConstPtr newParserDeck,
+                                    std::vector<const std::vector<double>*>& tensor,
+                                    std::array<int,9>&                     kmap)
+        {
+            PermeabilityKind kind = classifyPermeability(newParserDeck);
+            if (kind == Invalid) {
+                OPM_THROW(std::runtime_error, "Invalid set of permeability fields given.");
+            }
+            assert(tensor.size() == 1);
+            for (int i = 0; i < 9; ++i) { kmap[i] = 0; }
+
+            enum { xx, xy, xz,    // 0, 1, 2
+                   yx, yy, yz,    // 3, 4, 5
+                   zx, zy, zz };  // 6, 7, 8
+
+            // -----------------------------------------------------------
+            // 1st row: [kxx, kxy, kxz]
+            if (newParserDeck->hasKeyword("PERMX" )) {
+                kmap[xx] = tensor.size();
+                tensor.push_back(&newParserDeck->getKeyword("PERMX")->getSIDoubleData());
+
+                setScalarPermIfNeeded(kmap, xx, yy, zz);
+            }
+            if (newParserDeck->hasKeyword("PERMXY")) {
+                kmap[xy] = kmap[yx] = tensor.size();  // Enforce symmetry.
+                tensor.push_back(&newParserDeck->getKeyword("PERMXY")->getSIDoubleData());
+            }
+            if (newParserDeck->hasKeyword("PERMXZ")) {
+                kmap[xz] = kmap[zx] = tensor.size();  // Enforce symmetry.
+                tensor.push_back(&newParserDeck->getKeyword("PERMXZ")->getSIDoubleData());
+            }
+
+            // -----------------------------------------------------------
+            // 2nd row: [kyx, kyy, kyz]
+            if (newParserDeck->hasKeyword("PERMYX")) {
+                kmap[yx] = kmap[xy] = tensor.size();  // Enforce symmetry.
+                tensor.push_back(&newParserDeck->getKeyword("PERMYX")->getSIDoubleData());
+            }
+            if (newParserDeck->hasKeyword("PERMY" )) {
+                kmap[yy] = tensor.size();
+                tensor.push_back(&newParserDeck->getKeyword("PERMY")->getSIDoubleData());
+
+                setScalarPermIfNeeded(kmap, yy, zz, xx);
+            }
+            if (newParserDeck->hasKeyword("PERMYZ")) {
+                kmap[yz] = kmap[zy] = tensor.size();  // Enforce symmetry.
+                tensor.push_back(&newParserDeck->getKeyword("PERMYZ")->getSIDoubleData());
+            }
+
+            // -----------------------------------------------------------
+            // 3rd row: [kzx, kzy, kzz]
+            if (newParserDeck->hasKeyword("PERMZX")) {
+                kmap[zx] = kmap[xz] = tensor.size();  // Enforce symmetry.
+                tensor.push_back(&newParserDeck->getKeyword("PERMZX")->getSIDoubleData());
+            }
+            if (newParserDeck->hasKeyword("PERMZY")) {
+                kmap[zy] = kmap[yz] = tensor.size();  // Enforce symmetry.
+                tensor.push_back(&newParserDeck->getKeyword("PERMZY")->getSIDoubleData());
+            }
+            if (newParserDeck->hasKeyword("PERMZ" )) {
+                kmap[zz] = tensor.size();
+                tensor.push_back(&newParserDeck->getKeyword("PERMZ")->getSIDoubleData());
+
+                setScalarPermIfNeeded(kmap, zz, xx, yy);
+            }
+            return kind;
+        }
+
     } // anonymous namespace
 
 } // namespace Opm

opm/core/props/rock/RockFromDeck.hpp

@@ -22,6 +22,9 @@
 
 
 #include <opm/core/io/eclipse/EclipseGridParser.hpp>
+
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+
 #include <vector>
 
 struct UnstructuredGrid;
@@ -43,6 +46,14 @@ namespace Opm
         void init(const EclipseGridParser& deck,
                   const UnstructuredGrid& grid);
 
+        /// Initialize from deck and grid.
+        /// \param  newParserDeck Deck produced by the opm-parser code
+        /// \param  grid          Grid to which property object applies, needed for the
+        ///                       mapping from cell indices (typically from a processed grid)
+        ///                       to logical cartesian indices consistent with the deck.
+        void init(Opm::DeckConstPtr newParserDeck,
+                  const UnstructuredGrid& grid);
+
         /// \return   D, the number of spatial dimensions. Always 3 for deck input.
         int numDimensions() const
         {
@@ -72,9 +83,14 @@ namespace Opm
     private:
         void assignPorosity(const EclipseGridParser& parser,
                             const UnstructuredGrid& grid);
+        void assignPorosity(Opm::DeckConstPtr newParserDeck,
+                            const UnstructuredGrid& grid);
         void assignPermeability(const EclipseGridParser& parser,
                                 const UnstructuredGrid& grid,
                                 const double perm_threshold);
+        void assignPermeability(Opm::DeckConstPtr newParserDeck,
+                                const UnstructuredGrid& grid,
+                                double perm_threshold);
 
         std::vector<double> porosity_;
         std::vector<double> permeability_;

opm/core/props/satfunc/SatFuncBase.hpp

@@ -22,6 +22,11 @@
 #include <opm/core/utility/buildUniformMonotoneTable.hpp>
 #include <opm/core/utility/NonuniformTableLinear.hpp>
 #include <opm/core/props/BlackoilPhases.hpp>
+
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+#include <opm/parser/eclipse/Utility/SwofTable.hpp>
+#include <opm/parser/eclipse/Utility/SgofTable.hpp>
+
 #include <vector>
 
 namespace Opm
@@ -78,6 +83,10 @@ namespace Opm
                   const int table_num,
                   const PhaseUsage phase_usg,
                   const int samples);
+        void init(Opm::DeckConstPtr newParserDeck,
+                  const int table_num,
+                  const PhaseUsage phase_usg,
+                  const int samples);
         void updateSatHyst(const double* s,
                            const EPSTransforms* epst, 
                            const EPSTransforms* epst_hyst, 
@@ -248,6 +257,138 @@ namespace Opm
         }
     }
 
+    template <class TableType>
+    void SatFuncBase<TableType>::init(Opm::DeckConstPtr newParserDeck,
+                                      const int table_num,
+                                      const PhaseUsage phase_usg,
+                                      const int samples)
+    {
+        phase_usage = phase_usg;
+        double swco = 0.0;
+        double swmax = 1.0;
+        if (phase_usage.phase_used[Aqua]) {
+            Opm::SwofTable swof(newParserDeck->getKeyword("SWOF"), table_num);
+            const std::vector<double>& sw = swof.getSwColumn();
+            const std::vector<double>& krw = swof.getKrwColumn();
+            const std::vector<double>& krow = swof.getKrowColumn();
+            const std::vector<double>& pcow = swof.getPcowColumn();
+            if (krw.front() != 0.0 || krow.back() != 0.0) {
+                OPM_THROW(std::runtime_error, "Error SWOF data - non-zero krw(swco) and/or krow(1-sor)");
+            }
+
+            // Extend the tables with constant values such that the
+            // derivatives at the endpoints are zero
+            int n = sw.size();
+            std::vector<double> sw_ex(n+2);
+            std::vector<double> krw_ex(n+2);
+            std::vector<double> krow_ex(n+2);
+            std::vector<double> pcow_ex(n+2);
+
+            extendTable(sw,sw_ex,1);
+            extendTable(krw,krw_ex,0);
+            extendTable(krow,krow_ex,0);
+            extendTable(pcow,pcow_ex,0);
+
+            initializeTableType(krw_,sw_ex, krw_ex, samples);
+            initializeTableType(krow_,sw_ex, krow_ex, samples);
+            initializeTableType(pcow_,sw_ex, pcow_ex, samples);
+
+            krocw_ = krow[0]; // At connate water -> ecl. SWOF
+            swco = sw[0];
+            smin_[phase_usage.phase_pos[Aqua]] = sw[0];
+            swmax = sw.back();
+            smax_[phase_usage.phase_pos[Aqua]] = sw.back();
+
+            krwmax_ = krw.back();
+            kromax_ = krow.front();
+            swcr_ = swmax;
+            sowcr_ = 1.0 - swco;
+            krwr_ = krw.back();
+            krorw_ = krow.front();
+            for (std::vector<double>::size_type i=1; i<sw.size(); ++i) {
+                if (krw[i]> 0.0) {
+                   swcr_ = sw[i-1];
+                   krorw_ = krow[i-1];
+                   break;
+                }
+            }
+            for (std::vector<double>::size_type i=sw.size()-1; i>=1; --i) {
+                if (krow[i-1]> 0.0) {
+                   sowcr_ = 1.0 - sw[i];
+                   krwr_ = krw[i];
+                   break;
+                }
+            }
+        }
+        if (phase_usage.phase_used[Vapour]) {
+            Opm::SgofTable sgof(newParserDeck->getKeyword("SGOF"), table_num);
+            const std::vector<double>& sg = sgof.getSgColumn();
+            const std::vector<double>& krg = sgof.getKrgColumn();
+            const std::vector<double>& krog = sgof.getKrogColumn();
+            const std::vector<double>& pcog = sgof.getPcogColumn();
+
+            // Extend the tables with constant values such that the
+            // derivatives at the endpoints are zero
+            int n = sg.size();
+            std::vector<double> sg_ex(n+2);
+            std::vector<double> krg_ex(n+2);
+            std::vector<double> krog_ex(n+2);
+            std::vector<double> pcog_ex(n+2);
+
+            extendTable(sg,sg_ex,1);
+            extendTable(krg,krg_ex,0);
+            extendTable(krog,krog_ex,0);
+            extendTable(pcog,pcog_ex,0);
+
+            initializeTableType(krg_,sg_ex, krg_ex, samples);
+            initializeTableType(krog_,sg_ex, krog_ex, samples);
+            initializeTableType(pcog_,sg_ex, pcog_ex, samples);
+
+            smin_[phase_usage.phase_pos[Vapour]] = sg[0];
+            if (std::fabs(sg.back() + swco - 1.0) > 1e-3) {
+                OPM_THROW(std::runtime_error, "Gas maximum saturation in SGOF table = " << sg.back() <<
+                      ", should equal (1.0 - connate water sat) = " << (1.0 - swco));
+            }
+            smax_[phase_usage.phase_pos[Vapour]] = sg.back();
+            smin_[phase_usage.phase_pos[Vapour]] = sg.front();
+            krgmax_ = krg.back();
+
+            sgcr_ = sg.front();
+            sogcr_ = 1.0 - sg.back();
+            krgr_ = krg.back();
+            krorg_ = krg.front();
+            for (std::vector<double>::size_type i=1; i<sg.size(); ++i) {
+                if (krg[i]> 0.0) {
+                   sgcr_ = sg[i-1];
+                   krorg_ = krog[i-1];
+                   break;
+                }
+            }
+            for (std::vector<double>::size_type i=sg.size()-1; i>=1; --i) {
+                if (krog[i-1]> 0.0) {
+                   sogcr_ = 1.0 - sg[i];
+                   krgr_ = krg[i];
+                   break;
+                }
+            }
+
+        }
+
+        if (phase_usage.phase_used[Vapour] && phase_usage.phase_used[Aqua]) {
+            sowcr_ -= smin_[phase_usage.phase_pos[Vapour]];
+            sogcr_ -= smin_[phase_usage.phase_pos[Aqua]];
+            smin_[phase_usage.phase_pos[Liquid]] = 0.0;
+            smax_[phase_usage.phase_pos[Liquid]] = 1.0 - smin_[phase_usage.phase_pos[Aqua]]
+                                                       - smin_[phase_usage.phase_pos[Vapour]];  // First entry in SGOF-table supposed to be zero anyway ...
+        } else if (phase_usage.phase_used[Aqua]) {
+            smin_[phase_usage.phase_pos[Liquid]] = 1.0 - smax_[phase_usage.phase_pos[Aqua]];
+            smax_[phase_usage.phase_pos[Liquid]] = 1.0 - smin_[phase_usage.phase_pos[Aqua]];
+        } else if (phase_usage.phase_used[Vapour]) {
+            smin_[phase_usage.phase_pos[Liquid]] = 1.0 - smax_[phase_usage.phase_pos[Vapour]];
+            smax_[phase_usage.phase_pos[Liquid]] = 1.0 - smin_[phase_usage.phase_pos[Vapour]];
+        }
+    }
+
     template <class TableType>
     void SatFuncBase<TableType>::updateSatHyst(const double* s,
                                                const EPSTransforms* epst, 

opm/core/props/satfunc/SatFuncSimple.cpp

@@ -36,7 +36,7 @@ namespace Opm
     void SatFuncBase<NonuniformTableLinear<double> >::initializeTableType(NonuniformTableLinear<double> & table,
                                                                            const std::vector<double>& arg,
                                                                            const std::vector<double>& value,
-                                                                          const int /*samples*/)
+                                                                           const int samples)
     {
       table = NonuniformTableLinear<double>(arg, value);
     }
@@ -84,7 +84,7 @@ namespace Opm
         }
     }
 
-    double EPSTransforms::Transform::scaleSatDeriv(double s, double /*s_r*/, double /*s_cr*/, double /*s_max*/) const
+    double EPSTransforms::Transform::scaleSatDeriv(double s, double s_r, double s_cr, double s_max) const
     {
         if (doNotScale) {
             return 1.0;

opm/core/props/satfunc/SatFuncSimple.hpp

@@ -31,23 +31,23 @@ namespace Opm
         void evalKrDeriv(const double* s, double* kr, double* dkrds) const;
         void evalPc(const double* s, double* pc) const;
         void evalPcDeriv(const double* s, double* pc, double* dpcds) const;
-
+        
-        void evalKr(const double* /*s*/, double* /*kr*/, const EPSTransforms* /*epst*/) const
+        void evalKr(const double* s, double* kr, const EPSTransforms* epst) const 
         {OPM_THROW(std::runtime_error, "SatFuncSimple   --  need to be implemented ...");}
-        void evalKr(const double* /*s*/, double* /*kr*/, const EPSTransforms* /*epst*/, const EPSTransforms* /*epst_hyst*/, const SatHyst* /*sat_hyst*/) const
+        void evalKr(const double* s, double* kr, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const 
         {OPM_THROW(std::runtime_error, "SatFuncSimple   --  need to be implemented ...");}
         void evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst) const;
-        void evalKrDeriv(const double* /*s*/, double* /*kr*/, double* /*dkrds*/, const EPSTransforms* /*epst*/, const EPSTransforms* /*epst_hyst*/, const SatHyst* /*sat_hyst*/) const
+        void evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const 
         {OPM_THROW(std::runtime_error, "SatFuncSimple   --  need to be implemented ...");}
-        void evalPc(const double* /*s*/, double* /*pc*/, const EPSTransforms* /*epst*/) const
+        void evalPc(const double* s, double* pc, const EPSTransforms* epst) const 
         {OPM_THROW(std::runtime_error, "SatFuncSimple   --  need to be implemented ...");}
-        void evalPcDeriv(const double* /*s*/, double* /*pc*/, double* /*dpcds*/, const EPSTransforms* /*epst*/) const
+        void evalPcDeriv(const double* s, double* pc, double* dpcds, const EPSTransforms* epst) const 
         {OPM_THROW(std::runtime_error, "SatFuncSimple   --  need to be implemented ...");}
 
     private:
 
     };
-
+    
     typedef SatFuncSimple<UniformTableLinear<double> > SatFuncSimpleUniform;
     typedef SatFuncSimple<NonuniformTableLinear<double> > SatFuncSimpleNonuniform;
 
@@ -185,7 +185,7 @@ namespace Opm
             double dkrww = _dsdsw*epst->wat.scaleKrDeriv(s[wpos], this->krw_.derivative(_sw));
             double krg = epst->gas.scaleKr(s[gpos], this->krg_(_sg), this->krgr_);
             double dkrgg = _dsdsg*epst->gas.scaleKrDeriv(s[gpos], this->krg_.derivative(_sg));
-            // TODO Check the arguments to the krow- and krog-tables below...
+            // TODO Check the arguments to the krow- and krog-tables below...  
             double krow = epst->watoil.scaleKr(1.0-s[wpos]-s[gpos], this->krow_(1.0-_sow-this->smin_[gpos]), this->krorw_);   // ????
             double dkrow = _dsdsow*epst->watoil.scaleKrDeriv(1.0-s[wpos]-s[gpos], this->krow_.derivative(1.0-_sow-this->smin_[gpos])); // ????
             //double krog = epst->gasoil.scaleKr(this->krog_(1.0-_sog-this->smin_[wpos]), 1.0-s[wpos]-s[gpos], this->krorg_);   // ????

opm/core/props/satfunc/SatFuncStone2.hpp

@@ -31,25 +31,25 @@ namespace Opm
         void evalKrDeriv(const double* s, double* kr, double* dkrds) const;
         void evalPc(const double* s, double* pc) const;
         void evalPcDeriv(const double* s, double* pc, double* dpcds) const;
-
+               
-        void evalKr(const double* /*s*/, double* /*kr*/, const EPSTransforms* /*epst*/) const
+        void evalKr(const double* s, double* kr, const EPSTransforms* epst) const 
         {OPM_THROW(std::runtime_error, "SatFuncStone2   --  need to be implemented ...");}
-        void evalKr(const double* /*s*/, double* /*kr*/, const EPSTransforms* /*epst*/, const EPSTransforms* /*epst_hyst*/, const SatHyst* /*sat_hyst*/) const
+        void evalKr(const double* s, double* kr, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const 
         {OPM_THROW(std::runtime_error, "SatFuncStone2   --  need to be implemented ...");}
-        void evalKrDeriv(const double* /*s*/, double* /*kr*/, double* /*dkrds*/, const EPSTransforms* /*epst*/) const
+        void evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst) const
         {OPM_THROW(std::runtime_error, "SatFuncStone2   --  need to be implemented ...");}
-        void evalKrDeriv(const double* /*s*/, double* /*kr*/, double* /*dkrds*/, const EPSTransforms* /*epst*/, const EPSTransforms* /*epst_hyst*/, const SatHyst* /*sat_hyst*/) const
+        void evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const 
         {OPM_THROW(std::runtime_error, "SatFuncStone2   --  need to be implemented ...");}
-        void evalPc(const double* /*s*/, double* /*pc*/, const EPSTransforms* /*epst*/) const
+        void evalPc(const double* s, double* pc, const EPSTransforms* epst) const 
         {OPM_THROW(std::runtime_error, "SatFuncStone2   --  need to be implemented ...");}
-        void evalPcDeriv(const double* /*s*/, double* /*pc*/, double* /*dpcds*/, const EPSTransforms* /*epst*/) const
+        void evalPcDeriv(const double* s, double* pc, double* dpcds, const EPSTransforms* epst) const 
         {OPM_THROW(std::runtime_error, "SatFuncStone2   --  need to be implemented ...");}
 
 
     private:
 
     };
-
+    
     typedef SatFuncStone2<UniformTableLinear<double> > SatFuncStone2Uniform;
     typedef SatFuncStone2<NonuniformTableLinear<double> > SatFuncStone2Nonuniform;
 

opm/core/props/satfunc/SaturationPropsFromDeck.hpp

@@ -27,6 +27,9 @@
 #include <opm/core/props/satfunc/SatFuncStone2.hpp>
 #include <opm/core/props/satfunc/SatFuncSimple.hpp>
 #include <opm/core/props/satfunc/SatFuncGwseg.hpp>
+
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+
 #include <vector>
 
 struct UnstructuredGrid;
@@ -60,6 +63,17 @@ namespace Opm
                   const UnstructuredGrid& grid,
                   const int samples);
 
+        /// Initialize from deck and grid.
+        /// \param[in]  deck     Deck input parser
+        /// \param[in]  grid     Grid to which property object applies, needed for the
+        ///                      mapping from cell indices (typically from a processed grid)
+        ///                      to logical cartesian indices consistent with the deck.
+        /// \param[in]  samples  Number of uniform sample points for saturation tables.
+        /// NOTE: samples will only be used with the SatFuncSetUniform template argument.
+        void init(Opm::DeckConstPtr newParserDeck,
+                  const UnstructuredGrid& grid,
+                  const int samples);
+
         /// \return   P, the number of phases.
         int numPhases() const;
 
@@ -132,6 +146,14 @@ namespace Opm
                         const UnstructuredGrid& grid,
                         const std::string& keyword,
                         std::vector<double>& scaleparam);
+        void initEPS(Opm::DeckConstPtr newParserDeck,
+                     const UnstructuredGrid& grid);
+        void initEPSHyst(Opm::DeckConstPtr newParserDeck,
+                         const UnstructuredGrid& grid);
+        void initEPSKey(Opm::DeckConstPtr newParserDeck,
+                        const UnstructuredGrid& grid,
+                        const std::string& keyword,
+                        std::vector<double>& scaleparam);
         void initEPSParam(const int cell, 
                           EPSTransforms::Transform& data,
                           const bool oil,
@@ -149,6 +171,11 @@ namespace Opm
                           const std::vector<double>& s0,
                           const std::vector<double>& krsr,
                           const std::vector<double>& krmax);
+
+        bool columnIsMasked_(Opm::DeckConstPtr newParserDeck,
+                             const std::string &keywordName,
+                             int columnIdx)
+        { return newParserDeck->getKeyword(keywordName)->getRecord(0)->getItem(0)->getSIDouble(0) != -1.0; }
     };
 
 

opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp

@@ -26,6 +26,11 @@
 #include <opm/core/props/phaseUsageFromDeck.hpp>
 #include <opm/core/grid.h>
 
+#include <opm/parser/eclipse/Utility/EndscaleWrapper.hpp>
+#include <opm/parser/eclipse/Utility/ScalecrsWrapper.hpp>
+#include <opm/parser/eclipse/Utility/EnptvdTable.hpp>
+#include <opm/parser/eclipse/Utility/EnkrvdTable.hpp>
+
 #include <iostream>
 
 namespace Opm
@@ -140,6 +145,123 @@ namespace Opm
     }
 
 
+    /// Initialize from deck.
+    template <class SatFuncSet>
+    void SaturationPropsFromDeck<SatFuncSet>::init(Opm::DeckConstPtr newParserDeck,
+                                                   const UnstructuredGrid& grid,
+                                                   const int samples)
+    {
+        phase_usage_ = phaseUsageFromDeck(newParserDeck);
+
+        // Extract input data.
+        // Oil phase should be active.
+        if (!phase_usage_.phase_used[Liquid]) {
+            OPM_THROW(std::runtime_error, "SaturationPropsFromDeck::init()   --  oil phase must be active.");
+        }
+
+        // Obtain SATNUM, if it exists, and create cell_to_func_.
+        // Otherwise, let the cell_to_func_ mapping be just empty.
+        int satfuncs_expected = 1;
+        if (newParserDeck->hasKeyword("SATNUM")) {
+            const std::vector<int>& satnum = newParserDeck->getKeyword("SATNUM")->getIntData();
+            satfuncs_expected = *std::max_element(satnum.begin(), satnum.end());
+            const int num_cells = grid.number_of_cells;
+            cell_to_func_.resize(num_cells);
+            const int* gc = grid.global_cell;
+            for (int cell = 0; cell < num_cells; ++cell) {
+                const int newParserDeck_pos = (gc == NULL) ? cell : gc[cell];
+                cell_to_func_[cell] = satnum[newParserDeck_pos] - 1;
+            }
+        }
+
+        // Find number of tables, check for consistency.
+        enum { Uninitialized = -1 };
+        int num_tables = Uninitialized;
+        if (phase_usage_.phase_used[Aqua]) {
+            num_tables = newParserDeck->getKeyword("SWOF")->size();
+            if (num_tables < satfuncs_expected) {
+                OPM_THROW(std::runtime_error, "Found " << num_tables << " SWOF tables, SATNUM specifies at least " << satfuncs_expected);
+            }
+        }
+        if (phase_usage_.phase_used[Vapour]) {
+            int num_sgof_tables = newParserDeck->getKeyword("SGOF")->size();
+            if (num_sgof_tables < satfuncs_expected) {
+                OPM_THROW(std::runtime_error, "Found " << num_tables << " SGOF tables, SATNUM specifies at least " << satfuncs_expected);
+            }
+            if (num_tables == Uninitialized) {
+                num_tables = num_sgof_tables;
+            } else if (num_tables != num_sgof_tables) {
+                OPM_THROW(std::runtime_error, "Inconsistent number of tables in SWOF and SGOF.");
+            }
+        }
+
+        // Initialize tables.
+        satfuncset_.resize(num_tables);
+        for (int table = 0; table < num_tables; ++table) {
+            satfuncset_[table].init(newParserDeck, table, phase_usage_, samples);
+        }
+
+        // Saturation table scaling
+        do_hyst_ = false;
+        do_eps_ = false;
+        do_3pt_ = false;
+        if (newParserDeck->hasKeyword("ENDSCALE")) {
+            Opm::EndscaleWrapper endscale(newParserDeck->getKeyword("ENDSCALE"));
+            if (endscale.directionSwitch() != std::string("NODIR")) {
+                OPM_THROW(std::runtime_error,
+                          "SaturationPropsFromDeck::init()   --  ENDSCALE: "
+                          "Currently only 'NODIR' accepted.");
+            }
+            if (endscale.isReversible()) {
+                OPM_THROW(std::runtime_error,
+                          "SaturationPropsFromDeck::init()   --  ENDSCALE: "
+                          "Currently only 'REVERS' accepted.");
+            }
+            if (newParserDeck->hasKeyword("SCALECRS")) {
+                Opm::ScalecrsWrapper scalecrs(newParserDeck->getKeyword("SCALECRS"));
+                if (scalecrs.isEnabled()) {
+                    do_3pt_ = true;
+                }
+            }
+            do_eps_ = true;
+
+            initEPS(newParserDeck, grid);
+
+            // For now, a primitive detection of hysteresis. TODO: SATOPTS HYSTER/ and EHYSTR
+            do_hyst_ =
+                newParserDeck->hasKeyword("ISWL")
+                || newParserDeck->hasKeyword("ISWU")
+                || newParserDeck->hasKeyword("ISWCR")
+                || newParserDeck->hasKeyword("ISGL")
+                || newParserDeck->hasKeyword("ISGU")
+                || newParserDeck->hasKeyword("ISGCR")
+                ||  newParserDeck->hasKeyword("ISOWCR")
+                || newParserDeck->hasKeyword("ISOGCR");
+            if (do_hyst_) {
+                if (newParserDeck->hasKeyword("KRW")
+                    || newParserDeck->hasKeyword("KRG")
+                    || newParserDeck->hasKeyword("KRO")
+                    || newParserDeck->hasKeyword("KRWR")
+                    || newParserDeck->hasKeyword("KRGR")
+                    || newParserDeck->hasKeyword("KRORW")
+                    || newParserDeck->hasKeyword("KRORG")
+                    || newParserDeck->hasKeyword("IKRW")
+                    || newParserDeck->hasKeyword("IKRG")
+                    || newParserDeck->hasKeyword("IKRO")
+                    || newParserDeck->hasKeyword("IKRWR")
+                    || newParserDeck->hasKeyword("IKRGR")
+                    || newParserDeck->hasKeyword("IKRORW")
+                    || newParserDeck->hasKeyword("IKRORG") ) {
+                    OPM_THROW(std::runtime_error,
+                              "SaturationPropsFromDeck::init()   --  ENDSCALE: "
+                              "Currently hysteresis and relperm value scaling "
+                              "cannot be combined.");
+                }
+                initEPSHyst(newParserDeck, grid);
+            }
+        }
+    }
+
 
 
     /// \return   P, the number of phases.
@@ -383,6 +505,69 @@ namespace Opm
         }
     }
 
+    // Initialize saturation scaling parameters
+    template <class SatFuncSet>
+    void SaturationPropsFromDeck<SatFuncSet>::initEPS(Opm::DeckConstPtr newParserDeck,
+                                                      const UnstructuredGrid& grid)
+    {
+        std::vector<double> swl, swcr, swu, sgl, sgcr, sgu, sowcr, sogcr;
+        std::vector<double> krw, krg, kro, krwr, krgr, krorw, krorg;
+        // Initialize saturation scaling parameter
+        initEPSKey(newParserDeck, grid, std::string("SWL"),   swl);
+        initEPSKey(newParserDeck, grid, std::string("SWU"),   swu);
+        initEPSKey(newParserDeck, grid, std::string("SWCR"),  swcr);
+        initEPSKey(newParserDeck, grid, std::string("SGL"),   sgl);
+        initEPSKey(newParserDeck, grid, std::string("SGU"),   sgu);
+        initEPSKey(newParserDeck, grid, std::string("SGCR"),  sgcr);
+        initEPSKey(newParserDeck, grid, std::string("SOWCR"), sowcr);
+        initEPSKey(newParserDeck, grid, std::string("SOGCR"), sogcr);
+        initEPSKey(newParserDeck, grid, std::string("KRW"),   krw);
+        initEPSKey(newParserDeck, grid, std::string("KRG"),   krg);
+        initEPSKey(newParserDeck, grid, std::string("KRO"),   kro);
+        initEPSKey(newParserDeck, grid, std::string("KRWR"),  krwr);
+        initEPSKey(newParserDeck, grid, std::string("KRGR"),  krgr);
+        initEPSKey(newParserDeck, grid, std::string("KRORW"), krorw);
+        initEPSKey(newParserDeck, grid, std::string("KRORG"), krorg);
+
+        eps_transf_.resize(grid.number_of_cells);
+
+        const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua];
+        const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour];
+        const bool oilWater = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && !phase_usage_.phase_used[Vapour];
+        const bool oilGas = !phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour];
+        const bool threephase = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour];
+
+        for (int cell = 0; cell < grid.number_of_cells; ++cell) {
+            if (oilWater) {
+                // ### krw
+                initEPSParam(cell, eps_transf_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos],
+                  funcForCell(cell).sowcr_, -1.0, funcForCell(cell).krwr_, funcForCell(cell).krwmax_, swl, swcr, swu, sowcr, sgl, krwr, krw);
+                // ### krow
+                initEPSParam(cell, eps_transf_[cell].watoil, true, 0.0, funcForCell(cell).sowcr_, funcForCell(cell).smin_[wpos],
+                  funcForCell(cell).swcr_, -1.0, funcForCell(cell).krorw_, funcForCell(cell).kromax_, swl, sowcr, swl, swcr, sgl, krorw, kro);
+            } else if (oilGas) {
+                // ### krg
+                initEPSParam(cell, eps_transf_[cell].gas, false, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_, funcForCell(cell).smax_[gpos],
+                  funcForCell(cell).sogcr_, -1.0, funcForCell(cell).krgr_, funcForCell(cell).krgmax_, sgl, sgcr, sgu, sogcr, swl, krgr, krg);
+                // ### krog
+                initEPSParam(cell, eps_transf_[cell].gasoil, true, 0.0, funcForCell(cell).sogcr_, funcForCell(cell).smin_[gpos],
+                  funcForCell(cell).sgcr_, -1.0, funcForCell(cell).krorg_, funcForCell(cell).kromax_, sgl, sogcr, sgl, sgcr, swl, krorg, kro);
+            } else if (threephase) {
+                // ### krw
+                initEPSParam(cell, eps_transf_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos], funcForCell(cell).sowcr_,
+                  funcForCell(cell).smin_[gpos], funcForCell(cell).krwr_, funcForCell(cell).krwmax_, swl, swcr, swu, sowcr, sgl, krwr, krw);
+                // ### krow
+                initEPSParam(cell, eps_transf_[cell].watoil, true, 0.0, funcForCell(cell).sowcr_, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_,
+                  funcForCell(cell).smin_[gpos], funcForCell(cell).krorw_, funcForCell(cell).kromax_, swl, sowcr, swl, swcr, sgl, krorw, kro);
+                // ### krg
+                initEPSParam(cell, eps_transf_[cell].gas, false, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_, funcForCell(cell).smax_[gpos], funcForCell(cell).sogcr_,
+                  funcForCell(cell).smin_[wpos], funcForCell(cell).krgr_, funcForCell(cell).krgmax_, sgl, sgcr, sgu, sogcr, swl, krgr, krg);
+                // ### krog
+                initEPSParam(cell, eps_transf_[cell].gasoil, true, 0.0, funcForCell(cell).sogcr_, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_,
+                  funcForCell(cell).smin_[wpos], funcForCell(cell).krorg_, funcForCell(cell).kromax_, sgl, sogcr, sgl, sgcr, swl, krorg, kro);
+            }
+        }
+    }
 
     // Initialize hysteresis saturation scaling parameters
     template <class SatFuncSet>
@@ -450,6 +635,71 @@ namespace Opm
         }
     }
 
+    // Initialize hysteresis saturation scaling parameters
+    template <class SatFuncSet>
+    void SaturationPropsFromDeck<SatFuncSet>::initEPSHyst(Opm::DeckConstPtr newParserDeck,
+                                                          const UnstructuredGrid& grid)
+    {
+        std::vector<double> iswl, iswcr, iswu, isgl, isgcr, isgu, isowcr, isogcr;
+        std::vector<double> ikrw, ikrg, ikro, ikrwr, ikrgr, ikrorw, ikrorg;
+        // Initialize hysteresis saturation scaling parameters
+        initEPSKey(newParserDeck, grid, std::string("ISWL"),   iswl);
+        initEPSKey(newParserDeck, grid, std::string("ISWU"),   iswu);
+        initEPSKey(newParserDeck, grid, std::string("ISWCR"),  iswcr);
+        initEPSKey(newParserDeck, grid, std::string("ISGL"),   isgl);
+        initEPSKey(newParserDeck, grid, std::string("ISGU"),   isgu);
+        initEPSKey(newParserDeck, grid, std::string("ISGCR"),  isgcr);
+        initEPSKey(newParserDeck, grid, std::string("ISOWCR"), isowcr);
+        initEPSKey(newParserDeck, grid, std::string("ISOGCR"), isogcr);
+        initEPSKey(newParserDeck, grid, std::string("IKRW"),   ikrw);
+        initEPSKey(newParserDeck, grid, std::string("IKRG"),   ikrg);
+        initEPSKey(newParserDeck, grid, std::string("IKRO"),   ikro);
+        initEPSKey(newParserDeck, grid, std::string("IKRWR"),  ikrwr);
+        initEPSKey(newParserDeck, grid, std::string("IKRGR"),  ikrgr);
+        initEPSKey(newParserDeck, grid, std::string("IKRORW"), ikrorw);
+        initEPSKey(newParserDeck, grid, std::string("IKRORG"), ikrorg);
+
+        eps_transf_hyst_.resize(grid.number_of_cells);
+        sat_hyst_.resize(grid.number_of_cells);
+
+        const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua];
+        const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour];
+        const bool oilWater = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && !phase_usage_.phase_used[Vapour];
+        const bool oilGas = !phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour];
+        const bool threephase = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour];
+
+        for (int cell = 0; cell < grid.number_of_cells; ++cell) {
+            if (oilWater) {
+                // ### krw
+                initEPSParam(cell, eps_transf_hyst_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos],
+                  funcForCell(cell).sowcr_, -1.0, funcForCell(cell).krwr_, funcForCell(cell).krwmax_, iswl, iswcr, iswu, isowcr, isgl, ikrwr, ikrw);
+                // ### krow
+                initEPSParam(cell, eps_transf_hyst_[cell].watoil, true, 0.0, funcForCell(cell).sowcr_, funcForCell(cell).smin_[wpos],
+                  funcForCell(cell).swcr_, -1.0, funcForCell(cell).krorw_, funcForCell(cell).kromax_, iswl, isowcr, iswl, iswcr, isgl, ikrorw, ikro);
+            } else if (oilGas) {
+                // ### krg
+                initEPSParam(cell, eps_transf_hyst_[cell].gas, false, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_, funcForCell(cell).smax_[gpos],
+                  funcForCell(cell).sogcr_, -1.0, funcForCell(cell).krgr_, funcForCell(cell).krgmax_, isgl, isgcr, isgu, isogcr, iswl, ikrgr, ikrg);
+                // ### krog
+                initEPSParam(cell, eps_transf_hyst_[cell].gasoil, true, 0.0, funcForCell(cell).sogcr_, funcForCell(cell).smin_[gpos],
+                  funcForCell(cell).sgcr_, -1.0, funcForCell(cell).krorg_, funcForCell(cell).kromax_, isgl, isogcr, isgl, isgcr, iswl, ikrorg, ikro);
+            } else if (threephase) {
+                // ### krw
+                initEPSParam(cell, eps_transf_hyst_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos], funcForCell(cell).sowcr_,
+                  funcForCell(cell).smin_[gpos], funcForCell(cell).krwr_, funcForCell(cell).krwmax_, iswl, iswcr, iswu, isowcr, isgl, ikrwr, ikrw);
+                // ### krow
+                initEPSParam(cell, eps_transf_hyst_[cell].watoil, true, 0.0, funcForCell(cell).sowcr_, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_,
+                  funcForCell(cell).smin_[gpos], funcForCell(cell).krorw_, funcForCell(cell).kromax_, iswl, isowcr, iswl, iswcr, isgl, ikrorw, ikro);
+                // ### krg
+                initEPSParam(cell, eps_transf_hyst_[cell].gas, false, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_, funcForCell(cell).smax_[gpos], funcForCell(cell).sogcr_,
+                  funcForCell(cell).smin_[wpos], funcForCell(cell).krgr_, funcForCell(cell).krgmax_, isgl, isgcr, isgu, isogcr, iswl, ikrgr, ikrg);
+                // ### krog
+                initEPSParam(cell, eps_transf_hyst_[cell].gasoil, true, 0.0, funcForCell(cell).sogcr_, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_,
+                  funcForCell(cell).smin_[wpos], funcForCell(cell).krorg_, funcForCell(cell).kromax_, isgl, isogcr, isgl, isgcr, iswl, ikrorg, ikro);
+            }
+        }
+    }
+
     // Initialize saturation scaling parameter
     template <class SatFuncSet>
     void SaturationPropsFromDeck<SatFuncSet>::initEPSKey(const EclipseGridParser& deck,
@@ -624,6 +874,213 @@ namespace Opm
         }
     }
 
+    // Initialize saturation scaling parameter
+    template <class SatFuncSet>
+    void SaturationPropsFromDeck<SatFuncSet>::initEPSKey(Opm::DeckConstPtr newParserDeck,
+                                                         const UnstructuredGrid& grid,
+                                                         const std::string& keyword,
+                                                         std::vector<double>& scaleparam)
+    { 
+        const bool useAqua = phase_usage_.phase_used[Aqua];
+        const bool useLiquid = phase_usage_.phase_used[Liquid];
+        const bool useVapour = phase_usage_.phase_used[Vapour];
+        bool useKeyword = newParserDeck->hasKeyword(keyword);
+        bool hasENPTVD = newParserDeck->hasKeyword("ENPTVD");
+        bool hasENKRVD = newParserDeck->hasKeyword("ENKRVD");
+        int itab = 0;
+        std::vector<std::vector<std::vector<double> > > table_dummy;
+        std::vector<std::vector<std::vector<double> > >& table = table_dummy;
+
+        // Active keyword assigned default values for each cell (in case of possible box-wise assignment)
+        int phase_pos_aqua = phase_usage_.phase_pos[BlackoilPhases::Aqua];
+        int phase_pos_vapour = phase_usage_.phase_pos[BlackoilPhases::Vapour];
+        if ((keyword[0] == 'S' && (useKeyword || hasENPTVD)) || (keyword[1] == 'S' && useKeyword) ) {
+            if (keyword == std::string("SWL") || keyword == std::string("ISWL") ) {
+                if (useAqua && (useKeyword || columnIsMasked_(newParserDeck, "ENPTVD", 0))) {
+                    itab = 1;
+                    scaleparam.resize(grid.number_of_cells);
+                    for (int i=0; i<grid.number_of_cells; ++i)
+                        scaleparam[i] = funcForCell(i).smin_[phase_pos_aqua];
+                }
+            } else if (keyword == std::string("SWCR") || keyword == std::string("ISWCR") ) {
+                if (useAqua && (useKeyword || columnIsMasked_(newParserDeck, "ENPTVD", 1))) {
+                    itab = 2;
+                    scaleparam.resize(grid.number_of_cells);
+                    for (int i=0; i<grid.number_of_cells; ++i)
+                        scaleparam[i] = funcForCell(i).swcr_;
+                }
+            } else if (keyword == std::string("SWU") || keyword == std::string("ISWU") ) {
+                if (useAqua && (useKeyword || columnIsMasked_(newParserDeck, "ENPTVD", 2))) {
+                    itab = 3;
+                    scaleparam.resize(grid.number_of_cells);
+                    for (int i=0; i<grid.number_of_cells; ++i)
+                        scaleparam[i] = funcForCell(i).smax_[phase_pos_aqua];
+                }
+            } else if (keyword == std::string("SGL") || keyword == std::string("ISGL") ) {
+                if (useVapour && (useKeyword || columnIsMasked_(newParserDeck, "ENPTVD", 3))) {
+                    itab = 4;
+                    scaleparam.resize(grid.number_of_cells);
+                    for (int i=0; i<grid.number_of_cells; ++i)
+                        scaleparam[i] = funcForCell(i).smin_[phase_pos_vapour];
+                }
+            } else if (keyword == std::string("SGCR") || keyword == std::string("ISGCR") ) {
+                if (useVapour && (useKeyword || columnIsMasked_(newParserDeck, "ENPTVD", 4))) {
+                    itab = 5;
+                    scaleparam.resize(grid.number_of_cells);
+                    for (int i=0; i<grid.number_of_cells; ++i)
+                        scaleparam[i] = funcForCell(i).sgcr_;
+                }
+            } else if (keyword == std::string("SGU") || keyword == std::string("ISGU") ) {
+                if (useVapour && (useKeyword || columnIsMasked_(newParserDeck, "ENPTVD", 5))) {
+                    itab = 6;
+                    scaleparam.resize(grid.number_of_cells);
+                    for (int i=0; i<grid.number_of_cells; ++i)
+                        scaleparam[i] = funcForCell(i).smax_[phase_pos_vapour];
+                }
+            } else if (keyword == std::string("SOWCR") || keyword == std::string("ISOWCR") ) {
+                if (useAqua && (useKeyword || columnIsMasked_(newParserDeck, "ENPTVD", 6))) {
+                    itab = 7;
+                    scaleparam.resize(grid.number_of_cells);
+                    for (int i=0; i<grid.number_of_cells; ++i)
+                        scaleparam[i] = funcForCell(i).sowcr_;
+                }
+            } else if (keyword == std::string("SOGCR") || keyword == std::string("ISOGCR") ) {
+                if (useVapour && (useKeyword || columnIsMasked_(newParserDeck, "ENPTVD", 7))) {
+                    itab = 8;
+                    scaleparam.resize(grid.number_of_cells);
+                    for (int i=0; i<grid.number_of_cells; ++i)
+                        scaleparam[i] = funcForCell(i).sogcr_;
+                }
+            } else {
+                OPM_THROW(std::runtime_error, " -- unknown keyword: '" << keyword << "'");
+            }
+            if (!useKeyword && itab > 0) {
+                Opm::EnptvdTable enptvd(newParserDeck->getKeyword("ENPTVD"));
+                table.resize(1); // only one region supported so far
+                for (unsigned i = 0; i < table.size(); ++i) {
+                    table[i].resize(9);
+
+                    for (int k = 0; k < enptvd.numRows(); ++k) {
+                        table[i][0] = enptvd.getDepthColumn();
+                        table[i][1] = enptvd.getSwcoColumn();
+                        table[i][2] = enptvd.getSwcritColumn();
+                        table[i][3] = enptvd.getSwmaxColumn();
+                        table[i][4] = enptvd.getSgcoColumn();
+                        table[i][5] = enptvd.getSgcritColumn();
+                        table[i][6] = enptvd.getSgmaxColumn();
+                        table[i][7] = enptvd.getSowcritColumn();
+                        table[i][8] = enptvd.getSogcritColumn();
+                    }
+                }
+            }
+        } else if ((keyword[0] == 'K' && (useKeyword || hasENKRVD)) || (keyword[1] == 'K' && useKeyword) ) {
+            if (keyword == std::string("KRW") || keyword == std::string("IKRW") ) {
+                if (useAqua && (useKeyword || columnIsMasked_(newParserDeck, "ENKRVD", 0))) {
+                    itab = 1;
+                    scaleparam.resize(grid.number_of_cells);
+                    for (int i=0; i<grid.number_of_cells; ++i)
+                        scaleparam[i] = funcForCell(i).krwmax_;
+                }
+            } else if (keyword == std::string("KRG") || keyword == std::string("IKRG") ) {
+                if (useVapour && (useKeyword || columnIsMasked_(newParserDeck, "ENKRVD", 1))) {
+                    itab = 2;
+                    scaleparam.resize(grid.number_of_cells);
+                    for (int i=0; i<grid.number_of_cells; ++i)
+                        scaleparam[i] = funcForCell(i).krgmax_;
+                }
+                if (useLiquid && (useKeyword || columnIsMasked_(newParserDeck, "ENKRVD", 2))) {
+                    itab = 3;
+                    scaleparam.resize(grid.number_of_cells);
+                    for (int i=0; i<grid.number_of_cells; ++i)
+                        scaleparam[i] = funcForCell(i).kromax_;
+                }
+            } else if (keyword == std::string("KRWR") || keyword == std::string("IKRWR") ) {
+                if (useAqua && (useKeyword || columnIsMasked_(newParserDeck, "ENKRVD", 3))) {
+                    itab = 4;
+                    scaleparam.resize(grid.number_of_cells);
+                    for (int i=0; i<grid.number_of_cells; ++i)
+                        scaleparam[i] = funcForCell(i).krwr_;
+                }
+            } else if (keyword == std::string("KRGR") || keyword == std::string("IKRGR") ) {
+                if (useVapour && (useKeyword || columnIsMasked_(newParserDeck, "ENKRVD", 4))) {
+                    itab = 5;
+                    scaleparam.resize(grid.number_of_cells);
+                    for (int i=0; i<grid.number_of_cells; ++i)
+                        scaleparam[i] = funcForCell(i).krgr_;
+                }
+            } else if (keyword == std::string("KRORW") || keyword == std::string("IKRORW") ) {
+                if (useAqua && (useKeyword || columnIsMasked_(newParserDeck, "ENKRVD", 5))) {
+                    itab = 6;
+                    scaleparam.resize(grid.number_of_cells);
+                    for (int i=0; i<grid.number_of_cells; ++i)
+                        scaleparam[i] = funcForCell(i).krorw_;
+                }
+            } else if (keyword == std::string("KRORG") || keyword == std::string("IKRORG") ) {
+                if (useVapour && (useKeyword || columnIsMasked_(newParserDeck, "ENKRVD", 6))) {
+                    itab = 7;
+                    scaleparam.resize(grid.number_of_cells);
+                    for (int i=0; i<grid.number_of_cells; ++i)
+                        scaleparam[i] = funcForCell(i).krorg_;
+                }
+            } else {
+                OPM_THROW(std::runtime_error, " -- unknown keyword: '" << keyword << "'");
+            }
+            if (!useKeyword && itab > 0) {
+                Opm::EnkrvdTable enkrvd(newParserDeck->getKeyword("ENKRVD"));
+                table.resize(1); // only one region supported so far
+                for (unsigned i = 0; i < table.size(); ++i) {
+                    table[i].resize(8);
+
+                    for (int k = 0; k < enkrvd.numRows(); ++k) {
+                        table[i][0] = enkrvd.getDepthColumn();
+                        table[i][1] = enkrvd.getKrwmaxColumn();
+                        table[i][2] = enkrvd.getKrgmaxColumn();
+                        table[i][3] = enkrvd.getKromaxColumn();
+                        table[i][4] = enkrvd.getKrwcritColumn();
+                        table[i][5] = enkrvd.getKrgcritColumn();
+                        table[i][6] = enkrvd.getKrocritgColumn();
+                        table[i][7] = enkrvd.getKrocritwColumn();
+                    }
+                }
+            }
+        }
+
+        if (scaleparam.empty()) {
+            return;
+        } else if (useKeyword) {
+            // Keyword values from deck
+            std::cout << "--- Scaling parameter '" << keyword << "' assigned." << std::endl;
+            const int* gc = grid.global_cell;
+            const std::vector<double>& val = newParserDeck->getKeyword(keyword)->getSIDoubleData();
+            for (int c = 0; c < int(scaleparam.size()); ++c) {
+                const int deck_pos = (gc == NULL) ? c : gc[c];
+                scaleparam[c] = val[deck_pos];
+            }
+        } else {
+            // TODO for new parser
+            /*
+            std::cout << "--- Scaling parameter '" << keyword << "' assigned via ";
+            if (keyword[0] == 'S')
+                newParserDeck.getENPTVD().write(std::cout);
+            else
+                newParserDeck.getENKRVD().write(std::cout);
+            */
+            const double* cc = grid.cell_centroids;
+            const int dim = grid.dimensions;
+            for (int cell = 0; cell < grid.number_of_cells; ++cell) {
+                int jtab = cell_to_func_.empty() ? 0 : cell_to_func_[cell];
+                if (table[itab][jtab][0] != -1.0) {
+                    std::vector<double>& depth = table[0][jtab];
+                    std::vector<double>& val = table[itab][jtab];
+                    double zc = cc[dim*cell+dim-1];
+                    if (zc >= depth.front() && zc <= depth.back()) { //don't want extrap outside depth interval
+                        scaleparam[cell] = linearInterpolation(depth, val, zc);
+                    }
+                }
+            }
+        }
+    }
+
     // Saturation scaling
     template <class SatFuncSet>
     void SaturationPropsFromDeck<SatFuncSet>::initEPSParam(const int cell,

opm/core/simulator/SimulatorTimer.cpp

@@ -57,10 +57,21 @@ namespace Opm
         start_date_ = deck.getStartDate();
     }
 
+    /// Use the SimulatorTimer as a shim around opm-parser's Opm::TimeMap
+    void SimulatorTimer::init(Opm::TimeMapConstPtr timeMap,
+                              int timeStepIdx)
+    {
+        timeMap_ = timeMap;
+        current_step_ = timeStepIdx;
+    }
+
     /// Total number of steps.
     int SimulatorTimer::numSteps() const
     {
-        return timesteps_.size();
+        if (timeMap_)
+            return timeMap_->numTimesteps();
+        else
+            return timesteps_.size();
     }
 
     /// Current step number.
@@ -72,13 +83,14 @@ namespace Opm
     /// Set current step number.
     void SimulatorTimer::setCurrentStepNum(int step)
     {
-        if (current_step_ < 0 || current_step_ > int(timesteps_.size())) {
+        if (current_step_ < 0 || current_step_ > int(numSteps())) {
             // Note that we do allow current_step_ == timesteps_.size(),
             // that is the done() state.
             OPM_THROW(std::runtime_error, "Trying to set invalid step number: " << step);
         }
         current_step_ = step;
-        current_time_ = std::accumulate(timesteps_.begin(), timesteps_.begin() + step, 0.0);
+        if (timeMap_)
+            current_time_ = std::accumulate(timesteps_.begin(), timesteps_.begin() + step, 0.0);
     }
 
 
@@ -86,25 +98,37 @@ namespace Opm
     double SimulatorTimer::currentStepLength() const
     {
         assert(!done());
-        return timesteps_[current_step_];
+        if (timeMap_)
+            return timeMap_->getTimeStepLength(current_step_);
+        else
+            return timesteps_[current_step_];
     }
 
     double SimulatorTimer::stepLengthTaken() const
     {
         assert(current_step_ > 0);
-        return timesteps_[current_step_ - 1];
+        if (timeMap_)
+            return timeMap_->getTimeStepLength(current_step_ - 1);
+        else
+            return timesteps_[current_step_ - 1];
     }
 
     /// Current time.
     double SimulatorTimer::currentTime() const
     {
-        return current_time_;
+        if (timeMap_)
+            return timeMap_->getTimePassedUntil(current_step_);
+        else
+            return current_time_;
     }
 
 
     boost::posix_time::ptime SimulatorTimer::currentDateTime() const
     {
-      return boost::posix_time::ptime(start_date_) + boost::posix_time::seconds( (int) current_time_ );
+        if (timeMap_)
+            return timeMap_->getStartTime(current_step_);
+        else
+            return boost::posix_time::ptime(start_date_) + boost::posix_time::seconds( (int) current_time_ );
     }
 
 
@@ -112,7 +136,10 @@ namespace Opm
     /// Total time.
     double SimulatorTimer::totalTime() const
     {
-        return total_time_;
+        if (timeMap_)
+            return timeMap_->getTotalTime();
+        else
+            return total_time_;
     }
 
     /// Set total time.
@@ -121,7 +148,13 @@ namespace Opm
     /// access to later timesteps.
     void SimulatorTimer::setTotalTime(double time)
     {
-        total_time_ = time;
+        if (timeMap_) {
+            // well, what can we do if we use opm-parser's TimeMap?
+            OPM_THROW(std::logic_error,
+                      "Not implemented: SimulatorTimer::setTotalTime() if using a TimeMap.");
+        }
+        else
+            total_time_ = time;
     }
 
     /// Print a report with current and total time etc.
@@ -138,7 +171,8 @@ namespace Opm
     SimulatorTimer& SimulatorTimer::operator++()
     {
         assert(!done());
-        current_time_ += timesteps_[current_step_];
+        if (!timeMap_)
+            current_time_ += timesteps_[current_step_];
         ++current_step_;
         return *this;
     }
@@ -146,7 +180,10 @@ namespace Opm
     /// Return true if op++() has been called numSteps() times.
     bool SimulatorTimer::done() const
     {
-        return int(timesteps_.size()) == current_step_;
+        if (timeMap_)
+            return int(timeMap_->numTimesteps()) <= current_step_;
+        else
+            return int(timesteps_.size()) == current_step_;
     }
 
 

opm/core/simulator/SimulatorTimer.hpp

@@ -20,6 +20,8 @@
 #ifndef OPM_SIMULATORTIMER_HEADER_INCLUDED
 #define OPM_SIMULATORTIMER_HEADER_INCLUDED
 
+#include <opm/parser/eclipse/EclipseState/Schedule/TimeMap.hpp>
+
 #include <iosfwd>
 #include <vector>
 #include <boost/date_time/gregorian/gregorian.hpp>
@@ -47,6 +49,10 @@ namespace Opm
         /// Note that DATES are folded into TSTEP by the parser.
         void init(const EclipseGridParser& deck);
 
+        /// Use the SimulatorTimer as a shim around opm-parser's Opm::TimeMap
+        void init(TimeMapConstPtr timeMap,
+                  int timeStepIdx = 0);
+
         /// Total number of steps.
         int numSteps() const;
 
@@ -99,6 +105,7 @@ namespace Opm
         bool done() const;
 
     private:
+        Opm::TimeMapConstPtr timeMap_;
         std::vector<double> timesteps_;
         int current_step_;
         double current_time_;

opm/core/simulator/initState_impl.hpp

@@ -32,6 +32,8 @@
 #include <opm/core/props/phaseUsageFromDeck.hpp>
 #include <opm/core/utility/miscUtilitiesBlackoil.hpp>
 
+#include <opm/parser/eclipse/Utility/EquilWrapper.hpp>
+
 #include <iostream>
 #include <cmath>
 
@@ -471,6 +473,103 @@ namespace Opm
     }
 
 
+    /// Initialize a state from input deck.
+    template <class Props, class State>
+    void initStateFromDeck(const UnstructuredGrid& grid,
+                           const Props& props,
+                           Opm::DeckConstPtr newParserDeck,
+                           const double gravity,
+                           State& state)
+    {
+        const int num_phases = props.numPhases();
+        const PhaseUsage pu = phaseUsageFromDeck(newParserDeck);
+        if (num_phases != pu.num_phases) {
+            OPM_THROW(std::runtime_error, "initStateFromDeck():  user specified property object with " << num_phases << " phases, "
+                  "found " << pu.num_phases << " phases in deck.");
+        }
+        state.init(grid, num_phases);
+        if (newParserDeck->hasKeyword("EQUIL")) {
+            if (num_phases != 2) {
+                OPM_THROW(std::runtime_error, "initStateFromDeck(): EQUIL-based init currently handling only two-phase scenarios.");
+            }
+            if (pu.phase_used[BlackoilPhases::Vapour]) {
+                OPM_THROW(std::runtime_error, "initStateFromDeck(): EQUIL-based init currently handling only oil-water scenario (no gas).");
+            }
+            // Set saturations depending on oil-water contact.
+            EquilWrapper equil(newParserDeck->getKeyword("EQUIL"));
+            if (equil.numRegions() != 1) {
+                OPM_THROW(std::runtime_error, "initStateFromDeck(): No region support yet.");
+            }
+            const double woc = equil.waterOilContactDepth(0);
+            initWaterOilContact(grid, props, woc, WaterBelow, state);
+            // Set pressure depending on densities and depths.
+            const double datum_z = equil.datumDepth(0);
+            const double datum_p = equil.datumDepthPressure(0);
+            initHydrostaticPressure(grid, props, woc, gravity, datum_z, datum_p, state);
+        } else if (newParserDeck->hasKeyword("PRESSURE")) {
+            // Set saturations from SWAT/SGAS, pressure from PRESSURE.
+            std::vector<double>& s = state.saturation();
+            std::vector<double>& p = state.pressure();
+            const std::vector<double>& p_deck = newParserDeck->getKeyword("PRESSURE")->getSIDoubleData();
+            const int num_cells = grid.number_of_cells;
+            if (num_phases == 2) {
+                if (!pu.phase_used[BlackoilPhases::Aqua]) {
+                    // oil-gas: we require SGAS
+                    if (!newParserDeck->hasKeyword("SGAS")) {
+                        OPM_THROW(std::runtime_error, "initStateFromDeck(): missing SGAS keyword in 2-phase init");
+                    }
+                    const std::vector<double>& sg_deck = newParserDeck->getKeyword("SGAS")->getSIDoubleData();
+                    const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
+                    const int opos = pu.phase_pos[BlackoilPhases::Liquid];
+                    for (int c = 0; c < num_cells; ++c) {
+                        int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
+                        s[2*c + gpos] = sg_deck[c_deck];
+                        s[2*c + opos] = 1.0 - sg_deck[c_deck];
+                        p[c] = p_deck[c_deck];
+                    }
+                } else {
+                    // water-oil or water-gas: we require SWAT
+                    if (!newParserDeck->hasKeyword("SWAT")) {
+                        OPM_THROW(std::runtime_error, "initStateFromDeck(): missing SWAT keyword in 2-phase init");
+                    }
+                    const std::vector<double>& sw_deck = newParserDeck->getKeyword("SWAT")->getSIDoubleData();
+                    const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
+                    const int nwpos = (wpos + 1) % 2;
+                    for (int c = 0; c < num_cells; ++c) {
+                        int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
+                        s[2*c + wpos] = sw_deck[c_deck];
+                        s[2*c + nwpos] = 1.0 - sw_deck[c_deck];
+                        p[c] = p_deck[c_deck];
+                    }
+                }
+            } else if (num_phases == 3) {
+                const bool has_swat_sgas = newParserDeck->hasKeyword("SWAT") && newParserDeck->hasKeyword("SGAS");
+                if (!has_swat_sgas) {
+                    OPM_THROW(std::runtime_error, "initStateFromDeck(): missing SGAS or SWAT keyword in 3-phase init.");
+                }
+                const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
+                const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
+                const int opos = pu.phase_pos[BlackoilPhases::Liquid];
+                const std::vector<double>& sw_deck = newParserDeck->getKeyword("SWAT")->getSIDoubleData();
+                const std::vector<double>& sg_deck = newParserDeck->getKeyword("SGAS")->getSIDoubleData();
+                for (int c = 0; c < num_cells; ++c) {
+                    int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
+                    s[3*c + wpos] = sw_deck[c_deck];
+                    s[3*c + opos] = 1.0 - (sw_deck[c_deck] + sg_deck[c_deck]);
+                    s[3*c + gpos] = sg_deck[c_deck];
+                    p[c] = p_deck[c_deck];
+                 }
+            } else {
+                OPM_THROW(std::runtime_error, "initStateFromDeck(): init with SWAT etc. only available with 2 or 3 phases.");
+            }
+        } else {
+            OPM_THROW(std::runtime_error, "initStateFromDeck(): we must either have EQUIL, or PRESSURE and SWAT/SOIL/SGAS.");
+        }
+
+        // Finally, init face pressures.
+        initFacePressure(grid, state);
+    }
+
     /// Initialize a state from input deck.
     template <class Props, class State>
     void initStateFromDeck(const UnstructuredGrid& grid,
@@ -568,10 +667,6 @@ namespace Opm
         initFacePressure(grid, state);
     }
 
-
-
-
-
     /// Initialize surface volume from pressure and saturation by z = As.
     /// Here  saturation is used as an intial guess for z in the
     /// computation of A.
@@ -770,6 +865,39 @@ namespace Opm
         }
     }
 
+    /// Initialize a blackoil state from input deck.
+    template <class Props, class State>
+    void initBlackoilStateFromDeck(const UnstructuredGrid& grid,
+                                   const Props& props,
+                                   Opm::DeckConstPtr newParserDeck,
+                                   const double gravity,
+                                   State& state)
+    {
+        initStateFromDeck(grid, props, newParserDeck, gravity, state);
+        if (newParserDeck->hasKeyword("RS")) {
+            const std::vector<double>& rs_deck = newParserDeck->getKeyword("RS")->getSIDoubleData();
+            const int num_cells = grid.number_of_cells;
+            for (int c = 0; c < num_cells; ++c) {
+                int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
+                state.gasoilratio()[c] = rs_deck[c_deck];
+            }
+            initBlackoilSurfvolUsingRSorRV(grid, props, state);
+            computeSaturation(props,state);
+        } else if (newParserDeck->hasKeyword("RV")){
+            const std::vector<double>& rv_deck = newParserDeck->getKeyword("RV")->getSIDoubleData();
+            const int num_cells = grid.number_of_cells;
+            for (int c = 0; c < num_cells; ++c) {
+                int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
+                state.rv()[c] = rv_deck[c_deck];
+            }
+            initBlackoilSurfvolUsingRSorRV(grid, props, state);
+            computeSaturation(props,state);
+        }
+        else {
+            OPM_THROW(std::runtime_error, "Temporarily, we require the RS or the RV field.");
+        }
+    }
+
 
 } // namespace Opm
 

opm/core/well_controls.h

@@ -34,10 +34,10 @@ enum WellControlType  {
     SURFACE_RATE      /**< Well constrained by surface volume flow rate */
 };
 
-struct WellControls;
+struct WellControls; 
 
 bool 
-well_controls_equal(const struct WellControls *ctrls1, const struct WellControls *ctrls2);
+well_controls_equal(const struct WellControls *ctrls1, const struct WellControls *ctrls2 , bool verbose);
 
 struct WellControls * 
 well_controls_create(void);

opm/core/wells.h

@@ -263,7 +263,7 @@ struct Wells *
 clone_wells(const struct Wells *W);
 
 bool
-wells_equal(const struct Wells *W1, const struct Wells *W2);
+wells_equal(const struct Wells *W1, const struct Wells *W2 , bool verbose);
 
 
 

opm/core/wells/WellCollection.cpp

@@ -19,10 +19,66 @@
 
 #include "config.h"
 #include <opm/core/wells/WellCollection.hpp>
+
+#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Group.hpp>
+
+#include <boost/lexical_cast.hpp>
+
 #include <memory>
 
 namespace Opm
 {
+    void WellCollection::addField(GroupConstPtr fieldGroup, size_t timeStep, const PhaseUsage& phaseUsage) {
+        WellsGroupInterface* fieldNode = findNode(fieldGroup->name());
+        if (fieldNode) {
+            OPM_THROW(std::runtime_error, "Trying to add FIELD node, but this already exists. Can only have one FIELD node.");
+        }
+
+        roots_.push_back(createGroupWellsGroup(fieldGroup, timeStep, phaseUsage));
+    }
+
+    void WellCollection::addGroup(GroupConstPtr groupChild, std::string parent_name,
+                                  size_t timeStep, const PhaseUsage& phaseUsage) {
+        WellsGroupInterface* parent = findNode(parent_name);
+        if (!parent) {
+            OPM_THROW(std::runtime_error, "Trying to add child group to group named " << parent_name << ", but this does not exist in the WellCollection.");
+        }
+
+        if (findNode(groupChild->name())) {
+            OPM_THROW(std::runtime_error, "Trying to add child group named " << groupChild->name() << ", but this group is already in the WellCollection.");
+
+        }
+
+        std::shared_ptr<WellsGroupInterface> child = createGroupWellsGroup(groupChild, timeStep, phaseUsage);
+
+        WellsGroup* parent_as_group = static_cast<WellsGroup*> (parent);
+        if (!parent_as_group) {
+            OPM_THROW(std::runtime_error, "Trying to add child group to group named " << parent->name() << ", but it's not a group.");
+        }
+        parent_as_group->addChild(child);
+        child->setParent(parent);
+    }
+
+    void WellCollection::addWell(WellConstPtr wellChild, size_t timeStep, const PhaseUsage& phaseUsage) {
+        WellsGroupInterface* parent = findNode(wellChild->getGroupName(timeStep));
+        if (!parent) {
+            OPM_THROW(std::runtime_error, "Trying to add well " << wellChild->name() << " Step: " << boost::lexical_cast<std::string>(timeStep) << " to group named " << wellChild->getGroupName(timeStep) << ", but this group does not exist in the WellCollection.");
+        }
+
+        std::shared_ptr<WellsGroupInterface> child = createWellWellsGroup(wellChild, timeStep, phaseUsage);
+
+        WellsGroup* parent_as_group = static_cast<WellsGroup*> (parent);
+        if (!parent_as_group) {
+            OPM_THROW(std::runtime_error, "Trying to add well to group named " << wellChild->getGroupName(timeStep) << ", but it's not a group.");
+        }
+        parent_as_group->addChild(child);
+
+        leaf_nodes_.push_back(static_cast<WellNode*>(child.get()));
+
+        child->setParent(parent);
+    }
+
 
     void WellCollection::addChild(const std::string& child_name,
                                   const std::string& parent_name,
@@ -33,6 +89,7 @@ namespace Opm
             roots_.push_back(createWellsGroup(parent_name, deck));
             parent = roots_[roots_.size() - 1].get();
         }
+
         std::shared_ptr<WellsGroupInterface> child;
 
         for (size_t i = 0; i < roots_.size(); ++i) {
@@ -47,6 +104,7 @@ namespace Opm
                 break;
             }
         }
+
         if (!child.get()) {
             child = createWellsGroup(child_name, deck);
         }
@@ -65,7 +123,6 @@ namespace Opm
     }
 
 
-
     const std::vector<WellNode*>& WellCollection::getLeafNodes() const {
         return leaf_nodes_;
     }

opm/core/wells/WellCollection.hpp

@@ -23,10 +23,15 @@
 #define	OPM_WELLCOLLECTION_HPP
 
 #include <vector>
+#include <memory>
+
 #include <opm/core/wells/WellsGroup.hpp>
 #include <opm/core/grid.h>
 #include <opm/core/io/eclipse/EclipseGridParser.hpp>
-#include <memory>
+#include <opm/core/props/phaseUsageFromDeck.hpp>
+
+#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Group.hpp>
 
 namespace Opm
 {
@@ -34,6 +39,14 @@ namespace Opm
     class WellCollection
     {
     public:
+
+        void addField(GroupConstPtr fieldGroup, size_t timeStep, const PhaseUsage& phaseUsage);
+
+        void addWell(WellConstPtr wellChild, size_t timeStep, const PhaseUsage& phaseUsage);
+
+        void addGroup(GroupConstPtr groupChild, std::string parent_name,
+                      size_t timeStep, const PhaseUsage& phaseUsage);
+
         /// Adds and creates if necessary the child to the collection
         /// and appends it to parent's children. Also adds and creates the parent
         /// if necessary.

opm/core/wells/WellsGroup.cpp

@@ -74,7 +74,7 @@ namespace Opm
     {
         return parent_;
     }
-    const std::string& WellsGroupInterface::name()
+    const std::string& WellsGroupInterface::name() const
     {
         return name_;
     }
@@ -1141,4 +1141,54 @@ namespace Opm
 
         return return_value;
     }
+
+    std::shared_ptr<WellsGroupInterface> createGroupWellsGroup(GroupConstPtr group, size_t timeStep, const PhaseUsage& phase_usage )
+    {
+        InjectionSpecification injection_specification;
+        ProductionSpecification production_specification;
+        if (group->isInjectionGroup(timeStep)) {
+            injection_specification.injector_type_ = toInjectorType(Phase::PhaseEnum2String(group->getInjectionPhase(timeStep)));
+            injection_specification.control_mode_ = toInjectionControlMode(GroupInjection::ControlEnum2String(group->getInjectionControlMode(timeStep)));
+            injection_specification.surface_flow_max_rate_ = group->getSurfaceMaxRate(timeStep);
+            injection_specification.reservoir_flow_max_rate_ = group->getReservoirMaxRate(timeStep);
+            injection_specification.reinjection_fraction_target_ = group->getTargetReinjectFraction(timeStep);
+            injection_specification.voidage_replacment_fraction_ = group->getTargetVoidReplacementFraction(timeStep);
+        }
+        else if (group->isProductionGroup(timeStep)) {
+            production_specification.oil_max_rate_ = group->getOilTargetRate(timeStep);
+            production_specification.control_mode_ = toProductionControlMode(GroupProduction::ControlEnum2String(group->getProductionControlMode(timeStep)));
+            production_specification.water_max_rate_ = group->getWaterTargetRate(timeStep);
+            production_specification.gas_max_rate_ = group->getGasTargetRate(timeStep);
+            production_specification.liquid_max_rate_ = group->getLiquidTargetRate(timeStep);
+            production_specification.procedure_ = toProductionProcedure(GroupProductionExceedLimit::ActionEnum2String(group->getProductionExceedLimitAction(timeStep)));
+            production_specification.reservoir_flow_max_rate_ = group->getReservoirMaxRate(timeStep);
+        }
+
+        std::shared_ptr<WellsGroupInterface> wells_group(new WellsGroup(group->name(), production_specification, injection_specification, phase_usage));
+        return wells_group;
+    }
+
+    std::shared_ptr<WellsGroupInterface> createWellWellsGroup(WellConstPtr well, size_t timeStep, const PhaseUsage& phase_usage )
+    {
+        InjectionSpecification injection_specification;
+        ProductionSpecification production_specification;
+        if (well->isInjector(timeStep)) {
+            injection_specification.BHP_limit_ = well->getBHPLimit(timeStep);
+            injection_specification.injector_type_ = toInjectorType(WellInjector::Type2String(well->getInjectorType(timeStep)));
+            injection_specification.control_mode_ = toInjectionControlMode(WellInjector::ControlMode2String(well->getInjectorControlMode(timeStep)));
+            injection_specification.surface_flow_max_rate_ = well->getSurfaceInjectionRate(timeStep);
+            injection_specification.reservoir_flow_max_rate_ = well->getReservoirInjectionRate(timeStep);
+            production_specification.guide_rate_ = 0.0; // We know we're not a producer
+        }
+        else if (well->isProducer(timeStep)) {
+            production_specification.BHP_limit_ = well->getBHPLimit(timeStep);
+            production_specification.reservoir_flow_max_rate_ = well->getResVRate(timeStep);
+            production_specification.oil_max_rate_ = well->getOilRate(timeStep);
+            production_specification.control_mode_ = toProductionControlMode(WellProducer::ControlMode2String(well->getProducerControlMode(timeStep)));
+            production_specification.water_max_rate_ = well->getWaterRate(timeStep);
+            injection_specification.guide_rate_ = 0.0; // we know we're not an injector
+        }
+        std::shared_ptr<WellsGroupInterface> wells_group(new WellNode(well->name(), production_specification, injection_specification, phase_usage));
+        return wells_group;
+    }
 }

opm/core/wells/WellsGroup.hpp

@@ -25,6 +25,10 @@
 #include <opm/core/io/eclipse/EclipseGridParser.hpp>
 #include <opm/core/grid.h>
 #include <opm/core/props/BlackoilPhases.hpp>
+
+#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Group.hpp>
+
 #include <string>
 #include <memory>
 
@@ -58,7 +62,7 @@ namespace Opm
         virtual ~WellsGroupInterface();
 
         /// The unique identifier for the well or well group.
-        const std::string& name();
+        const std::string& name() const;
 
         /// Production specifications for the well or well group.
         const ProductionSpecification& prodSpec() const;
@@ -403,9 +407,21 @@ namespace Opm
     /// \param[in] name the name of the wells group.
     /// \param[in] deck the deck from which to fetch information.
     std::shared_ptr<WellsGroupInterface> createWellsGroup(const std::string& name,
-                                                               const EclipseGridParser& deck);
+                                                          const EclipseGridParser& deck);
 
+    /// Creates the WellsGroupInterface for the given well
+    /// \param[in] well the Well to construct object for
+    /// \param[in] timeStep the time step in question
+    /// \param[in] the phase usage
+    std::shared_ptr<WellsGroupInterface> createWellWellsGroup(WellConstPtr well, size_t timeStep,
+                                                              const PhaseUsage& phase_usage );
 
+    /// Creates the WellsGroupInterface for the given Group
+    /// \param[in] group the Group to construct object for
+    /// \param[in] timeStep the time step in question
+    /// \param[in] the phase usage
+    std::shared_ptr<WellsGroupInterface> createGroupWellsGroup(GroupConstPtr group, size_t timeStep,
+                                                               const PhaseUsage& phase_usage );
 }
 #endif	/* OPM_WELLSGROUP_HPP */
 

opm/core/wells/WellsManager.cpp

@@ -18,6 +18,8 @@
 */
 
 #include "config.h"
+
+
 #include <opm/core/wells/WellsManager.hpp>
 #include <opm/core/io/eclipse/EclipseGridParser.hpp>
 #include <opm/core/grid.h>
@@ -28,6 +30,8 @@
 #include <opm/core/wells/WellCollection.hpp>
 #include <opm/core/props/phaseUsageFromDeck.hpp>
 
+#include <opm/parser/eclipse/EclipseState/Schedule/ScheduleEnums.hpp>
+
 #include <array>
 #include <algorithm>
 #include <cassert>
@@ -43,23 +47,8 @@
 namespace
 {
 
-    struct WellData
-    {
-        WellType type;
-        // WellControlType control;
-        // double target;
-        double reference_bhp_depth;
-        // Opm::InjectionSpecification::InjectorType injected_phase;
-        int welspecsline;
-    };
 
 
-    struct PerfData
-    {
-        int cell;
-        double well_index;
-    };
-
     namespace ProductionControl
     {
         enum Mode { ORAT, WRAT, GRAT,
@@ -101,6 +90,33 @@ namespace
                       << control << " in input file");
             }
         }
+
+
+        Mode mode(Opm::WellProducer::ControlModeEnum controlMode)
+        {
+            switch( controlMode ) {
+            case Opm::WellProducer::ORAT:
+                return ORAT;
+            case Opm::WellProducer::WRAT:
+                return WRAT;
+            case Opm::WellProducer::GRAT:
+                return GRAT;
+            case Opm::WellProducer::LRAT:
+                return LRAT;
+            case Opm::WellProducer::CRAT:
+                return CRAT;
+            case Opm::WellProducer::RESV:
+                return RESV;
+            case Opm::WellProducer::BHP:
+                return BHP;
+            case Opm::WellProducer::THP:
+                return THP;
+            case Opm::WellProducer::GRUP:
+                return GRUP;
+            default:
+                throw std::invalid_argument("unhandled enum value");
+            }
+        }
     } // namespace ProductionControl
 
 
@@ -140,6 +156,25 @@ namespace
                       << control << " in input file");
             }
         }
+
+        Mode mode(Opm::WellInjector::ControlModeEnum controlMode)
+        {
+            switch ( controlMode  ) {
+            case Opm::WellInjector::GRUP:
+                return GRUP;
+            case Opm::WellInjector::RESV:
+                return RESV;
+            case Opm::WellInjector::RATE:
+                return RATE;
+            case Opm::WellInjector::THP:
+                return THP;
+            case Opm::WellInjector::BHP:
+                return BHP;
+            default:
+                throw std::invalid_argument("unhandled enum value");
+            }
+        }
+
     } // namespace InjectionControl
 
 
@@ -233,6 +268,85 @@ namespace Opm
     {
     }
 
+    /// Construct wells from deck.
+    WellsManager::WellsManager(const Opm::EclipseStateConstPtr eclipseState,
+                               const size_t timeStep,
+                               const UnstructuredGrid& grid,
+                               const double* permeability)
+        : w_(0)
+    {
+        if (grid.dimensions != 3) {
+            OPM_THROW(std::runtime_error, "We cannot initialize wells from a deck unless the corresponding grid is 3-dimensional.");
+        }
+
+        if (eclipseState->getSchedule()->numWells() == 0) {
+            OPM_MESSAGE("No wells specified in Schedule section, initializing no wells");
+            return;
+        }
+
+        std::map<int,int> cartesian_to_compressed;
+        setupCompressedToCartesian(grid, cartesian_to_compressed);
+
+        // Obtain phase usage data.
+        PhaseUsage pu = phaseUsageFromDeck(eclipseState);
+
+        // These data structures will be filled in this constructor,
+        // then used to initialize the Wells struct.
+        std::vector<std::string> well_names;
+        std::vector<WellData> well_data;
+
+
+        // For easy lookup:
+        std::map<std::string, int> well_names_to_index;
+
+        ScheduleConstPtr schedule = eclipseState->getSchedule();
+        std::vector<WellConstPtr> wells = schedule->getWells(timeStep);
+
+        well_names.reserve(wells.size());
+        well_data.reserve(wells.size());
+
+        createWellsFromSpecs(wells, timeStep, grid, well_names, well_data, well_names_to_index, pu, cartesian_to_compressed, permeability);
+        setupWellControls(wells, timeStep, well_names, pu);
+
+        {
+            GroupTreeNodeConstPtr fieldNode = eclipseState->getSchedule()->getGroupTree(timeStep)->getNode("FIELD");
+            GroupConstPtr fieldGroup = eclipseState->getSchedule()->getGroup(fieldNode->name());
+            well_collection_.addField(fieldGroup, timeStep, pu);
+            addChildGroups(fieldNode, eclipseState->getSchedule(), timeStep, pu);
+        }
+
+        for (auto wellIter = wells.begin(); wellIter != wells.end(); ++wellIter ) {
+            well_collection_.addWell((*wellIter), timeStep, pu);
+        }
+
+        well_collection_.setWellsPointer(w_);
+        well_collection_.applyGroupControls();
+
+
+        setupGuideRates(wells, timeStep, well_data, well_names_to_index);
+
+        // Debug output.
+#define EXTRA_OUTPUT
+#ifdef EXTRA_OUTPUT
+        /*
+        std::cout << "\t WELL DATA" << std::endl;
+        for(int i = 0; i< num_wells; ++i) {
+            std::cout << i << ": " << well_data[i].type << "  "
+                      << well_data[i].control << "  " << well_data[i].target
+                      << std::endl;
+        }
+
+        std::cout << "\n\t PERF DATA" << std::endl;
+        for(int i=0; i< int(wellperf_data.size()); ++i) {
+            for(int j=0; j< int(wellperf_data[i].size()); ++j) {
+                std::cout << i << ": " << wellperf_data[i][j].cell << "  "
+                          << wellperf_data[i][j].well_index << std::endl;
+            }
+        }
+        */
+#endif
+    }
+
 
 
     /// Construct wells from deck.
@@ -879,4 +993,404 @@ namespace Opm
         well_collection_.applyExplicitReinjectionControls(well_reservoirrates_phase, well_surfacerates_phase);
     }
 
+    void WellsManager::setupCompressedToCartesian(const UnstructuredGrid& grid, std::map<int,int>& cartesian_to_compressed ) {
+        // global_cell is a map from compressed cells to Cartesian grid cells.
+        // We must make the inverse lookup.
+        const int* global_cell = grid.global_cell;
+
+        if (global_cell) {
+            for (int i = 0; i < grid.number_of_cells; ++i) {
+                cartesian_to_compressed.insert(std::make_pair(global_cell[i], i));
+            }
+        }
+        else {
+            for (int i = 0; i < grid.number_of_cells; ++i) {
+                cartesian_to_compressed.insert(std::make_pair(i, i));
+            }
+        }
+
+    }
+
+    void WellsManager::createWellsFromSpecs(std::vector<WellConstPtr>& wells, size_t timeStep,
+                                            const UnstructuredGrid& grid,
+                                            std::vector<std::string>& well_names,
+                                            std::vector<WellData>& well_data,
+                                            std::map<std::string, int>& well_names_to_index,
+                                            const PhaseUsage& phaseUsage,
+                                            std::map<int,int> cartesian_to_compressed,
+                                            const double* permeability)
+    {
+        std::vector<std::vector<PerfData> > wellperf_data;
+        wellperf_data.resize(wells.size());
+
+        int well_index = 0;
+        for (auto wellIter= wells.begin(); wellIter != wells.end(); ++wellIter) {
+            WellConstPtr well = (*wellIter);
+            {   // WELSPECS handling
+                well_names_to_index[well->name()] = well_index;
+                well_names.push_back(well->name());
+                {
+                    WellData wd;
+                    // If negative (defaulted), set refdepth to a marker
+                    // value, will be changed after getting perforation
+                    // data to the centroid of the cell of the top well
+                    // perforation.
+                    wd.reference_bhp_depth = (well->getRefDepth() < 0.0) ? -1e100 : well->getRefDepth();
+                    wd.welspecsline = -1;
+                    if (well->isInjector( timeStep ))
+                        wd.type = INJECTOR;
+                    else
+                        wd.type = PRODUCER;
+                    well_data.push_back(wd);
+                }
+            }
+
+            {   // COMPDAT handling
+                CompletionSetConstPtr completionSet = well->getCompletions(timeStep);
+                for (size_t c=0; c<completionSet->size(); c++) {
+                    CompletionConstPtr completion = completionSet->get(c);
+                    int i = completion->getI();
+                    int j = completion->getJ();
+                    int k = completion->getK();
+
+                    const int* cpgdim = grid.cartdims;
+                    int cart_grid_indx = i + cpgdim[0]*(j + cpgdim[1]*k);
+                    std::map<int, int>::const_iterator cgit = cartesian_to_compressed.find(cart_grid_indx);
+                    if (cgit == cartesian_to_compressed.end()) {
+                        OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << i << ' ' << j << ' '
+                                  << k << " not found in grid (well = " << well->name() << ')');
+                    }
+                    int cell = cgit->second;
+                    PerfData pd;
+                    pd.cell = cell;
+                    if (completion->getCF() > 0.0) {
+                        pd.well_index = completion->getCF();
+                    } else {
+                        double radius = 0.5*completion->getDiameter();
+                        if (radius <= 0.0) {
+                            radius = 0.5*unit::feet;
+                            OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius);
+                        }
+                        std::array<double, 3> cubical = getCubeDim(grid, cell);
+                        const double* cell_perm = &permeability[grid.dimensions*grid.dimensions*cell];
+                        pd.well_index = computeWellIndex(radius, cubical, cell_perm, completion->getDiameter());
+                    }
+                    wellperf_data[well_index].push_back(pd);
+                }
+            }
+            well_index++;
+        }
+
+        // Set up reference depths that were defaulted. Count perfs.
+
+        const int num_wells = well_data.size();
+
+        int num_perfs = 0;
+        assert(grid.dimensions == 3);
+        for (int w = 0; w < num_wells; ++w) {
+            num_perfs += wellperf_data[w].size();
+            if (well_data[w].reference_bhp_depth < 0.0) {
+                // It was defaulted. Set reference depth to minimum perforation depth.
+                double min_depth = 1e100;
+                int num_wperfs = wellperf_data[w].size();
+                for (int perf = 0; perf < num_wperfs; ++perf) {
+                    double depth = grid.cell_centroids[3*wellperf_data[w][perf].cell + 2];
+                    min_depth = std::min(min_depth, depth);
+                }
+                well_data[w].reference_bhp_depth = min_depth;
+            }
+        }
+
+        // Create the well data structures.
+        w_ = create_wells(phaseUsage.num_phases, num_wells, num_perfs);
+        if (!w_) {
+            OPM_THROW(std::runtime_error, "Failed creating Wells struct.");
+        }
+
+
+        // Add wells.
+        for (int w = 0; w < num_wells; ++w) {
+            const int w_num_perf = wellperf_data[w].size();
+            std::vector<int> perf_cells(w_num_perf);
+            std::vector<double> perf_prodind(w_num_perf);
+            for (int perf = 0; perf < w_num_perf; ++perf) {
+                perf_cells[perf] = wellperf_data[w][perf].cell;
+                perf_prodind[perf] = wellperf_data[w][perf].well_index;
+            }
+            const double* comp_frac = NULL;
+            // We initialize all wells with a null component fraction,
+            // and must (for injection wells) overwrite it later.
+            int ok = add_well(well_data[w].type, well_data[w].reference_bhp_depth, w_num_perf,
+                              comp_frac, &perf_cells[0], &perf_prodind[0], well_names[w].c_str(), w_);
+            if (!ok) {
+                OPM_THROW(std::runtime_error, "Failed adding well " << well_names[w] << " to Wells data structure.");
+            }
+        }
+
+    }
+
+    void WellsManager::setupWellControls(std::vector<WellConstPtr>& wells, size_t timeStep,
+                                         std::vector<std::string>& well_names, const PhaseUsage& phaseUsage) {
+        int well_index = 0;
+        for (auto wellIter= wells.begin(); wellIter != wells.end(); ++wellIter) {
+            WellConstPtr well = (*wellIter);
+
+            if ( !( well->getStatus( timeStep ) == WellCommon::SHUT || well->getStatus( timeStep ) == WellCommon::OPEN) ) {
+                OPM_THROW(std::runtime_error, "Currently we do not support well status " << WellCommon::Status2String(well->getStatus( timeStep )));
+            }
+
+            if (well->isInjector(timeStep)) {
+                clear_well_controls(well_index, w_);
+                int ok = 1;
+                int control_pos[5] = { -1, -1, -1, -1, -1 };
+
+                if (well->hasInjectionControl(timeStep , WellInjector::RATE)) {
+                    control_pos[InjectionControl::RATE] = well_controls_get_num(w_->ctrls[well_index]);
+                    double distr[3] = { 0.0, 0.0, 0.0 };
+                    WellInjector::TypeEnum injectorType = well->getInjectorType(timeStep);
+
+                    if (injectorType == WellInjector::TypeEnum::WATER) {
+                        distr[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0;
+                    } else if (injectorType == WellInjector::TypeEnum::OIL) {
+                        distr[phaseUsage.phase_pos[BlackoilPhases::Liquid]] = 1.0;
+                    } else if (injectorType == WellInjector::TypeEnum::GAS) {
+                        distr[phaseUsage.phase_pos[BlackoilPhases::Vapour]] = 1.0;
+                    }
+
+                    ok = append_well_controls(SURFACE_RATE,
+                                              well->getSurfaceInjectionRate( timeStep ) ,
+                                              distr,
+                                              well_index,
+                                              w_);
+                }
+
+                if (ok && well->hasInjectionControl(timeStep , WellInjector::RESV)) {
+                    control_pos[InjectionControl::RESV] = well_controls_get_num(w_->ctrls[well_index]);
+                    double distr[3] = { 0.0, 0.0, 0.0 };
+                    WellInjector::TypeEnum injectorType = well->getInjectorType(timeStep);
+
+                    if (injectorType == WellInjector::TypeEnum::WATER) {
+                        distr[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0;
+                    } else if (injectorType == WellInjector::TypeEnum::OIL) {
+                        distr[phaseUsage.phase_pos[BlackoilPhases::Liquid]] = 1.0;
+                    } else if (injectorType == WellInjector::TypeEnum::GAS) {
+                        distr[phaseUsage.phase_pos[BlackoilPhases::Vapour]] = 1.0;
+                    }
+
+                    ok = append_well_controls(RESERVOIR_RATE,
+                                              well->getReservoirInjectionRate( timeStep ),
+                                              distr,
+                                              well_index,
+                                              w_);
+                }
+
+                if (ok && well->hasInjectionControl(timeStep , WellInjector::BHP)) {
+                    control_pos[InjectionControl::BHP] = well_controls_get_num(w_->ctrls[well_index]);
+                    ok = append_well_controls(BHP,
+                                              well->getBHPLimit(timeStep),
+                                              NULL,
+                                              well_index,
+                                              w_);
+                }
+
+                if (ok && well->hasInjectionControl(timeStep , WellInjector::THP)) {
+                    OPM_THROW(std::runtime_error, "We cannot handle THP limit for well " << well_names[well_index]);
+                }
+
+
+                if (!ok) {
+                    OPM_THROW(std::runtime_error, "Failure occured appending controls for well " << well_names[well_index]);
+                }
+
+
+                {
+                    InjectionControl::Mode mode = InjectionControl::mode( well->getInjectorControlMode(timeStep) );
+                    int cpos = control_pos[mode];
+                    if (cpos == -1 && mode != InjectionControl::GRUP) {
+                        OPM_THROW(std::runtime_error, "Control not specified in well " << well_names[well_index]);
+                    }
+
+                    // We need to check if the well is shut or not
+                    if (well->getStatus( timeStep ) == WellCommon::SHUT) {
+                        cpos = ~cpos;
+                    }
+                    set_current_control(well_index, cpos, w_);
+                }
+
+
+                // Set well component fraction.
+                double cf[3] = { 0.0, 0.0, 0.0 };
+                if (well->getInjectorType(timeStep) == WellInjector::WATER) {
+                    if (!phaseUsage.phase_used[BlackoilPhases::Aqua]) {
+                        OPM_THROW(std::runtime_error, "Water phase not used, yet found water-injecting well.");
+                    }
+                    cf[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0;
+                } else if (well->getInjectorType(timeStep) == WellInjector::OIL) {
+                    if (!phaseUsage.phase_used[BlackoilPhases::Liquid]) {
+                        OPM_THROW(std::runtime_error, "Oil phase not used, yet found oil-injecting well.");
+                    }
+                    cf[phaseUsage.phase_pos[BlackoilPhases::Liquid]] = 1.0;
+                } else if (well->getInjectorType(timeStep) == WellInjector::GAS) {
+                    if (!phaseUsage.phase_used[BlackoilPhases::Vapour]) {
+                        OPM_THROW(std::runtime_error, "Gas phase not used, yet found gas-injecting well.");
+                    }
+                    cf[phaseUsage.phase_pos[BlackoilPhases::Vapour]] = 1.0;
+                }
+                std::copy(cf, cf + phaseUsage.num_phases, w_->comp_frac + well_index*phaseUsage.num_phases);
+
+            }
+
+            if (well->isProducer(timeStep)) {
+                // Add all controls that are present in well.
+                // First we must clear existing controls, in case the
+                // current WCONPROD line is modifying earlier controls.
+                clear_well_controls(well_index, w_);
+                int control_pos[9] = { -1, -1, -1, -1, -1, -1, -1, -1, -1 };
+                int ok = 1;
+                if (ok && well->hasProductionControl(timeStep , WellProducer::ORAT)) {
+                    if (!phaseUsage.phase_used[BlackoilPhases::Liquid]) {
+                        OPM_THROW(std::runtime_error, "Oil phase not active and ORAT control specified.");
+                    }
+
+                    control_pos[ProductionControl::ORAT] = well_controls_get_num(w_->ctrls[well_index]);
+                    double distr[3] = { 0.0, 0.0, 0.0 };
+                    distr[phaseUsage.phase_pos[BlackoilPhases::Liquid]] = 1.0;
+                    ok = append_well_controls(SURFACE_RATE,
+                                              -well->getOilRate( timeStep ),
+                                              distr,
+                                               well_index,
+                                              w_);
+                }
+
+                if (ok && well->hasProductionControl(timeStep , WellProducer::WRAT)) {
+                    if (!phaseUsage.phase_used[BlackoilPhases::Aqua]) {
+                        OPM_THROW(std::runtime_error, "Water phase not active and WRAT control specified.");
+                    }
+                    control_pos[ProductionControl::WRAT] = well_controls_get_num(w_->ctrls[well_index]);
+                    double distr[3] = { 0.0, 0.0, 0.0 };
+                    distr[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0;
+                    ok = append_well_controls(SURFACE_RATE,
+                                              -well->getWaterRate(timeStep),
+                                              distr,
+                                              well_index,
+                                              w_);
+                }
+
+                if (ok && well->hasProductionControl(timeStep , WellProducer::GRAT)) {
+                    if (!phaseUsage.phase_used[BlackoilPhases::Vapour]) {
+                        OPM_THROW(std::runtime_error, "Gas phase not active and GRAT control specified.");
+                    }
+                    control_pos[ProductionControl::GRAT] = well_controls_get_num(w_->ctrls[well_index]);
+                    double distr[3] = { 0.0, 0.0, 0.0 };
+                    distr[phaseUsage.phase_pos[BlackoilPhases::Vapour]] = 1.0;
+                    ok = append_well_controls(SURFACE_RATE,
+                                              -well->getGasRate( timeStep ),
+                                              distr,
+                                              well_index,
+                                              w_);
+                }
+
+                if (ok && well->hasProductionControl(timeStep , WellProducer::LRAT)) {
+                    if (!phaseUsage.phase_used[BlackoilPhases::Aqua]) {
+                        OPM_THROW(std::runtime_error, "Water phase not active and LRAT control specified.");
+                    }
+                    if (!phaseUsage.phase_used[BlackoilPhases::Liquid]) {
+                        OPM_THROW(std::runtime_error, "Oil phase not active and LRAT control specified.");
+                    }
+                    control_pos[ProductionControl::LRAT] = well_controls_get_num(w_->ctrls[well_index]);
+                    double distr[3] = { 0.0, 0.0, 0.0 };
+                    distr[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0;
+                    distr[phaseUsage.phase_pos[BlackoilPhases::Liquid]] = 1.0;
+                    ok = append_well_controls(SURFACE_RATE,
+                                              -well->getLiquidRate(timeStep),
+                                              distr,
+                                              well_index,
+                                              w_);
+                }
+
+                if (ok && well->hasProductionControl(timeStep , WellProducer::RESV)) {
+                    control_pos[ProductionControl::RESV] = well_controls_get_num(w_->ctrls[well_index]);
+                    double distr[3] = { 1.0, 1.0, 1.0 };
+                    ok = append_well_controls(RESERVOIR_RATE,
+                                              -well->getResVRate(timeStep),
+                                              distr,
+                                              well_index,
+                                              w_);
+                }
+
+                if (ok && well->hasProductionControl(timeStep , WellProducer::BHP)) {
+                    control_pos[ProductionControl::BHP] = well_controls_get_num(w_->ctrls[well_index]);
+                    ok = append_well_controls(BHP,
+                                              well->getBHPLimit( timeStep ) ,
+                                              NULL,
+                                              well_index,
+                                              w_);
+                }
+
+                if (ok && well->hasProductionControl(timeStep , WellProducer::THP)) {
+                    OPM_THROW(std::runtime_error, "We cannot handle THP limit for well " << well_names[well_index]);
+                }
+
+                if (!ok) {
+                    OPM_THROW(std::runtime_error, "Failure occured appending controls for well " << well_names[well_index]);
+                }
+
+                ProductionControl::Mode mode = ProductionControl::mode(well->getProducerControlMode(timeStep));
+                int cpos = control_pos[mode];
+                if (cpos == -1 && mode != ProductionControl::GRUP) {
+                    OPM_THROW(std::runtime_error, "Control mode type " << mode << " not present in well " << well_names[well_index]);
+                }
+                // If it's shut, we complement the cpos
+                if (well->getStatus(timeStep) == WellCommon::SHUT) {
+                    cpos = ~cpos; // So we can easily retrieve the cpos later
+                }
+                set_current_control(well_index, cpos, w_);
+            }
+            well_index++;
+        }
+
+    }
+
+    void WellsManager::addChildGroups(GroupTreeNodeConstPtr parentNode, ScheduleConstPtr schedule, size_t timeStep, const PhaseUsage& phaseUsage) {
+        for (auto childIter = parentNode->begin(); childIter != parentNode->end(); ++childIter) {
+            GroupTreeNodeConstPtr childNode = (*childIter).second;
+            well_collection_.addGroup(schedule->getGroup(childNode->name()), parentNode->name(), timeStep, phaseUsage);
+            addChildGroups(childNode, schedule, timeStep, phaseUsage);
+        }
+    }
+
+
+
+
+    void WellsManager::setupGuideRates(std::vector<WellConstPtr>& wells, const size_t timeStep, std::vector<WellData>& well_data, std::map<std::string, int>& well_names_to_index)
+    {
+        for (auto wellIter = wells.begin(); wellIter != wells.end(); ++wellIter ) {
+            WellConstPtr well = *wellIter;
+            const int wix = well_names_to_index[well->name()];
+            WellNode& wellnode = *well_collection_.getLeafNodes()[wix];
+
+            if (well->getGuideRatePhase(timeStep) != GuideRate::UNDEFINED) {
+                if (well_data[wix].type == PRODUCER) {
+                    wellnode.prodSpec().guide_rate_ = well->getGuideRate(timeStep);
+                    if (well->getGuideRatePhase(timeStep) == GuideRate::OIL) {
+                        wellnode.prodSpec().guide_rate_type_ = ProductionSpecification::OIL;
+                    } else {
+                        OPM_THROW(std::runtime_error, "Guide rate type " << GuideRate::GuideRatePhaseEnum2String(well->getGuideRatePhase(timeStep)) << " specified for producer "
+                                  << well->name() << " in WGRUPCON, cannot handle.");
+                    }
+                } else if (well_data[wix].type == INJECTOR) {
+                    wellnode.injSpec().guide_rate_ = well->getGuideRate(timeStep);
+                    if (well->getGuideRatePhase(timeStep) == GuideRate::RAT) {
+                        wellnode.injSpec().guide_rate_type_ = InjectionSpecification::RAT;
+                    } else {
+                        OPM_THROW(std::runtime_error, "Guide rate type " << GuideRate::GuideRatePhaseEnum2String(well->getGuideRatePhase(timeStep)) << " specified for injector "
+                                  << well->name() << " in WGRUPCON, cannot handle.");
+                    }
+                } else {
+                    OPM_THROW(std::runtime_error, "Unknown well type " << well_data[wix].type << " for well " << well->name());
+                }
+            }
+        }
+    }
+
 } // namespace Opm

opm/core/wells/WellsManager.hpp

@@ -21,8 +21,11 @@
 #define OPM_WELLSMANAGER_HEADER_INCLUDED
 
 
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+
 #include <opm/core/wells/WellCollection.hpp>
 #include <opm/core/wells/WellsGroup.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/GroupTree.hpp>
 
 struct Wells;
 struct UnstructuredGrid;
@@ -32,7 +35,22 @@ namespace Opm
 {
 
     class EclipseGridParser;
+    struct WellData
+    {
+        WellType type;
+        // WellControlType control;
+        // double target;
+        double reference_bhp_depth;
+        // Opm::InjectionSpecification::InjectorType injected_phase;
+        int welspecsline;
+    };
 
+
+    struct PerfData
+    {
+        int cell;
+        double well_index;
+    };
     /// This class manages a Wells struct in the sense that it
     /// encapsulates creation and destruction of the wells
     /// data structure.
@@ -40,34 +58,40 @@ namespace Opm
     class WellsManager
     {
     public:
-	/// Default constructor -- no wells.
+        /// Default constructor -- no wells.
-    WellsManager();
+        WellsManager();
 
-    /// Construct from existing wells object.
+        /// Construct from existing wells object.
-    /// WellsManager is not properly initialised in the sense that the logic to
+        /// WellsManager is not properly initialised in the sense that the logic to
-    /// manage control switching does not exist.
+        /// manage control switching does not exist.
-    ///
+        ///
-    /// @param[in] W Existing wells object.
+        /// @param[in] W Existing wells object.
-    WellsManager(struct Wells* W);
+        WellsManager(struct Wells* W);
 
-	/// Construct from input deck and grid.
+        /// Construct from input deck and grid.
-	/// The permeability argument may be zero if the input contain
+        /// The permeability argument may be zero if the input contain
-	/// well productivity indices, otherwise it must be given in
+        /// well productivity indices, otherwise it must be given in
-	/// order to approximate these by the Peaceman formula.
+        /// order to approximate these by the Peaceman formula.
-	WellsManager(const Opm::EclipseGridParser& deck,
+        WellsManager(const Opm::EclipseGridParser& deck,
-		     const UnstructuredGrid& grid,
+                     const UnstructuredGrid& grid,
-		     const double* permeability);
+                     const double* permeability);
 
-	/// Destructor.
+
-	~WellsManager();
+        WellsManager(const Opm::EclipseStateConstPtr eclipseState,
+                     const size_t timeStep,
+                     const UnstructuredGrid& grid,
+                     const double* permeability);
+
+        /// Destructor.
+        ~WellsManager();
 
         /// Does the "deck" define any wells?
         bool empty() const;
 
-	/// Access the managed Wells.
+        /// Access the managed Wells.
-	/// The method is named similarly to c_str() in std::string,
+        /// The method is named similarly to c_str() in std::string,
-	/// to make it clear that we are returning a C-compatible struct.
+        /// to make it clear that we are returning a C-compatible struct.
-	const Wells* c_wells() const;
+        const Wells* c_wells() const;
 
         /// Access the well group hierarchy.
         const WellCollection& wellCollection() const;
@@ -108,18 +132,31 @@ namespace Opm
         void applyExplicitReinjectionControls(const std::vector<double>& well_reservoirrates_phase,
                                               const std::vector<double>& well_surfacerates_phase);
 
+
     private:
-	// Disable copying and assignment.
+        // Disable copying and assignment.
-	WellsManager(const WellsManager& other);
+        WellsManager(const WellsManager& other);
-	WellsManager& operator=(const WellsManager& other);
+        WellsManager& operator=(const WellsManager& other);
+        static void setupCompressedToCartesian(const UnstructuredGrid& grid, std::map<int,int>& cartesian_to_compressed );
+        void setupWellControls(std::vector<WellConstPtr>& wells, size_t timeStep,
+                               std::vector<std::string>& well_names, const PhaseUsage& phaseUsage);
 
-	// Data
+        void createWellsFromSpecs( std::vector<WellConstPtr>& wells, size_t timeStep,
-	Wells* w_;
+                                   const UnstructuredGrid& grid,
+                                   std::vector<std::string>& well_names,
+                                   std::vector<WellData>& well_data,
+                                   std::map<std::string, int> & well_names_to_index,
+                                   const PhaseUsage& phaseUsage,
+                                   const std::map<int,int> cartesian_to_compressed,
+                                   const double* permeability);
+
+        void addChildGroups(GroupTreeNodeConstPtr parentNode, ScheduleConstPtr schedule, size_t timeStep, const PhaseUsage& phaseUsage);
+        void setupGuideRates(std::vector<WellConstPtr>& wells, const size_t timeStep, std::vector<WellData>& well_data, std::map<std::string, int>& well_names_to_index);
+
+
+        // Data
+        Wells* w_;
         WellCollection well_collection_;
-
-
-
-
     };
 
 } // namespace Opm

opm/core/wells/well_controls.c

@@ -293,19 +293,38 @@ well_controls_add_new(enum WellControlType type , double target , const double *
 
 
 bool
-well_controls_equal(const struct WellControls *ctrls1, const struct WellControls *ctrls2)
+well_controls_equal(const struct WellControls *ctrls1, const struct WellControls *ctrls2 , bool verbose)
 /* ---------------------------------------------------------------------- */
 {
     bool are_equal = true;
-    are_equal = (ctrls1->num == ctrls2->num);
+
-    are_equal = are_equal && (ctrls1->number_of_phases == ctrls2->number_of_phases);
+    if (ctrls1->num !=  ctrls2->num) {
+        are_equal = false;
+        if (verbose)
+            printf("ctrls1->num:%d    ctrls2->num:%d \n",ctrls1->num , ctrls2->num);
+    }
+
+    if (ctrls1->number_of_phases !=  ctrls2->number_of_phases) {
+        are_equal = false;
+        if (verbose)
+            printf("ctrls1->number_of_phases:%d    ctrls2->number_of_phases:%d \n",ctrls1->number_of_phases , ctrls2->number_of_phases);
+    }
+
     if (!are_equal) {
         return are_equal;
     }
 
-    are_equal = are_equal && (memcmp(ctrls1->type, ctrls2->type, ctrls1->num * sizeof *ctrls1->type ) == 0);
+    if (memcmp(ctrls1->type, ctrls2->type, ctrls1->num * sizeof *ctrls1->type ) != 0) {
-    are_equal = are_equal && (memcmp(ctrls1->target, ctrls2->target, ctrls1->num * sizeof *ctrls1->target ) == 0);
+        are_equal = false;
-    are_equal = are_equal && (memcmp(ctrls1->distr, ctrls2->distr, ctrls1->num * ctrls1->number_of_phases * sizeof *ctrls1->distr ) == 0);
+        if (verbose) 
+            printf("The ->type vectors are different \n");
+    }
+
+    if (memcmp(ctrls1->target, ctrls2->target, ctrls1->num * sizeof *ctrls1->target ) != 0) {
+        are_equal = false;
+        if (verbose) 
+            printf("The ->target vectors are different \n");
+    }
     are_equal = are_equal && (ctrls1->cpty == ctrls2->cpty);
 
     return are_equal;

opm/core/wells/wells.c

@@ -541,7 +541,7 @@ clone_wells(const struct Wells *W)
 
 /* ---------------------------------------------------------------------- */
 bool
-wells_equal(const struct Wells *W1, const struct Wells *W2)
+wells_equal(const struct Wells *W1, const struct Wells *W2 , bool verbose)
 /* ---------------------------------------------------------------------- */
 {
     bool are_equal = true;
@@ -567,10 +567,25 @@ wells_equal(const struct Wells *W1, const struct Wells *W2)
                 are_equal = are_equal && (strcmp(W1->name[i], W2->name[i]) == 0);
             else 
                 are_equal = are_equal && (W1->name[i] == W2->name[i]);  
+
+            if (verbose && !are_equal) 
+                printf("Well name[%d] %s and %s are different \n",  i , W1->name[i] ,  W2->name[i]);
+        }
+        if (W1->type[i] != W2->type[i]) {
+            are_equal = false;
+            if (verbose)
+                printf("Well->type[%d] different %d %d \n",i , W1->type[i] , W2->type[i] );
+        }
+        if (W1->depth_ref[i] != W2->depth_ref[i]) {
+            are_equal = false;
+            if (verbose)
+                printf("Well->depth_ref[%d] different %g %g \n",i , W1->depth_ref[i] , W2->depth_ref[i] );
+        }
+        if (!well_controls_equal(W1->ctrls[i], W2->ctrls[i],verbose)) {
+            are_equal = false;
+            if (verbose)
+                printf("Well controls are different for well[%d]:%s \n",i,W1->name[i]);
         }
-        are_equal = are_equal && (W1->type[i] == W2->type[i]);
-        are_equal = are_equal && (W1->depth_ref[i] == W2->depth_ref[i]);
-        are_equal = are_equal && (well_controls_equal(W1->ctrls[i], W2->ctrls[i]));
     }
 
 
@@ -581,10 +596,16 @@ wells_equal(const struct Wells *W1, const struct Wells *W2)
         are_equal = are_equal && (mgmt1->well_cpty == mgmt2->well_cpty);
     }
 
-    are_equal = are_equal && (memcmp(W1->comp_frac, W2->comp_frac, W1->number_of_wells * W1->number_of_phases * sizeof *W1->comp_frac ) == 0);
+    if (memcmp(W1->comp_frac, W2->comp_frac, W1->number_of_wells * W1->number_of_phases * sizeof *W1->comp_frac ) != 0) {
-    are_equal = are_equal && (memcmp(W1->well_connpos, W2->well_connpos, (1 + W1->number_of_wells) * sizeof *W1->well_connpos ) == 0);
+        are_equal = false;
-    if (!are_equal) {
+        if (verbose)
-        return are_equal;
+            printf("Component fractions different \n");
+    }
+    
+    if (memcmp(W1->well_connpos, W2->well_connpos, (1 + W1->number_of_wells) * sizeof *W1->well_connpos ) != 0) {
+        are_equal = false;
+        if (verbose)
+            printf("perforation position map difference \n");
     }
 
     {

tests/testBlackoilState1.DATA

@@ -10,5 +10,9 @@ DYV
 DZV 
 10*10 /
 
+ACTNUM
+   1 998*2 3 /
+
+
 DEPTHZ
 121*2000 /

tests/test_parser.cpp

@@ -0,0 +1,54 @@
+/*
+  Copyright 2012 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 <http://www.gnu.org/licenses/>.
+*/
+
+#include <config.h>
+
+#if HAVE_DYNAMIC_BOOST_TEST
+#define BOOST_TEST_DYN_LINK
+#endif
+
+#define NVERBOSE  // Suppress own messages when throw()ing
+
+#define BOOST_TEST_MODULE OPM-ParserTest
+#include <boost/test/unit_test.hpp>
+
+#include <opm/parser/eclipse/Parser/Parser.hpp>
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+#include <opm/parser/eclipse/Deck/DeckDoubleItem.hpp>
+
+#include <string>
+#include <iostream>
+#include <vector>
+#include <memory>
+
+BOOST_AUTO_TEST_CASE(CreateParser)
+{
+    const std::string filename1 = "testBlackoilState1.DATA";
+    Opm::ParserPtr parser(new Opm::Parser() );
+    Opm::DeckConstPtr deck = parser->parseFile( filename1 );
+
+    BOOST_CHECK_EQUAL( 6U , deck->size() );
+    Opm::DeckItemConstPtr actnum = deck->getKeyword("ACTNUM")->getRecord(0)->getItem(0);
+    const std::vector<int>& actnum_data = actnum->getIntData();
+
+    BOOST_CHECK_EQUAL( 1000U , actnum->size() );
+    BOOST_CHECK_EQUAL( 1, actnum_data[0] );
+    BOOST_CHECK_EQUAL( 2, actnum_data[400] );
+    BOOST_CHECK_EQUAL( 3, actnum_data[999] );
+}

tests/test_wellcollection.cpp

@@ -0,0 +1,91 @@
+/*
+  Copyright 2012 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 <http://www.gnu.org/licenses/>.
+*/
+
+#include <config.h>
+
+#if HAVE_DYNAMIC_BOOST_TEST
+#define BOOST_TEST_DYN_LINK
+#endif
+
+#define NVERBOSE  // Suppress own messages when throw()ing
+
+#define BOOST_TEST_MODULE WellCollectionTest
+#include <boost/test/unit_test.hpp>
+#include <opm/core/wells/WellCollection.hpp>
+#include <opm/parser/eclipse/Parser/Parser.hpp>
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Group.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/GroupTreeNode.hpp>
+
+using namespace Opm;
+
+BOOST_AUTO_TEST_CASE(AddWellsAndGroupToCollection) {
+    ParserPtr parser(new Parser());
+    boost::filesystem::path scheduleFile("wells_group.data");
+    DeckConstPtr deck =  parser->parseFile(scheduleFile.string());
+    EclipseStateConstPtr eclipseState(new EclipseState(deck));
+    PhaseUsage pu = phaseUsageFromDeck(eclipseState);
+
+    GroupTreeNodePtr field=eclipseState->getSchedule()->getGroupTree(2)->getNode("FIELD");
+    GroupTreeNodePtr g1=eclipseState->getSchedule()->getGroupTree(2)->getNode("G1");
+    GroupTreeNodePtr g2=eclipseState->getSchedule()->getGroupTree(2)->getNode("G2");
+
+    WellCollection collection;
+
+    // Add groups to WellCollection
+    GroupConstPtr fieldGroup =  eclipseState->getSchedule()->getGroup(field->name());
+    collection.addField(fieldGroup, 2, pu);
+
+    for (auto iter = field->begin(); iter != field->end(); ++iter) {
+        GroupConstPtr childGroupNode = eclipseState->getSchedule()->getGroup((*iter).second->name());
+        collection.addGroup(childGroupNode, fieldGroup->name(), 2, pu);
+    }
+
+    GroupConstPtr g1Group =  eclipseState->getSchedule()->getGroup(g1->name());
+    for (auto iter = g1->begin(); iter != g1->end(); ++iter) {
+        GroupConstPtr childGroupNode = eclipseState->getSchedule()->getGroup((*iter).second->name());
+        collection.addGroup(childGroupNode, g1Group->name(), 2, pu);
+    }
+
+
+    GroupConstPtr g2Group =  eclipseState->getSchedule()->getGroup(g2->name());
+    for (auto iter = g2->begin(); iter != g2->end(); ++iter) {
+        GroupConstPtr childGroupNode = eclipseState->getSchedule()->getGroup((*iter).second->name());
+        collection.addGroup(childGroupNode, g2Group->name(), 2, pu);
+    }
+
+    BOOST_CHECK_EQUAL("FIELD", collection.findNode("FIELD")->name());
+    BOOST_CHECK_EQUAL("FIELD", collection.findNode("G1")->getParent()->name());
+    BOOST_CHECK_EQUAL("FIELD", collection.findNode("G2")->getParent()->name());
+
+    // Add wells to WellCollection
+    WellCollection wellCollection;
+    std::vector<WellConstPtr> wells = eclipseState->getSchedule()->getWells();
+    for (size_t i=0; i<wells.size(); i++) {
+        collection.addWell(wells[i], 2, pu);
+    }
+
+    BOOST_CHECK_EQUAL("G1", collection.findNode("INJ1")->getParent()->name());
+    BOOST_CHECK_EQUAL("G1", collection.findNode("INJ2")->getParent()->name());
+    BOOST_CHECK_EQUAL("G2", collection.findNode("PROD1")->getParent()->name());
+    BOOST_CHECK_EQUAL("G2", collection.findNode("PROD2")->getParent()->name());
+}
+

tests/test_wells.cpp

@@ -204,7 +204,7 @@ BOOST_AUTO_TEST_CASE(Copy)
                         BOOST_CHECK_EQUAL( dist1[p] , dist2[p]);
                 }
             }
-            BOOST_CHECK( well_controls_equal( c1 , c2 ));
+            BOOST_CHECK( well_controls_equal( c1 , c2 , false) );
         }
     }
 }
@@ -219,7 +219,7 @@ BOOST_AUTO_TEST_CASE(Equals_WellsEqual_ReturnsTrue) {
     std::shared_ptr<Wells> W2(create_wells(nphases, nwells, nperfs),
                                 destroy_wells);
 
-    BOOST_CHECK(wells_equal(W1.get(), W2.get()));
+    BOOST_CHECK(wells_equal(W1.get(), W2.get() , false));
 }
 
 
@@ -232,5 +232,5 @@ BOOST_AUTO_TEST_CASE(Equals_WellsDiffer_ReturnsFalse) {
     std::shared_ptr<Wells> W2(create_wells(nphases, 3, nperfs),
                                 destroy_wells);
 
-    BOOST_CHECK(!wells_equal(W1.get(), W2.get()));
+    BOOST_CHECK(!wells_equal(W1.get(), W2.get() , false ));
 }

tests/test_wellsgroup.cpp

@@ -0,0 +1,108 @@
+/*
+  Copyright 2014 Statoil.
+
+  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 <http://www.gnu.org/licenses/>.
+*/
+
+#include <config.h>
+
+#if HAVE_DYNAMIC_BOOST_TEST
+#define BOOST_TEST_DYN_LINK
+#endif
+
+#define NVERBOSE  // Suppress own messages when throw()ing
+
+#define BOOST_TEST_MODULE WellsGroupTest
+
+#include <memory>
+#include <vector>
+
+#include <boost/test/unit_test.hpp>
+#include <boost/date_time/posix_time/posix_time.hpp>
+
+#include <opm/core/props/phaseUsageFromDeck.hpp>
+#include <opm/core/wells/WellsGroup.hpp>
+
+#include <opm/parser/eclipse/Parser/Parser.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Group.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/GroupTreeNode.hpp>
+
+#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+
+using namespace Opm;
+
+BOOST_AUTO_TEST_CASE(ConstructGroupFromWell) {
+    ParserPtr parser(new Parser());
+    boost::filesystem::path scheduleFile("wells_group.data");
+    DeckConstPtr deck =  parser->parseFile(scheduleFile.string());
+    EclipseStateConstPtr eclipseState(new EclipseState(deck));
+    PhaseUsage pu = phaseUsageFromDeck(eclipseState);
+
+    std::vector<WellConstPtr> wells = eclipseState->getSchedule()->getWells();
+
+    for (size_t i=0; i<wells.size(); i++) {
+        WellConstPtr well = wells[i];
+        std::shared_ptr<WellsGroupInterface> wellsGroup = createWellWellsGroup(well, 2, pu);
+        BOOST_CHECK_EQUAL(well->name(), wellsGroup->name());
+        if (well->isInjector(2)) {
+            BOOST_CHECK_EQUAL(well->getSurfaceInjectionRate(2), wellsGroup->injSpec().surface_flow_max_rate_);
+            BOOST_CHECK_EQUAL(well->getBHPLimit(2), wellsGroup->injSpec().BHP_limit_);
+            BOOST_CHECK_EQUAL(well->getReservoirInjectionRate(2), wellsGroup->injSpec().reservoir_flow_max_rate_);
+            BOOST_CHECK_EQUAL(0.0, wellsGroup->prodSpec().guide_rate_);
+        }
+        if (well->isProducer(2)) {
+            BOOST_CHECK_EQUAL(well->getResVRate(2), wellsGroup->prodSpec().reservoir_flow_max_rate_);
+            BOOST_CHECK_EQUAL(well->getBHPLimit(2), wellsGroup->prodSpec().BHP_limit_);
+            BOOST_CHECK_EQUAL(well->getOilRate(2), wellsGroup->prodSpec().oil_max_rate_);
+            BOOST_CHECK_EQUAL(well->getWaterRate(2), wellsGroup->prodSpec().water_max_rate_);
+            BOOST_CHECK_EQUAL(0.0, wellsGroup->injSpec().guide_rate_);
+        }
+    }
+}
+
+
+BOOST_AUTO_TEST_CASE(ConstructGroupFromGroup) {
+    ParserPtr parser(new Parser());
+    boost::filesystem::path scheduleFile("wells_group.data");
+    DeckConstPtr deck =  parser->parseFile(scheduleFile.string());
+    EclipseStateConstPtr eclipseState(new EclipseState(deck));
+    PhaseUsage pu = phaseUsageFromDeck(eclipseState);
+
+    std::vector<GroupTreeNodeConstPtr> nodes = eclipseState->getSchedule()->getGroupTree(2)->getNodes();
+
+    for (size_t i=0; i<nodes.size(); i++) {
+        GroupConstPtr group = eclipseState->getSchedule()->getGroup(nodes[i]->name());
+        std::shared_ptr<WellsGroupInterface> wellsGroup = createGroupWellsGroup(group, 2, pu);
+        BOOST_CHECK_EQUAL(group->name(), wellsGroup->name());
+        if (group->isInjectionGroup(2)) {
+            BOOST_CHECK_EQUAL(group->getSurfaceMaxRate(2), wellsGroup->injSpec().surface_flow_max_rate_);
+            BOOST_CHECK_EQUAL(group->getReservoirMaxRate(2), wellsGroup->injSpec().reservoir_flow_max_rate_);
+            BOOST_CHECK_EQUAL(group->getTargetReinjectFraction(2), wellsGroup->injSpec().reinjection_fraction_target_);
+            BOOST_CHECK_EQUAL(group->getTargetVoidReplacementFraction(2), wellsGroup->injSpec().voidage_replacment_fraction_);
+        }
+        if (group->isProductionGroup(2)) {
+            BOOST_CHECK_EQUAL(group->getReservoirMaxRate(2), wellsGroup->prodSpec().reservoir_flow_max_rate_);
+            BOOST_CHECK_EQUAL(group->getGasTargetRate(2), wellsGroup->prodSpec().gas_max_rate_);
+            BOOST_CHECK_EQUAL(group->getOilTargetRate(2), wellsGroup->prodSpec().oil_max_rate_);
+            BOOST_CHECK_EQUAL(group->getWaterTargetRate(2), wellsGroup->prodSpec().water_max_rate_);
+            BOOST_CHECK_EQUAL(group->getLiquidTargetRate(2), wellsGroup->prodSpec().liquid_max_rate_);
+        }
+    }
+}
+
+
+

tests/test_wellsmanager.cpp

@@ -27,6 +27,10 @@
 
 #define BOOST_TEST_MODULE WellsManagerTests
 #include <boost/test/unit_test.hpp>
+
+#include <opm/parser/eclipse/Parser/Parser.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
+
 #include <opm/core/wells/WellsManager.hpp>
 #include <opm/core/wells.h>
 #include <opm/core/well_controls.h>
@@ -57,7 +61,7 @@ void wells_static_check(const Wells* wells) {
 }
 
 
-/* 
+/*
    The number of controls is determined by looking at which elements
    have been given explicit - non-default - values in the WCONxxxx
    keyword. Is that at all interesting?
@@ -67,7 +71,7 @@ void wells_static_check(const Wells* wells) {
 void check_controls_epoch0( struct WellControls ** ctrls) {
     // The injector
     {
-        const struct WellControls * ctrls0 = ctrls[0];        
+        const struct WellControls * ctrls0 = ctrls[0];
         BOOST_CHECK_EQUAL( 3 , well_controls_get_num(ctrls0));   // The number of controls for the injector == 3??
 
         BOOST_CHECK_EQUAL( SURFACE_RATE   , well_controls_iget_type(ctrls0 , 0) );
@@ -83,7 +87,7 @@ void check_controls_epoch0( struct WellControls ** ctrls) {
         BOOST_CHECK_EQUAL( 0 , well_controls_get_current(ctrls0) );
         
         // The phase distribution in the active target
-        { 
+        {
              const double * distr = well_controls_iget_distr( ctrls0 , 0 );
              BOOST_CHECK_EQUAL( 0 , distr[0] );  // Water
              BOOST_CHECK_EQUAL( 0 , distr[1] );  // Oil
@@ -106,7 +110,7 @@ void check_controls_epoch0( struct WellControls ** ctrls) {
         BOOST_CHECK_EQUAL( 0 , well_controls_get_current(ctrls1));
 
         // The phase distribution in the active target
-       { 
+       {
             const double * distr = well_controls_iget_distr( ctrls1 , 0 );
             BOOST_CHECK_EQUAL( 0 , distr[0] );  // Water
             BOOST_CHECK_EQUAL( 1 , distr[1] );  // Oil
@@ -121,7 +125,7 @@ void check_controls_epoch0( struct WellControls ** ctrls) {
 void check_controls_epoch1( struct WellControls ** ctrls) {
     // The injector
     {
-        const struct WellControls * ctrls0 = ctrls[0];        
+        const struct WellControls * ctrls0 = ctrls[0];
         BOOST_CHECK_EQUAL( 3 , well_controls_get_num(ctrls0));   // The number of controls for the injector == 3??
 
         BOOST_CHECK_EQUAL( SURFACE_RATE   , well_controls_iget_type(ctrls0 , 0 ));
@@ -136,7 +140,7 @@ void check_controls_epoch1( struct WellControls ** ctrls) {
         // Which control is active
         BOOST_CHECK_EQUAL( 1 , well_controls_get_current(ctrls0));
 
-        { 
+        {
             const double * distr = well_controls_iget_distr( ctrls0 , 1 );
             BOOST_CHECK_EQUAL( 1 , distr[0] );  // Water
             BOOST_CHECK_EQUAL( 0 , distr[1] );  // Oil
@@ -160,7 +164,7 @@ void check_controls_epoch1( struct WellControls ** ctrls) {
         // Which control is active
         BOOST_CHECK_EQUAL( 1 , well_controls_get_current(ctrls1) );
 
-        { 
+        {
             const double * distr = well_controls_iget_distr( ctrls1 , 1 );
             BOOST_CHECK_EQUAL( 1 , distr[0] );  // Water
             BOOST_CHECK_EQUAL( 1 , distr[1] );  // Oil
@@ -170,8 +174,6 @@ void check_controls_epoch1( struct WellControls ** ctrls) {
 }
 
 
-
-
 BOOST_AUTO_TEST_CASE(Constructor_Works) {
     Opm::EclipseGridParser Deck("wells_manager_data.data");
     Opm::GridManager gridManager(Deck);
@@ -190,11 +192,86 @@ BOOST_AUTO_TEST_CASE(Constructor_Works) {
         Opm::WellsManager wellsManager(Deck, *gridManager.c_grid(), NULL);
         const Wells* wells = wellsManager.c_wells();
         
-        wells_static_check( wells );    
+        wells_static_check( wells );
         check_controls_epoch1( wells->ctrls );
     }
 }
 
+BOOST_AUTO_TEST_CASE(New_Constructor_Works) {
+
+    Opm::ParserPtr parser(new Opm::Parser());
+    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(parser->parseFile("wells_manager_data.data")));
+
+    Opm::EclipseGridParser Deck("wells_manager_data.data");
+    Opm::GridManager gridManager(Deck);
+
+    Deck.setCurrentEpoch(0);
+    {
+        Opm::WellsManager wellsManager(eclipseState, 0, *gridManager.c_grid(), NULL);
+        Opm::WellsManager oldWellsManager(Deck, *gridManager.c_grid(), NULL);
+
+        std::cout << "Checking new well structure, epoch 0" << std::endl;
+        wells_static_check( wellsManager.c_wells() );
+
+        std::cout << "Checking old well structure, epoch 0" << std::endl;
+        wells_static_check( oldWellsManager.c_wells() );
+
+        check_controls_epoch0( wellsManager.c_wells()->ctrls );
+
+        BOOST_CHECK(wells_equal(wellsManager.c_wells(), oldWellsManager.c_wells() , false));
+    }
+
+    Deck.setCurrentEpoch(1);
+    {
+        Opm::WellsManager wellsManager(eclipseState, 1, *gridManager.c_grid(), NULL);
+        Opm::WellsManager oldWellsManager(Deck, *gridManager.c_grid(), NULL);
+
+        std::cout << "Checking new well structure, epoch 1" << std::endl;
+        wells_static_check( wellsManager.c_wells() );
+
+        std::cout << "Checking old well structure, epoch 1" << std::endl;
+        wells_static_check( oldWellsManager.c_wells() );
+
+        check_controls_epoch1( wellsManager.c_wells()->ctrls );
+
+        BOOST_CHECK(wells_equal( wellsManager.c_wells(), oldWellsManager.c_wells(),false));
+    }
+}
+
+
+BOOST_AUTO_TEST_CASE(New_Constructor_Works_ExpandedData) {
+
+    Opm::ParserPtr parser(new Opm::Parser());
+    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(parser->parseFile("wells_manager_data_expanded.data")));
+
+    Opm::EclipseGridParser Deck("wells_manager_data_expanded.data");
+    Opm::GridManager gridManager(Deck);
+
+    Deck.setCurrentEpoch(0);
+    {
+        Opm::WellsManager wellsManager(eclipseState, 0, *gridManager.c_grid(), NULL);
+        Opm::WellsManager oldWellsManager(Deck, *gridManager.c_grid(), NULL);
+
+        BOOST_CHECK(wells_equal(wellsManager.c_wells(), oldWellsManager.c_wells(),false));
+    }
+
+    Deck.setCurrentEpoch(1);
+    {
+        Opm::WellsManager wellsManager(eclipseState, 1, *gridManager.c_grid(), NULL);
+        Opm::WellsManager oldWellsManager(Deck, *gridManager.c_grid(), NULL);
+
+        BOOST_CHECK(wells_equal( wellsManager.c_wells(), oldWellsManager.c_wells(), true));
+    }
+
+    Deck.setCurrentEpoch(2);
+    {
+        Opm::WellsManager wellsManager(eclipseState, 2, *gridManager.c_grid(), NULL);
+        Opm::WellsManager oldWellsManager(Deck, *gridManager.c_grid(), NULL);
+
+        BOOST_CHECK(wells_equal( wellsManager.c_wells(), oldWellsManager.c_wells(), true));
+    }
+}
+
 
 
 BOOST_AUTO_TEST_CASE(WellsEqual) {
@@ -207,8 +284,8 @@ BOOST_AUTO_TEST_CASE(WellsEqual) {
     Deck.setCurrentEpoch(1);
     Opm::WellsManager wellsManager1(Deck, *gridManager.c_grid(), NULL);
 
-    BOOST_CHECK(  wells_equal( wellsManager0.c_wells() , wellsManager0.c_wells()) ); 
+    BOOST_CHECK(  wells_equal( wellsManager0.c_wells() , wellsManager0.c_wells(),false) );
-    BOOST_CHECK( !wells_equal( wellsManager0.c_wells() , wellsManager1.c_wells()) ); 
+    BOOST_CHECK( !wells_equal( wellsManager0.c_wells() , wellsManager1.c_wells(),false) );
 }
 
 
@@ -222,15 +299,30 @@ BOOST_AUTO_TEST_CASE(ControlsEqual) {
     Deck.setCurrentEpoch(1);
     Opm::WellsManager wellsManager1(Deck, *gridManager.c_grid(), NULL);
 
-    BOOST_CHECK(  well_controls_equal( wellsManager0.c_wells()->ctrls[0] , wellsManager0.c_wells()->ctrls[0]));
+    BOOST_CHECK(  well_controls_equal( wellsManager0.c_wells()->ctrls[0] , wellsManager0.c_wells()->ctrls[0] , false));
-    BOOST_CHECK(  well_controls_equal( wellsManager0.c_wells()->ctrls[1] , wellsManager0.c_wells()->ctrls[1]));
+    BOOST_CHECK(  well_controls_equal( wellsManager0.c_wells()->ctrls[1] , wellsManager0.c_wells()->ctrls[1] , false));
-    BOOST_CHECK(  well_controls_equal( wellsManager1.c_wells()->ctrls[0] , wellsManager1.c_wells()->ctrls[0]));
+    BOOST_CHECK(  well_controls_equal( wellsManager1.c_wells()->ctrls[0] , wellsManager1.c_wells()->ctrls[0] , false));
-    BOOST_CHECK(  well_controls_equal( wellsManager1.c_wells()->ctrls[1] , wellsManager1.c_wells()->ctrls[1]));
+    BOOST_CHECK(  well_controls_equal( wellsManager1.c_wells()->ctrls[1] , wellsManager1.c_wells()->ctrls[1] , false));
 
-    BOOST_CHECK(  !well_controls_equal( wellsManager0.c_wells()->ctrls[0] , wellsManager0.c_wells()->ctrls[1]));
+    BOOST_CHECK(  !well_controls_equal( wellsManager0.c_wells()->ctrls[0] , wellsManager0.c_wells()->ctrls[1] , false));
-    BOOST_CHECK(  !well_controls_equal( wellsManager0.c_wells()->ctrls[1] , wellsManager0.c_wells()->ctrls[0]));
+    BOOST_CHECK(  !well_controls_equal( wellsManager0.c_wells()->ctrls[1] , wellsManager0.c_wells()->ctrls[0] , false));
-    BOOST_CHECK(  !well_controls_equal( wellsManager1.c_wells()->ctrls[0] , wellsManager0.c_wells()->ctrls[0]));
+    BOOST_CHECK(  !well_controls_equal( wellsManager1.c_wells()->ctrls[0] , wellsManager0.c_wells()->ctrls[0] , false));
-    BOOST_CHECK(  !well_controls_equal( wellsManager1.c_wells()->ctrls[1] , wellsManager0.c_wells()->ctrls[1]));
+    BOOST_CHECK(  !well_controls_equal( wellsManager1.c_wells()->ctrls[1] , wellsManager0.c_wells()->ctrls[1] , false));
 }
 
 
+
+BOOST_AUTO_TEST_CASE(WellHasSTOP_ExceptionIsThrown) {
+    Opm::EclipseGridParser Deck("wells_manager_data_wellSTOP.data");
+    Opm::GridManager gridManager(Deck);
+
+    Opm::ParserPtr parser(new Opm::Parser());
+    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(parser->parseFile("wells_manager_data_wellSTOP.data")));
+
+    Deck.setCurrentEpoch(0);
+
+    BOOST_CHECK_THROW( new Opm::WellsManager(eclipseState, 0, *gridManager.c_grid(), NULL), std::runtime_error );
+}
+
+
+

tests/wells_group.data

@@ -0,0 +1,68 @@
+OIL
+GAS
+WATER
+
+DIMENS
+   10 10  5  /
+
+GRID
+
+DXV
+10*1000.0 /
+
+DYV
+10*1000.0 /
+
+DZV
+10.0 20.0 30.0 10.0 5.0 /
+
+DEPTHZ
+121*2000
+/
+
+SCHEDULE
+
+GRUPTREE
+ 'G1' 'FIELD' /
+ 'G2' 'FIELD' /
+/
+
+
+WELSPECS
+    'INJ1' 'G1'    1  1    8335 'GAS'  /
+    'PROD1' 'G2'   10 10    8400 'OIL'  /
+/
+
+TSTEP
+  14.0 /
+/
+
+WELSPECS
+    'INJ2' 'G1'    1  1    8335 'GAS'  /
+    'PROD2' 'G2'   10 10    8400 'OIL'  /
+/
+
+GCONINJE
+'G1' GAS RATE 30000 /
+/
+
+GCONPROD
+'G2' ORAT 10000 /
+/
+
+WCONINJE
+   'INJ1' 'WATER' 'OPEN' 'RESV' 10 20 40 /
+   'INJ2' 'WATER' 'OPEN' 'RESV' 10 20 40 /
+/
+
+WCONPROD
+   'PROD1' 'OPEN' 'RESV' 999  3* 123 100 /
+   'PROD2' 'OPEN' 'RESV' 999  3* 123 100 /
+/
+
+TSTEP
+  3 /
+/
+
+
+END

tests/wells_manager_data.data

@@ -17,8 +17,8 @@ DZV
 10.0 20.0 30.0 10.0 5.0 /
 
 DEPTHZ
-121*2000
+121*2000 /
-/
+
 
 SCHEDULE
 
@@ -34,11 +34,11 @@ COMPDAT
 
 WCONPROD
      'PROD1' 'OPEN' 'ORAT' 20000  4* 1000 /
-  /
+/
 
 WCONINJE
      'INJ1' 'GAS' 'OPEN' 'RATE' 100 200 400 /
-  /
+/
 
 
 DATES
@@ -53,10 +53,26 @@ WCONINJE
    'INJ1' 'WATER' 'OPEN' 'RESV' 10 20 40 /
 /
 
-
+END
 
 TSTEP
-14.0 
+  14.0 /
- /
+/
+
+
+WELSPECS
+    'TEST1' 'G1'    1  1    8335 'GAS'  /
+    'TEST2' 'G2'   10 10    8400 'OIL'  /
+/
+
+GRUPTREE
+ 'G1' 'SHIP' /
+ 'G2' 'SHIP' /
+/
+
+TSTEP
+  3 /
+/
+
 
 END

tests/wells_manager_data_expanded.data

@@ -0,0 +1,86 @@
+OIL
+GAS
+WATER
+
+DIMENS
+   10 10  5  /
+
+GRID
+
+DXV
+10*1000.0 /
+
+DYV
+10*1000.0 /
+
+DZV
+10.0 20.0 30.0 10.0 5.0 /
+
+-- The DEPTHZ keyword is only here to satisfy the old parser; content might
+-- completely bogus.
+DEPTHZ
+121*2000 /
+
+
+SCHEDULE
+
+WELSPECS
+    'INJ1' 'G'    1  1    8335 'GAS'  /
+    'PROD1' 'G'   10 10    8400 'OIL'  /
+/
+
+COMPDAT
+    'INJ1'   1  1 1  2 'OPEN' 1   10.6092   0.5  /
+    'INJ1'   1  1 3  5 'OPEN' 1   12.6092   0.5  /
+    'INJ1'   2  2 1  1 'OPEN' 1   14.6092   0.5  /
+    'PROD1'  10 3 3  3 'OPEN' 0   10.6092   0.5  /
+/
+
+WCONPROD
+     'PROD1' 'OPEN' 'ORAT' 20000  4* 1000 /
+/
+
+WCONINJE
+     'INJ1' 'GAS' 'OPEN' 'RATE' 100 200 400 /
+/
+
+
+DATES
+ 1 'FEB' 2000 /
+/
+
+WELSPECS
+    'INJ1' 'G3'    1  1    8335 'GAS'  /
+    'QNJ2' 'G3'    1  1    8335 'GAS'  /
+/
+
+
+COMPDAT
+    'QNJ2'   3  4 1  2 'OPEN' 1   10.6092   0.5  /
+    'QNJ2'   4  4 3  5 'OPEN' 1   12.6092   0.5  /
+/
+
+WCONPROD
+   'PROD1' 'OPEN' 'RESV' 999  3* 123 100 /
+/
+
+WCONINJE
+   'INJ1' 'WATER' 'OPEN' 'RESV' 10 20 40 /
+   'QNJ2' 'WATER' 'OPEN' 'RESV' 7 33 39 /
+/
+
+
+
+TSTEP
+14.0 /
+/
+
+WELOPEN
+   'INJ1' 'SHUT' 5* /
+/
+
+TSTEP
+14.0 /
+/
+
+END

tests/wells_manager_data_wellSTOP.data

@@ -0,0 +1,42 @@
+OIL
+GAS
+WATER
+
+DIMENS
+   10 10  5  /
+
+GRID
+
+DXV
+10*1000.0 /
+
+DYV
+10*1000.0 /
+
+DZV
+10.0 20.0 30.0 10.0 5.0 /
+
+DEPTHZ
+121*2000
+/
+
+SCHEDULE
+
+WELSPECS
+    'INJ1' 'G'    1  1    8335 'GAS'  /
+    'PROD1' 'G'   10 10    8400 'OIL'  /
+/
+
+COMPDAT
+    'INJ1'   1  1 1  1 'OPEN' 1   10.6092   0.5  /
+    'PROD1'  10 3 3  3 'OPEN' 0   10.6092   0.5  /
+/
+
+WELOPEN
+   'INJ1' 'STOP' 5* /
+/
+
+TSTEP
+  10 /
+
+END