Merge pull request #454 from jalvestad/restart-group+mswell

Restart group+mswell

CMakeLists_files.cmake

@@ -135,6 +135,9 @@ if(ENABLE_ECL_INPUT)
 endif()
 if(ENABLE_ECL_OUTPUT)
   list( APPEND MAIN_SOURCE_FILES
+          src/opm/output/eclipse/AggregateConnectionData.cpp
+          src/opm/output/eclipse/AggregateGroupData.cpp
+          src/opm/output/eclipse/AggregateMSWData.cpp
           src/opm/output/eclipse/AggregateWellData.cpp
           src/opm/output/eclipse/CreateDoubHead.cpp
           src/opm/output/eclipse/CreateInteHead.cpp
@@ -144,7 +147,9 @@ if(ENABLE_ECL_OUTPUT)
           src/opm/output/eclipse/EclipseGridInspector.cpp
           src/opm/output/eclipse/EclipseIO.cpp
           src/opm/output/eclipse/InteHEAD.cpp
+          src/opm/output/eclipse/libECLRestart.cpp
           src/opm/output/eclipse/LinearisedOutputTable.cpp
+          src/opm/output/eclipse/LoadRestart.cpp
           src/opm/output/eclipse/LogiHEAD.cpp
           src/opm/output/eclipse/RestartIO.cpp
           src/opm/output/eclipse/Summary.cpp
@@ -481,8 +486,12 @@ if(ENABLE_ECL_OUTPUT)
         opm/output/data/Cells.hpp
         opm/output/data/Solution.hpp
         opm/output/data/Wells.hpp
+        opm/output/eclipse/VectorItems/connection.hpp
         opm/output/eclipse/VectorItems/intehead.hpp
         opm/output/eclipse/VectorItems/well.hpp
+        opm/output/eclipse/AggregateGroupData.hpp
+        opm/output/eclipse/AggregateConnectionData.hpp
+        opm/output/eclipse/AggregateMSWData.hpp
         opm/output/eclipse/AggregateWellData.hpp
         opm/output/eclipse/CharArrayNullTerm.hpp
         opm/output/eclipse/DoubHEAD.hpp
@@ -490,8 +499,9 @@ if(ENABLE_ECL_OUTPUT)
         opm/output/eclipse/EclipseIO.hpp
         opm/output/eclipse/EclipseIOUtil.hpp
         opm/output/eclipse/InteHEAD.hpp
-        opm/output/eclipse/LogiHEAD.hpp
+        opm/output/eclipse/libECLRestart.hpp
         opm/output/eclipse/LinearisedOutputTable.hpp
+        opm/output/eclipse/LogiHEAD.hpp
         opm/output/eclipse/RegionCache.hpp
         opm/output/eclipse/RestartIO.hpp
         opm/output/eclipse/RestartValue.hpp

opm/output/eclipse/AggregateConnectionData.hpp

@@ -0,0 +1,75 @@
+/*
+  Copyright (c) 2018 Equinor ASA
+
+  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/>.
+*/
+
+#ifndef OPM_AGGREGATE_CONNECTION_DATA_HPP
+#define OPM_AGGREGATE_CONNECTION_DATA_HPP
+
+#include <opm/output/eclipse/WindowedArray.hpp>
+
+#include <cstddef>
+#include <string>
+#include <vector>
+
+namespace Opm {
+    class EclipseGrid;
+    class Schedule;
+    class UnitSystem;
+} // Opm
+
+namespace Opm { namespace data {
+    class WellRates;
+}}
+
+namespace Opm { namespace RestartIO { namespace Helpers {
+
+    class AggregateConnectionData
+    {
+    public:
+        explicit AggregateConnectionData(const std::vector<int>& inteHead);
+
+        void captureDeclaredConnData(const Opm::Schedule&        sched,
+                                     const Opm::EclipseGrid&     grid,
+                                     const Opm::UnitSystem&      units,
+                                     const Opm::data::WellRates& xw,
+                                     const std::size_t           sim_step);
+
+        const std::vector<int>& getIConn() const
+        {
+            return this->iConn_.data();
+        }
+
+        const std::vector<float>& getSConn() const
+        {
+            return this->sConn_.data();
+        }
+
+        const std::vector<double>& getXConn() const
+        {
+            return this->xConn_.data();
+        }
+
+    private:
+        WindowedMatrix<int> iConn_;
+        WindowedMatrix<float> sConn_;
+        WindowedMatrix<double> xConn_;
+    };
+
+}}} // Opm::RestartIO::Helpers
+
+#endif // OPM_AGGREGATE_CONNECTION_DATA_HPP

opm/output/eclipse/AggregateGroupData.hpp

@@ -0,0 +1,159 @@
+/*
+  Copyright (c) 2018 Statoil ASA
+
+  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/>.
+*/
+
+#ifndef OPM_AGGREGATE_GROUP_DATA_HPP
+#define OPM_AGGREGATE_GROUP_DATA_HPP
+
+#include <opm/output/eclipse/CharArrayNullTerm.hpp>
+#include <opm/output/eclipse/WindowedArray.hpp>
+
+#include <cstddef>
+#include <string>
+#include <vector>
+#include <map>
+
+namespace Opm {
+    class Schedule;
+    class SummaryState;
+    class Group;
+} // Opm
+
+namespace Opm { namespace RestartIO { namespace Helpers {
+
+    class groupMaps {
+    public:
+	const std::map <size_t, const Opm::Group*>& indexGroupMap() const;
+	const std::map <const std::string, size_t>& groupNameIndexMap() const;
+
+	void currentGrpTreeNameSeqIndMap(const Opm::Schedule&                        sched,
+	                                 const size_t                                simStep,
+	                                 const std::map<const std::string , size_t>& GnIMap,
+	                                 const std::map<size_t, const Opm::Group*>&  IGMap);
+
+    private:
+	std::map <size_t, const Opm::Group*> m_indexGroupMap;
+	std::map <const std::string, size_t> m_groupNameIndexMap;
+    };
+
+    class AggregateGroupData
+    {
+    public:
+	explicit AggregateGroupData(const std::vector<int>& inteHead);
+
+	void captureDeclaredGroupData(const Opm::Schedule&                 sched,
+				      const std::vector<std::string>&      restart_group_keys,
+				      const std::vector<std::string>&      restart_field_keys,
+				      const std::map<std::string, size_t>& groupKeyToIndex,
+				      const std::map<std::string, size_t>& fieldKeyToIndex,
+				      const bool                           ecl_compatible_rst,
+				      const std::size_t                    simStep,
+				      const Opm::SummaryState&             sumState,
+				      const std::vector<int>&              inteHead);
+
+        const std::vector<int>& getIGroup() const
+        {
+            return this->iGroup_.data();
+        }
+
+        const std::vector<float>& getSGroup() const
+        {
+            return this->sGroup_.data();
+        }
+
+        const std::vector<double>& getXGroup() const
+        {
+            return this->xGroup_.data();
+        }
+
+        const std::vector<CharArrayNullTerm<8>>& getZGroup() const
+        {
+            return this->zGroup_.data();
+        }
+
+	const std::vector<std::string> restart_group_keys = {"GOPP", "GWPP", "GOPR", "GWPR", "GGPR",
+							     "GVPR", "GWIR", "GGIR", "GWCT", "GGOR",
+							     "GOPT", "GWPT", "GGPT", "GVPT", "GWIT",
+							     "GGIT"};
+
+	const std::vector<std::string> restart_field_keys = {"FOPP", "FWPP", "FOPR", "FWPR", "FGPR",
+							     "FVPR", "FWIR", "FGIR", "FWCT", "FGOR",
+							     "FOPT", "FWPT", "FGPT", "FVPT", "FWIT",
+							     "FGIT"};
+
+	const std::map<std::string, size_t> groupKeyToIndex = {
+	  {"GOPR",  0},
+	  {"GWPR",  1},
+	  {"GGPR",  2},
+	  {"GVPR",  3},
+	  {"GWIR",  5},
+	  {"GGIR",  6},
+	  {"GWCT",  8},
+	  {"GGOR",  9},
+	  {"GOPT", 10},
+	  {"GWPT", 11},
+	  {"GGPT", 12},
+	  {"GVPT", 13},
+	  {"GWIT", 15},
+	  {"GGIT", 16},
+	  {"GOPP", 22},
+	  {"GWPP", 23},
+	};
+
+	const std::map<std::string, size_t> fieldKeyToIndex = {
+	  {"FOPR",  0},
+	  {"FWPR",  1},
+	  {"FGPR",  2},
+	  {"FVPR",  3},
+	  {"FWIR",  5},
+	  {"FGIR",  6},
+	  {"FWCT",  8},
+	  {"FGOR",  9},
+	  {"FOPT", 10},
+	  {"FWPT", 11},
+	  {"FGPT", 12},
+	  {"FVPT", 13},
+	  {"FWIT", 15},
+	  {"FGIT", 16},
+	  {"FOPP", 22},
+	  {"FWPP", 23},
+	};
+
+    private:
+        /// Aggregate 'IWEL' array (Integer) for all wells.
+        WindowedArray<int> iGroup_;
+
+        /// Aggregate 'SWEL' array (Real) for all wells.
+        WindowedArray<float> sGroup_;
+
+        /// Aggregate 'XWEL' array (Double Precision) for all wells.
+        WindowedArray<double> xGroup_;
+
+        /// Aggregate 'ZWEL' array (Character) for all wells.
+        WindowedArray<CharArrayNullTerm<8>> zGroup_;
+
+        /// Maximum number of wells in a group.
+        int nWGMax_;
+
+        /// Maximum number of groups
+        int nGMaxz_;
+    };
+
+}}} // Opm::RestartIO::Helpers
+
+#endif // OPM_AGGREGATE_WELL_DATA_HPP

opm/output/eclipse/AggregateMSWData.hpp

@@ -0,0 +1,100 @@
+/*
+  Copyright (c) 2018 Statoil ASA
+
+  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/>.
+*/
+
+#ifndef OPM_AGGREGATE_MSW_DATA_HPP
+#define OPM_AGGREGATE_MSW_DATA_HPP
+
+#include <opm/output/data/Wells.hpp>
+#include <opm/output/eclipse/CharArrayNullTerm.hpp>
+#include <opm/output/eclipse/WindowedArray.hpp>
+
+#include <string>
+#include <vector>
+
+namespace Opm {
+    class Phases;
+    class Schedule;
+    class EclipseGrid;
+    class UnitSystem;
+} // Opm
+
+namespace Opm { namespace RestartIO { namespace Helpers {
+
+    struct BranchSegmentPar {
+      int outletS;
+      int noSegInBranch;
+      int firstSeg;
+      int lastSeg;
+      int branch;
+    };
+
+    class AggregateMSWData
+    {
+    public:
+        explicit AggregateMSWData(const std::vector<int>& inteHead);
+
+        void captureDeclaredMSWData(const Opm::Schedule& sched,
+                                     const std::size_t    rptStep,
+				     const Opm::UnitSystem& units,
+				     const std::vector<int>& inteHead,
+				     const Opm::EclipseGrid&  grid);
+
+        /// Retrieve Integer Multisegment well data Array.
+        const std::vector<int>& getISeg() const
+        {
+            return this->iSeg_.data();
+        }
+
+        /// Retrieve Double precision segment data Array.
+        const std::vector<double>& getRSeg() const
+        {
+            return this->rSeg_.data();
+        }
+
+        /// Retrieve Integer multisegment well data Array for lateral branches  (ILBS)
+        const std::vector<int>& getILBs() const
+        {
+            return this->iLBS_.data();
+        }
+
+        /// Retrieve Integer multisegment well data Array for lateral branches (ILBR)
+        const std::vector<int>& getILBr() const
+        {
+            return this->iLBR_.data();
+        }
+
+
+    private:
+        /// Aggregate 'ISEG' array (Integer) for all multisegment wells
+        WindowedArray<int> iSeg_;
+
+        /// Aggregate 'RSEG' array (Double) for all multisegment wells
+        WindowedArray<double> rSeg_;
+
+        /// Aggregate 'ILBS' array (Integer) for all multisegment wells
+        WindowedArray<int> iLBS_;
+
+        /// Aggregate 'ILBR' array (Integer) for all multisegment wells
+        WindowedArray<int> iLBR_;
+
+    };
+
+}}} // Opm::RestartIO::Helpers
+
+#endif // OPM_AGGREGATE_WELL_DATA_HPP

opm/output/eclipse/AggregateWellData.hpp

@@ -44,12 +44,15 @@ namespace Opm { namespace RestartIO { namespace Helpers {
     public:
         explicit AggregateWellData(const std::vector<int>& inteHead);
 
-        void captureDeclaredWellData(const Opm::Schedule&   sched,
+	void captureDeclaredWellData(const Schedule&   	sched,
-                                     const Opm::UnitSystem& units,
+                        const UnitSystem& 		units,
-                                     const std::size_t      sim_step);
+                        const std::size_t 		sim_step,
+			const ::Opm::SummaryState&  	smry,
+			const std::vector<int>& 	inteHead);
 
         void captureDynamicWellData(const Opm::Schedule&        sched,
                                     const std::size_t           sim_step,
+				    const bool ecl_compatible_rst,
                                     const Opm::data::WellRates& xw,
                                     const Opm::SummaryState&    smry);
 
@@ -77,6 +80,8 @@ namespace Opm { namespace RestartIO { namespace Helpers {
             return this->zWell_.data();
         }
 
+        
+        
     private:
         /// Aggregate 'IWEL' array (Integer) for all wells.
         WindowedArray<int> iWell_;

opm/output/eclipse/EclipseIO.hpp

@@ -35,6 +35,7 @@
 #include <opm/output/data/Solution.hpp>
 #include <opm/output/data/Wells.hpp>
 #include <opm/output/eclipse/RestartValue.hpp>
+#include <opm/output/eclipse/SummaryState.hpp>
 
 namespace Opm {
 
@@ -177,7 +178,7 @@ public:
                         const std::map<std::string, double>& single_summary_values,
                         const std::map<std::string, std::vector<double>>& region_summary_values,
                         const std::map<std::pair<std::string, int>, double>& block_summary_values,
-                        bool write_double = false);
+                        const bool write_double = false);
 
 
     /*

opm/output/eclipse/InteHEAD.hpp

@@ -85,6 +85,10 @@ namespace Opm { namespace RestartIO {
             int mxwsit;
             int mxwpit;
         };
+	
+	struct Group {
+	  int ngroups;
+	};
 
         InteHEAD();
         ~InteHEAD() = default;
@@ -112,6 +116,7 @@ namespace Opm { namespace RestartIO {
         InteHEAD& variousParam(const int version, const int iprog);
         InteHEAD& wellSegDimensions(const WellSegDims& wsdim);
         InteHEAD& regionDimensions(const RegDims& rdim);
+	InteHEAD& ngroups(const Group& gr);
 
         const std::vector<int>& data() const
         {

opm/output/eclipse/RestartIO.hpp

@@ -46,6 +46,7 @@ class EclipseGrid;
 class EclipseState;
 class Phases;
 class Schedule;
+class SummaryState;
 
 namespace RestartIO {
 
@@ -72,6 +73,17 @@ namespace RestartIO {
    the report step argument '99'.
 */
 
+/*void save(const std::string& filename,
+          int report_step,
+          double seconds_elapsed,
+          data::Solution cells,
+          data::Wells wells,
+          const EclipseState& es,
+          const EclipseGrid& grid,
+          const Schedule& schedule,
+          std::map<std::string, std::vector<double>> extra_data = {},
+          bool write_double = false);
+*/
 void save(const std::string& filename,
           int report_step,
           double seconds_elapsed,
@@ -79,9 +91,9 @@ void save(const std::string& filename,
           const EclipseState& es,
           const EclipseGrid& grid,
           const Schedule& schedule,
+          const SummaryState& sumState,
           bool write_double = false);
 
-
 RestartValue load( const std::string& filename,
                    int report_step,
                    const std::vector<RestartKey>& solution_keys,

opm/output/eclipse/RestartValue.hpp

@@ -71,6 +71,9 @@ namespace Opm {
         void addExtra(const std::string& key, UnitSystem::measure dimension, std::vector<double> data);
         void addExtra(const std::string& key, std::vector<double> data);
         const std::vector<double>& getExtra(const std::string& key) const;
+
+        void convertFromSI(const UnitSystem& units);
+        void convertToSI(const UnitSystem& units);
     };
 
 }

opm/output/eclipse/VectorItems/connection.hpp

@@ -0,0 +1,74 @@
+/*
+  Copyright (c) 2018 Equinor ASA
+
+  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/>.
+*/
+
+#ifndef OPM_OUTPUT_ECLIPSE_VECTOR_CONNECTION_HPP
+#define OPM_OUTPUT_ECLIPSE_VECTOR_CONNECTION_HPP
+
+#include <vector>
+
+namespace Opm { namespace RestartIO { namespace Helpers { namespace VectorItems {
+    namespace IConn {
+        enum index : std::vector<int>::size_type {
+            SeqIndex   =  0, // Connection sequence index
+            CellI      =  1, // I-location (1-based cell index) of connection
+            CellJ      =  2, // J-location (1-based cell index) of connection
+            CellK      =  3, // K-location (1-based cell index) of connection
+            ConnStat   =  5, // Connection status.
+                             //    > 0 => open, shut otherwise
+
+            Drainage   =  6, // Saturation function (table ID) for drainage
+            Imbibition =  9, // Saturation function (table ID) for imbibition
+
+            ComplNum   = 12, // Completion ID (1-based)
+            ConnDir    = 13, // Penetration direction (1:X, 2:Y, 3:Z)
+            Segment    = 14, // Segment ID of connection
+                             //   0 for regular connections, > 0 for MSW.
+        };
+    } // IConn
+
+    namespace SConn {
+        enum index : std::vector<float>::size_type {
+            ConnTrans    =  0, // Connection transmissibility factor
+            Depth        =  1, // Connection centre depth
+            Diameter     =  2, // Connection diameter
+
+            EffectiveKH  =  3, // Effective Kh product of connection
+
+            item12       = 11, // Unknown
+
+            SegDistEnd   = 20, // Distance to end of connection in segment
+            SegDistStart = 21, // Distance to start of connection in segment
+
+            item30       = 29, // Unknown
+            item31       = 30, // Unknown
+        };
+    } // SConn
+
+    namespace XConn {
+        enum index : std::vector<double>::size_type {
+            OilRate    =  0,  // Surface flow rate (oil)
+            WaterRate  =  1,  // Surface flow rate (water)
+            GasRate    =  2,  // Surface Flow rate (gas)
+
+            ResVRate   = 49,  // Reservoir voidage rate
+        };
+    } // XConn
+}}}} // Opm::RestartIO::Helpers::VectorItems
+
+#endif // OPM_OUTPUT_ECLIPSE_VECTOR_CONNECTION_HPP

opm/output/eclipse/WriteRestartHelpers.hpp

@@ -42,7 +42,7 @@ namespace Opm {
 namespace Opm { namespace RestartIO { namespace Helpers {
 
     const double UNIMPLEMENTED_VALUE = 1e-100; // placeholder for values not yet available 
-            
+
     std::vector<double>
     createDoubHead(const EclipseState& es,
                    const Schedule&     sched,
@@ -57,12 +57,12 @@ namespace Opm { namespace RestartIO { namespace Helpers {
                    const Schedule&     sched,
                    const double        simTime,
                    const int           num_solver_steps,
-                   const int           lookup_step,   // The integer index used to look up dynamic properties, e.g. the number of well.
+                   const int           lookup_step);   // The integer index used to look up dynamic properties, e.g. the number of well.
-                   const int           report_step);  // The integer number this INTEHEAD keyword will be saved to, typically report_step = lookup_step + 1.
 
     std::vector<bool>
     createLogiHead(const EclipseState& es);
 
+
     std::vector<int> serialize_ICON(int lookup_step, // The integer index used to look up dynamic properties, e.g. the number of well.
                                     int ncwmax,      // Max number of completions per well, should be entry 17 from createInteHead.
                                     int niconz,      // Number of elements per completion in ICON, should be entry 32 from createInteHead.

opm/parser/eclipse/EclipseState/Runspec.hpp

@@ -75,7 +75,12 @@ public:
         return this->nGMax;
     }
 
+        int maxWellsInField() const
+    {
+        return this->nWMax;
+    }
 private:
+    int nWMax  { 0 };
     int nCWMax { 0 };
     int nWGMax { 0 };
     int nGMax  { 0 };

opm/parser/eclipse/EclipseState/Schedule/Connection.hpp

@@ -22,15 +22,16 @@
 #define COMPLETION_HPP_
 
 #include <array>
+#include <cstddef>
 #include <map>
 #include <memory>
 #include <string>
 #include <vector>
+
 #include <opm/parser/eclipse/Parser/ParseContext.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/ScheduleEnums.hpp>
 #include <opm/parser/eclipse/EclipseState/Util/Value.hpp>
 
-
 namespace Opm {
 
     class DeckKeyword;
@@ -46,7 +47,12 @@ namespace Opm {
                    double Kh,
                    double rw,
                    const int satTableId,
-                   const WellCompletion::DirectionEnum direction);
+                   const WellCompletion::DirectionEnum direction,
+		   const std::size_t seqIndex,
+		   const double segDistStart,
+		   const double segDistEnd,
+		   const bool defaultSatTabId
+		  );
 
 
         bool attachedToSegment() const;
@@ -68,7 +74,16 @@ namespace Opm {
         void setState(WellCompletion::StateEnum state);
         void setComplnum(int compnum);
         void scaleWellPi(double wellPi);
-        void updateSegment(int segment_number, double center_depth);
+        void updateSegment(int segment_number, double center_depth, std::size_t seqIndex);
+	const std::size_t& getSeqIndex() const;
+	const bool& getDefaultSatTabId() const;
+	const std::size_t& getCompSegSeqIndex() const;
+	void setCompSegSeqIndex(std::size_t index);
+	void setDefaultSatTabId(bool id);
+	const double& getSegDistStart() const;
+	const double& getSegDistEnd() const;
+	void setSegDistStart(const double& distStart);
+	void setSegDistEnd(const double& distEnd);
 
         bool operator==( const Connection& ) const;
         bool operator!=( const Connection& ) const;
@@ -83,6 +98,11 @@ namespace Opm {
         double m_rw;
 
         std::array<int,3> ijk;
+	std::size_t m_seqIndex;
+	double m_segDistStart;
+	double m_segDistEnd;
+	bool m_defaultSatTabId;
+	std::size_t m_compSeg_seqIndex=0;
 
         // related segment number
         // -1 means the completion is not related to segment
@@ -91,7 +111,5 @@ namespace Opm {
     };
 }
 
-
-
 #endif /* COMPLETION_HPP_ */
 

opm/parser/eclipse/EclipseState/Schedule/Group.hpp

@@ -24,6 +24,7 @@
 #include <memory>
 #include <set>
 #include <string>
+#include <stddef.h>
 
 #include <opm/parser/eclipse/EclipseState/Schedule/ScheduleEnums.hpp>
 #include <opm/parser/eclipse/EclipseState/Runspec.hpp>
@@ -66,9 +67,10 @@ namespace Opm {
 
     class Group {
     public:
-        Group(const std::string& name, const TimeMap& timeMap , size_t creationTimeStep);
+        Group(const std::string& name, const size_t& seqIndex, const TimeMap& timeMap , size_t creationTimeStep);
         bool hasBeenDefined(size_t timeStep) const;
         const std::string& name() const;
+	const size_t& seqIndex() const;
         bool isProductionGroup(size_t timeStep) const;
         bool isInjectionGroup(size_t timeStep) const;
         void setProductionGroup(size_t timeStep, bool isProductionGroup);
@@ -133,6 +135,7 @@ namespace Opm {
     private:
         size_t m_creationTimeStep;
         std::string m_name;
+	size_t m_seqIndex;
         GroupInjection::InjectionData m_injection;
         GroupProduction::ProductionData m_production;
         DynamicState< std::set< std::string > > m_wells;

opm/parser/eclipse/EclipseState/Schedule/GroupTree.hpp

@@ -22,17 +22,21 @@
 
 #include <string>
 #include <vector>
+#include <map>
 
 namespace Opm {
 
 class GroupTree {
     public:
-        void update( const std::string& name );
+        void update( const std::string& name);
-        void update( const std::string& name, const std::string& parent );
+        void update( const std::string& name, const std::string& parent);
+	void updateSeqIndex( const std::string& name, const std::string& other_parent);
         bool exists( const std::string& group ) const;
         const std::string& parent( const std::string& name ) const;
         std::vector< std::string > children( const std::string& parent ) const;
-
+	const std::map<const std::string , size_t>& nameSeqIndMap() const;
+	const std::map<size_t, const std::string >& seqIndNameMap() const;
+	size_t groupTreeSize();
         bool operator==( const GroupTree& ) const;
         bool operator!=( const GroupTree& ) const;
 
@@ -53,6 +57,8 @@ class GroupTree {
         std::vector< group > groups = { group { "FIELD", "" } };
         friend bool operator<( const std::string&, const group& );
         std::vector< group >::iterator find( const std::string& );
+	std::map<const std::string , size_t> m_nameSeqIndMap;
+	std::map<size_t, const std::string > m_seqIndNameMap;
 };
 
 }

opm/parser/eclipse/EclipseState/Schedule/MSW/updatingConnectionsWithSegments.hpp

@@ -24,9 +24,11 @@
 
 #include <opm/parser/eclipse/EclipseState/Schedule/WellConnections.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/MSW/WellSegments.hpp>
+#include <opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp>
 
 namespace Opm {
-    WellConnections * newConnectionsWithSegments(const DeckKeyword& compsegs, const WellConnections& input_connections, const WellSegments& segments);
+    WellConnections * newConnectionsWithSegments(const DeckKeyword& compsegs, const WellConnections& input_connections, 
+						 const WellSegments& segments, const EclipseGrid& grid, std::size_t& totNC);
 }
 
 #endif

opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp

@@ -94,8 +94,11 @@ namespace Opm
 
           is an inefficient way to get all the wells defined at time
           't'.
-        */
+        */ 
+	//std::vector< const Group& > getChildGroups(const std::string& group_name, size_t timeStep) const;
+        std::vector< const Group* > getChildGroups(const std::string& group_name, size_t timeStep) const;
         std::vector< const Well* > getWells(const std::string& group, size_t timeStep) const;
+	std::vector< const Well* > getChildWells(const std::string& group_name, size_t timeStep) const;
         std::vector< const Well* > getWellsMatching( const std::string& ) const;
         const OilVaporizationProperties& getOilVaporizationProperties(size_t timestep) const;
 
@@ -103,9 +106,11 @@ namespace Opm
 
         const GroupTree& getGroupTree(size_t t) const;
         size_t numGroups() const;
+	size_t numGroups(size_t timeStep) const;
         bool hasGroup(const std::string& groupName) const;
         const Group& getGroup(const std::string& groupName) const;
         std::vector< const Group* > getGroups() const;
+	std::vector< const Group* > getGroups(size_t timeStep) const;
         const Tuning& getTuning() const;
         const MessageLimits& getMessageLimits() const;
         void invalidNamePattern (const std::string& namePattern, const ParseContext& parseContext, const DeckKeyword& keyword) const;
@@ -159,7 +164,7 @@ namespace Opm
         void handleCOMPDAT( const DeckKeyword& keyword,  size_t currentStep, const EclipseGrid& grid, const Eclipse3DProperties& eclipseProperties, const ParseContext& parseContext);
         void handleCOMPLUMP( const DeckKeyword& keyword,  size_t currentStep );
         void handleWELSEGS( const DeckKeyword& keyword, size_t currentStep);
-        void handleCOMPSEGS( const DeckKeyword& keyword, size_t currentStep);
+        void handleCOMPSEGS( const DeckKeyword& keyword, size_t currentStep, const EclipseGrid& grid);
         void handleWCONINJE( const SCHEDULESection&,  const DeckKeyword& keyword, size_t currentStep, const ParseContext& parseContext);
         void handleWPOLYMER( const DeckKeyword& keyword, size_t currentStep, const ParseContext& parseContext);
         void handleWSOLVENT( const DeckKeyword& keyword, size_t currentStep, const ParseContext& parseContext);

opm/parser/eclipse/EclipseState/Schedule/Well.hpp

@@ -35,6 +35,7 @@
 #include <opm/parser/eclipse/EclipseState/Schedule/WellProductionProperties.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/MSW/WellSegments.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/ScheduleEnums.hpp>
+#include <opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp>
 
 namespace Opm {
 
@@ -49,13 +50,15 @@ namespace Opm {
 
     class Well {
     public:
-        Well(const std::string& name, int headI,
+        Well(const std::string& name, const size_t& seqIndex, int headI,
              int headJ, double refDepth, Phase preferredPhase,
              const TimeMap& timeMap, size_t creationTimeStep,
              WellCompletion::CompletionOrderEnum completionOrdering = WellCompletion::TRACK,
              bool allowCrossFlow = true, bool automaticShutIn = true);
         const std::string& name() const;
-
+	const size_t& seqIndex() const;
+	const std::size_t getTotNoConn() const;
+	void setTotNoConn(std::size_t noConn); 
         bool hasBeenDefined(size_t timeStep) const;
         const std::string getGroupName(size_t timeStep) const;
         void setGroupName(size_t timeStep , const std::string& groupName);
@@ -184,7 +187,7 @@ namespace Opm {
         bool hasEvent(uint64_t eventMask, size_t reportStep) const;
         void handleCOMPLUMP(const DeckRecord& record, size_t time_step);
         void handleCOMPDAT(size_t time_step, const DeckRecord& record, const EclipseGrid& grid, const Eclipse3DProperties& eclipseProperties);
-        void handleCOMPSEGS(const DeckKeyword& keyword, size_t time_step);
+        void handleCOMPSEGS(const DeckKeyword& keyword, const EclipseGrid& grid, size_t time_step);
         void handleWELOPEN(const DeckRecord& record, size_t time_step, WellCompletion::StateEnum status);
         void handleWPIMULT(const DeckRecord& record, size_t time_step);
         void handleWELSEGS(const DeckKeyword& keyword, size_t time_step);
@@ -198,7 +201,9 @@ namespace Opm {
     private:
         size_t m_creationTimeStep;
         std::string m_name;
-
+	std::size_t m_seqIndex;
+	std::size_t m_totNoConn=0;
+	
         DynamicState< WellCommon::StatusEnum > m_status;
 
         DynamicState< int > m_isAvailableForGroupControl;

opm/parser/eclipse/EclipseState/Schedule/WellConnections.hpp

@@ -39,8 +39,12 @@ namespace Opm {
                            double Kh,
                            double rw,
                            const int satTableId,
-                           const WellCompletion::DirectionEnum direction = WellCompletion::DirectionEnum::Z);
+                           const WellCompletion::DirectionEnum direction = WellCompletion::DirectionEnum::Z,
-        void loadCOMPDAT(const DeckRecord& record, const EclipseGrid& grid, const Eclipse3DProperties& eclipseProperties);
+			   const std::size_t seqIndex = 0,
+			   const double segDistStart= 0.0,
+			   const double segDistEnd= 0.0,
+			   const bool defaultSatTabId = true);
+        void loadCOMPDAT(const DeckRecord& record, const EclipseGrid& grid, const Eclipse3DProperties& eclipseProperties, std::size_t& totNC);
 
         using const_iterator = std::vector< Connection >::const_iterator;
 
@@ -53,6 +57,9 @@ namespace Opm {
 
         const_iterator begin() const { return this->m_connections.begin(); }
         const_iterator end() const { return this->m_connections.end(); }
+        
+        std::size_t totNoConn() const { return this->m_connections.size(); }
+        
         void filter(const EclipseGrid& grid);
         bool allConnectionsShut() const;
         /// Order connections irrespective of input order.
@@ -81,7 +88,11 @@ namespace Opm {
                            double Kh,
                            double rw,
                            const int satTableId,
-                           const WellCompletion::DirectionEnum direction = WellCompletion::DirectionEnum::Z);
+                           const WellCompletion::DirectionEnum direction = WellCompletion::DirectionEnum::Z,
+			   const std::size_t seqIndex = 0,
+			   const double segDistStart= 0.0,
+			   const double segDistEnd= 0.0,
+			   const bool defaultSatTabId = true);
 
         std::vector< Connection > m_connections;
         size_t findClosestConnection(int oi, int oj, double oz, size_t start_pos);

src/opm/output/eclipse/AggregateConnectionData.cpp

@@ -0,0 +1,326 @@
+/*
+  Copyright 2018 Statoil ASA
+
+  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 <opm/output/eclipse/AggregateConnectionData.hpp>
+
+#include <opm/output/eclipse/VectorItems/connection.hpp>
+#include <opm/output/eclipse/VectorItems/intehead.hpp>
+
+#include <opm/output/data/Wells.hpp>
+
+#include <opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/ScheduleEnums.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
+
+#include <opm/parser/eclipse/Units/UnitSystem.hpp>
+
+#include <cassert>
+#include <cstddef>
+#include <iostream>
+#include <cmath>
+
+namespace VI = Opm::RestartIO::Helpers::VectorItems;
+
+// #####################################################################
+// Class Opm::RestartIO::Helpers::AggregateConnectionData
+// ---------------------------------------------------------------------
+
+namespace {
+    std::size_t numWells(const std::vector<int>& inteHead)
+    {
+        return inteHead[VI::intehead::NWELLS];
+    }
+
+    std::size_t maxNumConn(const std::vector<int>& inteHead)
+    {
+        return inteHead[VI::intehead::NCWMAX];
+    }
+
+    std::map <std::size_t, const Opm::Connection*>  mapSeqIndexToConnection(const Opm::WellConnections& conns)
+    {
+	// make seqIndex to Connection map
+	std::map <std::size_t, const Opm::Connection*> seqIndConnMap;
+	for (const auto & conn : conns) {
+	    std::size_t sI = conn.getSeqIndex();
+	    seqIndConnMap.insert(std::make_pair(sI, &conn));
+	}
+	return seqIndConnMap;
+    }
+
+    std::map <std::size_t, const Opm::Connection*>  mapCompSegSeqIndexToConnection(const Opm::WellConnections& conns)
+    {
+	// make CompSegSeqIndex to Connection map
+	std::map <std::size_t, const Opm::Connection*> cs_seqIndConnMap;
+	for (const auto & conn : conns) {
+	    std::size_t sI = conn.getCompSegSeqIndex();
+	    cs_seqIndConnMap.insert(std::make_pair(sI, &conn));
+	}
+	return cs_seqIndConnMap;
+    }
+
+
+    template <class ConnOp>
+    void connectionLoop(const std::vector<const Opm::Well*>& wells,
+                        const Opm::EclipseGrid&              grid,
+                        const std::size_t                    sim_step,
+                        ConnOp&&                             connOp)
+    {
+        for (auto nWell = wells.size(), wellID = 0*nWell;
+             wellID < nWell; ++wellID)
+        {
+            const auto* well = wells[wellID];
+
+            if (well == nullptr) { continue; }
+            
+            const auto& conns = well->getActiveConnections(sim_step, grid);
+	    const int niSI = static_cast<int>(well->getTotNoConn());
+	    std::map <std::size_t, const Opm::Connection*> sIToConn;
+
+	    //Branch according to MSW well or not and 
+	    //sort active connections according to appropriate seqIndex
+	    if (well->isMultiSegment(sim_step)) {
+		//sort connections according to input sequence in COMPSEGS
+		sIToConn = mapCompSegSeqIndexToConnection(conns);
+	    } else 
+	    {
+		//sort connections according to input sequence in COMPDAT
+		sIToConn = mapSeqIndexToConnection(conns);
+	    }
+	    std::vector<const Opm::Connection*> connSI;
+	    for (int iSI = 0; iSI < niSI; iSI++) {
+		const auto searchSI = sIToConn.find(static_cast<std::size_t>(iSI));
+		if (searchSI != sIToConn.end()) {		  
+		  connSI.push_back(searchSI->second);
+		}
+	    }
+	    for (auto nConn = connSI.size(), connID = 0*nConn;
+                 connID < nConn; ++connID)
+            {
+                connOp(*well, wellID, *(connSI[connID]), connID);
+            }
+        }
+    }
+
+    namespace IConn {
+        std::size_t entriesPerConn(const std::vector<int>& inteHead)
+        {
+            return inteHead[VI::intehead::NICONZ];
+        }
+
+        Opm::RestartIO::Helpers::WindowedMatrix<int>
+        allocate(const std::vector<int>& inteHead)
+        {
+            using WM = Opm::RestartIO::Helpers::WindowedMatrix<int>;
+
+            return WM {
+                WM::NumRows   { numWells(inteHead) },
+                WM::NumCols   { maxNumConn(inteHead) },
+                WM::WindowSize{ entriesPerConn(inteHead) }
+            };
+        }
+
+        template <class IConnArray>
+        void staticContrib(const Opm::Connection& conn,
+                           const std::size_t      connID,
+                           IConnArray&            iConn)
+        {
+            using ConnState = ::Opm::WellCompletion::StateEnum;
+            using Ix = ::Opm::RestartIO::Helpers::VectorItems::IConn::index;
+
+            iConn[Ix::SeqIndex] = connID + 1;
+
+            iConn[Ix::CellI] = conn.getI() + 1;
+            iConn[Ix::CellJ] = conn.getJ() + 1;
+            iConn[Ix::CellK] = conn.getK() + 1;
+
+            iConn[Ix::ConnStat] = (conn.state() == ConnState::OPEN)
+                ? 1 : -1000;
+
+            iConn[Ix::Drainage] = conn.getDefaultSatTabId()
+                ? 0 : conn.satTableId();
+
+            // Don't support differing sat-func tables for
+            // draining and imbibition curves at connections.
+            iConn[Ix::Imbibition] = iConn[Ix::Drainage];
+
+            //iConn[Ix::ComplNum] = std::abs(conn.complnum);
+            iConn[Ix::ComplNum] = iConn[Ix::SeqIndex];
+
+            iConn[Ix::ConnDir] = conn.dir();
+            iConn[Ix::Segment] = conn.attachedToSegment()
+                ? conn.segment() : 0;
+        }
+    } // IConn
+
+    namespace SConn {
+        std::size_t entriesPerConn(const std::vector<int>& inteHead)
+        {
+            return inteHead[VI::intehead::NSCONZ];
+        }
+
+        Opm::RestartIO::Helpers::WindowedMatrix<float>
+        allocate(const std::vector<int>& inteHead)
+        {
+            using WM = Opm::RestartIO::Helpers::WindowedMatrix<float>;
+
+            return WM {
+                WM::NumRows   { numWells(inteHead) },
+                WM::NumCols   { maxNumConn(inteHead) },
+                WM::WindowSize{ entriesPerConn(inteHead) }
+            };
+        }
+
+        template <class SConnArray>
+        void staticContrib(const Opm::Connection& conn,
+                           const Opm::UnitSystem& units,
+                           SConnArray&            sConn)
+        {
+            using M  = ::Opm::UnitSystem::measure;
+            using Ix = ::Opm::RestartIO::Helpers::VectorItems::SConn::index;
+
+            auto scprop = [&units](const M u, const double x) -> float
+            {
+                return static_cast<float>(units.from_si(u, x));
+            };
+
+	    sConn[Ix::ConnTrans] =
+                        scprop(M::transmissibility, conn.CF());
+
+            sConn[Ix::Depth]    = scprop(M::length, conn.depth());
+            sConn[Ix::Diameter] = scprop(M::length, 2*conn.rw());
+
+	    sConn[Ix::EffectiveKH] =
+                        scprop(M::effective_Kh, conn.Kh());
+
+            sConn[Ix::item12] = sConn[Ix::ConnTrans];
+
+            sConn[Ix::SegDistEnd]   = scprop(M::length, conn.getSegDistEnd());
+            sConn[Ix::SegDistStart] = scprop(M::length, conn.getSegDistStart());
+
+            sConn[Ix::item30] = -1.0e+20f;
+            sConn[Ix::item31] = -1.0e+20f;
+        }
+    } // SConn
+
+    namespace XConn {
+        std::size_t entriesPerConn(const std::vector<int>& inteHead)
+        {
+            return inteHead[VI::intehead::NXCONZ];
+        }
+
+        Opm::RestartIO::Helpers::WindowedMatrix<double>
+        allocate(const std::vector<int>& inteHead)
+        {
+            using WM = Opm::RestartIO::Helpers::WindowedMatrix<double>;
+
+            return WM {
+                WM::NumRows   { numWells(inteHead) },
+                WM::NumCols   { maxNumConn(inteHead) },
+                WM::WindowSize{ entriesPerConn(inteHead) }
+            };
+        }
+
+        template <class XConnArray>
+        void dynamicContrib(const Opm::data::Connection& x,
+                            const Opm::UnitSystem&       units,
+                            XConnArray&                  xConn)
+        {
+            using M  = ::Opm::UnitSystem::measure;
+            using Ix = ::Opm::RestartIO::Helpers::VectorItems::XConn::index;
+            using R  = ::Opm::data::Rates::opt;
+
+            // Note flow rate sign.  Treat production rates as positive.
+            const auto& Q = x.rates;
+
+            if (Q.has(R::oil)) {
+                xConn[Ix::OilRate] =
+                    - units.from_si(M::liquid_surface_rate, Q.get(R::oil));
+            }
+
+            if (Q.has(R::wat)) {
+                xConn[Ix::WaterRate] =
+                    - units.from_si(M::liquid_surface_rate, Q.get(R::wat));
+            }
+
+            if (Q.has(R::gas)) {
+                xConn[Ix::WaterRate] =
+                    - units.from_si(M::gas_surface_rate, Q.get(R::gas));
+            }
+
+            xConn[Ix::ResVRate] = 0.0;
+
+            if (Q.has(R::reservoir_oil)) {
+                xConn[Ix::ResVRate] -=
+                    units.from_si(M::rate, Q.get(R::reservoir_oil));
+            }
+
+            if (Q.has(R::reservoir_water)) {
+                xConn[Ix::ResVRate] -=
+                    units.from_si(M::rate, Q.get(R::reservoir_water));
+            }
+
+            if (Q.has(R::reservoir_gas)) {
+                xConn[Ix::ResVRate] -=
+                    units.from_si(M::rate, Q.get(R::reservoir_gas));
+            }
+        }
+    } // XConn
+} // Anonymous
+
+Opm::RestartIO::Helpers::AggregateConnectionData::
+AggregateConnectionData(const std::vector<int>& inteHead)
+    : iConn_(IConn::allocate(inteHead))
+    , sConn_(SConn::allocate(inteHead))
+    , xConn_(XConn::allocate(inteHead))
+{}
+
+// ---------------------------------------------------------------------
+
+void
+Opm::RestartIO::Helpers::AggregateConnectionData::
+captureDeclaredConnData(const Schedule&        sched,
+                        const EclipseGrid&     grid,
+                        const UnitSystem&      units,
+                        const data::WellRates& xw,
+                        const std::size_t      sim_step)
+{
+    const auto& wells = sched.getWells(sim_step);
+
+    connectionLoop(wells, grid, sim_step, [&units, &xw, this]
+        (const Well&       well, const std::size_t wellID,
+         const Connection& conn, const std::size_t connID) -> void
+    {
+        auto ic = this->iConn_(wellID, connID);
+        auto sc = this->sConn_(wellID, connID);
+
+        IConn::staticContrib(conn, connID, ic);
+        SConn::staticContrib(conn, units, sc);
+
+        auto xi = xw.find(well.name());
+        if ((xi != xw.end()) &&
+            (connID < xi->second.connections.size()))
+        {
+            auto xc = this->xConn_(wellID, connID);
+
+            XConn::dynamicContrib(xi->second.connections[connID],
+                                  units, xc);
+        }
+    });
+}

src/opm/output/eclipse/AggregateGroupData.cpp

@@ -0,0 +1,624 @@
+/*
+  Copyright 2018 Statoil ASA
+
+  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 <opm/output/eclipse/AggregateGroupData.hpp>
+
+#include <opm/output/eclipse/SummaryState.hpp>
+#include <opm/output/eclipse/WriteRestartHelpers.hpp>
+
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+#include <opm/parser/eclipse/EclipseState/Runspec.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/ScheduleEnums.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/GroupTree.hpp>
+
+#include <algorithm>
+#include <cstddef>
+#include <cstring>
+#include <string>
+#include <iostream>
+
+// #####################################################################
+// Class Opm::RestartIO::Helpers::AggregateGroupData
+// ---------------------------------------------------------------------
+
+
+namespace {
+
+    int nigrpz(const std::vector<int>& inteHead)
+    {
+        // INTEHEAD(36) = NIGRPZ
+        return inteHead[36];
+    }
+
+    int nsgrpz(const std::vector<int>& inteHead)
+    {
+        return inteHead[37];
+    }
+
+    int nxgrpz(const std::vector<int>& inteHead)
+    {
+        return inteHead[38];
+    }
+
+    int nzgrpz(const std::vector<int>& inteHead)
+    {
+        return inteHead[39];
+    }
+
+    // maximum number of groups
+    int ngmaxz(const std::vector<int>& inteHead)
+    {
+        return inteHead[20];
+    }
+
+    // maximum number of wells in any group
+    int nwgmax(const std::vector<int>& inteHead)
+    {
+        return inteHead[19];
+    }
+
+    std::string trim(const std::string& s)
+    {
+        const auto b = s.find_first_not_of(" \t");
+        if (b == std::string::npos) {
+            // No blanks.  Return unchanged.
+            return s;
+        }
+
+        auto e = s.find_last_not_of(" \t");
+        if (e == std::string::npos) {
+            // No trailing blanks.  Return [b, end)
+            return s.substr(b, std::string::npos);
+        }
+
+        // Remove leading/trailing blanks.
+        return s.substr(b, e - b + 1);
+    }
+
+    std::vector<std::string>
+    groupNames(const std::vector<const Opm::Group*>& groups)
+    {
+        auto gnms = std::vector<std::string>{};
+        gnms.reserve(groups.size());
+
+        for (const auto* group : groups) {
+            gnms.push_back(trim(group->name()));
+        }
+
+        return gnms;
+    }
+
+    const int groupType(const Opm::Schedule& sched,
+		      const Opm::Group& group,
+		      const std::size_t simStep)
+    {
+	const std::string& groupName = group.name();
+	if (!sched.hasGroup(groupName))
+            throw std::invalid_argument("No such group: " + groupName);
+        {
+            if (group.hasBeenDefined( simStep )) {
+                const auto& groupTree = sched.getGroupTree( simStep );
+                const auto& childGroups = groupTree.children( groupName );
+
+                if (childGroups.size()) {
+		    return 1;
+                }
+                else {
+		    return 0;
+		}
+            }
+            return 0;
+        }
+    }
+
+    template <typename GroupOp>
+    void groupLoop(const std::vector<const Opm::Group*>& groups,
+                  GroupOp&&                             groupOp)
+    {
+        auto groupID = std::size_t{0};
+        for (const auto* group : groups) {
+            groupID += 1;
+
+            if (group == nullptr) { continue; }
+
+            groupOp(*group, groupID - 1);
+        }
+    }
+
+    std::map <size_t, const Opm::Group*>  currentGroupMapIndexGroup(const Opm::Schedule& sched, const size_t simStep, const std::vector<int>& inteHead)
+    {
+	const auto& groups = sched.getGroups(simStep);
+	// make group index for current report step
+	std::map <size_t, const Opm::Group*> indexGroupMap;
+	for (const auto* group : groups) {
+	    int ind = (group->name() == "FIELD")
+	    ? ngmaxz(inteHead)-1 : group->seqIndex()-1;
+	    const std::pair<size_t, const Opm::Group*> groupPair = std::make_pair(static_cast<size_t>(ind), group); 
+	    indexGroupMap.insert(groupPair);
+	}
+	return indexGroupMap;
+    }
+
+    std::map <const std::string, size_t>  currentGroupMapNameIndex(const Opm::Schedule& sched, const size_t simStep, const std::vector<int>& inteHead)
+    {
+	const auto& groups = sched.getGroups(simStep);
+	// make group name to index map for the current time step
+	std::map <const std::string, size_t> groupIndexMap;
+	for (const auto* group : groups) {
+	    int ind = (group->name() == "FIELD")
+                    ? ngmaxz(inteHead)-1 : group->seqIndex()-1;
+	    std::pair<const std::string, size_t> groupPair = std::make_pair(group->name(), ind); 
+	    groupIndexMap.insert(groupPair);
+	}
+	return groupIndexMap;
+    }
+
+    const std::map <std::string, size_t>  currentWellIndex(const Opm::Schedule& sched, const size_t simStep)
+    {
+	const auto& wells = sched.getWells(simStep);
+	// make group index for current report step
+	std::map <std::string, size_t> wellIndexMap;
+	for (const auto* well : wells) {
+	    std::pair<std::string, size_t> wellPair = std::make_pair(well->name(), well->seqIndex()); 
+	    wellIndexMap.insert(wellPair);
+	}
+	return wellIndexMap;
+    }
+
+    const int currentGroupLevel(const Opm::Schedule& sched, const Opm::Group& group, const size_t simStep)
+    {
+	int level = 0;
+	const std::vector< const Opm::Group* >  groups = sched.getGroups(simStep);
+      	const std::string& groupName = group.name();
+	if (!sched.hasGroup(groupName))
+            throw std::invalid_argument("No such group: " + groupName);
+        {
+            if (group.hasBeenDefined( simStep )) {
+                const auto& groupTree = sched.getGroupTree( simStep );
+		//find group level - field level is 0
+		std::string tstGrpName = groupName;
+		while (((tstGrpName.size())>0) && (!(tstGrpName=="FIELD"))) {
+		    std::string curParent = groupTree.parent(tstGrpName);
+		    level+=1;
+		    tstGrpName = curParent;
+		}
+		return level;
+	    }
+	    else {
+		throw std::invalid_argument("actual group has not been defined at report time: " + simStep );
+	    }
+	}
+	return level;
+    }
+
+
+    namespace IGrp {
+        std::size_t entriesPerGroup(const std::vector<int>& inteHead)
+        {
+            // INTEHEAD[36] = NIGRPZ
+            return inteHead[36];
+        }
+
+        Opm::RestartIO::Helpers::WindowedArray<int>
+        allocate(const std::vector<int>& inteHead)
+        {
+            using WV = Opm::RestartIO::Helpers::WindowedArray<int>;
+
+            return WV {
+                WV::NumWindows{ ngmaxz(inteHead) },
+                WV::WindowSize{ entriesPerGroup(inteHead) }
+            };
+        }
+
+        template <class IGrpArray>
+        void staticContrib(const Opm::Schedule&    sched,
+			   const Opm::Group&       group,
+			   const int               nwgmax,
+			   const int               ngmaxz,
+			   const std::size_t       simStep,
+			   IGrpArray&              iGrp,
+			   const std::vector<int>& inteHead)
+        {
+	  // find the number of wells or child groups belonging to a group and store in 
+	  // location nwgmax +1 in the iGrp array
+
+	    const auto childGroups = sched.getChildGroups(group.name(), simStep);
+	    const auto childWells  = sched.getChildWells(group.name(), simStep);
+	    const auto groupMapNameIndex =  currentGroupMapNameIndex(sched, simStep, inteHead);
+	    const auto mapIndexGroup = currentGroupMapIndexGroup(sched, simStep, inteHead);
+	    if ((childGroups.size() != 0) && (childWells.size()!=0))
+	      throw std::invalid_argument("group has both wells and child groups" + group.name());
+            int igrpCount = 0;
+	    if (childWells.size() != 0) {
+		//group has child wells
+		for ( auto it = childWells.begin() ; it != childWells.end(); it++) {
+		    iGrp[igrpCount] = (*it)->seqIndex()+1;
+		    igrpCount+=1;
+		}
+	    }
+	    else if (childGroups.size() != 0) {
+		//group has child groups
+		//The field group always has seqIndex = 0 because it is always defined first
+	        //Hence the all groups except the Field group uses the seqIndex assigned
+		std::vector<size_t> childGroupIndex;
+		Opm::RestartIO::Helpers::groupMaps groupMap;
+		groupMap.currentGrpTreeNameSeqIndMap(sched, simStep, groupMapNameIndex,mapIndexGroup);
+		const auto indGroupMap = groupMap.indexGroupMap();
+		const auto gNameIndMap = groupMap.groupNameIndexMap();
+		for (auto* grp : childGroups) {
+		  auto groupName = grp->name();
+		  auto searchGTName = gNameIndMap.find(groupName);
+		  if (searchGTName != gNameIndMap.end()) {
+		      childGroupIndex.push_back(searchGTName->second);
+		  }
+		  else {
+		      throw std::invalid_argument( "Invalid group name" );
+		  }
+		}
+		std::sort(childGroupIndex.begin(), childGroupIndex.end());
+		for (auto groupTreeIndex : childGroupIndex) {
+		    auto searchGTIterator = indGroupMap.find(groupTreeIndex);
+		    if (searchGTIterator != indGroupMap.end()) {
+			auto gname = (searchGTIterator->second)->name();
+			auto gSeqNoItr = groupMapNameIndex.find(gname);
+			if (gSeqNoItr != groupMapNameIndex.end()) {
+			    iGrp[igrpCount] = (gSeqNoItr->second) + 1;
+			    igrpCount+=1;
+			}
+			else {
+			    std::cout << "AggregateGroupData, ChildGroups - Group name: groupMapNameIndex: " << gSeqNoItr->first << std::endl;
+			    throw std::invalid_argument( "Invalid group name" );
+			}
+		    }
+		    else {
+			std::cout << "AggregateGroupData, ChildGroups - GroupTreeIndex: " << groupTreeIndex << std::endl;
+			throw std::invalid_argument( "Invalid GroupTree index" );
+		    }
+		}
+	    }
+
+	    //assign the number of child wells or child groups to
+	    // location nwgmax
+	    iGrp[nwgmax] =  (childGroups.size() == 0)
+                    ? childWells.size() : childGroups.size();
+
+	    // find the group type (well group (type 0) or node group (type 1) and store the type in  
+	    // location nwgmax + 26
+	    const auto grpType = groupType(sched, group, simStep);
+	    iGrp[nwgmax+26] = grpType;
+
+	    //find group level ("FIELD" is level 0) and store the level in
+	    //location nwgmax + 27
+	    const auto grpLevel = currentGroupLevel(sched, group, simStep);
+	    iGrp[nwgmax+27] = grpLevel;
+	    
+	    // set values for group probably connected to GCONPROD settings
+	    //
+	    if (group.name() != "FIELD")
+	    {
+		iGrp[nwgmax+ 5] = -1;
+		iGrp[nwgmax+12] = -1;
+		iGrp[nwgmax+17] = -1;
+		iGrp[nwgmax+22] = -1;
+
+		//assign values to group number (according to group sequence)
+		iGrp[nwgmax+88] = group.seqIndex();
+		iGrp[nwgmax+89] = group.seqIndex();
+		iGrp[nwgmax+95] = group.seqIndex();
+		iGrp[nwgmax+96] = group.seqIndex();
+	    }
+	    else
+	    {
+		//assign values to group number (according to group sequence)
+		iGrp[nwgmax+88] = ngmaxz;
+		iGrp[nwgmax+89] = ngmaxz;
+		iGrp[nwgmax+95] = ngmaxz;
+		iGrp[nwgmax+96] = ngmaxz;
+	    }
+	    //find parent group and store index of parent group in
+	    //location nwgmax + 28
+	    const auto& groupTree = sched.getGroupTree( simStep );
+	    const std::string& parent = groupTree.parent(group.name());
+	    if (group.name() == "FIELD")
+	      iGrp[nwgmax+28] = 0;
+	    else {
+		if (parent.size() == 0)
+		    throw std::invalid_argument("parent group does not exist" + group.name());
+		const auto itr = groupMapNameIndex.find(parent);
+		iGrp[nwgmax+28] = (itr->second)+1;
+	    }
+        }
+    } // Igrp
+
+    namespace SGrp {
+        std::size_t entriesPerGroup(const std::vector<int>& inteHead)
+        {
+            // INTEHEAD[37] = NSGRPZ
+            return inteHead[37];
+        }
+
+        Opm::RestartIO::Helpers::WindowedArray<float>
+        allocate(const std::vector<int>& inteHead)
+        {
+            using WV = Opm::RestartIO::Helpers::WindowedArray<float>;
+
+            return WV {
+                WV::NumWindows{ ngmaxz(inteHead) },
+                WV::WindowSize{ entriesPerGroup(inteHead) }
+            };
+        }
+
+	template <class SGrpArray>
+	void staticContrib(SGrpArray& sGrp)
+	{
+		const auto dflt   = -1.0e+20f;
+		const auto dflt_2 = -2.0e+20f;
+		const auto infty  =  1.0e+20f;
+		const auto zero   =  0.0f;
+		const auto one    =  1.0f;
+
+		const auto init = std::vector<float> { // 112 Items (0..111)
+		    // 0     1      2      3      4
+		    infty, infty, dflt , infty,  zero ,     //   0..  4  ( 0)
+		    zero , infty, infty, infty , infty,     //   5..  9  ( 1)
+		    infty, infty, infty, infty , dflt ,     //  10.. 14  ( 2)
+		    infty, infty, infty, infty , dflt ,     //  15.. 19  ( 3)
+		    infty, infty, infty, infty , dflt ,     //  20.. 24  ( 4)
+		    zero , zero , zero , dflt_2, zero ,     //  24.. 29  ( 5)
+		    zero , zero , zero , zero  , zero ,     //  30.. 34  ( 6)
+		    infty ,zero , zero , zero  , infty,     //  35.. 39  ( 7)
+		    zero , zero , zero , zero  , zero ,     //  40.. 44  ( 8)
+		    zero , zero , zero , zero  , zero ,     //  45.. 49  ( 9)
+		    zero , infty, infty, infty , infty,     //  50.. 54  (10)
+		    infty, infty, infty, infty , infty,     //  55.. 59  (11)
+		    infty, infty, infty, infty , infty,     //  60.. 64  (12)
+		    infty, infty, infty, infty , zero ,     //  65.. 69  (13)
+		    zero , zero , zero , zero  , zero ,     //  70.. 74  (14)
+		    zero , zero , zero , zero  , infty,     //  75.. 79  (15)
+		    infty, zero , infty, zero  , zero ,     //  80.. 84  (16)
+		    zero , zero , zero , zero  , zero ,     //  85.. 89  (17)
+		    zero , zero , one  , zero  , zero ,     //  90.. 94  (18)
+		    zero , zero , zero , zero  , zero ,     //  95.. 99  (19)
+		    zero , zero , zero , zero  , zero ,     // 100..104  (20)
+		    zero , zero , zero , zero  , zero ,     // 105..109  (21)
+		    zero , zero                             // 110..111  (22)
+		};
+
+		const auto sz = static_cast<
+		    decltype(init.size())>(sGrp.size());
+
+		auto b = std::begin(init);
+		auto e = b + std::min(init.size(), sz);
+
+		std::copy(b, e, std::begin(sGrp));
+	}
+    } // SGrp
+
+    namespace XGrp {
+        std::size_t entriesPerGroup(const std::vector<int>& inteHead)
+        {
+            // INTEHEAD[38] = NXGRPZ
+            return inteHead[38];
+        }
+
+        Opm::RestartIO::Helpers::WindowedArray<double>
+        allocate(const std::vector<int>& inteHead)
+        {
+            using WV = Opm::RestartIO::Helpers::WindowedArray<double>;
+
+            return WV {
+                WV::NumWindows{ ngmaxz(inteHead) },
+                WV::WindowSize{ entriesPerGroup(inteHead) }
+            };
+        }
+
+        std::vector<std::string>
+        filter_cumulative(const bool     ecl_compatible_rst,
+                          const std::vector<std::string>& keys)
+        {
+            if (ecl_compatible_rst) {
+                // User wants ECLIPSE-compatible output.  Write all vectors.
+                return keys;
+            }
+
+            auto ret = std::vector<std::string>{};
+            ret.reserve(keys.size());
+
+            for (const auto& key : keys) {
+                if ((key[3] == 'T') && ((key[2] == 'I') || (key[2] == 'P'))) {
+                    // Don't write cumulative quantities in case of 
+                    continue;
+                }
+
+                ret.push_back(key);
+            }
+
+            return ret;
+        }
+
+        // here define the dynamic group quantities to be written to the restart file
+        template <class XGrpArray>
+        void dynamicContrib(const std::vector<std::string>&      restart_group_keys,
+			    const std::vector<std::string>&      restart_field_keys,
+			    const std::map<std::string, size_t>& groupKeyToIndex,
+			    const std::map<std::string, size_t>& fieldKeyToIndex,
+			    const Opm::Group&                    group,
+			    const Opm::SummaryState&             sumState,
+			    const bool                           ecl_compatible_rst,
+			    XGrpArray&                           xGrp)
+        {
+	  std::string groupName = group.name();
+	  const std::vector<std::string>& keys = (groupName == "FIELD")
+	  ? restart_field_keys : restart_group_keys;
+	  const std::map<std::string, size_t>& keyToIndex = (groupName == "FIELD")
+	  ? fieldKeyToIndex : groupKeyToIndex;
+
+	  for (const auto key : filter_cumulative(ecl_compatible_rst,keys)) {
+	      std::string compKey = (groupName == "FIELD")
+	      ? key : key + ":" + groupName;
+
+	      if (sumState.has(compKey)) {
+		  double keyValue = sumState.get(compKey);
+		  const auto itr = keyToIndex.find(key);
+		  xGrp[itr->second] = keyValue;
+	      }
+	      /*else {
+		  std::cout << "AggregateGroupData_compkey: " << compKey << std::endl;
+		  std::cout << "AggregateGroupData, empty " << std::endl;
+		  //throw std::invalid_argument("No such keyword: " + compKey);
+	    }*/
+	  }
+	}
+    }
+
+    namespace ZGrp {
+        std::size_t entriesPerGroup(const std::vector<int>& inteHead)
+        {
+            // INTEHEAD[39] = NZGRPZ
+            return inteHead[39];
+        }
+
+        Opm::RestartIO::Helpers::WindowedArray<
+            Opm::RestartIO::Helpers::CharArrayNullTerm<8>
+        >
+        allocate(const std::vector<int>& inteHead)
+        {
+            using WV = Opm::RestartIO::Helpers::WindowedArray<
+                Opm::RestartIO::Helpers::CharArrayNullTerm<8>
+            >;
+
+            return WV {
+                WV::NumWindows{ ngmaxz(inteHead) },
+                WV::WindowSize{ entriesPerGroup(inteHead) }
+            };
+        }
+    }
+} // Anonymous
+
+void
+Opm::RestartIO::Helpers::groupMaps::
+currentGrpTreeNameSeqIndMap(const Opm::Schedule&                        sched,
+                            const size_t                                simStep,
+                            const std::map<const std::string , size_t>& GnIMap,
+                            const std::map<size_t, const Opm::Group*>&  IGMap)
+    {
+	const auto& grpTreeNSIMap = (sched.getGroupTree(simStep)).nameSeqIndMap();
+	const auto& grpTreeSINMap = (sched.getGroupTree(simStep)).seqIndNameMap();
+	// make group index for current report step
+	size_t index = 0;
+	for (auto itr = IGMap.begin(); itr != IGMap.end(); itr++) {
+	    auto name = (itr->second)->name();
+	    auto srchGN = grpTreeNSIMap.find(name);
+	    if (srchGN != grpTreeNSIMap.end()) {
+		//come here if group is in gruptree
+		this->m_indexGroupMap.insert(std::make_pair(srchGN->second, itr->second));
+		this->m_groupNameIndexMap.insert(std::make_pair(srchGN->first, srchGN->second));
+	    }
+	    else {
+	    //come here if group is not in gruptree to put in from global group list
+		this->m_indexGroupMap.insert(std::make_pair(itr->first, itr->second));
+		this->m_groupNameIndexMap.insert(std::make_pair(name, itr->first));
+	    }
+	}
+    }
+
+const std::map <size_t, const Opm::Group*>&
+Opm::RestartIO::Helpers::groupMaps::indexGroupMap() const
+{
+	return this->m_indexGroupMap;
+}
+
+const std::map <const std::string, size_t>&
+Opm::RestartIO::Helpers::groupMaps::groupNameIndexMap() const
+{
+	return this->m_groupNameIndexMap;
+}
+
+// =====================================================================
+
+Opm::RestartIO::Helpers::AggregateGroupData::
+AggregateGroupData(const std::vector<int>& inteHead)
+    : iGroup_ (IGrp::allocate(inteHead))
+    , sGroup_ (SGrp::allocate(inteHead))
+    , xGroup_ (XGrp::allocate(inteHead))
+    , zGroup_ (ZGrp::allocate(inteHead))
+    , nWGMax_ (nwgmax(inteHead))
+    , nGMaxz_ (ngmaxz(inteHead))
+{}
+
+// ---------------------------------------------------------------------
+
+void
+Opm::RestartIO::Helpers::AggregateGroupData::
+captureDeclaredGroupData(const Opm::Schedule&                 sched,
+			 const std::vector<std::string>&      restart_group_keys,
+			 const std::vector<std::string>&      restart_field_keys,
+			 const std::map<std::string, size_t>& groupKeyToIndex,
+			 const std::map<std::string, size_t>& fieldKeyToIndex,
+			 const bool                           ecl_compatible_rst,
+			 const std::size_t                    simStep,
+			 const Opm::SummaryState&             sumState,
+			 const std::vector<int>&              inteHead)
+{
+    const auto indexGroupMap = currentGroupMapIndexGroup(sched, simStep, inteHead);
+    const auto nameIndexMap = currentGroupMapNameIndex(sched, simStep, inteHead);
+
+    std::vector<const Opm::Group*> curGroups;
+    curGroups.resize(ngmaxz(inteHead));
+    //
+    //initialize curgroups
+    for (auto* cg : curGroups) cg = nullptr;
+
+    auto it = indexGroupMap.begin();
+    while (it != indexGroupMap.end())
+	{
+	    curGroups[static_cast<int>(it->first)] = it->second;
+	    it++;
+	}
+    {
+	groupLoop(curGroups, [sched, simStep, inteHead, this]
+            (const Group& group, const std::size_t groupID) -> void
+	    {
+		auto ig = this->iGroup_[groupID];
+		IGrp::staticContrib(sched, group, this->nWGMax_, this->nGMaxz_, simStep, ig, inteHead);
+	    });
+    }
+
+    // Define Static Contributions to SGrp Array.
+    groupLoop(curGroups,
+        [this](const Group& group, const std::size_t groupID) -> void
+    {
+        auto sw = this->sGroup_[groupID];
+        SGrp::staticContrib(sw);
+    });
+
+    // Define DynamicContributions to XGrp Array.
+    groupLoop(curGroups,
+        [restart_group_keys, restart_field_keys, groupKeyToIndex, fieldKeyToIndex, ecl_compatible_rst, sumState, this]
+	(const Group& group, const std::size_t groupID) -> void
+    {
+        auto xg = this->xGroup_[groupID];
+        XGrp::dynamicContrib( restart_group_keys, restart_field_keys, groupKeyToIndex, fieldKeyToIndex, group, sumState, ecl_compatible_rst, xg);
+    });
+
+}
+
+// ---------------------------------------------------------------------

src/opm/output/eclipse/AggregateMSWData.cpp

@@ -0,0 +1,638 @@
+/*
+  Copyright 2018 Statoil ASA
+
+  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 <opm/output/eclipse/AggregateMSWData.hpp>
+
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+
+#include <opm/parser/eclipse/EclipseState/Runspec.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/ScheduleEnums.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Connection.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/WellConnections.hpp>
+
+#include <opm/parser/eclipse/Units/UnitSystem.hpp>
+
+//#include <opm/output/data/Wells.hpp>
+
+#include <algorithm>
+#include <cstddef>
+#include <cstring>
+#include <string>
+#include <iostream>
+#include <limits>
+
+// #####################################################################
+// Class Opm::RestartIO::Helpers::AggregateMSWData
+// ---------------------------------------------------------------------
+
+namespace {
+
+    std::size_t nummsw(const std::vector<int>& inteHead)
+    {
+        // inteHead(174) = NSEGWL
+        return inteHead[174];
+    }
+
+    std::size_t nswlmx(const std::vector<int>& inteHead)
+    {
+        // inteHead(175) = NSWLMX
+        return inteHead[175];
+    }
+
+    std::size_t nisegz(const std::vector<int>& inteHead)
+    {
+        // inteHead(178) = NISEGZ
+        return inteHead[178];
+    }
+
+    std::size_t nrsegz(const std::vector<int>& inteHead)
+    {
+        // inteHead(179) = NRSEGZ
+        return inteHead[179];
+    }
+
+    std::size_t nilbrz(const std::vector<int>& inteHead)
+    {
+        // inteHead(180) = NILBRZ
+        return inteHead[180];
+    }
+
+    Opm::RestartIO::Helpers::BranchSegmentPar
+    getBranchSegmentParam(const Opm::WellSegments& segSet, const int branch)
+    {
+	int noSegInBranch = 0;
+	int firstSeg = -1;
+	int lastSeg  = -1;
+	int outletS = 0;
+	for (std::size_t segNo = 1; segNo <= segSet.size(); segNo++) {
+	    auto segInd = segSet.segmentNumberToIndex(segNo);
+	    auto i_branch = segSet[segInd].branchNumber();
+	    auto i_outS = segSet[segInd].outletSegment();
+	    if (i_branch == branch) {
+		noSegInBranch +=1;
+		if (firstSeg < 0) {
+		    firstSeg = segNo;
+		    outletS = (branch > 1) ? i_outS : 0;
+		}
+		lastSeg = segNo;
+	    }
+	}
+
+	return {
+	  outletS,
+	  noSegInBranch,
+	  firstSeg,
+	  lastSeg,
+	  branch
+	};
+    }
+
+    std::vector<std::size_t> SegmentSetBranches(const Opm::WellSegments& segSet) {
+	std::vector<std::size_t> branches;
+	for (std::size_t segNo = 1; segNo <= segSet.size(); segNo++) {
+	    auto segInd = segSet.segmentNumberToIndex(segNo);
+	    auto i_branch = segSet[segInd].branchNumber();
+	    if (std::find(branches.begin(), branches.end(), i_branch) == branches.end()) { 
+		branches.push_back(i_branch);
+	    }
+	}
+	return branches;
+    }
+
+    int firstSegmentInBranch(const Opm::WellSegments& segSet, const int branch) {
+	int firstSegNo = 0;
+	std::size_t segNo = 0;
+	while ((segNo <= segSet.size()) && (firstSegNo == 0)) {
+	    segNo+=1;
+	    auto segInd = segSet.segmentNumberToIndex(segNo);
+	    auto i_branch = segSet[segInd].branchNumber();
+	    if (branch == i_branch) {
+		firstSegNo = segNo;
+	    }
+	}
+	return firstSegNo;
+    }
+
+     int noConnectionsSegment(const Opm::WellConnections& compSet,
+                             const Opm::WellSegments&    segSet,
+                             const std::size_t              segIndex)
+    {
+	auto segNumber  = segSet[segIndex].segmentNumber();
+	int noConnections = 0;
+	for (auto it : compSet) {
+	    auto cSegment = it.segment();
+	    if (segNumber == cSegment) {
+		noConnections+=1;
+	    }
+	}
+
+	return noConnections;
+    }
+    
+    int sumConnectionsSegment(const Opm::WellConnections& compSet,
+                                 const Opm::WellSegments&    segSet,
+                                 const std::size_t              segIndex)
+    {
+      // This function returns (for a given segment) the sum of number of connections for each segment
+      // with lower segment index than the currnet segment
+      // If the segment contains no connections, the number returned is zero.
+	int sumConn = 0;
+	if (noConnectionsSegment(compSet, segSet, segIndex) > 0) {
+	// add up the number of connections for å segments with lower segment index than current segment
+	    	    auto segNumber  = segSet[segIndex].segmentNumber();
+	    for (std::size_t i_seg = 1; i_seg <= segNumber; i_seg++) {
+		auto  ind = segSet.segmentNumberToIndex(i_seg);
+		int addCon = (ind == static_cast<int>(segIndex)) ? 1 : noConnectionsSegment(compSet, segSet, ind);
+		sumConn += addCon;
+	    }
+	}
+	return sumConn;
+    }
+
+    
+    std::vector<std::size_t>
+    inflowSegmentsIndex(const Opm::WellSegments& segSet, std::size_t segIndex) {
+	auto segNumber  = segSet[segIndex].segmentNumber();
+	std::vector<std::size_t> inFlowSegNum;
+	for (std::size_t ind = 0; ind < segSet.size(); ind++) {
+	    auto i_outletSeg = segSet[ind].outletSegment();
+	    if (segNumber == i_outletSeg) {
+		inFlowSegNum.push_back(ind);
+	    }
+	}
+	return inFlowSegNum;
+    }
+
+    int noInFlowBranches(const Opm::WellSegments& segSet, std::size_t segIndex) {
+	auto segNumber  = segSet[segIndex].segmentNumber();
+	auto branch     = segSet[segIndex].branchNumber();
+	int noIFBr = 0;
+	for (std::size_t ind = 0; ind < segSet.size(); ind++) {
+	    auto o_segNum = segSet[ind].outletSegment();
+	    auto i_branch = segSet[ind].branchNumber();
+	    if ((segNumber == o_segNum) && (branch != i_branch)){
+		noIFBr+=1;
+	    }
+	}
+	return noIFBr;
+    }
+    //find the number of inflow branches (different from the current branch)
+    int sumNoInFlowBranches(const Opm::WellSegments& segSet, std::size_t segIndex) {
+	int sumIFB = 0;
+	auto segBranch     = segSet[segIndex].branchNumber();
+	const auto iSInd = inflowSegmentsIndex(segSet, segIndex);
+	for (auto ind : iSInd) {
+	    auto inflowBranch = segSet[ind].branchNumber();
+	    // if inflow segment belongs to different branch add contribution
+	    if (segBranch != inflowBranch) {
+		sumIFB+=1;
+		// search recursively down this branch to find more inflow branches
+		sumIFB += sumNoInFlowBranches(segSet, ind);
+	    }
+	}
+	return sumIFB;
+    }
+    
+
+     std::vector<std::size_t> segmentOrder(const Opm::WellSegments& segSet, const std::size_t segIndex) {
+	std::vector<std::size_t> ordSegNumber;
+	std::vector<std::size_t> tempOrdVect;
+	std::vector<std::size_t> segIndCB;
+	// Store "heel" segment since that will not always be at the end of the list
+	segIndCB.push_back(segIndex);
+	std::size_t sInd = segIndex;
+	std::size_t newSInd = segIndex;
+	auto origBranchNo = segSet[segIndex].branchNumber();
+	bool endOrigBranch = true;
+	//std::cout << "SegmentOrder-segIndex:"  << segIndex  << " origBranchno: " <<  origBranchNo  << std::endl;
+	// loop down branch to find all segments in branch and number from "toe" to "heel"
+	while (newSInd < segSet.size()) {
+	       endOrigBranch = true;
+		const auto iSInd = inflowSegmentsIndex(segSet, newSInd);
+		//std::cout << " SO- inflowSegmentsIndex:" <<   std::endl;
+		for (auto isi : iSInd )  {
+			auto inflowBranch = segSet[isi].branchNumber();
+			if (origBranchNo == inflowBranch) {
+			    endOrigBranch = false;
+			}
+			//std::cout << " SO- isi:" <<  isi  << std::endl;
+		}
+		if (iSInd.size() > 0) {
+		    for (auto ind : iSInd) {
+			auto inflowBranch = segSet[ind].branchNumber();
+			if (origBranchNo == inflowBranch) {
+			    // if inflow segment belongs to same branch add contribution
+			    segIndCB.insert(segIndCB.begin(), ind);
+			    // search recursively down this branch to find more inflow branches
+			    newSInd = ind;
+			    //std::cout << "SO-ind-loop: origB=iflowB - ind:" <<  ind  << std::endl;
+			}
+			else {
+			    // if inflow segment belongs to different branch, start new search
+			    //std::cout << "SO-ind-loop: origB!=iflowB - ind:" <<  ind  << std::endl;
+			    auto nSOrd = segmentOrder(segSet, ind);
+			    // copy the segments already found and indexed into the total ordered segment vector
+			    for (std::size_t indOS = 0; indOS < nSOrd.size(); indOS++) {
+				ordSegNumber.push_back(nSOrd[indOS]);
+				//std::cout << "SO-ind-loop: origB!=iflowB - indOS:" <<  indOS  << " nSOrd[indOS] " << nSOrd[indOS]  << std::endl;
+			    }
+			    if (endOrigBranch) {
+			       newSInd = segSet.size();
+			    }
+			    // increment the local branch sequence number counter
+			}
+		    }
+		}
+		if (endOrigBranch || (iSInd.size()==0)) {
+		    // have come to toe of current branch - store segment indicies of current branch
+		    //std::cout << "SO-Toe of current branch - newSInd :" <<  newSInd  << std::endl;
+		    for (std::size_t indOS = 0; indOS < segIndCB.size(); indOS++) {
+			ordSegNumber.push_back(segIndCB[indOS]);
+			//std::cout << "SO  end CB - indOS:" <<  indOS  << " segIndCB[indOS] " << segIndCB[indOS]  << std::endl;
+		    }
+		    // set new index to exit while loop
+		    newSInd = segSet.size();
+		}
+	    //}
+	}
+
+	if (origBranchNo == 1) {
+	    // make the vector of ordered segments
+	    //std::cout << "SO-OrBr=1 -ordSegNumber.size():" <<  ordSegNumber.size()  << std::endl;
+	    tempOrdVect.resize(ordSegNumber.size());
+	    for (std::size_t ov_ind = 0; ov_ind < ordSegNumber.size(); ov_ind++) {
+		tempOrdVect[ordSegNumber[ov_ind]] = ov_ind+1;
+		//std::cout << "SO_OrBr=1- ov_ind:" <<  ov_ind  << " ordSegNumber[ov_ind] " << ordSegNumber[ov_ind]  << std::endl;
+	    }
+	    return tempOrdVect;
+	} else {
+	    return ordSegNumber;
+	}
+    }
+
+    
+    int inflowSegmentCurBranch(const Opm::WellSegments& segSet, std::size_t segIndex) {
+	auto branch = segSet[segIndex].branchNumber();
+	auto segNumber  = segSet[segIndex].segmentNumber();
+	int inFlowSegNum = 0;
+	for (std::size_t ind = 0; ind < segSet.size(); ind++) {
+	    auto i_segNum = segSet[ind].segmentNumber();
+	    auto i_branch = segSet[ind].branchNumber();
+	    auto i_outFlowSeg = segSet[ind].outletSegment();
+	    if ((branch == i_branch) && (segNumber == i_outFlowSeg)) {
+		if (inFlowSegNum == 0) {
+		    inFlowSegNum = i_segNum;
+		}
+		else {
+		  std::cout << "Non-unique inflow segment in same branch, Well: " << segSet.wellName() << std::endl;
+		  std::cout <<  "Segment number: " << segNumber << std::endl;
+		  std::cout <<  "Branch number: " << branch << std::endl;
+		  std::cout <<  "Inflow segment number 1: " << inFlowSegNum << std::endl;
+		  std::cout <<  "Inflow segment number 2: " << segSet[ind].segmentNumber() << std::endl;
+		  throw std::invalid_argument("Non-unique inflow segment in same branch, Well " + segSet.wellName());
+		}
+	    }
+	}
+	return inFlowSegNum;
+    }
+
+    template <typename MSWOp>
+    void MSWLoop(const std::vector<const Opm::Well*>& wells,
+                  MSWOp&&                             mswOp)
+    {
+        auto mswID = std::size_t{0};
+        for (const auto* well : wells) {
+            mswID += 1;
+
+            if (well == nullptr) { continue; }
+
+            mswOp(*well, mswID - 1);
+        }
+    }
+
+    namespace ISeg {
+        std::size_t entriesPerMSW(const std::vector<int>& inteHead)
+        {
+            // inteHead(176) = NSEGMX
+            // inteHead(178) = NISEGZ
+            return inteHead[176] * inteHead[178];
+        }
+
+        Opm::RestartIO::Helpers::WindowedArray<int>
+        allocate(const std::vector<int>& inteHead)
+        {
+            using WV = Opm::RestartIO::Helpers::WindowedArray<int>;
+
+            return WV {
+                WV::NumWindows{ nswlmx(inteHead) },
+                WV::WindowSize{ entriesPerMSW(inteHead) }
+            };
+        }
+
+        template <class ISegArray>
+        void staticContrib(const Opm::Well&        well,
+                           const std::size_t       rptStep,
+                           const std::vector<int>& inteHead,
+			   const Opm::EclipseGrid&  grid,
+                           ISegArray&              iSeg
+			  )
+        {
+            if (well.isMultiSegment(rptStep)) {
+		//loop over segment set and print out information
+		auto welSegSet = well.getWellSegments(rptStep);
+		auto completionSet = well.getConnections(rptStep);
+		auto noElmSeg = nisegz(inteHead);
+		std::size_t segmentInd = 0;
+		auto orderedSegmentNo = segmentOrder(welSegSet, segmentInd);
+		for (std::size_t ind_seg = 1; ind_seg <= welSegSet.size(); ind_seg++) {
+		    auto ind = welSegSet.segmentNumberToIndex(ind_seg);
+		    auto iS = (ind_seg-1)*noElmSeg;
+		    iSeg[iS + 0] = orderedSegmentNo[ind];
+		    iSeg[iS + 1] = welSegSet[ind].outletSegment();
+		    iSeg[iS + 2] = inflowSegmentCurBranch(welSegSet, ind);
+		    iSeg[iS + 3] = welSegSet[ind].branchNumber();
+		    iSeg[iS + 4] = noInFlowBranches(welSegSet, ind);
+		    iSeg[iS + 5] = sumNoInFlowBranches(welSegSet, ind);
+		    iSeg[iS + 6] = noConnectionsSegment(completionSet, welSegSet, ind);
+		    iSeg[iS + 7] = sumConnectionsSegment(completionSet, welSegSet, ind);
+		    iSeg[iS + 8] = iSeg[iS+0];
+		}
+	    }
+	    else {
+	      throw std::invalid_argument("No such multisegment well: " + well.name());
+	    }
+        }
+    } // ISeg
+
+    namespace RSeg {
+        std::size_t entriesPerMSW(const std::vector<int>& inteHead)
+        {
+            // inteHead(176) = NSEGMX
+            // inteHead(179) = NRSEGZ
+            return inteHead[176] * inteHead[179];
+        }
+
+        Opm::RestartIO::Helpers::WindowedArray<double>
+        allocate(const std::vector<int>& inteHead)
+        {
+            using WV = Opm::RestartIO::Helpers::WindowedArray<double>;
+
+            return WV {
+                WV::NumWindows{ nswlmx(inteHead) },
+                WV::WindowSize{ entriesPerMSW(inteHead) }
+            };
+        }
+
+        template <class RSegArray>
+        void staticContrib(const Opm::Well&        well,
+                           const std::size_t       rptStep,
+                           const std::vector<int>& inteHead,
+			   const Opm::UnitSystem& units,
+                           RSegArray&              rSeg)
+        {
+	    if (well.isMultiSegment(rptStep)) {
+		using M = ::Opm::UnitSystem::measure;
+		//loop over segment set and print out information
+		auto welSegSet = well.getWellSegments(rptStep);
+		auto completionSet = well.getCompletions(rptStep);
+		auto noElmSeg = nrsegz(inteHead);
+		    //treat the top segment individually
+		    int ind_seg = 1;
+		    auto ind = welSegSet.segmentNumberToIndex(ind_seg);
+		    rSeg[0] = units.from_si(M::length, welSegSet.lengthTopSegment());
+		    rSeg[1] = units.from_si(M::length, welSegSet.depthTopSegment());
+		    rSeg[5] = units.from_si(M::volume, welSegSet.volumeTopSegment());
+		    rSeg[6] = rSeg[0];
+		    rSeg[7] = rSeg[1];
+		    // set item ind + 10 to 0.5 based on tests on E100
+		    rSeg[10] = 0.5;
+
+		    //  segment pressure (to be added!!)
+		    rSeg[ 39] = 0;
+
+		    //Default values
+		    rSeg[ 39] = 1.0;
+
+		    rSeg[105] = 1.0;
+		    rSeg[106] = 1.0;
+		    rSeg[107] = 1.0;
+		    rSeg[108] = 1.0;
+		    rSeg[109] = 1.0;
+		    rSeg[110] = 1.0;
+
+		//Treat subsequent segments
+		for (std::size_t ind_seg = 2; ind_seg <= welSegSet.size(); ind_seg++) {
+		    // set the elements of the rSeg array
+		    auto ind = welSegSet.segmentNumberToIndex(ind_seg);
+		    auto outSeg = welSegSet[ind].outletSegment();
+		    auto ind_ofs = welSegSet.segmentNumberToIndex(outSeg);
+		    auto iS = (ind_seg-1)*noElmSeg;
+		    rSeg[iS +   0] = units.from_si(M::length, (welSegSet[ind].totalLength() - welSegSet[ind_ofs].totalLength()));
+		    rSeg[iS +   1] = units.from_si(M::length, (welSegSet[ind].depth() - welSegSet[ind_ofs].depth()));
+		    rSeg[iS +   2] = units.from_si(M::length, (welSegSet[ind].internalDiameter()));
+		    rSeg[iS +   3] = units.from_si(M::length, (welSegSet[ind].roughness()));
+		    rSeg[iS +   4] = units.from_si(M::length, (welSegSet[ind].crossArea()));
+		    //repeat unit conversion to convert for area instead of length
+		    rSeg[iS +   4] = units.from_si(M::length, rSeg[iS +   4]);
+		    rSeg[iS +   5] = units.from_si(M::volume, (welSegSet[ind].volume()));
+		    rSeg[iS +   6] = units.from_si(M::length, (welSegSet[ind].totalLength()));
+		    rSeg[iS +   7] = units.from_si(M::length, (welSegSet[ind].depth()));
+
+		    // set item ind + 10 to 0.5 based on tests on E100
+		    rSeg[10] = 0.5;
+		    //  segment pressure (to be added!!)
+		    rSeg[iS +  11] = 0;
+
+		    //Default values
+		    rSeg[iS +  39] = 1.0;
+
+		    rSeg[iS + 105] = 1.0;
+		    rSeg[iS + 106] = 1.0;
+		    rSeg[iS + 107] = 1.0;
+		    rSeg[iS + 108] = 1.0;
+		    rSeg[iS + 109] = 1.0;
+		    rSeg[iS + 110] = 1.0;
+		}
+	    }
+	    else {
+	      throw std::invalid_argument("No such multisegment well: " + well.name());
+	    }
+        }
+    } // RSeg
+
+    namespace ILBS {
+        std::size_t entriesPerMSW(const std::vector<int>& inteHead)
+        {
+            // inteHead(177) = NLBRMX
+            return inteHead[177];
+        }
+
+        Opm::RestartIO::Helpers::WindowedArray<int>
+        allocate(const std::vector<int>& inteHead)
+        {
+            using WV = Opm::RestartIO::Helpers::WindowedArray<int>;
+
+            return WV {
+                WV::NumWindows{ nswlmx(inteHead) },
+                WV::WindowSize{ entriesPerMSW(inteHead) }
+            };
+        }
+
+        template <class ILBSArray>
+        void staticContrib(const Opm::Well&  well,
+                           const std::size_t rptStep,
+                           ILBSArray&        iLBS)
+        {
+	    if (well.isMultiSegment(rptStep)) {
+		//
+		auto welSegSet = well.getWellSegments(rptStep);
+		auto branches = SegmentSetBranches(welSegSet);
+		for (auto it = branches.begin()+1; it != branches.end(); it++){
+		    iLBS[*it-2] = firstSegmentInBranch(welSegSet, *it);
+		}
+	    }
+	    else {
+	      throw std::invalid_argument("No such multisegment well: " + well.name());
+	    }
+        }
+    } // ILBS
+
+    namespace ILBR {
+        std::size_t entriesPerMSW(const std::vector<int>& inteHead)
+        {
+            // inteHead(177) = NLBRMX
+            // inteHead(180) = NILBRZ
+            return inteHead[177] * inteHead[180];
+        }
+
+        Opm::RestartIO::Helpers::WindowedArray<int>
+        allocate(const std::vector<int>& inteHead)
+        {
+            using WV = Opm::RestartIO::Helpers::WindowedArray<int>;
+
+            return WV {
+                WV::NumWindows{ nswlmx(inteHead) },
+                WV::WindowSize{ entriesPerMSW(inteHead) }
+            };
+        }
+
+        template <class ILBRArray>
+        void staticContrib(const Opm::Well&  well,
+                           const std::size_t rptStep,
+			   const std::vector<int>& inteHead,
+                           ILBRArray&        iLBR)
+        {
+	    if (well.isMultiSegment(rptStep)) {
+		//
+		auto welSegSet = well.getWellSegments(rptStep);
+		auto branches = SegmentSetBranches(welSegSet);
+		auto noElmBranch = nilbrz(inteHead);
+		for (auto it = branches.begin(); it != branches.end(); it++){
+		    auto iB = (*it-1)*noElmBranch;
+		    auto branchParam = getBranchSegmentParam(welSegSet, *it);
+		    iLBR[iB  ] = branchParam.outletS;
+		    iLBR[iB+1] = branchParam.noSegInBranch;
+		    iLBR[iB+2] = branchParam.firstSeg;
+		    iLBR[iB+3] = branchParam.lastSeg;
+		    iLBR[iB+4] = branchParam.branch - 1;
+		}
+	    }
+	    else {
+	      throw std::invalid_argument("No such multisegment well: " + well.name());
+	    }
+        }
+    } // ILBR
+
+} // Anonymous
+
+// =====================================================================
+
+Opm::RestartIO::Helpers::AggregateMSWData::
+AggregateMSWData(const std::vector<int>& inteHead)
+    : iSeg_ (ISeg::allocate(inteHead))
+    , rSeg_ (RSeg::allocate(inteHead))
+    , iLBS_ (ILBS::allocate(inteHead))
+    , iLBR_ (ILBR::allocate(inteHead))
+{}
+
+// ---------------------------------------------------------------------
+
+void
+Opm::RestartIO::Helpers::AggregateMSWData::
+captureDeclaredMSWData(const Schedule&         sched,
+                       const std::size_t       rptStep,
+		       const Opm::UnitSystem& units,
+                       const std::vector<int>& inteHead,
+		       const Opm::EclipseGrid&  grid)
+{
+    const auto& wells = sched.getWells(rptStep);
+    auto msw = std::vector<const Opm::Well*>{};
+
+    //msw.reserve(wells.size());
+    for (const auto well : wells) {
+	if (well->isMultiSegment(rptStep)) msw.push_back(well);
+    }
+
+    // Extract Contributions to ISeg Array
+    {
+        MSWLoop(msw, [rptStep, inteHead, grid, this]
+            (const Well& well, const std::size_t mswID) -> void
+        {
+            auto imsw = this->iSeg_[mswID];
+
+            ISeg::staticContrib(well, rptStep, inteHead, grid, imsw);
+        });
+    }
+
+    // Extract Contributions to RSeg Array
+    {
+        MSWLoop(msw, [&units, rptStep, inteHead, this]
+            (const Well& well, const std::size_t mswID) -> void
+        {
+            auto rmsw = this->rSeg_[mswID];
+
+            RSeg::staticContrib(well, rptStep, inteHead, units, rmsw);
+        });
+    }
+
+    // Extract Contributions to ILBS Array
+    {
+        MSWLoop(msw, [rptStep, this]
+            (const Well& well, const std::size_t mswID) -> void
+        {
+            auto ilbs_msw = this->iLBS_[mswID];
+
+            ILBS::staticContrib(well, rptStep, ilbs_msw);
+        });
+    }
+
+    // Extract Contributions to ILBR Array
+    {
+        MSWLoop(msw, [rptStep, inteHead, this]
+            (const Well& well, const std::size_t mswID) -> void
+        {
+            auto ilbr_msw = this->iLBR_[mswID];
+
+            ILBR::staticContrib(well, rptStep, inteHead, ilbr_msw);
+        });
+    }
+}

src/opm/output/eclipse/AggregateWellData.cpp

@@ -31,6 +31,7 @@
 #include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/ScheduleEnums.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
+//#include <opm/output/data/Wells.hpp>
 
 #include <opm/parser/eclipse/Units/UnitSystem.hpp>
 #include <opm/parser/eclipse/Units/Units.hpp>
@@ -119,7 +120,37 @@ namespace {
             };
         }
 
-        int groupIndex(const std::string&              grpName,
+	std::map <const std::string, size_t>  currentGroupMapNameIndex(const Opm::Schedule& sched, const size_t simStep, const std::vector<int>& inteHead)
+	{
+	    const auto& groups = sched.getGroups(simStep);
+	    // make group name to index map for the current time step
+	    std::map <const std::string, size_t> groupIndexMap;
+	    for (const auto* group : groups) {
+		int ind = (group->name() == "FIELD")
+			? inteHead[VI::intehead::NGMAXZ]-1 : group->seqIndex()-1;
+		std::pair<const std::string, size_t> groupPair = std::make_pair(group->name(), ind); 
+		groupIndexMap.insert(groupPair);
+	    }
+	    return groupIndexMap;
+      }
+    
+      int groupIndex(const std::string&              grpName,
+                     const std::map <const std::string, size_t>&  currentGroupMapNameIndex)
+      {
+	    int ind = 0;
+            auto searchGTName = currentGroupMapNameIndex.find(grpName);
+	    if (searchGTName != currentGroupMapNameIndex.end()) 
+	    {
+		ind = searchGTName->second + 1;	
+	    }
+	    else 
+	    {
+		std::cout << "group Name: " << grpName << std::endl;
+		throw std::invalid_argument( "Invalid group name" );
+	    }
+	return ind;
+      }
+        /*int groupIndex(const std::string&              grpName,
                        const std::vector<std::string>& groupNames,
                        const int                       maxGroups)
         {
@@ -141,7 +172,7 @@ namespace {
 
             // One-based indexing.
             return std::distance(b, i) + 1;
-        }
+        }  */
 
         int wellType(const Opm::Well&  well,
                      const std::size_t sim_step)
@@ -264,7 +295,8 @@ namespace {
         template <class IWellArray>
         void staticContrib(const Opm::Well&                well,
                            const std::size_t               msWellID,
-                           const std::vector<std::string>& groupNames,
+                           const std::map <const std::string, size_t>&  GroupMapNameInd,
+			   /*const std::vector<std::string>& groupNames,*/
                            const int                       maxGroups,
                            const std::size_t               sim_step,
                            IWellArray&                     iWell)
@@ -279,16 +311,26 @@ namespace {
                 const auto& conn = well.getConnections(sim_step);
 
                 iWell[Ix::NConn]  = static_cast<int>(conn.size());
-                iWell[Ix::FirstK] = (iWell[Ix::NConn] == 0)
+		
+		if (well.isMultiSegment(sim_step))  
+		{
+		    // Set top and bottom connections to zero for multi segment wells
+		  iWell[Ix::FirstK] = 0;
+		  iWell[Ix::LastK] = 0;
+		}
+		else
+		{
+		  iWell[Ix::FirstK] = (iWell[Ix::NConn] == 0)
                     ? 0 : conn.get(0).getK() + 1;
 
-                iWell[Ix::LastK] = (iWell[Ix::NConn] == 0)
+		  iWell[Ix::LastK] = (iWell[Ix::NConn] == 0)
-                    ? 0 : conn.get(conn.size() - 1).getK() + 1;
+		    ? 0 : conn.get(conn.size() - 1).getK() + 1;
+		}
             }
 
             iWell[Ix::Group] =
-                groupIndex(trim(well.getGroupName(sim_step)),
+	     groupIndex(trim(well.getGroupName(sim_step)),
-                           groupNames, maxGroups);
+                GroupMapNameInd);
 
             iWell[Ix::WType]  = wellType  (well, sim_step);
             iWell[Ix::WCtrl]  = ctrlMode  (well, sim_step);
@@ -433,6 +475,7 @@ namespace {
         void staticContrib(const Opm::Well&       well,
                            const Opm::UnitSystem& units,
                            const std::size_t      sim_step,
+			   const ::Opm::SummaryState& smry,
                            SWellArray&            sWell)
         {
             using Ix = ::Opm::RestartIO::Helpers::VectorItems::SWell::index;
@@ -442,58 +485,77 @@ namespace {
             {
                 return static_cast<float>(units.from_si(u, x));
             };
-
+	    
             assignDefaultSWell(sWell);
 
             if (well.isProducer(sim_step)) {
                 const auto& pp = well.getProductionProperties(sim_step);
 
-                using PP = ::Opm::WellProducer::ControlModeEnum;
+                if (pp.OilRate != 0.0) {
-
-                if (pp.hasProductionControl(PP::ORAT)) {
                     sWell[Ix::OilRateTarget] =
                         swprop(M::liquid_surface_rate, pp.OilRate);
-                }
+		}
 
-                if (pp.hasProductionControl(PP::WRAT)) {
+                if (pp.WaterRate != 0.0) {
                     sWell[Ix::WatRateTarget] =
                         swprop(M::liquid_surface_rate, pp.WaterRate);
                 }
 
-                if (pp.hasProductionControl(PP::GRAT)) {
+                if (pp.GasRate != 0.0) {
                     sWell[Ix::GasRateTarget] =
                         swprop(M::gas_surface_rate, pp.GasRate);
                 }
 
-                if (pp.hasProductionControl(PP::LRAT)) {
+                if (pp.LiquidRate != 0.0) {
                     sWell[Ix::LiqRateTarget] =
                         swprop(M::liquid_surface_rate, pp.LiquidRate);
                 }
-                else if (pp.hasProductionControl(PP::ORAT) &&
+                else  {
-                         pp.hasProductionControl(PP::WRAT))
-                {
                     sWell[Ix::LiqRateTarget] =
                         swprop(M::liquid_surface_rate, pp.OilRate + pp.WaterRate);
                 }
 
-                if (pp.hasProductionControl(PP::RESV)) {
+                if (pp.ResVRate != 0.0) {
                     sWell[Ix::ResVRateTarget] =
                         swprop(M::rate, pp.ResVRate);
                 }
+		else if (smry.has("WVPR:" + well.name())) {
+                    // Write out summary voidage production rate if
+                    // target/limit is not set
+                    sWell[Ix::ResVRateTarget] =
+                        static_cast<float>(smry.get("WVPR:" + well.name()));
+                }
 
-                if (pp.hasProductionControl(PP::THP)) {
+                if (pp.THPLimit != 0.0) {
                     sWell[Ix::THPTarget] =
                         swprop(M::pressure, pp.THPLimit);
                 }
-
+           
-                sWell[Ix::BHPTarget] = pp.hasProductionControl(PP::BHP)
+		sWell[Ix::BHPTarget] = pp.BHPLimit != 0.0
                     ? swprop(M::pressure, pp.BHPLimit)
-                    : swprop(M::pressure, 100.0e3*::Opm::unit::psia);
+                    : swprop(M::pressure, 1.0*::Opm::unit::atm);
             }
             else if (well.isInjector(sim_step)) {
                 const auto& ip = well.getInjectionProperties(sim_step);
 
                 using IP = ::Opm::WellInjector::ControlModeEnum;
+		using IT = ::Opm::WellInjector::TypeEnum;
+
+		if (ip.hasInjectionControl(IP::RATE)) {
+		    if (ip.injectorType == IT::OIL) {
+			sWell[Ix::OilRateTarget] = swprop(M::liquid_surface_rate, ip.surfaceInjectionRate);
+		    }
+		    if (ip.injectorType == IT::WATER) {
+			sWell[Ix::WatRateTarget] = swprop(M::liquid_surface_rate, ip.surfaceInjectionRate);
+		    }
+		    if (ip.injectorType == IT::GAS) {
+			sWell[Ix::GasRateTarget] = swprop(M::gas_surface_rate, ip.surfaceInjectionRate);
+		    }
+                }
+
+		if (ip.hasInjectionControl(IP::RESV)) {
+		    sWell[Ix::ResVRateTarget] = swprop(M::rate, ip.reservoirInjectionRate);
+                }		
 
                 if (ip.hasInjectionControl(IP::THP)) {
                     sWell[Ix::THPTarget] = swprop(M::pressure, ip.THPLimit);
@@ -501,7 +563,7 @@ namespace {
 
                 sWell[Ix::BHPTarget] = ip.hasInjectionControl(IP::BHP)
                     ? swprop(M::pressure, ip.BHPLimit)
-                    : swprop(M::pressure, 1.0*::Opm::unit::atm);
+                    : swprop(M::pressure, 1.0E05*::Opm::unit::psia);
             }
 
             sWell[Ix::DatumDepth] =
@@ -548,13 +610,16 @@ namespace {
         template <class XWellArray>
         void assignProducer(const std::string&         well,
                             const ::Opm::SummaryState& smry,
+			    const bool  ecl_compatible_rst, 
                             XWellArray&                xWell)
         {
             using Ix = ::Opm::RestartIO::Helpers::VectorItems::XWell::index;
 
             auto get = [&smry, &well](const std::string& vector)
             {
-                return smry.get(vector + ':' + well);
+                const auto key = vector + ':' + well;
+
+                return smry.has(key) ? smry.get(key) : 0.0;
             };
 
             xWell[Ix::OilPrRate] = get("WOPR");
@@ -570,11 +635,12 @@ namespace {
             xWell[Ix::WatCut]  = get("WWCT");
             xWell[Ix::GORatio] = get("WGOR");
 
-            xWell[Ix::OilPrTotal]  = get("WOPT");
+	    if (ecl_compatible_rst) {
-            xWell[Ix::WatPrTotal]  = get("WWPT");
+		xWell[Ix::OilPrTotal]  = get("WOPT");
-            xWell[Ix::GasPrTotal]  = get("WGPT");
+		xWell[Ix::WatPrTotal]  = get("WWPT");
-            xWell[Ix::VoidPrTotal] = get("WVPT");
+		xWell[Ix::GasPrTotal]  = get("WGPT");
-
+		xWell[Ix::VoidPrTotal] = get("WVPT");
+	    }
             // Not fully characterised.
             xWell[Ix::item37] = xWell[Ix::WatPrRate];
             xWell[Ix::item38] = xWell[Ix::GasPrRate];
@@ -583,13 +649,16 @@ namespace {
         template <class XWellArray>
         void assignWaterInjector(const std::string&         well,
                                  const ::Opm::SummaryState& smry,
+				 const bool  ecl_compatible_rst, 
                                  XWellArray&                xWell)
         {
             using Ix = ::Opm::RestartIO::Helpers::VectorItems::XWell::index;
 
             auto get = [&smry, &well](const std::string& vector)
             {
-                return smry.get(vector + ':' + well);
+                const auto key = vector + ':' + well;
+
+                return smry.has(key) ? smry.get(key) : 0.0;
             };
 
             // Injection rates reported as negative, cumulative
@@ -598,8 +667,10 @@ namespace {
             xWell[Ix::LiqPrRate] = xWell[Ix::WatPrRate];
 
             xWell[Ix::FlowBHP] = get("WBHP");
-
+  
-            xWell[Ix::WatInjTotal] = get("WWIT");
+	    if (ecl_compatible_rst) {
+		xWell[Ix::WatInjTotal] = get("WWIT");
+	    }
 
             xWell[Ix::item37] = xWell[Ix::WatPrRate];
             xWell[Ix::item82] = xWell[Ix::WatInjTotal];
@@ -610,13 +681,16 @@ namespace {
         template <class XWellArray>
         void assignGasInjector(const std::string&         well,
                                const ::Opm::SummaryState& smry,
+			       const bool  ecl_compatible_rst, 
                                XWellArray&                xWell)
         {
             using Ix = ::Opm::RestartIO::Helpers::VectorItems::XWell::index;
 
             auto get = [&smry, &well](const std::string& vector)
             {
-                return smry.get(vector + ':' + well);
+                const auto key = vector + ':' + well;
+
+                return smry.has(key) ? smry.get(key) : 0.0;
             };
 
             // Injection rates reported as negative production rates,
@@ -626,10 +700,13 @@ namespace {
 
             xWell[Ix::FlowBHP] = get("WBHP");
 
-            xWell[Ix::GasInjTotal] = get("WGIT");
+	    if (ecl_compatible_rst) {
+		xWell[Ix::GasInjTotal] = get("WGIT");
+	    }
 
-            xWell[Ix::GasFVF] = xWell[Ix::VoidPrRate]
+            xWell[Ix::GasFVF] = (std::abs(xWell[Ix::GasPrRate]) > 0.0)
-                              / xWell[Ix::GasPrRate];
+                ? xWell[Ix::VoidPrRate] / xWell[Ix::GasPrRate]
+                : 0.0;
 
             // Not fully characterised.
             xWell[Ix::item38] = xWell[Ix::GasPrRate];
@@ -642,10 +719,11 @@ namespace {
         void dynamicContrib(const ::Opm::Well&         well,
                             const ::Opm::SummaryState& smry,
                             const std::size_t          sim_step,
+			    const bool                 ecl_compatible_rst,
                             XWellArray&                xWell)
         {
             if (well.isProducer(sim_step)) {
-                assignProducer(well.name(), smry, xWell);
+                assignProducer(well.name(), smry, ecl_compatible_rst, xWell);
             }
             else if (well.isInjector(sim_step)) {
                 using IType = ::Opm::WellInjector::TypeEnum;
@@ -659,16 +737,16 @@ namespace {
                     break;
 
                 case IType::WATER:
-                    assignWaterInjector(well.name(), smry, xWell);
+                    assignWaterInjector(well.name(), smry, ecl_compatible_rst, xWell);
                     break;
 
                 case IType::GAS:
-                    assignGasInjector(well.name(), smry, xWell);
+                    assignGasInjector(well.name(), smry, ecl_compatible_rst, xWell);
                     break;
 
                 case IType::MULTI:
-                    assignWaterInjector(well.name(), smry, xWell);
+                    assignWaterInjector(well.name(), smry, ecl_compatible_rst, xWell);
-                    assignGasInjector  (well.name(), smry, xWell);
+                    assignGasInjector  (well.name(), smry, ecl_compatible_rst, xWell);
                     break;
                 }
             }
@@ -726,33 +804,36 @@ void
 Opm::RestartIO::Helpers::AggregateWellData::
 captureDeclaredWellData(const Schedule&   sched,
                         const UnitSystem& units,
-                        const std::size_t sim_step)
+                        const std::size_t sim_step,
+			const ::Opm::SummaryState&  smry,
+			const std::vector<int>& inteHead)
 {
     const auto& wells = sched.getWells(sim_step);
 
     // Static contributions to IWEL array.
     {
-        const auto grpNames = groupNames(sched.getGroups());
+        //const auto grpNames = groupNames(sched.getGroups());
+	const auto groupMapNameIndex = IWell::currentGroupMapNameIndex(sched, sim_step, inteHead);
         auto msWellID       = std::size_t{0};
 
-        wellLoop(wells, [&grpNames, &msWellID, sim_step, this]
+        wellLoop(wells, [&groupMapNameIndex, &msWellID, sim_step, this]
             (const Well& well, const std::size_t wellID) -> void
         {
             msWellID += well.isMultiSegment(sim_step);  // 1-based index.
             auto iw   = this->iWell_[wellID];
 
-            IWell::staticContrib(well, msWellID, grpNames,
+            IWell::staticContrib(well, msWellID, groupMapNameIndex,
                                  this->nWGMax_, sim_step, iw);
         });
     }
 
     // Static contributions to SWEL array.
-    wellLoop(wells, [&units, sim_step, this]
+    wellLoop(wells, [&units, sim_step, &smry, this]
         (const Well& well, const std::size_t wellID) -> void
     {
         auto sw = this->sWell_[wellID];
 
-        SWell::staticContrib(well, units, sim_step, sw);
+        SWell::staticContrib(well, units, sim_step, smry, sw);
     });
 
     // Static contributions to XWEL array.
@@ -780,6 +861,7 @@ void
 Opm::RestartIO::Helpers::AggregateWellData::
 captureDynamicWellData(const Schedule&             sched,
                        const std::size_t           sim_step,
+		       const bool ecl_compatible_rst,
                        const Opm::data::WellRates& xw,
                        const ::Opm::SummaryState&  smry)
 {
@@ -801,11 +883,11 @@ captureDynamicWellData(const Schedule&             sched,
     });
 
     // Dynamic contributions to XWEL array.
-    wellLoop(wells, [this, sim_step, &smry]
+    wellLoop(wells, [this, sim_step, ecl_compatible_rst, &smry]
         (const Well& well, const std::size_t wellID) -> void
     {
         auto xw = this->xWell_[wellID];
 
-        XWell::dynamicContrib(well, smry, sim_step, xw);
+        XWell::dynamicContrib(well, smry, sim_step, ecl_compatible_rst, xw);
     });
 }

src/opm/output/eclipse/CreateInteHead.cpp

@@ -145,7 +145,8 @@ namespace {
                    const ::Opm::Schedule& sched,
                    const std::size_t      lookup_step)
     {
-        const auto& wsd = rspec.wellSegmentDimensions();
+	const auto& wd = rspec.wellDimensions();
+	const auto& wsd = rspec.wellSegmentDimensions();
 
         const auto& sched_wells = sched.getWells(lookup_step);
 
@@ -160,7 +161,7 @@ namespace {
         const auto nsegmx = wsd.maxSegmentsPerWell();
         const auto nlbrmx = wsd.maxLateralBranchesPerWell();
         const auto nisegz = 22;  // Number of entries per segment in ISEG.
-        const auto nrsegz = 140; // Number of entries per segment in RSEG array.
+        const auto nrsegz = 146; // Number of entries per segment in RSEG array.  (Eclipse v.2017)
         const auto nilbrz = 10;  // Number of entries per segment in ILBR array.
 
         return {
@@ -192,6 +193,18 @@ namespace {
             static_cast<int>(nplmix),
         };
     }
+    
+    
+    Opm::RestartIO::InteHEAD::Group
+    getNoGroups(const ::Opm::Schedule& sched,
+		const std::size_t      step)
+    {
+        const auto ngroups = sched.numGroups(step)-1;
+
+        return {
+	    ngroups
+	};
+    }
 } // Anonymous
 
 // #####################################################################
@@ -205,8 +218,8 @@ createInteHead(const EclipseState& es,
                const Schedule&     sched,
                const double        simTime,
                const int           num_solver_steps,
-               const int           lookup_step,
+               const int           lookup_step
-               const int           report_step)
+	      )
 {
     const auto& rspec = es.runspec();
     const auto& tdim  = es.getTableManager();
@@ -229,11 +242,12 @@ createInteHead(const EclipseState& es,
              // n{isx}aaqz: number of data elements per aquifer in {ISX}AAQ
              // n{isa}caqz: number of data elements per aquifer connection in {ISA}CAQ
         .params_NAAQZ       (1, 18, 24, 10, 7, 2, 4)
-        .stepParam          (num_solver_steps, report_step)
+        .stepParam          (num_solver_steps, lookup_step)
         .tuningParam        (getTuningPars(sched.getTuning(), lookup_step))
         .wellSegDimensions  (getWellSegDims(rspec, sched, lookup_step))
         .regionDimensions   (getRegDims(tdim, rdim))
-        .variousParam       (2014, 100) // Output should be compatible with Eclipse 100, 2014 version.
+	.ngroups(getNoGroups(sched, lookup_step))
+        .variousParam       (201702, 100)  // Output should be compatible with Eclipse 100, 2017.02 version.
         ;
 
     return ih.data();

src/opm/output/eclipse/DoubHEAD.cpp

@@ -375,7 +375,8 @@ namespace {
 // ---------------------------------------------------------------------
 
 Opm::RestartIO::DoubHEAD::DoubHEAD()
-    : data_(Index::NUMBER_OF_ITEMS, -1.0e20)
+    : data_(Index::NUMBER_OF_ITEMS, 0.0)
+    //: data_(Index::NUMBER_OF_ITEMS, -1.0e20)
 {
     // Numbers below have unknown usage, values have been determined by
     // experiments to be constant across a range of reference cases.

src/opm/output/eclipse/EclipseIO.cpp

@@ -20,8 +20,6 @@
   You should have received a copy of the GNU General Public License
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
-#include <unordered_map>
-
 #include "config.h"
 
 #include <opm/output/eclipse/EclipseIO.hpp>
@@ -38,12 +36,15 @@
 #include <opm/parser/eclipse/EclipseState/Schedule/ScheduleEnums.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
 #include <opm/parser/eclipse/Utility/Functional.hpp>
-#include <opm/output/eclipse/Summary.hpp>
+
-#include <opm/output/eclipse/Tables.hpp>
 #include <opm/output/eclipse/RestartIO.hpp>
+#include <opm/output/eclipse/Summary.hpp>
+#include <opm/output/eclipse/SummaryState.hpp>
+#include <opm/output/eclipse/Tables.hpp>
 
 #include <cstdlib>
 #include <memory>     // unique_ptr
+#include <unordered_map>
 #include <utility>    // move
 
 #include <ert/ecl/EclKW.hpp>
@@ -422,7 +423,7 @@ void EclipseIO::writeTimeStep(int report_step,
                               const std::map<std::string, double>& single_summary_values,
                               const std::map<std::string, std::vector<double> >& region_summary_values,
                               const std::map<std::pair<std::string, int>, double>& block_summary_values,
-                              bool write_double)
+                              const bool write_double)
  {
 
     if( !this->impl->output_enabled )
@@ -454,7 +455,6 @@ void EclipseIO::writeTimeStep(int report_step,
         this->impl->summary.write();
     }
 
-
     /*
       Current implementation will not write restart files for substep,
       but there is an unsupported option to the RPTSCHED keyword which
@@ -468,7 +468,8 @@ void EclipseIO::writeTimeStep(int report_step,
                                                  report_step,
                                                  ioConfig.getFMTOUT() );
 
-        RestartIO::save( filename , report_step, secs_elapsed, value, es , grid , schedule, write_double);
+        RestartIO::save(filename, report_step, secs_elapsed, value, es, grid, schedule,
+                        this->impl->summary.get_restart_vectors(), write_double);
     }
 
 

src/opm/output/eclipse/InteHEAD.cpp

@@ -37,7 +37,7 @@ enum index : std::vector<int>::size_type {
   ih_015	=	15	,		//	0	0
   NWELLS	=	VI::intehead::NWELLS	,		//	NWELLS	39		NWELL = number of wells
   NCWMAX	=	VI::intehead::NCWMAX	,		//	NCWMAX	108	Weldims item2	NCWMAX = maximum number of completions per well
-  ih_018	=	18	,		//	NGRP?	0	Number of actual groups
+  NGRP	=	18	,		//	NGRP?	0	Number of actual groups
   NWGMAX	=	VI::intehead::NWGMAX	,		//	NWGMAX	0	maximum of weldims item3 or item4	NWGMAX = maximum number of wells in any well group
   NGMAXZ	=	VI::intehead::NGMAXZ	,		//	NGMAXZ	0	weldims item3 + 1	NGMAXZ = maximum number of groups in field
   ih_021	=	21	,		//	0	0
@@ -68,8 +68,8 @@ enum index : std::vector<int>::size_type {
   NSCAQZ	=	VI::intehead::NSCAQZ	,		//	2	0
   NACAQZ	=	VI::intehead::NACAQZ	,		//	4	0
   ih_048	=	48	,		//	0	0
-  ih_049	=	49	,		//	0	0
+  ih_049	=	49	,		//	1	// has been determined by testing
-  ih_050	=	50	,		//	0	0
+  ih_050	=	50	,		//	1	// has been determined by testing
   ih_051	=	51	,		//	0	0
   ih_052	=	52	,		//	0	0
   ih_053	=	53	,		//	0	0
@@ -198,7 +198,7 @@ enum index : std::vector<int>::size_type {
   NSEGMX	=	VI::intehead::NSEGMX	,		//	NSEGMX	0	Item 2 in WSEGDIMS keyword (runspec section)	NSEGMX = maximum number of segments per well
   NLBRMX	=	VI::intehead::NLBRMX	,		//	NLBRMX	0	Item 3 in WSEGDIMS keyword (runspec section)	NLBRMX = maximum number of lateral branches per well
   NISEGZ	=	VI::intehead::NISEGZ	,		//	22	0	22	NISEGZ = number of entries per segment in ISEG array
-  NRSEGZ	=	VI::intehead::NRSEGZ	,		//	140	0	140	NRSEGZ = number of entries per segment in RSEG array
+  NRSEGZ	=	VI::intehead::NRSEGZ	,		//	146	0	140	NRSEGZ = number of entries per segment in RSEG array
   NILBRZ	=	VI::intehead::NILBRZ	,		//	10		10	NILBRZ = number of entries per segment in ILBR array
   RSTSIZE	=	181	,		//	0
   ih_182	=	182	,		//	0
@@ -593,12 +593,10 @@ params_NAAQZ(const int ncamax,
 
 Opm::RestartIO::InteHEAD&
 Opm::RestartIO::InteHEAD::
-stepParam(const int num_solver_steps, const int report_step)
+stepParam(const int num_solver_steps, const int sim_step)
 {
     this -> data_[NUM_SOLVER_STEPS] = num_solver_steps;
-    this -> data_[REPORT_STEP] = report_step;
+    this -> data_[REPORT_STEP]      = sim_step + 1;
-
-
 
     return *this;
 }
@@ -622,14 +620,27 @@ Opm::RestartIO::InteHEAD::variousParam(const int version,
 {
     this->data_[VERSION] = version;
     this->data_[IPROG]   = iprog;
-    // ih_076: Usage unknown, experiments fails (zero determinant in well message) with too low numbers. 5 is highest observed across reference cases.
+
-    this->data_[ih_076]  = 5;
+    // ih_049: Usage unknown, value fixed across reference cases
+    this->data_[ih_049]  = 1;
+
+    // ih_050: Usage unknown, value fixed across reference cases
+    this->data_[ih_050]  = 1;
+
+    // ih_076: Usage unknown, experiments fails (zero determinant
+    //         in well message) with too low numbers.
+    //         5 is highest observed across reference cases.
+    this->data_[ih_076] = 5;
+
     // ih_101: Usage unknown, value fixed across reference cases.
-    this->data_[ih_101]  = 1;
+    this->data_[ih_101] = 1;
-    // ih_103: Usage unknown, value not fixed across reference cases, experiments generate warning with 0 but not with 1.
+
-    this->data_[ih_103]  = 1;
+    // ih_103: Usage unknown, value not fixed across reference cases,
+    //         experiments generate warning with 0 but not with 1.
+    this->data_[ih_103] = 1;
+
     // ih_200: Usage unknown, value fixed across reference cases.
-    this->data_[ih_200]  = 1;
+    this->data_[ih_200] = 1;
 
     return *this;
 }
@@ -657,6 +668,16 @@ Opm::RestartIO::InteHEAD::regionDimensions(const RegDims& rdim)
     return *this;
 }
 
+Opm::RestartIO::InteHEAD&
+Opm::RestartIO::InteHEAD::
+ngroups(const Group& gr)
+{
+    this -> data_[NGRP] = gr.ngroups;
+
+    return *this;
+}
+
+
 // =====================================================================
 // Free functions (calendar/time utilities)
 // =====================================================================

src/opm/output/eclipse/LoadRestart.cpp

@@ -0,0 +1,648 @@
+/*
+  Copyright (c) 2016 Statoil ASA
+  Copyright (c) 2013-2015 Andreas Lauser
+  Copyright (c) 2013 SINTEF ICT, Applied Mathematics.
+  Copyright (c) 2013 Uni Research AS
+  Copyright (c) 2015 IRIS AS
+
+  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 <opm/output/eclipse/RestartIO.hpp>
+
+#include <opm/output/eclipse/RestartValue.hpp>
+
+#include <opm/output/eclipse/VectorItems/connection.hpp>
+#include <opm/output/eclipse/VectorItems/intehead.hpp>
+#include <opm/output/eclipse/VectorItems/well.hpp>
+
+#include <opm/output/eclipse/libECLRestart.hpp>
+
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+#include <opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp>
+#include <opm/parser/eclipse/EclipseState/Runspec.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
+
+#include <algorithm>
+#include <cmath>
+#include <cstddef>
+#include <exception>
+#include <stdexcept>
+#include <string>
+#include <unordered_map>
+#include <utility>
+#include <vector>
+
+namespace VI = ::Opm::RestartIO::Helpers::VectorItems;
+
+class RestartFileView
+{
+public:
+    explicit RestartFileView(const std::string& filename,
+                             const int          report_step);
+
+    ~RestartFileView() = default;
+
+    RestartFileView(const RestartFileView& rhs) = delete;
+    RestartFileView(RestartFileView&& rhs);
+
+    RestartFileView& operator=(const RestartFileView& rhs) = delete;
+    RestartFileView& operator=(RestartFileView&& rhs);
+
+    std::size_t simStep() const
+    {
+        return this->sim_step_;
+    }
+
+    operator const Opm::RestartIO::ecl_file_view_type*() const
+    {
+        return this->step_view_;
+    }
+
+    const Opm::RestartIO::ecl_kw_type* getKeyword(const char* kw) const
+    {
+        namespace Load = Opm::RestartIO;
+
+        // Main grid only.  Does not handle/support LGR.
+        return Load::ecl_file_view_has_kw       (*this, kw)
+            ?  Load::ecl_file_view_iget_named_kw(*this, kw, 0)
+            :  nullptr;
+    }
+
+private:
+    using RstFile = Opm::RestartIO::ert_unique_ptr<
+        Opm::RestartIO::ecl_file_type,
+        Opm::RestartIO::ecl_file_close>;
+
+    std::size_t                         sim_step_;
+    RstFile                             rst_file_;
+    Opm::RestartIO::ecl_file_view_type* step_view_ = nullptr;
+
+    operator Opm::RestartIO::ecl_file_type*()
+    {
+        return this->rst_file_.get();
+    }
+
+    operator const Opm::RestartIO::ecl_file_type*() const
+    {
+        return this->rst_file_.get();
+    }
+};
+
+RestartFileView::RestartFileView(const std::string& filename,
+                                 const int          report_step)
+    : sim_step_(std::max(report_step - 1, 0))
+    , rst_file_(Opm::RestartIO::ecl_file_open(filename.c_str(), 0))
+{
+    namespace Load = Opm::RestartIO;
+
+    if (this->rst_file_ == nullptr) {
+        throw std::invalid_argument {
+            "Unable to open Restart File '" + filename
+            + "' at Report Step " + std::to_string(report_step)
+        };
+    }
+
+    this->step_view_ =
+        (Load::EclFiletype(filename) == Load::ECL_UNIFIED_RESTART_FILE)
+        ? Load::ecl_file_get_restart_view(*this, -1, report_step, -1, -1)
+        : Load::ecl_file_get_global_view (*this);    // Separate
+
+    if (this->step_view_ == nullptr) {
+        throw std::runtime_error {
+            "Unable to acquire restart information for report step "
+            + std::to_string(report_step)
+        };
+    }
+}
+
+RestartFileView::RestartFileView(RestartFileView&& rhs)
+    : sim_step_ (rhs.sim_step_)            // Scalar (size_t)
+    , rst_file_ (std::move(rhs.rst_file_))
+    , step_view_(rhs.step_view_)           // Pointer
+{}
+
+RestartFileView& RestartFileView::operator=(RestartFileView&& rhs)
+{
+    this->sim_step_  = rhs.sim_step_;            // Scalar (size_t)
+    this->rst_file_  = std::move(rhs.rst_file_);
+    this->step_view_ = rhs.step_view_;           // Pointer copy
+
+    return *this;
+}
+
+namespace {
+    void throwIfMissingRequired(const Opm::RestartKey& rst_key)
+    {
+        if (rst_key.required) {
+            throw std::runtime_error {
+                "Requisite restart vector '"
+                + rst_key.key +
+                "' is not available in restart file"
+            };
+        }
+    }
+
+    std::vector<double>
+    double_vector(const ::Opm::RestartIO::ecl_kw_type* ecl_kw)
+    {
+        namespace Load = ::Opm::RestartIO;
+
+        const auto size = static_cast<std::size_t>(
+            Load::ecl_kw_get_size(ecl_kw));
+
+        if (Load::ecl_type_get_type(Load::ecl_kw_get_data_type(ecl_kw)) == Load::ECL_DOUBLE_TYPE) {
+            const double* ecl_data = Load::ecl_kw_get_type_ptr<double>(ecl_kw, Load::ECL_DOUBLE_TYPE);
+
+            return { ecl_data , ecl_data + size };
+        }
+        else {
+            const float* ecl_data = Load::ecl_kw_get_type_ptr<float>(ecl_kw, Load::ECL_FLOAT_TYPE);
+
+            return { ecl_data , ecl_data + size };
+        }
+    }
+
+    Opm::data::Solution
+    restoreSOLUTION(const RestartFileView&              rst_view,
+                    const std::vector<Opm::RestartKey>& solution_keys,
+                    const int                           numcells)
+    {
+        Opm::data::Solution sol(/* init_si = */ false);
+
+        for (const auto& value : solution_keys) {
+            const auto& vector = value.key;
+            const auto* kw     = rst_view.getKeyword(vector.c_str());
+
+            if (kw == nullptr) {
+                throwIfMissingRequired(value);
+
+                // Requested vector not available, but caller does not
+                // actually require the vector for restart purposes.
+                // Skip this.
+                continue;
+            }
+
+            if (Opm::RestartIO::ecl_kw_get_size(kw) != numcells) {
+                throw std::runtime_error {
+                    "Restart file: Could not restore "
+                    + std::string(Opm::RestartIO::ecl_kw_get_header(kw))
+                    + ", mismatched number of cells"
+                };
+            }
+
+            sol.insert(vector, value.dim, double_vector(kw),
+                       Opm::data::TargetType::RESTART_SOLUTION);
+        }
+
+        return sol;
+    }
+
+    void restoreExtra(const RestartFileView&              rst_view,
+                      const std::vector<Opm::RestartKey>& extra_keys,
+                      const Opm::UnitSystem&              usys,
+                      Opm::RestartValue&                  rst_value)
+    {
+        for (const auto& extra : extra_keys) {
+            const auto& vector = extra.key;
+            const auto* kw     = rst_view.getKeyword(vector.c_str());
+
+            if (kw == nullptr) {
+                throwIfMissingRequired(extra);
+
+                // Requested vector not available, but caller does not
+                // actually require the vector for restart purposes.
+                // Skip this.
+                continue;
+            }
+
+            // Requisite vector available in result set.  Recover data.
+            rst_value.addExtra(vector, extra.dim, double_vector(kw));
+        }
+
+        for (auto& extra_value : rst_value.extra) {
+            const auto& restart_key = extra_value.first;
+            auto&       data        = extra_value.second;
+
+            usys.to_si(restart_key.dim, data);
+        }
+    }
+
+    void checkWellVectorSizes(const ::Opm::RestartIO::ecl_kw_type*      opm_xwel,
+                              const ::Opm::RestartIO::ecl_kw_type*      opm_iwel,
+                              const int                                 sim_step,
+                              const std::vector<Opm::data::Rates::opt>& phases,
+                              const std::vector<const Opm::Well*>&      sched_wells)
+    {
+        const auto expected_xwel_size =
+            std::accumulate(sched_wells.begin(), sched_wells.end(),
+                            std::size_t(0),
+                [&phases, sim_step](const std::size_t acc, const Opm::Well* w)
+                -> std::size_t
+            {
+                return acc
+                    + 2 + phases.size()
+                    + (w->getConnections(sim_step).size()
+                        * (phases.size() + Opm::data::Connection::restart_size));
+            });
+
+        if (static_cast<std::size_t>(::Opm::RestartIO::ecl_kw_get_size(opm_xwel)) != expected_xwel_size)
+        {
+            throw std::runtime_error {
+                "Mismatch between OPM_XWEL and deck; "
+                "OPM_XWEL size was " + std::to_string(::Opm::RestartIO::ecl_kw_get_size(opm_xwel)) +
+                ", expected " + std::to_string(expected_xwel_size)
+            };
+        }
+
+        if (::Opm::RestartIO::ecl_kw_get_size(opm_iwel) != int(sched_wells.size())) {
+            throw std::runtime_error {
+                "Mismatch between OPM_IWEL and deck; "
+                "OPM_IWEL size was " + std::to_string(::Opm::RestartIO::ecl_kw_get_size(opm_iwel)) +
+                ", expected " + std::to_string(sched_wells.size())
+            };
+        }
+    }
+
+    Opm::data::Wells
+    restore_wells_opm(const RestartFileView&     rst_view,
+                      const ::Opm::EclipseState& es,
+                      const ::Opm::EclipseGrid&  grid,
+                      const ::Opm::Schedule&     schedule)
+    {
+        namespace Load = ::Opm::RestartIO;
+
+        const auto* opm_iwel = rst_view.getKeyword("OPM_IWEL");
+        const auto* opm_xwel = rst_view.getKeyword("OPM_XWEL");
+
+        if ((opm_xwel == nullptr) || (opm_iwel == nullptr)) {
+            return {};
+        }
+
+        using rt = Opm::data::Rates::opt;
+
+        const auto& sched_wells = schedule.getWells(rst_view.simStep());
+        std::vector<rt> phases;
+        {
+            const auto& phase = es.runspec().phases();
+
+            if (phase.active(Opm::Phase::WATER)) { phases.push_back(rt::wat); }
+            if (phase.active(Opm::Phase::OIL))   { phases.push_back(rt::oil); }
+            if (phase.active(Opm::Phase::GAS))   { phases.push_back(rt::gas); }
+        }
+
+        checkWellVectorSizes(opm_xwel, opm_iwel,
+                             rst_view.simStep(),
+                             phases, sched_wells);
+
+        Opm::data::Wells wells;
+        const auto* opm_xwel_data = Load::ecl_kw_get_type_ptr<double>(opm_xwel, Load::ECL_DOUBLE_TYPE);
+        const auto* opm_iwel_data = Load::ecl_kw_get_type_ptr<int>(opm_iwel, Load::ECL_INT_TYPE);
+
+        for (const auto* sched_well : sched_wells) {
+            auto& well = wells[ sched_well->name() ];
+
+            well.bhp         = *opm_xwel_data++;
+            well.temperature = *opm_xwel_data++;
+            well.control     = *opm_iwel_data++;
+
+            for (const auto& phase : phases) {
+                well.rates.set(phase, *opm_xwel_data++);
+            }
+
+            for (const auto& sc : sched_well->getConnections(rst_view.simStep())) {
+                const auto i = sc.getI(), j = sc.getJ(), k = sc.getK();
+
+                if (!grid.cellActive(i, j, k) || sc.state() == Opm::WellCompletion::SHUT) {
+                    opm_xwel_data += Opm::data::Connection::restart_size + phases.size();
+                    continue;
+                }
+
+                const auto active_index = grid.activeIndex(i, j, k);
+
+                well.connections.emplace_back();
+                auto& connection = well.connections.back();
+
+                connection.index          = active_index;
+                connection.pressure       = *opm_xwel_data++;
+                connection.reservoir_rate = *opm_xwel_data++;
+
+                for (const auto& phase : phases) {
+                    connection.rates.set(phase, *opm_xwel_data++);
+                }
+            }
+        }
+
+        return wells;
+    }
+
+    template <typename T>
+    const T* getPtr(const ::Opm::RestartIO::ecl_kw_type* kw)
+    {
+        return (kw == nullptr) ? nullptr
+            : static_cast<const T*>(ecl_kw_iget_ptr(kw /* <- ADL */, 0));
+    }
+
+    int getInteHeadElem(const ::Opm::RestartIO::ecl_kw_type* intehead,
+                        const std::vector<int>::size_type    i)
+    {
+        return getPtr<int>(intehead)[i];
+    }
+
+    struct WellArrayDim
+    {
+        explicit WellArrayDim(const ::Opm::RestartIO::ecl_kw_type* intehead);
+
+        std::size_t maxConnPerWell;
+        std::size_t numIWelElem;
+        std::size_t numXWelElem;
+        std::size_t numIConElm;
+        std::size_t numXConElm;
+    };
+
+    WellArrayDim::WellArrayDim(const ::Opm::RestartIO::ecl_kw_type* intehead)
+        : maxConnPerWell(getInteHeadElem(intehead, VI::intehead::NCWMAX))
+        , numIWelElem   (getInteHeadElem(intehead, VI::intehead::NIWELZ))
+        , numXWelElem   (getInteHeadElem(intehead, VI::intehead::NXWELZ))
+        , numIConElm    (getInteHeadElem(intehead, VI::intehead::NICONZ))
+        , numXConElm    (getInteHeadElem(intehead, VI::intehead::NXCONZ))
+    {}
+
+    template <typename T>
+    boost::iterator_range<const T*>
+    getDataWindow(const T*          arr,
+                  const std::size_t windowSize,
+                  const std::size_t well,
+                  const std::size_t conn           = 0,
+                  const std::size_t maxConnPerWell = 1)
+    {
+        const auto  off   = windowSize * (conn + maxConnPerWell*well);
+        const auto* begin = arr   + off;
+        const auto* end   = begin + windowSize;
+
+        return { begin, end };
+    }
+
+    boost::iterator_range<const int*>
+    getIWelWindow(const int*          iwel,
+                  const WellArrayDim& wdim,
+                  const std::size_t   well)
+    {
+        return getDataWindow(iwel, wdim.numIWelElem, well);
+    }
+
+    boost::iterator_range<const double*>
+    getXWelWindow(const double*       xwel,
+                  const WellArrayDim& wdim,
+                  const std::size_t   well)
+    {
+        return getDataWindow(xwel, wdim.numXWelElem, well);
+    }
+
+    boost::iterator_range<const int*>
+    getIConWindow(const int*          icon,
+                  const WellArrayDim& wdim,
+                  const std::size_t   well,
+                  const std::size_t   conn)
+    {
+        return getDataWindow(icon, wdim.numIConElm, well,
+                             conn, wdim.maxConnPerWell);
+    }
+
+    boost::iterator_range<const double*>
+    getXConWindow(const double*       xcon,
+                  const WellArrayDim& wdim,
+                  const std::size_t   well,
+                  const std::size_t   conn)
+    {
+        return getDataWindow(xcon, wdim.numXConElm, well,
+                             conn, wdim.maxConnPerWell);
+    }
+
+    std::unordered_map<std::size_t, boost::iterator_range<const double*>::size_type>
+    seqID_to_resID(const WellArrayDim& wdim,
+                   const std::size_t   wellID,
+                   const std::size_t   nConn,
+                   const int*          icon_full)
+    {
+        using SizeT    = boost::iterator_range<const double*>::size_type;
+        auto  seqToRes = std::unordered_map<std::size_t, SizeT>{};
+
+        for (auto connID = 0*nConn; connID < nConn; ++connID) {
+            const auto icon =
+                getIConWindow(icon_full, wdim, wellID, connID);
+
+            seqToRes.emplace(icon[VI::IConn::index::SeqIndex] - 1, connID);
+        }
+
+        return seqToRes;
+    }
+
+    void restoreConnRates(const Opm::Well&        well,
+                          const std::size_t       wellID,
+                          const std::size_t       sim_step,
+                          const Opm::EclipseGrid& grid,
+                          const WellArrayDim&     wdim,
+                          const Opm::UnitSystem&  usys,
+                          const Opm::Phases&      phases,
+                          const int*              iwel_full,
+                          const int*              icon_full,
+                          const double*           xcon_full,
+                          Opm::data::Well&        xw)
+    {
+        using M = ::Opm::UnitSystem::measure;
+
+        const auto iwel  = getIWelWindow(iwel_full, wdim, wellID);
+        const auto nConn = static_cast<std::size_t>(
+            iwel[VI::IWell::index::NConn]);
+
+        xw.connections.resize(nConn, Opm::data::Connection{});
+
+        if ((icon_full == nullptr) || (xcon_full == nullptr)) {
+            // Result set does not provide certain pieces of
+            // information which are needed to reconstruct
+            // connection flow rates.  Nothing to do here.
+            return;
+        }
+
+        const auto oil = phases.active(Opm::Phase::OIL);
+        const auto gas = phases.active(Opm::Phase::GAS);
+        const auto wat = phases.active(Opm::Phase::WATER);
+
+        const auto conns      = well.getActiveConnections(sim_step, grid);
+        const auto seq_to_res =
+            seqID_to_resID(wdim, wellID, nConn, icon_full);
+
+        auto linConnID = std::size_t{0};
+        for (const auto& conn : conns) {
+            const auto connID = seq_to_res.at(conn.getSeqIndex());
+            const auto xcon   =
+                getXConWindow(xcon_full, wdim, wellID, connID);
+
+            auto& xc = xw.connections[linConnID++];
+
+            if (wat) {
+                xc.rates.set(Opm::data::Rates::opt::wat,
+                             - usys.to_si(M::liquid_surface_rate,
+                                          xcon[VI::XConn::index::WaterRate]));
+            }
+
+            if (oil) {
+                xc.rates.set(Opm::data::Rates::opt::oil,
+                             - usys.to_si(M::liquid_surface_rate,
+                                          xcon[VI::XConn::index::OilRate]));
+            }
+
+            if (gas) {
+                xc.rates.set(Opm::data::Rates::opt::gas,
+                             - usys.to_si(M::gas_surface_rate,
+                                          xcon[VI::XConn::index::GasRate]));
+            }
+        }
+    }
+
+    Opm::data::Well
+    restore_well(const Opm::Well&        well,
+                 const std::size_t       wellID,
+                 const std::size_t       sim_step,
+                 const Opm::EclipseGrid& grid,
+                 const WellArrayDim&     wdim,
+                 const Opm::UnitSystem&  usys,
+                 const Opm::Phases&      phases,
+                 const int*              iwel_full,
+                 const double*           xwel_full,
+                 const int*              icon_full,
+                 const double*           xcon_full)
+    {
+        if ((iwel_full == nullptr) || (xwel_full == nullptr)) {
+            // Result set does not provide well information.
+            // No wells?  In any case, nothing to do here.
+            return {};
+        }
+
+        using M = ::Opm::UnitSystem::measure;
+
+        const auto xwel = getXWelWindow(xwel_full, wdim, wellID);
+
+        const auto oil = phases.active(Opm::Phase::OIL);
+        const auto gas = phases.active(Opm::Phase::GAS);
+        const auto wat = phases.active(Opm::Phase::WATER);
+
+        auto xw = ::Opm::data::Well{};
+
+        // 1) Restore well rates (xw.rates)
+        if (wat) {
+            xw.rates.set(Opm::data::Rates::opt::wat,
+                         - usys.to_si(M::liquid_surface_rate,
+                                      xwel[VI::XWell::index::WatPrRate]));
+        }
+
+        if (oil) {
+            xw.rates.set(Opm::data::Rates::opt::oil,
+                         - usys.to_si(M::liquid_surface_rate,
+                                      xwel[VI::XWell::index::OilPrRate]));
+        }
+
+        if (gas) {
+            xw.rates.set(Opm::data::Rates::opt::gas,
+                         - usys.to_si(M::gas_surface_rate,
+                                      xwel[VI::XWell::index::GasPrRate]));
+        }
+
+        // 2) Restore other well quantities (really only xw.bhp)
+        xw.bhp = usys.to_si(M::pressure, xwel[VI::XWell::index::FlowBHP]);
+        xw.thp = xw.temperature = 0.0;
+
+        // 3) Restore connection flow rates (xw.connections[i].rates)
+        restoreConnRates(well, wellID, sim_step, grid, wdim, usys, phases,
+                         iwel_full, icon_full, xcon_full, xw);
+
+        return xw;
+    }
+
+    Opm::data::Wells
+    restore_wells_ecl(const RestartFileView&     rst_view,
+                      const ::Opm::EclipseState& es,
+                      const ::Opm::EclipseGrid&  grid,
+                      const ::Opm::Schedule&     schedule)
+    {
+        auto soln = ::Opm::data::Wells{};
+
+        const auto* intehead = rst_view.getKeyword("INTEHEAD");
+
+        if (intehead == nullptr) {
+            // Result set does not provide indexing information.
+            // Can't do anything here.
+            return soln;
+        }
+
+        const auto  wdim   = WellArrayDim{ intehead };
+        const auto& units  = es.getUnits();
+        const auto& phases = es.runspec().phases();
+
+        const auto* iwel_full = getPtr<int>   (rst_view.getKeyword("IWEL"));
+        const auto* xwel_full = getPtr<double>(rst_view.getKeyword("XWEL"));
+        const auto* icon_full = getPtr<int>   (rst_view.getKeyword("ICON"));
+        const auto* xcon_full = getPtr<double>(rst_view.getKeyword("XCON"));
+
+        const auto  sim_step = rst_view.simStep();
+        const auto& wells    = schedule.getWells(sim_step);
+        for (auto nWells = wells.size(), wellID = 0*nWells;
+             wellID < nWells; ++wellID)
+        {
+            const auto* well = wells[wellID];
+
+            soln[well->name()] =
+                restore_well(*well, wellID, sim_step, grid, wdim, units, phases,
+                             iwel_full, xwel_full, icon_full, xcon_full);
+        }
+
+        return soln;
+    }
+} // Anonymous namespace
+
+namespace Opm { namespace RestartIO  {
+
+    RestartValue
+    load(const std::string&             filename,
+         int                            report_step,
+         const std::vector<RestartKey>& solution_keys,
+         const EclipseState&            es,
+         const EclipseGrid&             grid,
+         const Schedule&                schedule,
+         const std::vector<RestartKey>& extra_keys)
+    {
+        const auto rst_view = RestartFileView{ filename, report_step };
+
+        auto xr = restoreSOLUTION(rst_view, solution_keys,
+                                  grid.getNumActive());
+
+        xr.convertToSI(es.getUnits());
+
+        auto xw = Opm::RestartIO::ecl_file_view_has_kw(rst_view, "OPM_XWEL")
+            ? restore_wells_opm(rst_view, es, grid, schedule)
+            : restore_wells_ecl(rst_view, es, grid, schedule);
+
+        auto rst_value = RestartValue{ std::move(xr), std::move(xw) };
+
+        if (! extra_keys.empty()) {
+            restoreExtra(rst_view, extra_keys, es.getUnits(), rst_value);
+        }
+
+        return rst_value;
+    }
+}} // Opm::RestartIO

src/opm/output/eclipse/LogiHEAD.cpp

@@ -148,18 +148,18 @@ variousParam(const bool e300_radial, const bool e100_radial, const int nswlmx)
     this -> data_[lh_001] = true;
     this -> data_[lh_003] = e300_radial;
     this -> data_[lh_004] = e100_radial;
-    this -> data_[lh_016] = true;
+    //this -> data_[lh_016] = true;
-    this -> data_[lh_018] = true;
+    //this -> data_[lh_018] = true;
-    this -> data_[lh_031] = true;
+    //this -> data_[lh_031] = true;
-    this -> data_[lh_044] = true;
+    //this -> data_[lh_044] = true;
     this -> data_[lh_075] = nswlmx >= 1; // True if MS Wells exist.
-    this -> data_[lh_076] = true;
+    //this -> data_[lh_076] = true;
-    this -> data_[lh_087] = true;
+    //this -> data_[lh_087] = true;
-    this -> data_[lh_099] = true;
+    //this -> data_[lh_099] = true;
-    this -> data_[lh_113] = true;
+    //this -> data_[lh_113] = true;
-    this -> data_[lh_114] = true;
+    //this -> data_[lh_114] = true;
-    this -> data_[lh_115] = true;
+    //this -> data_[lh_115] = true;
-    this -> data_[lh_117] = true;
+    //this -> data_[lh_117] = true;
 
     return *this;
 }

src/opm/output/eclipse/RestartValue.cpp

@@ -81,5 +81,25 @@ namespace Opm {
         this->addExtra(key, UnitSystem::measure::identity, std::move(data));
     }
 
+    void RestartValue::convertFromSI(const UnitSystem& units) {
+        this->solution.convertFromSI(units);
+        for (auto & extra_value : this->extra) {
+            const auto& restart_key = extra_value.first;
+            auto & data = extra_value.second;
+
+            units.from_si(restart_key.dim, data);
+        }
+    }
+
+    void RestartValue::convertToSI(const UnitSystem& units) {
+        this->solution.convertToSI(units);
+        for (auto & extra_value : this->extra) {
+            const auto& restart_key = extra_value.first;
+            auto & data = extra_value.second;
+
+            units.to_si(restart_key.dim, data);
+        }
+    }
+
 
 }

src/opm/parser/eclipse/EclipseState/Runspec.cpp

@@ -90,6 +90,7 @@ Welldims::Welldims(const Deck& deck)
         //
         // i.e., the maximum of item 1 and item 4 here.
         this->nGMax = wd.getItem("MAXGROUPS").get<int>(0);
+	this->nWMax = wd.getItem("MAXWELLS").get<int>(0);
     }
 }
 

src/opm/parser/eclipse/EclipseState/Schedule/Connection.cpp

@@ -44,7 +44,12 @@ namespace Opm {
                            double Kh,
                            double rw,
                            const int satTableId,
-                           const WellCompletion::DirectionEnum direction)
+                           const WellCompletion::DirectionEnum direction,
+			   const std::size_t seqIndex,
+			   const double segDistStart,
+			   const double segDistEnd,
+			   const bool defaultSatTabId
+			  )
         : direction(direction),
           center_depth(depth),
           open_state(stateArg),
@@ -53,9 +58,22 @@ namespace Opm {
           m_CF(CF),
           m_Kh(Kh),
           m_rw(rw),
-          ijk({i,j,k})
+          ijk({i,j,k}),
+          m_seqIndex(seqIndex),
+          m_segDistStart(segDistStart),
+          m_segDistEnd(segDistEnd),
+          m_defaultSatTabId(defaultSatTabId)
     {}
 
+    /*bool Connection::sameCoordinate(const Connection& other) const {
+        if ((m_i == other.m_i) &&
+            (m_j == other.m_j) &&
+            (m_k == other.m_k))
+            return true;
+        else
+            return false;
+    }*/
+
     bool Connection::sameCoordinate(const int i, const int j, const int k) const {
         if ((ijk[0] == i) && (ijk[1] == j) && (ijk[2] == k)) {
             return true;
@@ -81,11 +99,48 @@ namespace Opm {
     bool Connection::attachedToSegment() const {
         return (segment_number > 0);
     }
+    
+    const std::size_t& Connection::getSeqIndex() const {
+        return m_seqIndex;
+    }
+    
+    const bool& Connection::getDefaultSatTabId() const {
+        return m_defaultSatTabId;
+    }
+    
+    const std::size_t& Connection::getCompSegSeqIndex() const {
+        return m_compSeg_seqIndex;
+    }
 
     WellCompletion::DirectionEnum Connection::dir() const {
         return this->direction;
     }
 
+    const double& Connection::getSegDistStart() const {
+        return m_segDistStart;
+    }
+
+    const double& Connection::getSegDistEnd() const {
+        return m_segDistEnd;
+    }
+
+    
+    void Connection::setCompSegSeqIndex(std::size_t index) {
+        m_compSeg_seqIndex = index;
+    }
+    
+    void Connection::setDefaultSatTabId(bool id) {
+        m_defaultSatTabId = id;
+    }
+    
+    void Connection::setSegDistStart(const double& distStart) {
+        m_segDistStart = distStart;
+    }
+
+    void Connection::setSegDistEnd(const double& distEnd) {
+        m_segDistEnd = distEnd;
+    }
+    
     double Connection::depth() const {
         return this->center_depth;
     }
@@ -122,9 +177,10 @@ namespace Opm {
         this->open_state = state;
     }
 
-    void Connection::updateSegment(int segment_number, double center_depth) {
+  void Connection::updateSegment(int segment_number, double center_depth, std::size_t seqIndex) {
         this->segment_number = segment_number;
         this->center_depth = center_depth;
+        this->m_seqIndex = seqIndex;
     }
 
     int Connection::segment() const {
@@ -150,12 +206,11 @@ namespace Opm {
             && this->open_state == rhs.open_state
             && this->direction == rhs.direction
             && this->segment_number == rhs.segment_number
-            && this->center_depth == rhs.center_depth;
+            && this->center_depth == rhs.center_depth
+            && this->m_seqIndex == rhs.m_seqIndex;
     }
 
     bool Connection::operator!=( const Connection& rhs ) const {
         return !( *this == rhs );
     }
 }
-
-

src/opm/parser/eclipse/EclipseState/Schedule/Group.cpp

@@ -50,9 +50,10 @@ namespace Opm {
 
     /*****************************************************************/
 
-    Group::Group(const std::string& name_, const TimeMap& timeMap , size_t creationTimeStep) :
+    Group::Group(const std::string& name_, const size_t& seqIndex_, const TimeMap& timeMap , size_t creationTimeStep) :
         m_creationTimeStep( creationTimeStep ),
         m_name( name_ ),
+        m_seqIndex( seqIndex_),
         m_injection( timeMap ),
         m_production( timeMap ),
         m_wells( timeMap, {} ),
@@ -68,6 +69,9 @@ namespace Opm {
         return m_name;
     }
 
+    const size_t& Group::seqIndex() const {
+        return m_seqIndex;
+    }
 
     bool Group::hasBeenDefined(size_t timeStep) const {
         if (timeStep < m_creationTimeStep)

src/opm/parser/eclipse/EclipseState/Schedule/GroupTree.cpp

@@ -19,13 +19,14 @@
 
 #include <algorithm>
 #include <stdexcept>
+#include <iostream>
 
 #include <opm/parser/eclipse/EclipseState/Schedule/GroupTree.hpp>
 
 namespace Opm {
 
-void GroupTree::update( const std::string& name ) {
+void GroupTree::update( const std::string& name) {
-    this->update( name, "FIELD" );
+    this->update( name, "FIELD");
 }
 
 /*
@@ -35,8 +36,16 @@ void GroupTree::update( const std::string& name ) {
  * this grouptree class is just meta data for the actual group objects (stored
  * and represented elsewhere)
  */
+ 
+const std::map<const std::string , size_t>& GroupTree::nameSeqIndMap() const {
+  return m_nameSeqIndMap;
+}
 
-void GroupTree::update( const std::string& name, const std::string& other_parent ) {
+const std::map<size_t, const std::string >& GroupTree::seqIndNameMap() const {
+  return m_seqIndNameMap;
+}
+
+void GroupTree::update( const std::string& name, const std::string& other_parent) {
     if( name == "FIELD" )
         throw std::invalid_argument( "The FIELD group name is reserved." );
 
@@ -44,12 +53,10 @@ void GroupTree::update( const std::string& name, const std::string& other_parent
         throw std::invalid_argument( "Parent group must have a name." );
 
     auto root = this->find( other_parent );
-
+    if( root == this->groups.end() || root->name != other_parent ) 
-    if( root == this->groups.end() || root->name != other_parent )
         this->groups.insert( root, 1, group { other_parent, "FIELD" } );
 
     auto node = this->find( name );
-
     if( node == this->groups.end() || node->name != name ) {
         this->groups.insert( node, 1, group { name, other_parent } );
         return;
@@ -58,6 +65,30 @@ void GroupTree::update( const std::string& name, const std::string& other_parent
     node->parent = other_parent;
 }
 
+void GroupTree::updateSeqIndex( const std::string& name, const std::string& other_parent) {
+    if( name == "FIELD" )
+        throw std::invalid_argument( "The FIELD group name is reserved." );
+
+    if( other_parent.empty() )
+        throw std::invalid_argument( "Parent group must have a name." );
+
+    // add code to set an index that determine the sequence of the groups
+    // defined in the group tree
+    size_t index = this->m_nameSeqIndMap.size();
+    auto name_itr = this->m_nameSeqIndMap.find(name);
+    if (name_itr == this->m_nameSeqIndMap.end()) {
+	this->m_nameSeqIndMap.insert(std::make_pair(name, index));
+	this->m_seqIndNameMap.insert(std::make_pair(index, name));    	
+	index +=1;
+    }
+    auto parent_itr = this->m_nameSeqIndMap.find(other_parent);
+    if (parent_itr == this->m_nameSeqIndMap.end()) {
+	this->m_nameSeqIndMap.insert(std::make_pair(other_parent, index));
+	this->m_seqIndNameMap.insert(std::make_pair(index, other_parent));
+    }
+}
+
+
 bool GroupTree::exists( const std::string& name ) const {
     return std::binary_search( this->groups.begin(),
                                this->groups.end(),
@@ -78,10 +109,15 @@ std::vector< std::string > GroupTree::children( const std::string& other_parent
         throw std::out_of_range( "Node '" + other_parent + "' does not exist." );
 
     std::vector< std::string > kids;
-    for( const auto& node : this->groups ) {
+    /*    for( const auto& node : this->groups ) {
         if( node.parent != other_parent ) continue;
         kids.push_back( node.name );
     }
+*/ 
+    for( auto it = this->groups.begin(); it != this->groups.end(); it++ ) {
+        if( (*it).parent != other_parent ) continue;
+        kids.push_back( (*it).name );
+    }
 
     return kids;
 }

src/opm/parser/eclipse/EclipseState/Schedule/MSW/Compsegs.cpp

@@ -33,7 +33,7 @@ namespace Opm {
 
 
     Compsegs::Compsegs(int i_in, int j_in, int k_in, int branch_number_in, double distance_start_in, double distance_end_in,
-                       WellCompletion::DirectionEnum dir_in, double center_depth_in, int segment_number_in)
+                       WellCompletion::DirectionEnum dir_in, double center_depth_in, int segment_number_in, size_t seqIndex_in)
     : m_i(i_in),
       m_j(j_in),
       m_k(k_in),
@@ -42,11 +42,12 @@ namespace Opm {
       m_distance_end(distance_end_in),
       m_dir(dir_in),
       center_depth(center_depth_in),
-      segment_number(segment_number_in)
+      segment_number(segment_number_in),
+      m_seqIndex(seqIndex_in)
     {
     }
 
-    std::vector< Compsegs > Compsegs::compsegsFromCOMPSEGSKeyword( const DeckKeyword& compsegsKeyword ) {
+    std::vector< Compsegs > Compsegs::compsegsFromCOMPSEGSKeyword( const DeckKeyword& compsegsKeyword, const EclipseGrid& grid, std::size_t& totNC ) {
 
         // only handle the second record here
         // The first record here only contains the well name
@@ -121,12 +122,19 @@ namespace Opm {
             }
 
             if (!record.getItem<ParserKeywords::COMPSEGS::END_IJK>().hasValue(0)) { // only one compsegs
-                compsegs.emplace_back( I, J, K,
+		
-                                       branch,
+		if (grid.cellActive(I, J, K)) {
-                                       distance_start, distance_end,
+		    std::size_t seqIndex = compsegs.size();
-                                       direction,
+		    totNC = seqIndex;
-                                       center_depth,
+		    compsegs.emplace_back( I, J, K,
-                                       segment_number );
+					  branch,
+					  distance_start, distance_end,
+					  direction,
+					  center_depth,
+					  segment_number,
+					  seqIndex
+				      );
+		}
             } else { // a range is defined. genrate a range of Compsegs
                 throw std::runtime_error("entering COMPSEGS entries with a range is not supported yet!");
             }
@@ -231,15 +239,22 @@ namespace Opm {
     }
 
     void Compsegs::updateConnectionsWithSegment(const std::vector< Compsegs >& compsegs,
+						const EclipseGrid& grid,
                                                 WellConnections& connection_set) {
 
         for( const auto& compseg : compsegs ) {
             const int i = compseg.m_i;
             const int j = compseg.m_j;
             const int k = compseg.m_k;
+	    if (grid.cellActive(i, j, k)) {
+		Connection& connection = connection_set.getFromIJK( i, j, k );
+		connection.updateSegment(compseg.segment_number, compseg.center_depth,compseg.m_seqIndex);
 
-            Connection& connection = connection_set.getFromIJK( i, j, k );
+		//keep connection sequence number from input sequence
-            connection.updateSegment(compseg.segment_number, compseg.center_depth);
+		connection.setCompSegSeqIndex(compseg.m_seqIndex);
+		connection.setSegDistStart(compseg.m_distance_start);
+		connection.setSegDistEnd(compseg.m_distance_end);
+	    }
         }
 
         for (size_t ic = 0; ic < connection_set.size(); ++ic) {

src/opm/parser/eclipse/EclipseState/Schedule/MSW/Compsegs.hpp

@@ -25,12 +25,15 @@
 #include <string>
 
 #include <opm/parser/eclipse/EclipseState/Schedule/ScheduleEnums.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/WellConnections.hpp>
+#include <opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp>
 
 namespace Opm {
 
     class WellConnections;
     class DeckKeyword;
     class WellSegments;
+    class EclipseGrid;
 
     struct Compsegs {
         int m_i;
@@ -48,19 +51,21 @@ namespace Opm {
         // we do not handle thermal length for the moment
         // double m_thermal_length;
         int segment_number;
-
+	std::size_t m_seqIndex;
+		
         Compsegs(int i_in, int j_in, int k_in, int branch_number_in, double distance_start_in, double distance_end_in,
-                 WellCompletion::DirectionEnum dir_in, double center_depth_in, int segment_number_in);
+                 WellCompletion::DirectionEnum dir_in, double center_depth_in, int segment_number_in, std::size_t seqIndex_in);
 
         void calculateCenterDepthWithSegments(const WellSegments& segment_set);
 
-        static std::vector< Compsegs > compsegsFromCOMPSEGSKeyword( const DeckKeyword& compsegsKeyword );
+        static std::vector< Compsegs > compsegsFromCOMPSEGSKeyword( const DeckKeyword& compsegsKeyword, const EclipseGrid& grid, std::size_t& totNC );
 
         // get the segment number information and depth information based on the information from WellSegments
         static void processCOMPSEGS(std::vector< Compsegs >& compsegs, const WellSegments& segment_set );
 
         // update the segment related information for Connections
         static void updateConnectionsWithSegment(const std::vector< Compsegs >& compsegs,
+						 const EclipseGrid& grid,
                                                  WellConnections& connection_set);
 
     };

src/opm/parser/eclipse/EclipseState/Schedule/MSW/updatingConnectionsWithSegments.cpp

@@ -18,19 +18,23 @@
 */
 
 #include <opm/parser/eclipse/EclipseState/Schedule/MSW/updatingConnectionsWithSegments.hpp>
+#include <opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp>
 #include "Compsegs.hpp"
 
 
 namespace Opm {
     WellConnections * newConnectionsWithSegments(const DeckKeyword& compsegs,
                                                  const WellConnections& input_connections,
-                                                 const WellSegments& segment_set)
+                                                 const WellSegments& segment_set,
+						 const EclipseGrid& grid,
+						 std::size_t& totNC
+						)
     {
         WellConnections * new_connection_set = new WellConnections(input_connections);
 
-        std::vector<Compsegs> compsegs_vector = Compsegs::compsegsFromCOMPSEGSKeyword( compsegs );
+        std::vector<Compsegs> compsegs_vector = Compsegs::compsegsFromCOMPSEGSKeyword( compsegs, grid, totNC );
         Compsegs::processCOMPSEGS(compsegs_vector, segment_set);
-        Compsegs::updateConnectionsWithSegment(compsegs_vector, *new_connection_set);
+        Compsegs::updateConnectionsWithSegment(compsegs_vector, grid, *new_connection_set);
         return new_connection_set;
     }
 }

src/opm/parser/eclipse/EclipseState/Schedule/Schedule.cpp

@@ -20,6 +20,7 @@
 #include <string>
 #include <vector>
 #include <stdexcept>
+#include <iostream>
 
 #include <boost/algorithm/string.hpp>
 #include <boost/date_time/posix_time/posix_time.hpp>
@@ -106,6 +107,9 @@ namespace Opm {
     {}
 
 
+
+
+
     std::time_t Schedule::getStartTime() const {
         return this->posixStartTime( );
     }
@@ -205,7 +209,7 @@ namespace Opm {
             handleWELSEGS(keyword, currentStep);
 
         else if (keyword.name() == "COMPSEGS")
-            handleCOMPSEGS(keyword, currentStep);
+            handleCOMPSEGS(keyword, currentStep, grid);
 
         else if (keyword.name() == "WELOPEN")
             handleWELOPEN(keyword, currentStep, parseContext);
@@ -1348,11 +1352,11 @@ namespace Opm {
         well.handleWELSEGS(keyword, currentStep);WellSegments newSegmentset;
     }
 
-    void Schedule::handleCOMPSEGS( const DeckKeyword& keyword, size_t currentStep) {
+    void Schedule::handleCOMPSEGS( const DeckKeyword& keyword, size_t currentStep, const EclipseGrid& grid) {
         const auto& record1 = keyword.getRecord(0);
         const std::string& well_name = record1.getItem("WELL").getTrimmedString(0);
         auto& well = this->m_wells.get( well_name );
-        well.handleCOMPSEGS(keyword, currentStep);
+        well.handleCOMPSEGS(keyword, grid, currentStep);
     }
 
     void Schedule::handleWGRUPCON( const DeckKeyword& keyword, size_t currentStep) {
@@ -1382,6 +1386,7 @@ namespace Opm {
             const std::string& childName = record.getItem("CHILD_GROUP").getTrimmedString(0);
             const std::string& parentName = record.getItem("PARENT_GROUP").getTrimmedString(0);
             newTree.update(childName, parentName);
+	    newTree.updateSeqIndex(childName, parentName);
 
             if (!hasGroup(parentName))
                 addGroup( parentName , currentStep );
@@ -1473,8 +1478,10 @@ namespace Opm {
         if (automaticShutInStr == "STOP") {
             automaticShutIn = false;
         }
-
+	
-        Well well(wellName,
+	const size_t wseqIndex = m_wells.size(); 
+        
+        Well well(wellName, wseqIndex,
                   headI, headJ, refDepth,
                   preferredPhase, m_timeMap,
                   timeStep,
@@ -1530,6 +1537,51 @@ namespace Opm {
         }
     }
 
+    std::vector< const Group* > Schedule::getChildGroups(const std::string& group_name, size_t timeStep) const {
+        if (!hasGroup(group_name))
+            throw std::invalid_argument("No such group: " + group_name);
+        {
+            const auto& group = getGroup( group_name );
+	    std::vector<const Group*> child_groups;
+
+            if (group.hasBeenDefined( timeStep )) {
+                const GroupTree& group_tree = getGroupTree( timeStep );
+                const auto& ch_grps = group_tree.children( group_name );
+		//for (const std::string& group_name : ch_grps) {
+		for ( auto it = ch_grps.begin() ; it != ch_grps.end(); it++) {
+                        child_groups.push_back( &getGroup(*it));
+                    }
+	    }
+	    return child_groups;
+	}
+    }
+
+        std::vector< const Well* > Schedule::getChildWells(const std::string& group_name, size_t timeStep) const {
+        if (!hasGroup(group_name))
+            throw std::invalid_argument("No such group: " + group_name);
+        {
+            const auto& group = getGroup( group_name );
+            std::vector<const Well*> wells;
+
+            if (group.hasBeenDefined( timeStep )) {
+                const GroupTree& group_tree = getGroupTree( timeStep );
+                const auto& child_groups = group_tree.children( group_name );
+
+                if (!child_groups.size()) {
+                    //for (const auto& well_name : group.getWells( timeStep )) {
+		    const auto& ch_wells = group.getWells( timeStep );
+                    for (auto it= ch_wells.begin(); it != ch_wells.end(); it++) {		      
+                        wells.push_back( getWell( *it ));
+                    }
+                }
+            }
+            return wells;
+        }
+    }
+
+    
+    
+    
     std::vector< const Well* > Schedule::getWells(size_t timeStep) const {
         if (timeStep >= m_timeMap.size()) {
             throw std::invalid_argument("Timestep to large");
@@ -1598,13 +1650,17 @@ namespace Opm {
     }
 
     void Schedule::addGroup(const std::string& groupName, size_t timeStep) {
-        m_groups.insert( groupName, Group { groupName, m_timeMap, timeStep } );
+	const size_t gseqIndex = m_groups.size(); 
+        m_groups.insert( groupName, Group { groupName, gseqIndex, m_timeMap, timeStep } );
         m_events.addEvent( ScheduleEvents::NEW_GROUP , timeStep );
     }
 
     size_t Schedule::numGroups() const {
         return m_groups.size();
     }
+    size_t Schedule::numGroups(size_t timeStep) const {
+        return this->getGroups( timeStep ).size();
+    }
 
     bool Schedule::hasGroup(const std::string& groupName) const {
         return m_groups.hasKey(groupName);
@@ -1645,6 +1701,25 @@ namespace Opm {
         return groups;
     }
 
+    std::vector< const Group* > Schedule::getGroups(size_t timeStep) const {
+      
+	if (timeStep >= m_timeMap.size()) {
+            throw std::invalid_argument("Timestep to large");
+        }
+        
+        auto defined = [=]( const Group& g ) {
+            return g.hasBeenDefined( timeStep );
+        };
+	
+        std::vector< const Group* > groups;
+
+        for( const auto& group : m_groups ) {
+	    if( !defined( group ) ) continue;
+	    groups.push_back( &group );
+	}
+        return groups;
+    }
+	
     void Schedule::addWellToGroup( Group& newGroup, Well& well , size_t timeStep) {
         const std::string currentGroupName = well.getGroupName(timeStep);
         if (currentGroupName != "") {

src/opm/parser/eclipse/EclipseState/Schedule/Well.cpp

@@ -33,13 +33,14 @@
 
 namespace Opm {
 
-    Well::Well(const std::string& name_, int headI,
+    Well::Well(const std::string& name_, const size_t& seqIndex_, int headI,
                int headJ, double refDepth , Phase preferredPhase,
                const TimeMap& timeMap, size_t creationTimeStep,
                WellCompletion::CompletionOrderEnum completionOrdering,
                bool allowCrossFlow, bool automaticShutIn)
         : m_creationTimeStep( creationTimeStep ),
           m_name( name_ ),
+          m_seqIndex( seqIndex_),
           m_status( timeMap, WellCommon::SHUT ),
           m_isAvailableForGroupControl( timeMap, true ),
           m_guideRate( timeMap, -1.0 ),
@@ -74,6 +75,10 @@ namespace Opm {
         return m_name;
     }
 
+        const size_t& Well::seqIndex() const {
+        return m_seqIndex;
+    }
+
 
     void Well::switchToProducer( size_t timeStep) {
         WellInjectionProperties p = getInjectionPropertiesCopy(timeStep);
@@ -376,6 +381,14 @@ namespace Opm {
         return *m_completions.back();
     }
 
+    const std::size_t Well::getTotNoConn() const {
+        return this->m_totNoConn;
+    }
+
+    void Well::setTotNoConn(std::size_t noConn)  {
+        m_totNoConn = noConn;
+    }
+    
     const std::string Well::getGroupName(size_t time_step) const {
         return m_groupName.get(time_step);
     }
@@ -676,15 +689,23 @@ namespace Opm {
 
     void Well::handleCOMPDAT(size_t time_step, const DeckRecord& record, const EclipseGrid& grid, const Eclipse3DProperties& eclipseProperties) {
         WellConnections * connections = new WellConnections(this->getConnections(time_step));
-        connections->loadCOMPDAT(record, grid, eclipseProperties);
+	std::size_t totNC = 0;
+        connections->loadCOMPDAT(record, grid, eclipseProperties, totNC);
+	if (totNC > 0) {
+	    this->setTotNoConn(totNC+1);
+	}
         this->updateWellConnections(time_step, connections);
     }
 
 
-    void Well::handleCOMPSEGS(const DeckKeyword& keyword, size_t time_step) {
+    void Well::handleCOMPSEGS(const DeckKeyword& keyword, const EclipseGrid& grid, size_t time_step) {
         const auto& segment_set = this->getWellSegments(time_step);
         const auto& completion_set = this->getConnections( time_step );
-        WellConnections * new_connection_set = newConnectionsWithSegments(keyword, completion_set, segment_set);
+	std::size_t totNC = 0;
+        WellConnections * new_connection_set = newConnectionsWithSegments(keyword, completion_set, segment_set, grid, totNC);
+	if (totNC > 0) {
+	    this->setTotNoConn(totNC+1);
+	}
         this->updateWellConnections(time_step, new_connection_set);
     }
 

src/opm/parser/eclipse/EclipseState/Schedule/WellConnections.cpp

@@ -20,6 +20,7 @@
 #include <cassert>
 #include <cmath>
 #include <limits>
+#include <iostream>
 
 #include <opm/parser/eclipse/Units/Units.hpp>
 #include <opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp>
@@ -137,11 +138,15 @@ namespace {
                                         double Kh,
                                         double rw,
                                         const int satTableId,
-                                        const WellCompletion::DirectionEnum direction)
+                                        const WellCompletion::DirectionEnum direction,
+				        const std::size_t seqIndex,
+				        const double segDistStart,
+					const double segDistEnd,
+					const bool defaultSatTabId)
     {
         int conn_i = (i < 0) ? this->headI : i;
         int conn_j = (j < 0) ? this->headJ : j;
-        Connection conn(conn_i, conn_j, k, complnum, depth, state, CF, Kh, rw, satTableId, direction);
+        Connection conn(conn_i, conn_j, k, complnum, depth, state, CF, Kh, rw, satTableId, direction, seqIndex, segDistStart, segDistEnd, defaultSatTabId);
         this->add(conn);
     }
 
@@ -154,7 +159,11 @@ namespace {
                                         double Kh,
                                         double rw,
                                         const int satTableId,
-                                        const WellCompletion::DirectionEnum direction)
+                                        const WellCompletion::DirectionEnum direction,
+					const std::size_t seqIndex,
+					const double segDistStart,
+					const double segDistEnd,
+					const bool defaultSatTabId)
     {
         int complnum = -(this->m_connections.size() + 1);
         this->addConnection(i,
@@ -167,10 +176,14 @@ namespace {
                             Kh,
                             rw,
                             satTableId,
-                            direction);
+                            direction,
+			    seqIndex,
+			    segDistStart,
+			    segDistEnd,
+			    defaultSatTabId);
     }
 
-    void WellConnections::loadCOMPDAT(const DeckRecord& record, const EclipseGrid& grid, const Eclipse3DProperties& eclipseProperties) {
+    void WellConnections::loadCOMPDAT(const DeckRecord& record, const EclipseGrid& grid, const Eclipse3DProperties& eclipseProperties, std::size_t& totNC) {
         const auto& permx = eclipseProperties.getDoubleGridProperty("PERMX").getData();
         const auto& permy = eclipseProperties.getDoubleGridProperty("PERMY").getData();
         const auto& permz = eclipseProperties.getDoubleGridProperty("PERMZ").getData();
@@ -272,21 +285,27 @@ namespace {
             auto prev = std::find_if( this->m_connections.begin(),
                                       this->m_connections.end(),
                                       same_ijk );
-
+	    // Only add connection for active grid cells
-            if (prev == this->m_connections.end()) {
+            if (grid.cellActive(I, J, k)) {
-                this->addConnection(I,J,k,
+		if (prev == this->m_connections.end()) {
+		    std::size_t noConn = this->m_connections.size();
+		    totNC = noConn;
+		    this->addConnection(I,J,k,
                                     grid.getCellDepth( I,J,k ),
                                     state,
                                     CF,
                                     Kh,
                                     rw,
                                     satTableId,
-                                    direction );
+                                    direction,
-            } else {
+				    noConn, 0., 0., defaultSatTable);
-                // The complnum value carries over; the rest of the state is fully specified by
+		} 
-                // the current COMPDAT keyword.
+		else {
-                int complnum = prev->complnum();
+		    std::size_t noConn = prev->getSeqIndex();
-                *prev = Connection(I,J,k,
+		    // The complnum value carries over; the rest of the state is fully specified by
+		    // the current COMPDAT keyword.
+		    int complnum = prev->complnum();
+		    *prev = Connection(I,J,k,
                                    complnum,
                                    grid.getCellDepth(I,J,k),
                                    state,
@@ -294,14 +313,15 @@ namespace {
                                    Kh,
                                    rw,
                                    satTableId,
-                                   direction );
+                                   direction,
-            }
+				   noConn, 0., 0., defaultSatTable);
-        }
+		}
+	    }
+	}
     }
 
 
 
-
     size_t WellConnections::size() const {
         return m_connections.size();
     }

tests/FIRST_SIM.DATA

@@ -9,6 +9,15 @@ UNIFIN
 DIMENS
  10 10 10 /
 
+WELLDIMS
+-- Item 1: NWMAX  (Maximum number of wells in model)
+-- Item 2: NCWMAX (Maximum number of connections per well)
+-- Item 3: NGMAX  (Maximum number of groups in model--excluding FIELD)
+-- Item 4: NWGMAX (Maximum number of wells or child groups per group)
+-- NWMAX  NCWMAX  NGMAX  NWGMAX
+   6      3       1      6
+/
+
 GRID
 DXV
 10*0.25 /

tests/FIRST_SIM_THPRES.DATA

@@ -13,10 +13,18 @@ UNIFIN
 DIMENS
  10 10 10 /
 
+WELLDIMS
+-- Item 1: NWMAX  (Maximum number of wells in model)
+-- Item 2: NCWMAX (Maximum number of connections per well)
+-- Item 3: NGMAX  (Maximum number of groups in model--excluding FIELD)
+-- Item 4: NWGMAX (Maximum number of wells or child groups per group)
+-- NWMAX  NCWMAX  NGMAX  NWGMAX
+   6      3       1      6
+/
+
 EQLDIMS
    10 /
 
-
 GRID
 DXV
 10*0.25 /

tests/SUMMARY_EFF_FAC.DATA

@@ -12,6 +12,15 @@ DIMENS
 REGDIMS
   10 /
 
+WELLDIMS
+-- Item 1: NWMAX  (Maximum number of wells in model)
+-- Item 2: NCWMAX (Maximum number of connections per well)
+-- Item 3: NGMAX  (Maximum number of groups in model--excluding FIELD)
+-- Item 4: NWGMAX (Maximum number of wells or child groups per group)
+-- NWMAX  NCWMAX  NGMAX  NWGMAX
+   3      2       5      2
+/
+
 OIL
 GAS
 WATER

tests/parser/ConnectionTests.cpp

@@ -69,9 +69,9 @@ BOOST_AUTO_TEST_CASE(CreateWellConnectionsOK) {
 BOOST_AUTO_TEST_CASE(AddCompletionSizeCorrect) {
     Opm::WellCompletion::DirectionEnum dir = Opm::WellCompletion::DirectionEnum::Z;
     Opm::WellConnections completionSet;
-    Opm::Connection completion1( 10,10,10, 1, 0.0, Opm::WellCompletion::OPEN , 99.88, 355.113, 0.25, 0, dir);
+    Opm::Connection completion1( 10,10,10, 1, 0.0, Opm::WellCompletion::OPEN , 99.88, 355.113, 0.25, 0, dir,0, 0., 0., true);
-    Opm::Connection completion2( 10,10,11, 1, 0.0, Opm::WellCompletion::SHUT , 99.88, 355.113, 0.25, 0, dir);
+    Opm::Connection completion2( 10,10,11, 1, 0.0, Opm::WellCompletion::SHUT , 99.88, 355.113, 0.25, 0, dir,0, 0., 0., true);
-completionSet.add( completion1 );
+    completionSet.add( completion1 );
     BOOST_CHECK_EQUAL( 1U , completionSet.size() );
 
     completionSet.add( completion2 );
@@ -83,8 +83,8 @@ completionSet.add( completion1 );
 
 BOOST_AUTO_TEST_CASE(WellConnectionsGetOutOfRangeThrows) {
     Opm::WellCompletion::DirectionEnum dir = Opm::WellCompletion::DirectionEnum::Z;
-    Opm::Connection completion1( 10,10,10, 1, 0.0, Opm::WellCompletion::OPEN , 99.88, 355.113, 0.25, 0, dir);
+    Opm::Connection completion1( 10,10,10, 1, 0.0, Opm::WellCompletion::OPEN , 99.88, 355.113, 0.25, 0, dir,0,0., 0., true);
-    Opm::Connection completion2( 10,10,11, 1, 0.0, Opm::WellCompletion::SHUT , 99.88, 355.113, 0.25, 0, dir);
+    Opm::Connection completion2( 10,10,11, 1, 0.0, Opm::WellCompletion::SHUT , 99.88, 355.113, 0.25, 0, dir,0,0., 0., true);
     Opm::WellConnections completionSet;
     completionSet.add( completion1 );
     BOOST_CHECK_EQUAL( 1U , completionSet.size() );
@@ -103,9 +103,9 @@ BOOST_AUTO_TEST_CASE(AddCompletionCopy) {
     Opm::WellConnections completionSet;
     Opm::WellCompletion::DirectionEnum dir = Opm::WellCompletion::DirectionEnum::Z;
 
-    Opm::Connection completion1( 10,10,10, 1, 0.0, Opm::WellCompletion::OPEN , 99.88, 355.113, 0.25, 0, dir);
+    Opm::Connection completion1( 10,10,10, 1, 0.0, Opm::WellCompletion::OPEN , 99.88, 355.113, 0.25, 0, dir,0,0., 0., true);
-    Opm::Connection completion2( 10,10,11, 1, 0.0, Opm::WellCompletion::SHUT , 99.88, 355.113, 0.25, 0, dir);
+    Opm::Connection completion2( 10,10,11, 1, 0.0, Opm::WellCompletion::SHUT , 99.88, 355.113, 0.25, 0, dir,0,0., 0., true);
-    Opm::Connection completion3( 10,10,12, 1, 0.0, Opm::WellCompletion::SHUT , 99.88, 355.113, 0.25, 0, dir);
+    Opm::Connection completion3( 10,10,12, 1, 0.0, Opm::WellCompletion::SHUT , 99.88, 355.113, 0.25, 0, dir,0,0., 0., true);
 
     completionSet.add( completion1 );
     completionSet.add( completion2 );
@@ -125,9 +125,9 @@ BOOST_AUTO_TEST_CASE(ActiveCompletions) {
     Opm::EclipseGrid grid(10,20,20);
     Opm::WellCompletion::DirectionEnum dir = Opm::WellCompletion::DirectionEnum::Z;
     Opm::WellConnections completions;
-    Opm::Connection completion1( 0,0,0, 1, 0.0, Opm::WellCompletion::OPEN , 99.88, 355.113, 0.25, 0, dir);
+    Opm::Connection completion1( 0,0,0, 1, 0.0, Opm::WellCompletion::OPEN , 99.88, 355.113, 0.25, 0, dir,0,0., 0., true);
-    Opm::Connection completion2( 0,0,1, 1, 0.0, Opm::WellCompletion::SHUT , 99.88, 355.113, 0.25, 0, dir);
+    Opm::Connection completion2( 0,0,1, 1, 0.0, Opm::WellCompletion::SHUT , 99.88, 355.113, 0.25, 0, dir,0,0., 0., true);
-    Opm::Connection completion3( 0,0,2, 1, 0.0, Opm::WellCompletion::SHUT , 99.88, 355.113, 0.25, 0, dir);
+    Opm::Connection completion3( 0,0,2, 1, 0.0, Opm::WellCompletion::SHUT , 99.88, 355.113, 0.25, 0, dir,0,0., 0., true);
 
     completions.add( completion1 );
     completions.add( completion2 );
@@ -151,8 +151,9 @@ Opm::WellConnections loadCOMPDAT(const std::string& compdat_keyword) {
     Opm::Eclipse3DProperties props(deck, tables, grid );
     const auto& keyword = deck.getKeyword("COMPDAT", 0);
     Opm::WellConnections connections;
+    std::size_t totnc = 0;
     for (const auto& rec : keyword)
-        connections.loadCOMPDAT(rec, grid, props);
+        connections.loadCOMPDAT(rec, grid, props, totnc);
 
     return connections;
 }

tests/parser/GroupTests.cpp

@@ -49,14 +49,14 @@ static TimeMap createXDaysTimeMap(size_t numDays) {
 
 BOOST_AUTO_TEST_CASE(CreateGroup_CorrectNameAndDefaultValues) {
     auto timeMap = createXDaysTimeMap(10);
-    Opm::Group group("G1" , timeMap , 0);
+    Opm::Group group("G1" , 1, timeMap , 0);
     BOOST_CHECK_EQUAL( "G1" , group.name() );
 }
 
 
 BOOST_AUTO_TEST_CASE(CreateGroupCreateTimeOK) {
     auto timeMap = createXDaysTimeMap(10);
-    Opm::Group group("G1" , timeMap , 5);
+    Opm::Group group("G1" , 1, timeMap , 5);
     BOOST_CHECK_EQUAL( false, group.hasBeenDefined( 4 ));
     BOOST_CHECK_EQUAL( true, group.hasBeenDefined( 5 ));
     BOOST_CHECK_EQUAL( true, group.hasBeenDefined( 6 ));
@@ -66,8 +66,8 @@ BOOST_AUTO_TEST_CASE(CreateGroupCreateTimeOK) {
 
 BOOST_AUTO_TEST_CASE(CreateGroup_SetInjectorProducer_CorrectStatusSet) {
     auto timeMap = createXDaysTimeMap(10);
-    Opm::Group group1("IGROUP" , timeMap , 0);
+    Opm::Group group1("IGROUP" , 1, timeMap , 0);
-    Opm::Group group2("PGROUP" , timeMap , 0);
+    Opm::Group group2("PGROUP" , 2, timeMap , 0);
 
     group1.setProductionGroup(0, true);
     BOOST_CHECK(group1.isProductionGroup(1));
@@ -91,7 +91,7 @@ BOOST_AUTO_TEST_CASE(CreateGroup_SetInjectorProducer_CorrectStatusSet) {
 
 BOOST_AUTO_TEST_CASE(InjectRateOK) {
     auto timeMap = createXDaysTimeMap(10);
-    Opm::Group group("G1" , timeMap , 0);
+    Opm::Group group("G1" , 1, timeMap , 0);
     BOOST_CHECK_EQUAL( 0 , group.getInjectionRate( 0 ));
     group.setInjectionRate( 2 , 100 );
     BOOST_CHECK_EQUAL( 100 , group.getInjectionRate( 2 ));
@@ -101,7 +101,7 @@ BOOST_AUTO_TEST_CASE(InjectRateOK) {
 
 BOOST_AUTO_TEST_CASE(ControlModeOK) {
     auto timeMap = createXDaysTimeMap(10);
-    Opm::Group group("G1" , timeMap , 0);
+    Opm::Group group("G1" , 1, timeMap , 0);
     BOOST_CHECK_EQUAL( Opm::GroupInjection::NONE , group.getInjectionControlMode( 0 ));
     group.setInjectionControlMode( 2 , Opm::GroupInjection::RESV );
     BOOST_CHECK_EQUAL( Opm::GroupInjection::RESV , group.getInjectionControlMode( 2 ));
@@ -112,7 +112,7 @@ BOOST_AUTO_TEST_CASE(ControlModeOK) {
 
 BOOST_AUTO_TEST_CASE(GroupChangePhaseSameTimeThrows) {
     auto timeMap = createXDaysTimeMap(10);
-    Opm::Group group("G1" , timeMap , 0);
+    Opm::Group group("G1" , 1, timeMap , 0);
     BOOST_CHECK_EQUAL( Opm::Phase::WATER , group.getInjectionPhase( 0 )); // Default phase - assumed WATER
     group.setInjectionPhase( 5, Opm::Phase::WATER );
     group.setInjectionPhase( 5, Opm::Phase::WATER );
@@ -126,7 +126,7 @@ BOOST_AUTO_TEST_CASE(GroupChangePhaseSameTimeThrows) {
 
 BOOST_AUTO_TEST_CASE(GroupMiscInjection) {
     auto timeMap = createXDaysTimeMap(10);
-    Opm::Group group("G1" , timeMap , 0);
+    Opm::Group group("G1" , 1, timeMap , 0);
 
     group.setSurfaceMaxRate( 3 , 100 );
     BOOST_CHECK_EQUAL( 100 , group.getSurfaceMaxRate( 5 ));
@@ -145,7 +145,7 @@ BOOST_AUTO_TEST_CASE(GroupMiscInjection) {
 
 BOOST_AUTO_TEST_CASE(GroupDoesNotHaveWell) {
     auto timeMap = createXDaysTimeMap(10);
-    Opm::Group group("G1" , timeMap , 0);
+    Opm::Group group("G1" , 1, timeMap , 0);
 
     BOOST_CHECK_EQUAL(false , group.hasWell("NO", 2));
     BOOST_CHECK_EQUAL(0U , group.numWells(2));
@@ -155,9 +155,9 @@ BOOST_AUTO_TEST_CASE(GroupDoesNotHaveWell) {
 BOOST_AUTO_TEST_CASE(GroupAddWell) {
 
     auto timeMap = createXDaysTimeMap( 10 );
-    Opm::Group group("G1" , timeMap , 0);
+    Opm::Group group("G1" , 1, timeMap , 0);
-    auto well1 = std::make_shared< Well >("WELL1", 0, 0, 0.0, Opm::Phase::OIL, timeMap, 0);
+    auto well1 = std::make_shared< Well >("WELL1", 1, 0, 0, 0.0, Opm::Phase::OIL, timeMap, 0);
-    auto well2 = std::make_shared< Well >("WELL2", 0, 0, 0.0, Opm::Phase::OIL, timeMap, 0);
+    auto well2 = std::make_shared< Well >("WELL2", 2, 0, 0, 0.0, Opm::Phase::OIL, timeMap, 0);
 
     BOOST_CHECK_EQUAL(0U , group.numWells(2));
     group.addWell( 3 , well1.get() );
@@ -192,9 +192,9 @@ BOOST_AUTO_TEST_CASE(GroupAddWell) {
 BOOST_AUTO_TEST_CASE(GroupAddAndDelWell) {
 
     auto timeMap = createXDaysTimeMap( 10 );
-    Opm::Group group("G1" , timeMap , 0);
+    Opm::Group group("G1" , 1, timeMap , 0);
-    auto well1 = std::make_shared< Well >("WELL1", 0, 0, 0.0, Opm::Phase::OIL, timeMap, 0);
+    auto well1 = std::make_shared< Well >("WELL1", 1, 0, 0, 0.0, Opm::Phase::OIL, timeMap, 0);
-    auto well2 = std::make_shared< Well >("WELL2", 0, 0, 0.0, Opm::Phase::OIL, timeMap, 0);
+    auto well2 = std::make_shared< Well >("WELL2", 2, 0, 0, 0.0, Opm::Phase::OIL, timeMap, 0);
 
     BOOST_CHECK_EQUAL(0U , group.numWells(2));
     group.addWell( 3 , well1.get() );
@@ -225,9 +225,9 @@ BOOST_AUTO_TEST_CASE(GroupAddAndDelWell) {
 
 BOOST_AUTO_TEST_CASE(getWells) {
     auto timeMap = createXDaysTimeMap( 10 );
-    Opm::Group group("G1" , timeMap , 0);
+    Opm::Group group("G1" , 1, timeMap , 0);
-    auto well1 = std::make_shared< Well >("WELL1", 0, 0, 0.0, Opm::Phase::OIL, timeMap, 0);
+    auto well1 = std::make_shared< Well >("WELL1", 1, 0, 0, 0.0, Opm::Phase::OIL, timeMap, 0);
-    auto well2 = std::make_shared< Well >("WELL2", 0, 0, 0.0, Opm::Phase::OIL, timeMap, 0);
+    auto well2 = std::make_shared< Well >("WELL2", 2, 0, 0, 0.0, Opm::Phase::OIL, timeMap, 0);
 
     group.addWell( 2 , well1.get() );
     group.addWell( 3 , well1.get() );

tests/parser/MultisegmentWellTests.cpp

@@ -44,14 +44,15 @@
 BOOST_AUTO_TEST_CASE(MultisegmentWellTest) {
     Opm::WellCompletion::DirectionEnum dir = Opm::WellCompletion::DirectionEnum::Z;
     Opm::WellConnections connection_set;
-    connection_set.add(Opm::Connection( 19, 0, 0, 1, 0.0, Opm::WellCompletion::OPEN , 200, 17.29, 0.25, 0, dir) );
+    Opm::EclipseGrid grid(20,20,20);
-    connection_set.add(Opm::Connection( 19, 0, 1, 1, 0.0, Opm::WellCompletion::OPEN , 200, 17.29, 0.25, 0, dir) );
+    connection_set.add(Opm::Connection( 19, 0, 0, 1, 0.0, Opm::WellCompletion::OPEN , 200, 17.29, 0.25, 0, dir,0, 0., 0., true) );
-    connection_set.add(Opm::Connection( 19, 0, 2, 1, 0.0, Opm::WellCompletion::OPEN , 200, 17.29, 0.25, 0, dir) );
+    connection_set.add(Opm::Connection( 19, 0, 1, 1, 0.0, Opm::WellCompletion::OPEN , 200, 17.29, 0.25, 0, dir,0, 0., 0., true) );
+    connection_set.add(Opm::Connection( 19, 0, 2, 1, 0.0, Opm::WellCompletion::OPEN , 200, 17.29, 0.25, 0, dir,0, 0., 0., true) );
 
-    connection_set.add(Opm::Connection( 18, 0, 1, 1, 0.0, Opm::WellCompletion::OPEN , 200, 17.29, 0.25, 0,  Opm::WellCompletion::DirectionEnum::X) );
+    connection_set.add(Opm::Connection( 18, 0, 1, 1, 0.0, Opm::WellCompletion::OPEN , 200, 17.29, 0.25, 0,  Opm::WellCompletion::DirectionEnum::X,0, 0., 0., true) );
-    connection_set.add(Opm::Connection( 17, 0, 1, 1, 0.0, Opm::WellCompletion::OPEN , 200, 17.29, 0.25, 0,  Opm::WellCompletion::DirectionEnum::X) );
+    connection_set.add(Opm::Connection( 17, 0, 1, 1, 0.0, Opm::WellCompletion::OPEN , 200, 17.29, 0.25, 0,  Opm::WellCompletion::DirectionEnum::X,0, 0., 0., true) );
-    connection_set.add(Opm::Connection( 16, 0, 1, 1, 0.0, Opm::WellCompletion::OPEN , 200, 17.29, 0.25, 0,  Opm::WellCompletion::DirectionEnum::X) );
+    connection_set.add(Opm::Connection( 16, 0, 1, 1, 0.0, Opm::WellCompletion::OPEN , 200, 17.29, 0.25, 0,  Opm::WellCompletion::DirectionEnum::X,0, 0., 0., true) );
-    connection_set.add(Opm::Connection( 15, 0, 1, 1, 0.0, Opm::WellCompletion::OPEN , 200, 17.29, 0.25, 0,  Opm::WellCompletion::DirectionEnum::X) );
+    connection_set.add(Opm::Connection( 15, 0, 1, 1, 0.0, Opm::WellCompletion::OPEN , 200, 17.29, 0.25, 0,  Opm::WellCompletion::DirectionEnum::X,0, 0., 0., true) );
 
     BOOST_CHECK_EQUAL( 7U , connection_set.size() );
 
@@ -87,8 +88,8 @@ BOOST_AUTO_TEST_CASE(MultisegmentWellTest) {
     segment_set.loadWELSEGS(welsegs);
 
     BOOST_CHECK_EQUAL(6U, segment_set.size());
-
+    std::size_t totNC = 0;
-    const Opm::WellConnections * new_connection_set = Opm::newConnectionsWithSegments(compsegs, connection_set, segment_set);
+    const Opm::WellConnections * new_connection_set = Opm::newConnectionsWithSegments(compsegs, connection_set, segment_set, grid, totNC);
 
     BOOST_CHECK_EQUAL(7U, new_connection_set->size());
 

tests/parser/ParseContextTests.cpp

@@ -712,4 +712,4 @@ BOOST_AUTO_TEST_CASE( test_invalid_wtemplate_config ) {
 
         BOOST_CHECK_THROW( Schedule( deckUnSupported , grid , eclipseProperties, Phases(true, true, true) , parseContext), std::invalid_argument );
     }
-}
+}

tests/parser/WellTests.cpp

@@ -63,7 +63,7 @@ inline std::ostream& operator<<( std::ostream& stream, const Well& well ) {
 
 BOOST_AUTO_TEST_CASE(CreateWell_CorrectNameAndDefaultValues) {
     auto timeMap = createXDaysTimeMap(10);
-    Opm::Well well("WELL1" ,  0, 0, 0.0, Opm::Phase::OIL, timeMap , 0);
+    Opm::Well well("WELL1" ,  1, 0, 0, 0.0, Opm::Phase::OIL, timeMap , 0);
     BOOST_CHECK_EQUAL( "WELL1" , well.name() );
     BOOST_CHECK_EQUAL(0.0 , well.getProductionPropertiesCopy(5).OilRate);
 }
@@ -71,10 +71,10 @@ BOOST_AUTO_TEST_CASE(CreateWell_CorrectNameAndDefaultValues) {
 BOOST_AUTO_TEST_CASE(CreateWell_Equals) {
     auto timeMap = createXDaysTimeMap(10);
     auto timeMap2 = createXDaysTimeMap(11);
-    Opm::Well well1("WELL1" ,  0, 0, 0.0, Opm::Phase::OIL, timeMap , 0);
+    Opm::Well well1("WELL1" ,  1, 0, 0, 0.0, Opm::Phase::OIL, timeMap , 0);
-    Opm::Well well2("WELL1" ,  0, 0, 0.0, Opm::Phase::OIL, timeMap , 0);
+    Opm::Well well2("WELL1" ,  2, 0, 0, 0.0, Opm::Phase::OIL, timeMap , 0);
-    Opm::Well well3("WELL3" ,  0, 0, 0.0, Opm::Phase::OIL, timeMap , 0);
+    Opm::Well well3("WELL3" ,  3, 0, 0, 0.0, Opm::Phase::OIL, timeMap , 0);
-    Opm::Well well4("WELL3" ,  0, 0, 0.0, Opm::Phase::OIL, timeMap2 , 0);
+    Opm::Well well4("WELL3" ,  4, 0, 0, 0.0, Opm::Phase::OIL, timeMap2 , 0);
     BOOST_CHECK_EQUAL( well1, well1 );
     BOOST_CHECK_EQUAL( well2, well1 );
     BOOST_CHECK( well1 == well2 );
@@ -87,7 +87,7 @@ BOOST_AUTO_TEST_CASE(CreateWell_Equals) {
 
 BOOST_AUTO_TEST_CASE(CreateWell_GetProductionPropertiesShouldReturnSameObject) {
     auto timeMap = createXDaysTimeMap(10);
-    Opm::Well well("WELL1" ,  0, 0, 0.0, Opm::Phase::OIL, timeMap , 0);
+    Opm::Well well("WELL1" ,  1, 0, 0, 0.0, Opm::Phase::OIL, timeMap , 0);
 
     BOOST_CHECK_EQUAL(&(well.getProductionProperties(5)), &(well.getProductionProperties(5)));
     BOOST_CHECK_EQUAL(&(well.getProductionProperties(8)), &(well.getProductionProperties(8)));
@@ -96,7 +96,7 @@ BOOST_AUTO_TEST_CASE(CreateWell_GetProductionPropertiesShouldReturnSameObject) {
 
 BOOST_AUTO_TEST_CASE(CreateWell_GetInjectionPropertiesShouldReturnSameObject) {
     auto timeMap = createXDaysTimeMap(10);
-    Opm::Well well("WELL1" ,  0, 0, 0.0, Opm::Phase::WATER, timeMap , 0);
+    Opm::Well well("WELL1" ,  1, 0, 0, 0.0, Opm::Phase::WATER, timeMap , 0);
 
     BOOST_CHECK_EQUAL(&(well.getInjectionProperties(5)), &(well.getInjectionProperties(5)));
     BOOST_CHECK_EQUAL(&(well.getInjectionProperties(8)), &(well.getInjectionProperties(8)));
@@ -105,7 +105,7 @@ BOOST_AUTO_TEST_CASE(CreateWell_GetInjectionPropertiesShouldReturnSameObject) {
 
 BOOST_AUTO_TEST_CASE(CreateWellCreateTimeStepOK) {
     auto timeMap = createXDaysTimeMap(10);
-    Opm::Well well("WELL1" , 0, 0, 0.0, Opm::Phase::OIL, timeMap , 5);
+    Opm::Well well("WELL1" , 1, 0, 0, 0.0, Opm::Phase::OIL, timeMap , 5);
     BOOST_CHECK_EQUAL( false , well.hasBeenDefined(0) );
     BOOST_CHECK_EQUAL( false , well.hasBeenDefined(4) );
     BOOST_CHECK_EQUAL( true , well.hasBeenDefined(5) );
@@ -116,7 +116,7 @@ BOOST_AUTO_TEST_CASE(CreateWellCreateTimeStepOK) {
 
 BOOST_AUTO_TEST_CASE(setWellProductionProperties_PropertiesSetCorrect) {
     auto timeMap = createXDaysTimeMap(10);
-    Opm::Well well("WELL1" , 0, 0, 0.0, Opm::Phase::OIL, timeMap , 0);
+    Opm::Well well("WELL1" , 1, 0, 0, 0.0, Opm::Phase::OIL, timeMap , 0);
 
     BOOST_CHECK_EQUAL(0.0 , well.getProductionPropertiesCopy( 5 ).OilRate);
     Opm::WellProductionProperties props;
@@ -149,7 +149,7 @@ BOOST_AUTO_TEST_CASE(setWellProductionProperties_PropertiesSetCorrect) {
 
 BOOST_AUTO_TEST_CASE(setOilRate_RateSetCorrect) {
     auto timeMap = createXDaysTimeMap(10);
-    Opm::Well well("WELL1" , 0, 0, 0.0, Opm::Phase::OIL, timeMap , 0);
+    Opm::Well well("WELL1" , 1, 0, 0, 0.0, Opm::Phase::OIL, timeMap , 0);
 
     BOOST_CHECK_EQUAL(0.0 , well.getProductionPropertiesCopy(5).OilRate);
     Opm::WellProductionProperties props;
@@ -161,7 +161,7 @@ BOOST_AUTO_TEST_CASE(setOilRate_RateSetCorrect) {
 
 BOOST_AUTO_TEST_CASE(seLiquidRate_RateSetCorrect) {
     auto timeMap = createXDaysTimeMap(10);
-    Opm::Well well("WELL1" ,  0, 0, 0.0, Opm::Phase::OIL, timeMap , 0);
+    Opm::Well well("WELL1" ,  1, 0, 0, 0.0, Opm::Phase::OIL, timeMap , 0);
 
     BOOST_CHECK_EQUAL(0.0 , well.getProductionPropertiesCopy(5).LiquidRate);
     Opm::WellProductionProperties props;
@@ -174,7 +174,7 @@ BOOST_AUTO_TEST_CASE(seLiquidRate_RateSetCorrect) {
 
 BOOST_AUTO_TEST_CASE(setPredictionModeProduction_ModeSetCorrect) {
     auto timeMap = createXDaysTimeMap(10);
-    Opm::Well well("WELL1" , 0, 0, 0.0, Opm::Phase::OIL, timeMap , 0);
+    Opm::Well well("WELL1" , 1, 0, 0, 0.0, Opm::Phase::OIL, timeMap , 0);
 
     BOOST_CHECK_EQUAL( true, well.getProductionPropertiesCopy(5).predictionMode);
     Opm::WellProductionProperties props;
@@ -187,7 +187,7 @@ BOOST_AUTO_TEST_CASE(setPredictionModeProduction_ModeSetCorrect) {
 
 BOOST_AUTO_TEST_CASE(setpredictionModeInjection_ModeSetCorrect) {
     auto timeMap = createXDaysTimeMap(10);
-    Opm::Well well("WELL1" , 0, 0, 0.0, Opm::Phase::WATER, timeMap , 0);
+    Opm::Well well("WELL1" , 1, 0, 0, 0.0, Opm::Phase::WATER, timeMap , 0);
 
     BOOST_CHECK_EQUAL( true, well.getInjectionPropertiesCopy(5).predictionMode);
     Opm::WellInjectionProperties props;
@@ -333,28 +333,13 @@ BOOST_AUTO_TEST_CASE(WellCOMPDATtestINPUT) {
 
 BOOST_AUTO_TEST_CASE(NewWellZeroCompletions) {
     auto timeMap = createXDaysTimeMap(10);
-    Opm::Well well("WELL1" , 0, 0, 0.0, Opm::Phase::OIL, timeMap , 0);
+    Opm::Well well("WELL1" , 1, 0, 0, 0.0, Opm::Phase::OIL, timeMap , 0);
     BOOST_CHECK_EQUAL( 0U , well.getConnections( 0 ).size() );
 }
 
-// Helper function for CompletionOrder test.
-inline Opm::Connection connection( int i, int j, int k, int complnum = 1 ) {
-    return Opm::Connection { i, j, k,
-            complnum,
-            k*1.0,
-            Opm::WellCompletion::AUTO,
-            99.88,
-            17.29,
-            0.25,
-            0,
-            Opm::WellCompletion::DirectionEnum::Z };
-}
-
-
-
 BOOST_AUTO_TEST_CASE(setGasRate_RateSetCorrect) {
     auto timeMap = createXDaysTimeMap(10);
-    Opm::Well well("WELL1" , 0, 0, 0.0, Opm::Phase::GAS, timeMap , 0);
+    Opm::Well well("WELL1" , 1, 0, 0, 0.0, Opm::Phase::GAS, timeMap , 0);
 
     BOOST_CHECK_EQUAL(0.0 , well.getProductionPropertiesCopy(5).GasRate);
     Opm::WellProductionProperties properties;
@@ -368,7 +353,7 @@ BOOST_AUTO_TEST_CASE(setGasRate_RateSetCorrect) {
 
 BOOST_AUTO_TEST_CASE(setWaterRate_RateSetCorrect) {
     auto timeMap = createXDaysTimeMap(10);
-    Opm::Well well("WELL1" , 0, 0, 0.0, Opm::Phase::WATER, timeMap , 0);
+    Opm::Well well("WELL1" , 1, 0, 0, 0.0, Opm::Phase::WATER, timeMap , 0);
 
     BOOST_CHECK_EQUAL(0.0 , well.getProductionPropertiesCopy(5).WaterRate);
     Opm::WellProductionProperties properties;
@@ -381,7 +366,7 @@ BOOST_AUTO_TEST_CASE(setWaterRate_RateSetCorrect) {
 
 BOOST_AUTO_TEST_CASE(setSurfaceInjectionRate_RateSetCorrect) {
     auto timeMap = createXDaysTimeMap(10);
-    Opm::Well well("WELL1" , 0, 0, 0.0, Opm::Phase::WATER, timeMap , 0);
+    Opm::Well well("WELL1" , 1, 0, 0, 0.0, Opm::Phase::WATER, timeMap , 0);
 
     BOOST_CHECK_EQUAL(0.0 , well.getInjectionPropertiesCopy(5).surfaceInjectionRate);
     Opm::WellInjectionProperties props(well.getInjectionPropertiesCopy(5));
@@ -400,7 +385,7 @@ BOOST_AUTO_TEST_CASE(setSurfaceInjectionRate_RateSetCorrect) {
 
 BOOST_AUTO_TEST_CASE(setReservoirInjectionRate_RateSetCorrect) {
     auto timeMap = createXDaysTimeMap(10);
-    Opm::Well well("WELL1" , 0, 0, 0.0, Opm::Phase::WATER, timeMap , 0);
+    Opm::Well well("WELL1" , 1, 0, 0, 0.0, Opm::Phase::WATER, timeMap , 0);
 
     BOOST_CHECK_EQUAL(0.0 , well.getInjectionPropertiesCopy(5).reservoirInjectionRate);
     Opm::WellInjectionProperties properties(well.getInjectionPropertiesCopy(5));
@@ -415,7 +400,7 @@ BOOST_AUTO_TEST_CASE(isProducerCorrectlySet) {
     // HACK: This test checks correctly setting of isProducer/isInjector. This property depends on which of
     //       WellProductionProperties/WellInjectionProperties is set last, independent of actual values.
     auto timeMap = createXDaysTimeMap(10);
-    Opm::Well well("WELL1" , 0, 0, 0.0, Opm::Phase::OIL, timeMap ,0);
+    Opm::Well well("WELL1" , 1, 0, 0, 0.0, Opm::Phase::OIL, timeMap ,0);
 
     /* 1: Well is created as producer */
     BOOST_CHECK_EQUAL( false , well.isInjector(0));
@@ -471,7 +456,7 @@ BOOST_AUTO_TEST_CASE(isProducerCorrectlySet) {
 
 BOOST_AUTO_TEST_CASE(GroupnameCorretlySet) {
     auto timeMap = createXDaysTimeMap(10);
-    Opm::Well well("WELL1" ,  0, 0, 0.0, Opm::Phase::WATER, timeMap ,0);
+    Opm::Well well("WELL1" ,  1, 0, 0, 0.0, Opm::Phase::WATER, timeMap ,0);
 
     BOOST_CHECK_EQUAL("" , well.getGroupName(2));
 
@@ -485,7 +470,7 @@ BOOST_AUTO_TEST_CASE(GroupnameCorretlySet) {
 
 BOOST_AUTO_TEST_CASE(addWELSPECS_setData_dataSet) {
     auto timeMap = createXDaysTimeMap(10);
-    Opm::Well well("WELL1", 23, 42, 2334.32, Opm::Phase::WATER, timeMap, 3);
+    Opm::Well well("WELL1", 1, 23, 42, 2334.32, Opm::Phase::WATER, timeMap, 3);
 
     BOOST_CHECK(!well.hasBeenDefined(2));
     BOOST_CHECK(well.hasBeenDefined(3));
@@ -498,7 +483,7 @@ BOOST_AUTO_TEST_CASE(addWELSPECS_setData_dataSet) {
 
 BOOST_AUTO_TEST_CASE(XHPLimitDefault) {
     auto timeMap = createXDaysTimeMap(10);
-    Opm::Well well("WELL1", 1, 2, 2334.32, Opm::Phase::WATER, timeMap, 0);
+    Opm::Well well("WELL1", 1, 1, 2, 2334.32, Opm::Phase::WATER, timeMap, 0);
 
 
     Opm::WellProductionProperties productionProps(well.getProductionPropertiesCopy(1));
@@ -519,7 +504,7 @@ BOOST_AUTO_TEST_CASE(XHPLimitDefault) {
 
 BOOST_AUTO_TEST_CASE(InjectorType) {
     auto timeMap = createXDaysTimeMap(10);
-    Opm::Well well("WELL1", 1, 2, 2334.32, Opm::Phase::WATER, timeMap, 0);
+    Opm::Well well("WELL1", 1, 1, 2, 2334.32, Opm::Phase::WATER, timeMap, 0);
 
     Opm::WellInjectionProperties injectionProps(well.getInjectionPropertiesCopy(1));
     injectionProps.injectorType = Opm::WellInjector::WATER;
@@ -531,15 +516,13 @@ BOOST_AUTO_TEST_CASE(InjectorType) {
 
 
 
-
-
 /*****************************************************************/
 
 
 BOOST_AUTO_TEST_CASE(WellHaveProductionControlLimit) {
 
     auto timeMap = createXDaysTimeMap(20);
-    Opm::Well well("WELL1", 1, 2, 2334.32, Opm::Phase::OIL, timeMap, 0);
+    Opm::Well well("WELL1", 1, 1, 2, 2334.32, Opm::Phase::OIL, timeMap, 0);
 
 
     BOOST_CHECK( !well.getProductionPropertiesCopy(1).hasProductionControl( Opm::WellProducer::ORAT ));
@@ -589,7 +572,7 @@ BOOST_AUTO_TEST_CASE(WellHaveProductionControlLimit) {
 BOOST_AUTO_TEST_CASE(WellHaveInjectionControlLimit) {
 
     auto timeMap = createXDaysTimeMap(20);
-    Opm::Well well("WELL1", 1, 2, 2334.32, Opm::Phase::WATER, timeMap, 0);
+    Opm::Well well("WELL1", 1, 1, 2, 2334.32, Opm::Phase::WATER, timeMap, 0);
 
     BOOST_CHECK( !well.getInjectionPropertiesCopy(1).hasInjectionControl( Opm::WellInjector::RATE ));
     BOOST_CHECK( !well.getInjectionPropertiesCopy(1).hasInjectionControl( Opm::WellInjector::RESV ));
@@ -633,7 +616,7 @@ BOOST_AUTO_TEST_CASE(WellHaveInjectionControlLimit) {
 
 BOOST_AUTO_TEST_CASE(WellSetAvailableForGroupControl_ControlSet) {
     auto timeMap = createXDaysTimeMap(20);
-    Opm::Well well("WELL1", 1, 2, 2334.32, Opm::Phase::WATER, timeMap, 0);
+    Opm::Well well("WELL1", 1, 1, 2, 2334.32, Opm::Phase::WATER, timeMap, 0);
 
     BOOST_CHECK(well.isAvailableForGroupControl(10));
     well.setAvailableForGroupControl(12, false);
@@ -644,7 +627,7 @@ BOOST_AUTO_TEST_CASE(WellSetAvailableForGroupControl_ControlSet) {
 
 BOOST_AUTO_TEST_CASE(WellSetGuideRate_GuideRateSet) {
     auto timeMap = createXDaysTimeMap(20);
-    Opm::Well well("WELL1", 1, 2, 2334.32, Opm::Phase::WATER, timeMap, 0);
+    Opm::Well well("WELL1", 1, 1, 2, 2334.32, Opm::Phase::WATER, timeMap, 0);
 
     BOOST_CHECK_LT(well.getGuideRate(0), 0);
     well.setGuideRate(1, 32.2);
@@ -654,7 +637,7 @@ BOOST_AUTO_TEST_CASE(WellSetGuideRate_GuideRateSet) {
 
 BOOST_AUTO_TEST_CASE(WellGuideRatePhase_GuideRatePhaseSet) {
     auto timeMap = createXDaysTimeMap(20);
-    Opm::Well well("WELL1", 1, 2, 2334.32, Opm::Phase::WATER, timeMap, 0);
+    Opm::Well well("WELL1", 1, 1, 2, 2334.32, Opm::Phase::WATER, timeMap, 0);
     BOOST_CHECK_EQUAL(Opm::GuideRate::UNDEFINED, well.getGuideRatePhase(0));
     well.setGuideRatePhase(3, Opm::GuideRate::RAT);
     BOOST_CHECK_EQUAL(Opm::GuideRate::UNDEFINED, well.getGuideRatePhase(2));
@@ -663,7 +646,7 @@ BOOST_AUTO_TEST_CASE(WellGuideRatePhase_GuideRatePhaseSet) {
 
 BOOST_AUTO_TEST_CASE(WellEfficiencyFactorSet) {
     auto timeMap = createXDaysTimeMap(20);
-    Opm::Well well("WELL1", 1, 2, 2334.32, Opm::Phase::WATER, timeMap, 0);
+    Opm::Well well("WELL1", 1, 1, 2, 2334.32, Opm::Phase::WATER, timeMap, 0);
     BOOST_CHECK_EQUAL(1.0, well.getEfficiencyFactor(0));
     well.setEfficiencyFactor(3, 0.9);
     BOOST_CHECK_EQUAL(1.0, well.getEfficiencyFactor(0));
@@ -672,7 +655,7 @@ BOOST_AUTO_TEST_CASE(WellEfficiencyFactorSet) {
 
 BOOST_AUTO_TEST_CASE(WellSetScalingFactor_ScalingFactorSetSet) {
     auto timeMap = createXDaysTimeMap(20);
-    Opm::Well well("WELL1", 1, 2, 2334.32, Opm::Phase::WATER, timeMap, 0);
+    Opm::Well well("WELL1", 1, 1, 2, 2334.32, Opm::Phase::WATER, timeMap, 0);
     BOOST_CHECK_EQUAL(1.0, well.getGuideRateScalingFactor(0));
     well.setGuideRateScalingFactor(4, 0.6);
     BOOST_CHECK_EQUAL(1.0, well.getGuideRateScalingFactor(3));

tests/summary_deck.DATA

@@ -17,6 +17,15 @@ DIMENS
 REGDIMS
   10 /
 
+WELLDIMS
+-- Item 1: NWMAX  (Maximum number of wells in model)
+-- Item 2: NCWMAX (Maximum number of connections per well)
+-- Item 3: NGMAX  (Maximum number of groups in model--excluding FIELD)
+-- Item 4: NWGMAX (Maximum number of wells or child groups per group)
+-- NWMAX  NCWMAX  NGMAX  NWGMAX
+   5      2       3      2
+/
+
 OIL
 GAS
 WATER

tests/summary_deck_non_constant_porosity.DATA

@@ -21,6 +21,15 @@ OIL
 GAS
 WATER
 
+WELLDIMS
+-- Item 1: NWMAX  (Maximum number of wells in model)
+-- Item 2: NCWMAX (Maximum number of connections per well)
+-- Item 3: NGMAX  (Maximum number of groups in model--excluding FIELD)
+-- Item 4: NWGMAX (Maximum number of wells or child groups per group)
+-- NWMAX  NCWMAX  NGMAX  NWGMAX
+   4      2       3      2
+/
+
 GRID
 
 DX

tests/test_AggregateWellData.cpp

@@ -370,9 +370,10 @@ BOOST_AUTO_TEST_CASE (Declared_Well_Data)
 
     BOOST_CHECK_EQUAL(ih.nwells, MockIH::Sz{2});
 
+    const auto smry = sim_state();
     auto awd = Opm::RestartIO::Helpers::AggregateWellData{ih.value};
     awd.captureDeclaredWellData(simCase.sched,
-                                simCase.es.getUnits(), rptStep);
+                                simCase.es.getUnits(), rptStep, smry, ih.value);
 
     // IWEL (OP_1)
     {
@@ -427,15 +428,24 @@ BOOST_AUTO_TEST_CASE (Declared_Well_Data)
         const auto i0 = 0*ih.nswelz;
 
         const auto& swell = awd.getSWell();
-        BOOST_CHECK_CLOSE(swell[i0 + Ix::OilRateTarget] , 20.0e3f, 1.0e-7f);
+        BOOST_CHECK_CLOSE(swell[i0 + Ix::OilRateTarget], 20.0e3f, 1.0e-7f);
-        BOOST_CHECK_CLOSE(swell[i0 + Ix::WatRateTarget] , 1.0e20f, 1.0e-7f);
-        BOOST_CHECK_CLOSE(swell[i0 + Ix::GasRateTarget] , 1.0e20f, 1.0e-7f);
-        BOOST_CHECK_CLOSE(swell[i0 + Ix::LiqRateTarget] , 1.0e20f, 1.0e-7f);
-        BOOST_CHECK_CLOSE(swell[i0 + Ix::ResVRateTarget], 1.0e20f, 1.0e-7f);
-        BOOST_CHECK_CLOSE(swell[i0 + Ix::THPTarget]     , 1.0e20f, 1.0e-7f);
-        BOOST_CHECK_CLOSE(swell[i0 + Ix::BHPTarget]     , 1000.0f, 1.0e-7f);
 
-        BOOST_CHECK_CLOSE(swell[i0 + Ix::DatumDepth]    , 0.375f, 1.0e-7f);
+        // No WRAT limit
+        BOOST_CHECK_CLOSE(swell[i0 + Ix::WatRateTarget], 1.0e20f, 1.0e-7f);
+
+        // No GRAT limit
+        BOOST_CHECK_CLOSE(swell[i0 + Ix::GasRateTarget], 1.0e20f, 1.0e-7f);
+
+        // LRAT limit derived from ORAT + WRAT (= ORAT + 0.0)
+        BOOST_CHECK_CLOSE(swell[i0 + Ix::LiqRateTarget], 20.0e3f, 1.0e-7f);
+
+        // No direct limit, extract value from 'smry' (WVPR:OP_1)
+        BOOST_CHECK_CLOSE(swell[i0 + Ix::ResVRateTarget], 4.0f, 1.0e-7f);
+
+        // No THP limit
+        BOOST_CHECK_CLOSE(swell[i0 + Ix::THPTarget] , 1.0e20f, 1.0e-7f);
+        BOOST_CHECK_CLOSE(swell[i0 + Ix::BHPTarget] , 1000.0f, 1.0e-7f);
+        BOOST_CHECK_CLOSE(swell[i0 + Ix::DatumDepth],  0.375f, 1.0e-7f);
     }
 
     // SWEL (OP_2)
@@ -513,7 +523,7 @@ BOOST_AUTO_TEST_CASE (Dynamic_Well_Data_Step1)
     const auto smry = sim_state();
     auto awd = Opm::RestartIO::Helpers::AggregateWellData{ih.value};
 
-    awd.captureDynamicWellData(simCase.sched, rptStep, xw, smry);
+    awd.captureDynamicWellData(simCase.sched, rptStep, true, xw, smry);
 
     // IWEL (OP_1)
     {
@@ -609,7 +619,7 @@ BOOST_AUTO_TEST_CASE (Dynamic_Well_Data_Step2)
     const auto smry = sim_state();
     auto awd = Opm::RestartIO::Helpers::AggregateWellData{ih.value};
 
-    awd.captureDynamicWellData(simCase.sched, rptStep, xw, smry);
+    awd.captureDynamicWellData(simCase.sched, rptStep, true, xw, smry);
 
     // IWEL (OP_1) -- closed producer
     {

tests/test_InteHEAD.cpp

@@ -328,8 +328,8 @@ BOOST_AUTO_TEST_CASE(Time_and_report_step)
 
     const auto& v = ih.data();
 
-    BOOST_CHECK_EQUAL(v[67], 12); // TSTEP
+    BOOST_CHECK_EQUAL(v[67], 12);    // TSTEP
-    BOOST_CHECK_EQUAL(v[68],  2); // REP_STEP
+    BOOST_CHECK_EQUAL(v[68], 2 + 1); // REP_STEP (= sim_step + 1)
 }
 
 BOOST_AUTO_TEST_CASE(Tuning_param)
@@ -415,6 +415,17 @@ BOOST_AUTO_TEST_CASE(regionDimensions)
     BOOST_CHECK_EQUAL(v[99], nmfipr); // NMFIPR
 }
 
+BOOST_AUTO_TEST_CASE(ngroups)
+{
+    const auto ngroup  = 8;
+
+    const auto ih = Opm::RestartIO::InteHEAD{}
+        .ngroups({ ngroup });
+	
+    const auto& v = ih.data();
+
+    BOOST_CHECK_EQUAL(v[18], ngroup);  // NGRP
+}
 BOOST_AUTO_TEST_CASE(SimulationDate)
 {
     const auto input = std::string { R"(

tests/test_LogiHEAD.cpp

@@ -39,18 +39,7 @@ BOOST_AUTO_TEST_CASE(Radial_Settings_and_Init)
     BOOST_CHECK_EQUAL(v[  1], true);  //
     BOOST_CHECK_EQUAL(v[  3], false); // E300 Radial
     BOOST_CHECK_EQUAL(v[  4], true);  // E100 Radial
-    BOOST_CHECK_EQUAL(v[ 16], true); //
+    BOOST_CHECK_EQUAL(v[ 75], true);  // MS Well Simulation Case
-    BOOST_CHECK_EQUAL(v[ 18], true); //
-    BOOST_CHECK_EQUAL(v[ 31], true); //
-    BOOST_CHECK_EQUAL(v[ 44], true); //
-    BOOST_CHECK_EQUAL(v[ 75], true); // MS Well Simulation Case
-    BOOST_CHECK_EQUAL(v[ 76], true); //
-    BOOST_CHECK_EQUAL(v[ 87], true); //
-    BOOST_CHECK_EQUAL(v[ 99], true); //
-    BOOST_CHECK_EQUAL(v[113], true); //
-    BOOST_CHECK_EQUAL(v[114], true); //
-    BOOST_CHECK_EQUAL(v[115], true); //
-    BOOST_CHECK_EQUAL(v[117], true); //
 }
 
 BOOST_AUTO_TEST_SUITE_END()

tests/test_Restart.cpp

@@ -346,6 +346,66 @@ data::Solution mkSolution( int numCells ) {
     return sol;
 }
 
+Opm::SummaryState sim_state()
+{
+    auto state = Opm::SummaryState{};
+
+    state.add("WOPR:OP_1" ,    1.0);
+    state.add("WWPR:OP_1" ,    2.0);
+    state.add("WGPR:OP_1" ,    3.0);
+    state.add("WVPR:OP_1" ,    4.0);
+    state.add("WOPT:OP_1" ,   10.0);
+    state.add("WWPT:OP_1" ,   20.0);
+    state.add("WGPT:OP_1" ,   30.0);
+    state.add("WVPT:OP_1" ,   40.0);
+    state.add("WWIR:OP_1" ,    0.0);
+    state.add("WGIR:OP_1" ,    0.0);
+    state.add("WWIT:OP_1" ,    0.0);
+    state.add("WGIT:OP_1" ,    0.0);
+    state.add("WWCT:OP_1" ,    0.625);
+    state.add("WGOR:OP_1" ,  234.5);
+    state.add("WBHP:OP_1" ,  314.15);
+    state.add("WGVIR:OP_1",    0.0);
+    state.add("WWVIR:OP_1",    0.0);
+
+    state.add("WOPR:OP_2" ,    0.0);
+    state.add("WWPR:OP_2" ,    0.0);
+    state.add("WGPR:OP_2" ,    0.0);
+    state.add("WVPR:OP_2" ,    0.0);
+    state.add("WOPT:OP_2" ,    0.0);
+    state.add("WWPT:OP_2" ,    0.0);
+    state.add("WGPT:OP_2" ,    0.0);
+    state.add("WVPT:OP_2" ,    0.0);
+    state.add("WWIR:OP_2" ,  100.0);
+    state.add("WGIR:OP_2" ,  200.0);
+    state.add("WWIT:OP_2" , 1000.0);
+    state.add("WGIT:OP_2" , 2000.0);
+    state.add("WWCT:OP_2" ,    0.0);
+    state.add("WGOR:OP_2" ,    0.0);
+    state.add("WBHP:OP_2" ,  400.6);
+    state.add("WGVIR:OP_2", 1234.0);
+    state.add("WWVIR:OP_2", 4321.0);
+
+    state.add("WOPR:OP_3" ,   11.0);
+    state.add("WWPR:OP_3" ,   12.0);
+    state.add("WGPR:OP_3" ,   13.0);
+    state.add("WVPR:OP_3" ,   14.0);
+    state.add("WOPT:OP_3" ,  110.0);
+    state.add("WWPT:OP_3" ,  120.0);
+    state.add("WGPT:OP_3" ,  130.0);
+    state.add("WVPT:OP_3" ,  140.0);
+    state.add("WWIR:OP_3" ,    0.0);
+    state.add("WGIR:OP_3" ,    0.0);
+    state.add("WWIT:OP_3" ,    0.0);
+    state.add("WGIT:OP_3" ,    0.0);
+    state.add("WWCT:OP_3" ,    0.0625);
+    state.add("WGOR:OP_3" , 1234.5);
+    state.add("WBHP:OP_3" ,  314.15);
+    state.add("WGVIR:OP_3",    0.0);
+    state.add("WWVIR:OP_3",   43.21);
+
+    return state;
+}
 
 RestartValue first_sim(const EclipseState& es, EclipseIO& eclWriter, bool write_double) {
     const auto& grid = es.getInputGrid();
@@ -362,7 +422,7 @@ RestartValue first_sim(const EclipseState& es, EclipseIO& eclWriter, bool write_
                              false,
                              first_step - start_time,
                              restart_value,
-                             {}, {}, {}, write_double);
+			     {}, {}, {}, write_double);
 
     return restart_value;
 }
@@ -443,6 +503,7 @@ BOOST_AUTO_TEST_CASE(ECL_FORMATTED) {
         auto num_cells = setup.grid.getNumActive( );
         auto cells = mkSolution( num_cells );
         auto wells = mkWells();
+	auto sumState = sim_state();
         {
             RestartValue restart_value(cells, wells);
 
@@ -454,6 +515,7 @@ BOOST_AUTO_TEST_CASE(ECL_FORMATTED) {
                             setup.es,
                             setup.grid,
                             setup.schedule,
+			    sumState,
                             true);
 
             {
@@ -472,6 +534,7 @@ BOOST_AUTO_TEST_CASE(ECL_FORMATTED) {
                             setup.es,
                             setup.grid,
                             setup.schedule,
+			    sumState,
                             true);
             {
                 ecl_file_type * rst_file = ecl_file_open( "ECL_FILE.UNRST" , 0 );
@@ -540,13 +603,15 @@ BOOST_AUTO_TEST_CASE(WriteWrongSOlutionSize) {
         auto num_cells = setup.grid.getNumActive( ) + 1;
         auto cells = mkSolution( num_cells );
         auto wells = mkWells();
+	Opm::SummaryState sumState;
 
         BOOST_CHECK_THROW( RestartIO::save("FILE.UNRST", 1 ,
                                            100,
                                            RestartValue(cells, wells),
                                            setup.es,
                                            setup.grid ,
-                                           setup.schedule),
+                                           setup.schedule,
+					   sumState),
                                            std::runtime_error);
     }
     test_work_area_free(test_area);
@@ -583,14 +648,17 @@ BOOST_AUTO_TEST_CASE(ExtraData_content) {
         const auto& units = setup.es.getUnits();
         {
             RestartValue restart_value(cells, wells);
+            const auto sumState = sim_state();
 
             restart_value.addExtra("EXTRA", UnitSystem::measure::pressure, {10,1,2,3});
+
             RestartIO::save("FILE.UNRST", 1 ,
                             100,
                             restart_value,
                             setup.es,
                             setup.grid,
-                            setup.schedule);
+                            setup.schedule,
+                            sumState);
 
             {
                 ecl_file_type * f = ecl_file_open( "FILE.UNRST" , 0 );
@@ -600,8 +668,8 @@ BOOST_AUTO_TEST_CASE(ExtraData_content) {
                     BOOST_CHECK_EQUAL( ecl_kw_get_header( ex) , "EXTRA" );
                     BOOST_CHECK_EQUAL(  4 , ecl_kw_get_size( ex ));
 
-                    BOOST_CHECK_CLOSE( 10 ,                                                units.to_si( UnitSystem::measure::pressure, ecl_kw_iget_double( ex, 0 )), 0.00001);
+                    BOOST_CHECK_CLOSE( 10 , units.to_si( UnitSystem::measure::pressure, ecl_kw_iget_double( ex, 0 )), 0.00001);
-                    BOOST_CHECK_CLOSE( units.from_si( UnitSystem::measure::pressure, 3)  , ecl_kw_iget_double( ex, 3 ),                                              0.00001);
+                    BOOST_CHECK_CLOSE( units.from_si( UnitSystem::measure::pressure, 3), ecl_kw_iget_double( ex, 3 ), 0.00001);
                 }
                 ecl_file_close( f );
             }
@@ -609,11 +677,17 @@ BOOST_AUTO_TEST_CASE(ExtraData_content) {
             BOOST_CHECK_THROW( RestartIO::load( "FILE.UNRST" , 1 , {}, setup.es, setup.grid , setup.schedule,
                                                 {{"NOT-THIS", UnitSystem::measure::identity, true}}) , std::runtime_error );
             {
-              const auto rst_value = RestartIO::load( "FILE.UNRST" , 1 , { RestartKey("SWAT", UnitSystem::measure::identity),
+                const auto rst_value = RestartIO::load(
-                                                                           RestartKey("NO", UnitSystem::measure::identity, false)},
+                    "FILE.UNRST" , 1 ,
-                setup.es, setup.grid , setup.schedule,
+                    /* solution_keys = */ {
-                {{"EXTRA", UnitSystem::measure::pressure, true},
+                        RestartKey("SWAT", UnitSystem::measure::identity),
-                 {"EXTRA2", UnitSystem::measure::identity, false}});
+                        RestartKey("NO"  , UnitSystem::measure::identity, false)
+                    },
+                    setup.es, setup.grid , setup.schedule,
+                    /* extra_keys = */ {
+                        {"EXTRA" , UnitSystem::measure::pressure, true}  ,
+                        {"EXTRA2", UnitSystem::measure::identity, false}
+                    });
 
                 BOOST_CHECK(!rst_value.hasExtra("EXTRA2"));
                 BOOST_CHECK( rst_value.hasExtra("EXTRA"));
@@ -656,13 +730,16 @@ BOOST_AUTO_TEST_CASE(STORE_THPRES) {
             */
 
             restart_value.addExtra("THRESHPR", UnitSystem::measure::pressure, {0,1});
+            const auto sumState = sim_state();
+
             /* THPRES data has wrong size in extra container. */
             BOOST_CHECK_THROW( RestartIO::save("FILE.UNRST", 1 ,
                                                100,
                                                restart_value,
                                                setup.es,
                                                setup.grid,
-                                               setup.schedule), std::runtime_error);
+                                               setup.schedule,
+                                               sumState), std::runtime_error);
 
             int num_regions = setup.es.getTableManager().getEqldims().getNumEquilRegions();
             std::vector<double>  thpres(num_regions * num_regions, 78);
@@ -674,7 +751,7 @@ BOOST_AUTO_TEST_CASE(STORE_THPRES) {
                             restart_value2,
                             setup.es,
                             setup.grid,
-                            setup.schedule);
+                            setup.schedule, sumState);
 
             {
                 ecl_file_type * rst_file = ecl_file_open("FILE2.UNRST", 0);

tests/test_serialize_ICON.cpp

@@ -19,6 +19,10 @@
 
 #include <config.h>
 
+#include <iostream> // @@
+#include <algorithm> // @@
+#include <iterator> // @@ 
+
 #define BOOST_TEST_MODULE serialize_ICON_TEST
 #include <boost/test/unit_test.hpp>
 
@@ -84,6 +88,11 @@ BOOST_AUTO_TEST_CASE( serialize_icon_test )
             }
             w_offset += (ICONZ * ncwmax);
         }
-    }
-};
 
+        std::copy(icondata.begin(),
+                  icondata.end(),
+                  std::ostream_iterator<int>(std::cout, " "));
+        std::cout << std::endl;
+        BOOST_CHECK_EQUAL(1, 1);// @@@
+    }
+}

tests/testblackoilstate3.DATA

@@ -17,6 +17,15 @@ METRIC
 DIMENS
     10  10   10   /
 
+WELLDIMS
+-- Item 1: NWMAX  (Maximum number of wells in model)
+-- Item 2: NCWMAX (Maximum number of connections per well)
+-- Item 3: NGMAX  (Maximum number of groups in model--excluding FIELD)
+-- Item 4: NWGMAX (Maximum number of wells or child groups per group)
+-- NWMAX  NCWMAX  NGMAX  NWGMAX
+   9      10      2      6
+/
+
 GRID
 DXV 
 10*10 /