Merge remote-tracking branch 'upstream/master' into master-refactor-for-cpgrid-support

Manually resolved conflicts:
	opm/core/io/eclipse/EclipseWriter.cpp
	opm/core/io/eclipse/EclipseWriter.hpp
	opm/core/props/BlackoilPropertiesFromDeck.cpp
	opm/core/simulator/initState_impl.hpp
	opm/core/wells/WellsManager.cpp
	opm/core/wells/WellsManager.hpp

CMakeLists.txt

@@ -46,6 +46,7 @@ include (CMakeLists_files.cmake)
 
 macro (config_hook)
 #	opm_need_version_of ("dune-common")
+  opm_need_version_of ("dune-istl")
 	list (APPEND ${project}_CONFIG_IMPL_VARS
 		HAVE_DUNE_ISTL
 		)

CMakeLists_files.cmake

@@ -164,6 +164,7 @@ list (APPEND TEST_SOURCE_FILES
 	tests/test_wachspresscoord.cpp
 	tests/test_column_extract.cpp
 	tests/test_geom2d.cpp
+	tests/test_linearsolver.cpp
 	tests/test_param.cpp
 	tests/test_blackoilfluid.cpp
 	tests/test_shadow.cpp

examples/compute_tof.cpp

@@ -47,6 +47,9 @@
 #include <opm/core/tof/TofReorder.hpp>
 #include <opm/core/tof/TofDiscGalReorder.hpp>
 
+#include <opm/parser/eclipse/Parser/Parser.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
+
 #include <memory>
 #include <boost/filesystem.hpp>
 
@@ -111,13 +114,16 @@ try
     double gravity[3] = { 0.0 };
     if (use_deck) {
         std::string deck_filename = param.get<std::string>("deck_filename");
+        Opm::ParserPtr parser(new Opm::Parser());
+        Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(parser->parseFile(deck_filename)));
+
         deck.reset(new EclipseGridParser(deck_filename));
         // Grid init
         grid.reset(new GridManager(*deck));
         // Rock and fluid init
         props.reset(new IncompPropertiesFromDeck(*deck, *grid->c_grid()));
         // Wells init.
-        wells.reset(new Opm::WellsManager(*deck, *grid->c_grid(), props->permeability()));
+        wells.reset(new Opm::WellsManager(eclipseState , 0 , *grid->c_grid(), props->permeability()));
         // Gravity.
         gravity[2] = deck->hasField("NOGRAV") ? 0.0 : unit::gravity;
         // Init state variables (saturation and pressure).

examples/sim_2p_comp_reorder.cpp

@@ -44,6 +44,10 @@
 #include <opm/core/simulator/WellState.hpp>
 #include <opm/core/simulator/SimulatorCompressibleTwophase.hpp>
 
+#include <opm/parser/eclipse/Parser/Parser.hpp>
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+
+
 #include <memory>
 #include <boost/filesystem.hpp>
 
@@ -80,6 +84,7 @@ try
 
     // If we have a "deck_filename", grid and props will be read from that.
     bool use_deck = param.has("deck_filename");
+    EclipseStateConstPtr eclipseState;
     std::unique_ptr<EclipseGridParser> deck;
     std::unique_ptr<GridManager> grid;
     std::unique_ptr<BlackoilPropertiesInterface> props;
@@ -89,7 +94,10 @@ try
     // int max_well_control_iterations = 0;
     double gravity[3] = { 0.0 };
     if (use_deck) {
+        ParserPtr parser(new Opm::Parser());
         std::string deck_filename = param.get<std::string>("deck_filename");
+
+        eclipseState.reset( new EclipseState(parser->parseFile(deck_filename)));
         deck.reset(new EclipseGridParser(deck_filename));
         // Grid init
         grid.reset(new GridManager(*deck));
@@ -242,7 +250,7 @@ try
                       << simtimer.numSteps() - step << ")\n\n" << std::flush;
 
             // Create new wells, well_state
-            WellsManager wells(*deck, *grid->c_grid(), props->permeability());
+            WellsManager wells(eclipseState , epoch , *grid->c_grid(), props->permeability());
             // @@@ HACK: we should really make a new well state and
             // properly transfer old well state to it every epoch,
             // since number of wells may change etc.

examples/sim_2p_incomp.cpp

@@ -45,6 +45,10 @@
 #include <opm/core/simulator/WellState.hpp>
 #include <opm/core/simulator/SimulatorIncompTwophase.hpp>
 
+#include <opm/parser/eclipse/Parser/Parser.hpp>
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+
+
 #include <memory>
 #include <boost/filesystem.hpp>
 
@@ -94,6 +98,8 @@ try
 
     // If we have a "deck_filename", grid and props will be read from that.
     bool use_deck = param.has("deck_filename");
+    EclipseStateConstPtr eclipseState;
+
     std::unique_ptr<EclipseGridParser> deck;
     std::unique_ptr<GridManager> grid;
     std::unique_ptr<IncompPropertiesInterface> props;
@@ -103,8 +109,10 @@ try
     // int max_well_control_iterations = 0;
     double gravity[3] = { 0.0 };
     if (use_deck) {
+        ParserPtr parser(new Opm::Parser());
         std::string deck_filename = param.get<std::string>("deck_filename");
         deck.reset(new EclipseGridParser(deck_filename));
+        eclipseState.reset( new EclipseState(parser->parseFile(deck_filename)));
         // Grid init
         grid.reset(new GridManager(*deck));
         // Rock and fluid init
@@ -260,7 +268,7 @@ try
                       << simtimer.numSteps() - step << ")\n\n" << std::flush;
 
             // Create new wells, well_state
-            WellsManager wells(*deck, *grid->c_grid(), props->permeability());
+            WellsManager wells(eclipseState , epoch , *grid->c_grid(), props->permeability());
             // @@@ HACK: we should really make a new well state and
             // properly transfer old well state to it every epoch,
             // since number of wells may change etc.

examples/wells_example.cpp

@@ -19,6 +19,9 @@
 #include <opm/core/linalg/LinearSolverFactory.hpp>
 #include <opm/core/props/rock/RockCompressibility.hpp>
 
+#include <opm/parser/eclipse/Parser/Parser.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
+
 int main(int argc, char** argv)
 try
 {
@@ -39,9 +42,11 @@ try
     // Define rock and fluid properties
     IncompPropertiesFromDeck incomp_properties(parser, *grid.c_grid());
     RockCompressibility rock_comp(parser);
+    ParserPtr newParser(new Opm::Parser());
+    EclipseStateConstPtr eclipseState(new Opm::EclipseState(newParser->parseFile(file_name)));
 
     // Finally handle the wells
-    WellsManager wells(parser, *grid.c_grid(), incomp_properties.permeability());
+    WellsManager wells(eclipseState , 0 , *grid.c_grid(), incomp_properties.permeability());
 
     double gravity[3] = {0.0, 0.0, parameters.getDefault<double>("gravity", 0.0)};
     Opm::LinearSolverFactory linsolver(parameters);

opm/core/grid/GridManager.hpp

@@ -44,10 +44,10 @@ namespace Opm
     {
     public:
         /// Construct a 3d corner-point grid or tensor grid from a deck.
-        GridManager(const Opm::EclipseGridParser& deck);
+        explicit GridManager(const Opm::EclipseGridParser& deck);
 
         /// Construct a 3d corner-point grid or tensor grid from a deck.
-        GridManager(Opm::DeckConstPtr newParserDeck);
+        explicit GridManager(Opm::DeckConstPtr newParserDeck);
 
         /// Construct a 2d cartesian grid with cells of unit size.
         GridManager(int nx, int ny);
@@ -65,7 +65,7 @@ namespace Opm
         /// Construct a grid from an input file.
         /// The file format used is currently undocumented,
         /// and is therefore only suited for internal use.
-        GridManager(const std::string& input_filename);
+        explicit GridManager(const std::string& input_filename);
 
         /// Destructor.
         ~GridManager();

opm/core/io/OutputWriter.cpp

@@ -48,7 +48,7 @@ private:
 /// Psuedo-constructor, can appear in template
 template <typename Format> unique_ptr <OutputWriter>
 create (const ParameterGroup& params,
-        std::shared_ptr <const EclipseGridParser> parser,
+        std::shared_ptr <EclipseGridParser> parser,
         std::shared_ptr <const UnstructuredGrid> grid) {
     return unique_ptr <OutputWriter> (new Format (params, parser, grid));
 }
@@ -61,7 +61,7 @@ create (const ParameterGroup& params,
 /// to the list below!
 typedef map <const char*, unique_ptr <OutputWriter> (*)(
         const ParameterGroup&,
-        std::shared_ptr <const EclipseGridParser>,
+        std::shared_ptr <EclipseGridParser>,
         std::shared_ptr <const UnstructuredGrid>)> map_t;
 map_t FORMATS = {
     { "output_ecl", &create <EclipseWriter> },
@@ -71,7 +71,7 @@ map_t FORMATS = {
 
 unique_ptr <OutputWriter>
 OutputWriter::create (const ParameterGroup& params,
-                      std::shared_ptr <const EclipseGridParser> parser,
+                      std::shared_ptr <EclipseGridParser> parser,
                       std::shared_ptr <const UnstructuredGrid> grid) {
     // allocate a list which will be filled with writers. this list
     // is initially empty (no output).

opm/core/io/OutputWriter.hpp

@@ -108,7 +108,7 @@ public:
      */
     static std::unique_ptr <OutputWriter>
     create (const parameter::ParameterGroup& params,
-            std::shared_ptr <const EclipseGridParser> parser,
+            std::shared_ptr <EclipseGridParser> parser,
             std::shared_ptr <const UnstructuredGrid> grid);
 };
 

opm/core/io/eclipse/EclipseGridParser.cpp

@@ -126,6 +126,7 @@ namespace EclipseKeywords
           string("EQUIL"),    string("PVCDO"),    string("TSTEP"),
           string("PLYVISC"),  string("PLYROCK"),  string("PLYADS"),
           string("PLYMAX"),   string("TLMIXPAR"), string("WPOLYMER"),
+          string("PLYSHEAR"),
           string("GRUPTREE"), string("GCONINJE"), string("GCONPROD"),
           string("WGRUPCON"), string("ENDSCALE"),  string("SCALECRS"),
           string("ENPTVD"),   string("ENKRVD"),

opm/core/io/eclipse/EclipseGridParser.hpp

@@ -195,6 +195,7 @@ public:
     SPECIAL_FIELD(PLYROCK)
     SPECIAL_FIELD(PLYADS)
     SPECIAL_FIELD(PLYMAX)
+    SPECIAL_FIELD(PLYSHEAR)
     SPECIAL_FIELD(TLMIXPAR)
     SPECIAL_FIELD(WPOLYMER)
     SPECIAL_FIELD(GRUPTREE)

opm/core/io/eclipse/EclipseWriter.hpp

@@ -24,7 +24,10 @@
 #include <opm/core/io/OutputWriter.hpp>
 #include <opm/core/props/BlackoilPhases.hpp>
 
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+
 #include <string>
+#include <vector>
 #include <memory>  // std::unique_ptr
 
 struct UnstructuredGrid;
@@ -57,7 +60,7 @@ public:
      * binary files using ERT.
      */
     EclipseWriter(const parameter::ParameterGroup& params,
-                  std::shared_ptr <const EclipseGridParser> parser,
+                  std::shared_ptr <EclipseGridParser> parser,
                   std::shared_ptr <const UnstructuredGrid> grid);
 
     EclipseWriter (const parameter::ParameterGroup& params,
@@ -65,6 +68,25 @@ public:
                    int number_of_cells, const int* global_cell, const int* cart_dims,
                    int dimension);
 
+    /*!
+     * \brief Sets the common attributes required to write eclipse
+     *        binary files using ERT.
+     */
+    EclipseWriter(const parameter::ParameterGroup& params,
+                  Opm::DeckConstPtr newParserDeck,
+                  std::shared_ptr <const UnstructuredGrid> grid);
+
+    /*!
+     * \brief Sets the common attributes required to write eclipse
+     *        binary files using ERT.
+     */
+    EclipseWriter(const parameter::ParameterGroup& params,
+                  Opm::DeckConstPtr newParserDeck,
+                  int number_of_cells, const int* global_cell, const int* cart_dims,
+                  int dimension);
+
+
+                 
     /**
      * We need a destructor in the compilation unit to avoid the
      * EclipseSummary being a complete type here.
@@ -92,19 +114,26 @@ public:
 
 private:
     std::shared_ptr <const EclipseGridParser> parser_;
+    Opm::DeckConstPtr newParserDeck_;
     int number_of_cells_;
     int dimensions_;
     const int* cart_dims_;
     const int* global_cell_;
+    bool enableOutput_;
+    int outputInterval_;
+    int outputTimeStepIdx_;
     std::string outputDir_;
     std::string baseName_;
     PhaseUsage uses_;           // active phases in the input deck
     std::shared_ptr <EclipseSummary> summary_;
 
+    void activeToGlobalCellData_(std::vector<double> &globalCellsBuf,
+                                 const std::vector<double> &activeCellsBuf,
+                                 const std::vector<double> &inactiveCellsBuf) const;
+
     /// Write solution field variables (pressure and saturation)
-    void writeSolution (const SimulatorTimer& timer,
-                        const SimulatorState& reservoirState,
-                        const WellState& wellState);
+    void writeSolution_(const SimulatorTimer& timer,
+                        const SimulatorState& reservoirState);
 
     void init(const parameter::ParameterGroup& params);
 };

opm/core/io/eclipse/SpecialEclipseFields.hpp

@@ -2094,6 +2094,42 @@ struct PLYMAX : public SpecialBase
 };
 
 
+struct PLYSHEAR : public SpecialBase
+{
+    std::vector<double> water_velocity_;
+    std::vector<double> vrf_;
+
+    virtual std::string name() const {return std::string("PLYSHEAR");}
+
+    virtual void read(std::istream& is)
+    {
+        // Note. This function assumes that NTSFUN = 1, and reads only one table.
+        std::vector<double> plyshear;
+        readVectorData(is, plyshear);
+        for (int i=0; i<(int)plyshear.size(); i+=2) {
+            water_velocity_.push_back(plyshear[i]);
+            vrf_.push_back(plyshear[i+1]);
+        }
+    }
+
+    virtual void write(std::ostream& os) const
+    {
+        os << name() << '\n';
+        for (int i=0; i<(int)water_velocity_.size(); ++i) {
+            os << water_velocity_[i] << " " << vrf_[i] << '\n';
+        }
+        os << '\n';
+    }
+
+    virtual void convertToSI(const EclipseUnits& units)
+    {
+        for (int i=0; i<(int)water_velocity_.size(); ++i) {
+            water_velocity_[i] *= units.length / units.time;
+        }
+    }
+};
+
+
 struct TLMIXPAR : public SpecialBase
 {
     std::vector<double> tlmixpar_;

opm/core/linalg/LinearSolverIstl.cpp

@@ -30,6 +30,7 @@
 
 // TODO: clean up includes.
 #include <dune/common/deprecated.hh>
+#include <dune/common/version.hh>
 #include <dune/istl/bvector.hh>
 #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/operators.hh>
@@ -39,6 +40,10 @@
 #include <dune/istl/paamg/amg.hh>
 #include <dune/istl/paamg/kamg.hh>
 
+#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
+#include <dune/istl/paamg/fastamg.hh>
+#endif
+
 #include <stdexcept>
 #include <iostream>
 
@@ -61,7 +66,7 @@ namespace Opm
         solveCG_AMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity,
                     double prolongateFactor, int smoothsteps);
 
-#ifdef HAS_DUNE_FAST_AMG
+#if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
         LinearSolverInterface::LinearSolverReport
         solveKAMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity,
                   double prolongateFactor, int smoothsteps);
@@ -173,7 +178,7 @@ namespace Opm
                               linsolver_prolongate_factor_, linsolver_smooth_steps_);
             break;
         case KAMG:
-#ifdef HAS_DUNE_FAST_AMG
+#if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
             res = solveKAMG(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_,
                             linsolver_prolongate_factor_, linsolver_smooth_steps_);
 #else
@@ -181,7 +186,7 @@ namespace Opm
 #endif
             break;
         case FastAMG:
-#ifdef HAS_DUNE_FAST_AMG
+#if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
             res = solveFastAMG(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_,
                                linsolver_prolongate_factor_);
 #else
@@ -311,7 +316,7 @@ namespace Opm
     }
 
 
-#ifdef HAS_DUNE_FAST_AMG
+#if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
     LinearSolverInterface::LinearSolverReport
     solveKAMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity,
               double linsolver_prolongate_factor, int linsolver_smooth_steps)

opm/core/props/pvt/BlackoilPvtProperties.cpp

@@ -125,16 +125,17 @@ namespace Opm
         // Surface densities. Accounting for different orders in eclipse and our code.
         if (newParserDeck->hasKeyword("DENSITY")) {
             Opm::DeckKeywordConstPtr densityKeyword = newParserDeck->getKeyword("DENSITY");
-            const std::vector<double>& d = densityKeyword->getRecord(region_number_)->getItem(0)->getSIDoubleData();
-            enum { ECL_oil = 0, ECL_water = 1, ECL_gas = 2 };
+            if (phase_usage_.phase_used[Liquid]) {
+                densities_[phase_usage_.phase_pos[Liquid]]
+                    = densityKeyword->getRecord(region_number_)->getItem("OIL")->getSIDouble(0);
+            }
             if (phase_usage_.phase_used[Aqua]) {
-                densities_[phase_usage_.phase_pos[Aqua]]   = d[ECL_water];
+                densities_[phase_usage_.phase_pos[Aqua]]
+                    = densityKeyword->getRecord(region_number_)->getItem("WATER")->getSIDouble(0);
             }
             if (phase_usage_.phase_used[Vapour]) {
-                densities_[phase_usage_.phase_pos[Vapour]] = d[ECL_gas];
-            }
-            if (phase_usage_.phase_used[Liquid]) {
-                densities_[phase_usage_.phase_pos[Liquid]] = d[ECL_oil];
+                densities_[phase_usage_.phase_pos[Vapour]]
+                    = densityKeyword->getRecord(region_number_)->getItem("GAS")->getSIDouble(0);
             }
         } else {
             OPM_THROW(std::runtime_error, "Input is missing DENSITY\n");

opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp

@@ -240,7 +240,7 @@ namespace Opm
                           "SaturationPropsFromDeck::init()   --  ENDSCALE: "
                           "Currently only 'NODIR' accepted.");
             }
-            if (endscale.isReversible()) {
+            if (!endscale.isReversible()) {
                 OPM_THROW(std::runtime_error,
                           "SaturationPropsFromDeck::init()   --  ENDSCALE: "
                           "Currently only 'REVERS' accepted.");

opm/core/simulator/SimulatorCompressibleTwophase.cpp

@@ -498,13 +498,13 @@ namespace Opm
                       << std::endl;
             std::cout.precision(8);
 
-            watercut.push(timer.currentTime() + timer.currentStepLength(),
+            watercut.push(timer.simulationTimeElapsed() + timer.currentStepLength(),
                           produced[0]/(produced[0] + produced[1]),
                           tot_produced[0]/tot_porevol_init);
             if (wells_) {
                 wellreport.push(props_, *wells_,
                                 state.pressure(), state.surfacevol(), state.saturation(),
-                                timer.currentTime() + timer.currentStepLength(),
+                                timer.simulationTimeElapsed() + timer.currentStepLength(),
                                 well_state.bhp(), well_state.perfRates());
             }
             sreport.total_time =  step_timer.secsSinceStart();

opm/core/simulator/SimulatorIncompTwophase.cpp

@@ -578,12 +578,12 @@ namespace Opm
                     dynamic_cast<TransportSolverTwophaseReorder&>(*tsolver_)
                         .solveGravity(&initial_porevol[0], stepsize, state);
                 }
-                watercut.push(timer.currentTime() + timer.currentStepLength(),
+                watercut.push(timer.simulationTimeElapsed() + timer.currentStepLength(),
                               produced[0]/(produced[0] + produced[1]),
                               tot_produced[0]/tot_porevol_init);
                 if (wells_) {
                     wellreport.push(props_, *wells_, state.saturation(),
-                                    timer.currentTime() + timer.currentStepLength(),
+                                    timer.simulationTimeElapsed() + timer.currentStepLength(),
                                     well_state.bhp(), well_state.perfRates());
                 }
             }

opm/core/simulator/SimulatorOutput.cpp

@@ -30,7 +30,7 @@ using namespace Opm;
 
 SimulatorOutputBase::SimulatorOutputBase (
         const parameter::ParameterGroup& params,
-        std::shared_ptr <const EclipseGridParser> parser,
+        std::shared_ptr <EclipseGridParser> parser,
         std::shared_ptr <const UnstructuredGrid> grid,
         std::shared_ptr <const SimulatorTimer> timer,
         std::shared_ptr <const SimulatorState> state,
@@ -72,7 +72,7 @@ SimulatorOutputBase::operator std::function <void ()> () {
 
 void
 SimulatorOutputBase::writeOutput () {
-    const int this_time = timer_->currentTime ();
+    const int this_time = timer_->simulationTimeElapsed ();
 
     // if the simulator signals for timesteps that aren't reporting
     // times, then ignore them

opm/core/simulator/SimulatorOutput.hpp

@@ -53,7 +53,7 @@ protected:
      * need to pick them up from the object members.
      */
     SimulatorOutputBase (const parameter::ParameterGroup& p,
-                         std::shared_ptr <const EclipseGridParser> parser,
+                         std::shared_ptr <EclipseGridParser> parser,
                          std::shared_ptr <const UnstructuredGrid> grid,
                          std::shared_ptr <const SimulatorTimer> timer,
                          std::shared_ptr <const SimulatorState> state,
@@ -139,7 +139,7 @@ private:
 template <typename Simulator>
 struct SimulatorOutput : public SimulatorOutputBase {
 	SimulatorOutput (const parameter::ParameterGroup& params,
-                     std::shared_ptr <const EclipseGridParser> parser,
+                     std::shared_ptr <EclipseGridParser> parser,
                      std::shared_ptr <const UnstructuredGrid> grid,
                      std::shared_ptr <const SimulatorTimer> timer,
                      std::shared_ptr <const SimulatorState> state,
@@ -161,7 +161,7 @@ struct SimulatorOutput : public SimulatorOutputBase {
      * the arguments passed exceeds the lifetime of this object.
      */
     SimulatorOutput (const parameter::ParameterGroup& params,
-                     const EclipseGridParser& parser,
+                     EclipseGridParser& parser,
                      const UnstructuredGrid& grid,
                      const SimulatorTimer& timer,
                      const SimulatorState& state,

opm/core/simulator/SimulatorTimer.cpp

@@ -59,25 +59,50 @@ namespace Opm
 
     /// Use the SimulatorTimer as a shim around opm-parser's Opm::TimeMap
     void SimulatorTimer::init(Opm::TimeMapConstPtr timeMap,
-                              int timeStepIdx)
+                              size_t beginReportStepIdx,
+                              size_t endReportStepIdx)
     {
         timeMap_ = timeMap;
-        current_step_ = timeStepIdx;
+        current_step_ = 0;
+        beginReportStepIdx_ = beginReportStepIdx;
+        endReportStepIdx_ = std::min(timeMap_->numTimesteps() + 1, endReportStepIdx);
     }
 
     /// Total number of steps.
     int SimulatorTimer::numSteps() const
     {
         if (timeMap_)
-            return timeMap_->numTimesteps();
+            return endReportStepIdx_ - beginReportStepIdx_;
         else
             return timesteps_.size();
     }
 
+    /// Index of the first considered simulation episode
+    size_t SimulatorTimer::beginReportStepIndex() const
+    {
+        if (!timeMap_) {
+            OPM_THROW(std::runtime_error, "indexFirstEpisode() is only implemented "
+                      "for simulation timers which are based on Opm::TimeMap");
+        }
+
+        return beginReportStepIdx_;
+    }
+
+    /// Index of the last considered simulation episode
+    size_t SimulatorTimer::endReportStepIndex() const
+    {
+        if (!timeMap_) {
+            OPM_THROW(std::runtime_error, "indexLastEpisode() is only implemented "
+                      "for simulation timers which are based on Opm::TimeMap");
+        }
+
+        return endReportStepIdx_;
+    }
+
     /// Current step number.
     int SimulatorTimer::currentStepNum() const
     {
-        return current_step_;
+        return current_step_ + beginReportStepIdx_;
     }
 
     /// Set current step number.
@@ -99,7 +124,7 @@ namespace Opm
     {
         assert(!done());
         if (timeMap_)
-            return timeMap_->getTimeStepLength(current_step_);
+            return timeMap_->getTimeStepLength(beginReportStepIdx_ + current_step_);
         else
             return timesteps_[current_step_];
     }
@@ -108,25 +133,32 @@ namespace Opm
     {
         assert(current_step_ > 0);
         if (timeMap_)
-            return timeMap_->getTimeStepLength(current_step_ - 1);
+            return timeMap_->getTimeStepLength(beginReportStepIdx_ + current_step_ - 1);
         else
             return timesteps_[current_step_ - 1];
     }
 
-    /// Current time.
-    double SimulatorTimer::currentTime() const
+    /// time elapsed since the start of the simulation [s].
+    double SimulatorTimer::simulationTimeElapsed() const
     {
         if (timeMap_)
-            return timeMap_->getTimePassedUntil(current_step_);
+            return
+                timeMap_->getTimePassedUntil(beginReportStepIdx_ + current_step_);
         else
             return current_time_;
     }
 
+    /// time elapsed since the start of the POSIX epoch (Jan 1st, 1970) [s].
+    time_t SimulatorTimer::currentPosixTime() const
+    {
+        tm t = boost::posix_time::to_tm(currentDateTime());
+        return std::mktime(&t);
+    }
 
     boost::posix_time::ptime SimulatorTimer::currentDateTime() const
     {
         if (timeMap_)
-            return timeMap_->getStartTime(current_step_);
+            return timeMap_->getStartTime(beginReportStepIdx_ + current_step_);
         else
             return boost::posix_time::ptime(start_date_) + boost::posix_time::seconds( (int) current_time_ );
     }
@@ -137,7 +169,8 @@ namespace Opm
     double SimulatorTimer::totalTime() const
     {
         if (timeMap_)
-            return timeMap_->getTotalTime();
+            return
+                timeMap_->getTotalTime();
         else
             return total_time_;
     }
@@ -161,7 +194,7 @@ namespace Opm
     void SimulatorTimer::report(std::ostream& os) const
     {
         os << "\n\n---------------    Simulation step number " << currentStepNum() << "    ---------------"
-           << "\n      Current time (days)     " << Opm::unit::convert::to(currentTime(), Opm::unit::day)
+           << "\n      Current time (days)     " << Opm::unit::convert::to(simulationTimeElapsed(), Opm::unit::day)
            << "\n      Current stepsize (days) " << Opm::unit::convert::to(currentStepLength(), Opm::unit::day)
            << "\n      Total time (days)       " << Opm::unit::convert::to(totalTime(), Opm::unit::day)
            << "\n" << std::endl;
@@ -181,7 +214,7 @@ namespace Opm
     bool SimulatorTimer::done() const
     {
         if (timeMap_)
-            return int(timeMap_->numTimesteps()) <= current_step_;
+            return current_step_ > int(endReportStepIdx_ - beginReportStepIdx_ - 1);
         else
             return int(timesteps_.size()) == current_step_;
     }

opm/core/simulator/SimulatorTimer.hpp

@@ -51,11 +51,18 @@ namespace Opm
 
         /// Use the SimulatorTimer as a shim around opm-parser's Opm::TimeMap
         void init(TimeMapConstPtr timeMap,
-                  int timeStepIdx = 0);
+                  size_t beginReportStepIdx = 0,
+                  size_t endReportStepIdx = std::numeric_limits<size_t>::max());
 
         /// Total number of steps.
         int numSteps() const;
 
+        /// Index of the first report step considered
+        size_t beginReportStepIndex() const;
+
+        /// Index of the next-after-last report step to be considered
+        size_t endReportStepIndex() const;
+
         /// Current step number. This is the number of timesteps that
         /// has been completed from the start of the run. The time
         /// after initialization but before the simulation has started
@@ -79,8 +86,13 @@ namespace Opm
         /// it is an error to call stepLengthTaken().
         double stepLengthTaken () const;
 
-        /// Current time.
-        double currentTime() const;
+        /// Time elapsed since the start of the POSIX epoch (Jan 1st,
+        /// 1970) until the current time step begins [s].
+        time_t currentPosixTime() const;
+
+        /// Time elapsed since the start of the simulation until the
+        /// beginning of the current time step [s].
+        double simulationTimeElapsed() const;
 
         /// Return the current time as a posix time object.
         boost::posix_time::ptime currentDateTime() const;
@@ -107,6 +119,8 @@ namespace Opm
     private:
         Opm::TimeMapConstPtr timeMap_;
         std::vector<double> timesteps_;
+        size_t beginReportStepIdx_;
+        size_t endReportStepIdx_;
         int current_step_;
         double current_time_;
         double total_time_;

opm/core/simulator/WellState.hpp

@@ -23,6 +23,7 @@
 #include <opm/core/wells.h>
 #include <opm/core/well_controls.h>
 #include <vector>
+#include <cassert>
 
 namespace Opm
 {
@@ -45,26 +46,55 @@ namespace Opm
                 bhp_.resize(nw);
                 wellrates_.resize(nw * np, 0.0);
                 for (int w = 0; w < nw; ++w) {
+                    assert((wells->type[w] == INJECTOR) || (wells->type[w] == PRODUCER));
                     const WellControls* ctrl = wells->ctrls[w];
-                    // Initialize bhp to be target pressure if
-                    // bhp-controlled well, otherwise set to a little
-                    // above or below (depending on if the well is an
-                    // injector or producer) pressure in first perforation
-                    // cell.
-                    if (well_controls_well_is_shut(ctrl) || (well_controls_get_current_type(ctrl) != BHP)) {
-                        const int first_cell = wells->well_cells[wells->well_connpos[w]];
-                        const double safety_factor = (wells->type[w] == INJECTOR) ? 1.01 : 0.99;
-                        bhp_[w] = safety_factor*state.pressure()[first_cell];
-                    } else {
-                        bhp_[w] = well_controls_get_current_target( ctrl );
-                    }
-
-                    // Initialize well rates to match controls if type is SURFACE_RATE
-                    if (well_controls_well_is_open( ctrl ) || (well_controls_get_current_type(ctrl) == SURFACE_RATE)) {
-                        const double rate_target = well_controls_get_current_target(ctrl);
-                        const double * distr = well_controls_get_current_distr( ctrl );
+                    if (well_controls_well_is_shut(ctrl)) {
+                        // Shut well:
+                        // 1. Assign zero well rates.
                         for (int p = 0; p < np; ++p) {
-                            wellrates_[np*w + p] = rate_target * distr[p];
+                            wellrates_[np*w + p] = 0.0;
+                        }
+                        // 2. Assign bhp equal to bhp control, if
+                        //    applicable, otherwise assign equal to
+                        //    first perforation cell pressure.
+                        if (well_controls_get_current_type(ctrl) == BHP) {
+                            bhp_[w] = well_controls_get_current_target( ctrl );
+                        } else {
+                            const int first_cell = wells->well_cells[wells->well_connpos[w]];
+                            bhp_[w] = state.pressure()[first_cell];
+                        }
+                    } else {
+                        // Open well:
+                        // 1. Initialize well rates to match controls
+                        //    if type is SURFACE_RATE.  Otherwise, we
+                        //    cannot set the correct value here, so we
+                        //    assign a small rate with the correct
+                        //    sign so that any logic depending on that
+                        //    sign will work as expected.
+                        if (well_controls_get_current_type(ctrl) == SURFACE_RATE) {
+                            const double rate_target = well_controls_get_current_target(ctrl);
+                            const double * distr = well_controls_get_current_distr( ctrl );
+                            for (int p = 0; p < np; ++p) {
+                                wellrates_[np*w + p] = rate_target * distr[p];
+                            }
+                        } else {
+                            const double small_rate = 1e-14;
+                            const double sign = (wells->type[w] == INJECTOR) ? 1.0 : -1.0;
+                            for (int p = 0; p < np; ++p) {
+                                wellrates_[np*w + p] = small_rate * sign;
+                            }
+                        }
+                        // 2. Initialize bhp to be target pressure if
+                        //    bhp-controlled well, otherwise set to a
+                        //    little above or below (depending on if
+                        //    the well is an injector or producer)
+                        //    pressure in first perforation cell.
+                        if (well_controls_get_current_type(ctrl) == BHP) {
+                            bhp_[w] = well_controls_get_current_target( ctrl );
+                        } else {
+                            const int first_cell = wells->well_cells[wells->well_connpos[w]];
+                            const double safety_factor = (wells->type[w] == INJECTOR) ? 1.01 : 0.99;
+                            bhp_[w] = safety_factor*state.pressure()[first_cell];
                         }
                     }
                 }

opm/core/simulator/initState_impl.hpp

@@ -584,6 +584,11 @@ namespace Opm
                   "found " << pu.num_phases << " phases in deck.");
         }
         state.init(number_of_cells, number_of_faces, num_phases);
+        if (newParserDeck->hasKeyword("EQUIL") && newParserDeck->hasKeyword("PRESSURE")) {
+            OPM_THROW(std::runtime_error, "initStateFromDeck(): The deck must either specify the initial "
+                      "condition using the PRESSURE _or_ the EQUIL keyword (currently it has both)");
+        }
+
         if (newParserDeck->hasKeyword("EQUIL")) {
             if (num_phases != 2) {
                 OPM_THROW(std::runtime_error, "initStateFromDeck(): EQUIL-based init currently handling only two-phase scenarios.");
@@ -703,6 +708,10 @@ namespace Opm
             OPM_THROW(std::runtime_error, "initStateFromDeck():  user specified property object with " << num_phases << " phases, "
                   "found " << pu.num_phases << " phases in deck.");
         }
+        if (deck.hasField("EQUIL") && deck.hasField("PRESSURE")) {
+            OPM_THROW(std::runtime_error, "initStateFromDeck(): The deck must either specify the initial "
+                      "condition using the PRESSURE _or_ the EQUIL keyword (currently it has both)");
+        }
         state.init(number_of_cells, number_of_faces, num_phases);
         if (deck.hasField("EQUIL")) {
             if (num_phases != 2) {

opm/core/wells/WellsGroup.cpp

@@ -1171,19 +1171,21 @@ namespace Opm
         InjectionSpecification injection_specification;
         ProductionSpecification production_specification;
         if (well->isInjector(timeStep)) {
-            injection_specification.BHP_limit_ = well->getBHPLimit(timeStep);
-            injection_specification.injector_type_ = toInjectorType(WellInjector::Type2String(well->getInjectorType(timeStep)));
-            injection_specification.control_mode_ = toInjectionControlMode(WellInjector::ControlMode2String(well->getInjectorControlMode(timeStep)));
-            injection_specification.surface_flow_max_rate_ = well->getSurfaceInjectionRate(timeStep);
-            injection_specification.reservoir_flow_max_rate_ = well->getReservoirInjectionRate(timeStep);
+            const WellInjectionProperties& properties = well->getInjectionProperties(timeStep);
+            injection_specification.BHP_limit_ = properties.BHPLimit;
+            injection_specification.injector_type_ = toInjectorType(WellInjector::Type2String(properties.injectorType));
+            injection_specification.control_mode_ = toInjectionControlMode(WellInjector::ControlMode2String(properties.controlMode));
+            injection_specification.surface_flow_max_rate_ = properties.surfaceInjectionRate;
+            injection_specification.reservoir_flow_max_rate_ = properties.reservoirInjectionRate;
             production_specification.guide_rate_ = 0.0; // We know we're not a producer
         }
         else if (well->isProducer(timeStep)) {
-            production_specification.BHP_limit_ = well->getBHPLimit(timeStep);
-            production_specification.reservoir_flow_max_rate_ = well->getResVRate(timeStep);
-            production_specification.oil_max_rate_ = well->getOilRate(timeStep);
-            production_specification.control_mode_ = toProductionControlMode(WellProducer::ControlMode2String(well->getProducerControlMode(timeStep)));
-            production_specification.water_max_rate_ = well->getWaterRate(timeStep);
+            const WellProductionProperties& properties = well->getProductionProperties(timeStep);
+            production_specification.BHP_limit_ = properties.BHPLimit;
+            production_specification.reservoir_flow_max_rate_ = properties.ResVRate;
+            production_specification.oil_max_rate_ = properties.OilRate;
+            production_specification.control_mode_ = toProductionControlMode(WellProducer::ControlMode2String(properties.controlMode));
+            production_specification.water_max_rate_ = properties.WaterRate;
             injection_specification.guide_rate_ = 0.0; // we know we're not an injector
         }
         std::shared_ptr<WellsGroupInterface> wells_group(new WellNode(well->name(), production_specification, injection_specification, phase_usage));

opm/core/wells/WellsManager.cpp

@@ -231,113 +231,20 @@ namespace Opm
     {
     }
 
-
-    /// Construct from existing wells object.
-WellsManager::WellsManager(struct Wells* W, bool checkCellExistence)
-    : w_(clone_wells(W)), checkCellExistence_(checkCellExistence)
-    {
-    }
-
     /// Construct wells from deck.
     WellsManager::WellsManager(const Opm::EclipseStateConstPtr eclipseState,
                                const size_t timeStep,
                                const UnstructuredGrid& grid,
-                               const double* permeability,
-                               bool checkCellExistence)
-        : w_(0), checkCellExistence_(checkCellExistence)
+                               const double* permeability)
+        : w_(0)
     {
-        if (UgGridHelpers::dimensions(grid) != 3) {
-            OPM_THROW(std::runtime_error, "We cannot initialize wells from a deck unless the corresponding grid is 3-dimensional.");
-        }
-
-        if (eclipseState->getSchedule()->numWells() == 0) {
-            OPM_MESSAGE("No wells specified in Schedule section, initializing no wells");
-            return;
-        }
-
-        std::map<int,int> cartesian_to_compressed;
-        setupCompressedToCartesian(UgGridHelpers::globalCell(grid), 
-                                   UgGridHelpers::numCells(grid), cartesian_to_compressed);
-
-        // Obtain phase usage data.
-        PhaseUsage pu = phaseUsageFromDeck(eclipseState);
-
-        // These data structures will be filled in this constructor,
-        // then used to initialize the Wells struct.
-        std::vector<std::string> well_names;
-        std::vector<WellData> well_data;
-
-
-        // For easy lookup:
-        std::map<std::string, int> well_names_to_index;
-
-        ScheduleConstPtr schedule = eclipseState->getSchedule();
-        std::vector<WellConstPtr> wells = schedule->getWells(timeStep);
-
-        well_names.reserve(wells.size());
-        well_data.reserve(wells.size());
-
-        createWellsFromSpecs(wells, timeStep, UgGridHelpers::cell2Faces(grid),
-                             UgGridHelpers::cartDims(grid),
-                             UgGridHelpers::beginFaceCentroids(grid),
-                             UgGridHelpers::beginCellCentroids(grid),
-                             UgGridHelpers::dimensions(grid),
-                             well_names, well_data, well_names_to_index, pu, cartesian_to_compressed, permeability);
-        setupWellControls(wells, timeStep, well_names, pu);
-
-        {
-            GroupTreeNodeConstPtr fieldNode = eclipseState->getSchedule()->getGroupTree(timeStep)->getNode("FIELD");
-            GroupConstPtr fieldGroup = eclipseState->getSchedule()->getGroup(fieldNode->name());
-            well_collection_.addField(fieldGroup, timeStep, pu);
-            addChildGroups(fieldNode, eclipseState->getSchedule(), timeStep, pu);
-        }
-
-        for (auto wellIter = wells.begin(); wellIter != wells.end(); ++wellIter ) {
-            well_collection_.addWell((*wellIter), timeStep, pu);
-        }
-
-        well_collection_.setWellsPointer(w_);
-        well_collection_.applyGroupControls();
-
-
-        setupGuideRates(wells, timeStep, well_data, well_names_to_index);
-
-        // Debug output.
-#define EXTRA_OUTPUT
-#ifdef EXTRA_OUTPUT
-        /*
-        std::cout << "\t WELL DATA" << std::endl;
-        for(int i = 0; i< num_wells; ++i) {
-            std::cout << i << ": " << well_data[i].type << "  "
-                      << well_data[i].control << "  " << well_data[i].target
-                      << std::endl;
-        }
-
-        std::cout << "\n\t PERF DATA" << std::endl;
-        for(int i=0; i< int(wellperf_data.size()); ++i) {
-            for(int j=0; j< int(wellperf_data[i].size()); ++j) {
-                std::cout << i << ": " << wellperf_data[i][j].cell << "  "
-                          << wellperf_data[i][j].well_index << std::endl;
-            }
-        }
-        */
-#endif
-    }
-
-
-
-    /// Construct wells from deck.
-    WellsManager::WellsManager(const Opm::EclipseGridParser& deck,
-                               const UnstructuredGrid& grid,
-                               const double* permeability,
-                               bool checkCellExistence)
-        : w_(0), checkCellExistence_(checkCellExistence)
-    {
-        init(deck, grid.number_of_cells, grid.global_cell, grid.cartdims, grid.dimensions,
-             grid.cell_centroids, UgGridHelpers::cell2Faces(grid), grid.face_centroids,
+        init(eclipseState, timeStep, UgGridHelpers::numCells(grid),
+             UgGridHelpers::globalCell(grid), UgGridHelpers::cartDims(grid), 
+             UgGridHelpers::dimensions(grid), UgGridHelpers::beginCellCentroids(grid),
+             UgGridHelpers::cell2Faces(grid), UgGridHelpers::beginFaceCentroids(grid),
              permeability);
-    }
 
+    }
 
     /// Destructor.
     WellsManager::~WellsManager()
@@ -423,14 +330,15 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence)
             }
 
             if (well->isInjector(timeStep)) {
-                clear_well_controls(well_index, w_);
+                const WellInjectionProperties& injectionProperties = well->getInjectionProperties(timeStep);
                 int ok = 1;
                 int control_pos[5] = { -1, -1, -1, -1, -1 };
-
-                if (well->hasInjectionControl(timeStep , WellInjector::RATE)) {
+                                
+                clear_well_controls(well_index, w_);
+                if (injectionProperties.hasInjectionControl(WellInjector::RATE)) {
                     control_pos[WellsManagerDetail::InjectionControl::RATE] = well_controls_get_num(w_->ctrls[well_index]);
                     double distr[3] = { 0.0, 0.0, 0.0 };
-                    WellInjector::TypeEnum injectorType = well->getInjectorType(timeStep);
+                    WellInjector::TypeEnum injectorType = injectionProperties.injectorType;
 
                     if (injectorType == WellInjector::TypeEnum::WATER) {
                         distr[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0;
@@ -441,16 +349,16 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence)
                     }
 
                     ok = append_well_controls(SURFACE_RATE,
-                                              well->getSurfaceInjectionRate( timeStep ) ,
+                                              injectionProperties.surfaceInjectionRate, 
                                               distr,
                                               well_index,
                                               w_);
                 }
 
-                if (ok && well->hasInjectionControl(timeStep , WellInjector::RESV)) {
+                if (ok && injectionProperties.hasInjectionControl(WellInjector::RESV)) {
                     control_pos[WellsManagerDetail::InjectionControl::RESV] = well_controls_get_num(w_->ctrls[well_index]);
                     double distr[3] = { 0.0, 0.0, 0.0 };
-                    WellInjector::TypeEnum injectorType = well->getInjectorType(timeStep);
+                    WellInjector::TypeEnum injectorType = injectionProperties.injectorType;
 
                     if (injectorType == WellInjector::TypeEnum::WATER) {
                         distr[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0;
@@ -461,22 +369,23 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence)
                     }
 
                     ok = append_well_controls(RESERVOIR_RATE,
-                                              well->getReservoirInjectionRate( timeStep ),
+                                              injectionProperties.reservoirInjectionRate, 
                                               distr,
                                               well_index,
                                               w_);
                 }
 
-                if (ok && well->hasInjectionControl(timeStep , WellInjector::BHP)) {
+                if (ok && injectionProperties.hasInjectionControl(WellInjector::BHP)) {
+                    control_pos[WellsManagerDetail::InjectionControl::BHP] = well_controls_get_num(w_->ctrls[well_index]);
                     control_pos[WellsManagerDetail::InjectionControl::BHP] = well_controls_get_num(w_->ctrls[well_index]);
                     ok = append_well_controls(BHP,
-                                              well->getBHPLimit(timeStep),
+                                              injectionProperties.BHPLimit, 
                                               NULL,
                                               well_index,
                                               w_);
                 }
 
-                if (ok && well->hasInjectionControl(timeStep , WellInjector::THP)) {
+                if (ok && injectionProperties.hasInjectionControl(WellInjector::THP)) {
                     OPM_THROW(std::runtime_error, "We cannot handle THP limit for well " << well_names[well_index]);
                 }
 
@@ -487,7 +396,7 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence)
 
 
                 {
-                    WellsManagerDetail::InjectionControl::Mode mode = WellsManagerDetail::InjectionControl::mode( well->getInjectorControlMode(timeStep) );
+                    WellsManagerDetail::InjectionControl::Mode mode = WellsManagerDetail::InjectionControl::mode( injectionProperties.controlMode );
                     int cpos = control_pos[mode];
                     if (cpos == -1 && mode != WellsManagerDetail::InjectionControl::GRUP) {
                         OPM_THROW(std::runtime_error, "Control not specified in well " << well_names[well_index]);
@@ -495,7 +404,7 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence)
 
                     // We need to check if the well is shut or not
                     if (well->getStatus( timeStep ) == WellCommon::SHUT) {
-                        cpos = ~cpos;
+                        well_controls_shut_well( w_->ctrls[well_index] );
                     }
                     set_current_control(well_index, cpos, w_);
                 }
@@ -503,34 +412,39 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence)
 
                 // Set well component fraction.
                 double cf[3] = { 0.0, 0.0, 0.0 };
-                if (well->getInjectorType(timeStep) == WellInjector::WATER) {
-                    if (!phaseUsage.phase_used[BlackoilPhases::Aqua]) {
-                        OPM_THROW(std::runtime_error, "Water phase not used, yet found water-injecting well.");
-                    }
-                    cf[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0;
-                } else if (well->getInjectorType(timeStep) == WellInjector::OIL) {
-                    if (!phaseUsage.phase_used[BlackoilPhases::Liquid]) {
-                        OPM_THROW(std::runtime_error, "Oil phase not used, yet found oil-injecting well.");
-                    }
-                    cf[phaseUsage.phase_pos[BlackoilPhases::Liquid]] = 1.0;
-                } else if (well->getInjectorType(timeStep) == WellInjector::GAS) {
-                    if (!phaseUsage.phase_used[BlackoilPhases::Vapour]) {
-                        OPM_THROW(std::runtime_error, "Gas phase not used, yet found gas-injecting well.");
-                    }
-                    cf[phaseUsage.phase_pos[BlackoilPhases::Vapour]] = 1.0;
-                }
-                std::copy(cf, cf + phaseUsage.num_phases, w_->comp_frac + well_index*phaseUsage.num_phases);
+                {
+                    WellInjector::TypeEnum injectorType = injectionProperties.injectorType;
 
+                    if (injectorType == WellInjector::WATER) {
+                        if (!phaseUsage.phase_used[BlackoilPhases::Aqua]) {
+                            OPM_THROW(std::runtime_error, "Water phase not used, yet found water-injecting well.");
+                        }
+                        cf[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0;
+                    } else if (injectorType == WellInjector::OIL) {
+                        if (!phaseUsage.phase_used[BlackoilPhases::Liquid]) {
+                            OPM_THROW(std::runtime_error, "Oil phase not used, yet found oil-injecting well.");
+                        }
+                        cf[phaseUsage.phase_pos[BlackoilPhases::Liquid]] = 1.0;
+                    } else if (injectorType == WellInjector::GAS) {
+                        if (!phaseUsage.phase_used[BlackoilPhases::Vapour]) {
+                            OPM_THROW(std::runtime_error, "Gas phase not used, yet found gas-injecting well.");
+                        }
+                        cf[phaseUsage.phase_pos[BlackoilPhases::Vapour]] = 1.0;
+                    }
+                    std::copy(cf, cf + phaseUsage.num_phases, w_->comp_frac + well_index*phaseUsage.num_phases);
+                }
             }
 
             if (well->isProducer(timeStep)) {
                 // Add all controls that are present in well.
                 // First we must clear existing controls, in case the
                 // current WCONPROD line is modifying earlier controls.
-                clear_well_controls(well_index, w_);
+                const WellProductionProperties& productionProperties = well->getProductionProperties(timeStep);
                 int control_pos[9] = { -1, -1, -1, -1, -1, -1, -1, -1, -1 };
                 int ok = 1;
-                if (ok && well->hasProductionControl(timeStep , WellProducer::ORAT)) {
+                
+                clear_well_controls(well_index, w_);
+                if (ok && productionProperties.hasProductionControl(WellProducer::ORAT)) {
                     if (!phaseUsage.phase_used[BlackoilPhases::Liquid]) {
                         OPM_THROW(std::runtime_error, "Oil phase not active and ORAT control specified.");
                     }
@@ -539,13 +453,13 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence)
                     double distr[3] = { 0.0, 0.0, 0.0 };
                     distr[phaseUsage.phase_pos[BlackoilPhases::Liquid]] = 1.0;
                     ok = append_well_controls(SURFACE_RATE,
-                                              -well->getOilRate( timeStep ),
+                                              -productionProperties.OilRate, 
                                               distr,
                                                well_index,
                                               w_);
                 }
 
-                if (ok && well->hasProductionControl(timeStep , WellProducer::WRAT)) {
+                if (ok && productionProperties.hasProductionControl(WellProducer::WRAT)) {
                     if (!phaseUsage.phase_used[BlackoilPhases::Aqua]) {
                         OPM_THROW(std::runtime_error, "Water phase not active and WRAT control specified.");
                     }
@@ -553,13 +467,13 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence)
                     double distr[3] = { 0.0, 0.0, 0.0 };
                     distr[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0;
                     ok = append_well_controls(SURFACE_RATE,
-                                              -well->getWaterRate(timeStep),
+                                              -productionProperties.WaterRate, 
                                               distr,
                                               well_index,
                                               w_);
                 }
 
-                if (ok && well->hasProductionControl(timeStep , WellProducer::GRAT)) {
+                if (ok && productionProperties.hasProductionControl(WellProducer::GRAT)) {
                     if (!phaseUsage.phase_used[BlackoilPhases::Vapour]) {
                         OPM_THROW(std::runtime_error, "Gas phase not active and GRAT control specified.");
                     }
@@ -567,13 +481,13 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence)
                     double distr[3] = { 0.0, 0.0, 0.0 };
                     distr[phaseUsage.phase_pos[BlackoilPhases::Vapour]] = 1.0;
                     ok = append_well_controls(SURFACE_RATE,
-                                              -well->getGasRate( timeStep ),
+                                              -productionProperties.GasRate,
                                               distr,
                                               well_index,
                                               w_);
                 }
 
-                if (ok && well->hasProductionControl(timeStep , WellProducer::LRAT)) {
+                if (ok && productionProperties.hasProductionControl(WellProducer::LRAT)) {
                     if (!phaseUsage.phase_used[BlackoilPhases::Aqua]) {
                         OPM_THROW(std::runtime_error, "Water phase not active and LRAT control specified.");
                     }
@@ -585,32 +499,32 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence)
                     distr[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0;
                     distr[phaseUsage.phase_pos[BlackoilPhases::Liquid]] = 1.0;
                     ok = append_well_controls(SURFACE_RATE,
-                                              -well->getLiquidRate(timeStep),
+                                              -productionProperties.LiquidRate , 
                                               distr,
                                               well_index,
                                               w_);
                 }
 
-                if (ok && well->hasProductionControl(timeStep , WellProducer::RESV)) {
+                if (ok && productionProperties.hasProductionControl(WellProducer::RESV)) {
                     control_pos[WellsManagerDetail::ProductionControl::RESV] = well_controls_get_num(w_->ctrls[well_index]);
                     double distr[3] = { 1.0, 1.0, 1.0 };
                     ok = append_well_controls(RESERVOIR_RATE,
-                                              -well->getResVRate(timeStep),
+                                              -productionProperties.ResVRate , 
                                               distr,
                                               well_index,
                                               w_);
                 }
 
-                if (ok && well->hasProductionControl(timeStep , WellProducer::BHP)) {
+                if (ok && productionProperties.hasProductionControl(WellProducer::BHP)) {
                     control_pos[WellsManagerDetail::ProductionControl::BHP] = well_controls_get_num(w_->ctrls[well_index]);
                     ok = append_well_controls(BHP,
-                                              well->getBHPLimit( timeStep ) ,
+                                              productionProperties.BHPLimit , 
                                               NULL,
                                               well_index,
                                               w_);
                 }
 
-                if (ok && well->hasProductionControl(timeStep , WellProducer::THP)) {
+                if (ok && productionProperties.hasProductionControl(WellProducer::THP)) {
                     OPM_THROW(std::runtime_error, "We cannot handle THP limit for well " << well_names[well_index]);
                 }
 
@@ -618,14 +532,16 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence)
                     OPM_THROW(std::runtime_error, "Failure occured appending controls for well " << well_names[well_index]);
                 }
 
-                WellsManagerDetail::ProductionControl::Mode mode = WellsManagerDetail::ProductionControl::mode(well->getProducerControlMode(timeStep));
+                WellsManagerDetail::ProductionControl::Mode mode = WellsManagerDetail::ProductionControl::mode(productionProperties.controlMode);
                 int cpos = control_pos[mode];
                 if (cpos == -1 && mode != WellsManagerDetail::ProductionControl::GRUP) {
                     OPM_THROW(std::runtime_error, "Control mode type " << mode << " not present in well " << well_names[well_index]);
                 }
                 // If it's shut, we complement the cpos
                 if (well->getStatus(timeStep) == WellCommon::SHUT) {
-                    cpos = ~cpos; // So we can easily retrieve the cpos later
+                    well_controls_shut_well( w_->ctrls[well_index] );
+                } else if (cpos == -1 && mode != WellsManagerDetail::ProductionControl::GRUP) {
+                    OPM_THROW(std::runtime_error, "Control mode type " << mode << " not present in well " << well_names[well_index]);
                 }
                 set_current_control(well_index, cpos, w_);
             }

opm/core/wells/WellsManager.hpp

@@ -66,20 +66,15 @@ namespace Opm
         /// manage control switching does not exist.
         ///
         /// @param[in] W Existing wells object.
-        WellsManager(struct Wells* W, bool checkCellExistence=true);
+        explicit WellsManager(struct Wells* W, bool checkCellExistence=true);
 
         /// Construct from input deck and grid.
         /// The permeability argument may be zero if the input contain
         /// well productivity indices, otherwise it must be given in
         /// order to approximate these by the Peaceman formula.
-        WellsManager(const Opm::EclipseGridParser& deck,
-                     const UnstructuredGrid& grid,
-                     const double* permeability,
-                     bool checkCellExistence=true);
-
-
         template<class CC, class F2C, class FC>
-        WellsManager(const Opm::EclipseGridParser& deck,
+        WellsManager(const Opm::EclipseStateConstPtr eclipseState,
+                     const size_t timeStep,
                      int num_cells,
                      const int* global_cell,
                      const int* cart_dims,
@@ -87,15 +82,12 @@ namespace Opm
                      CC begin_cell_centroids,
                      const F2C& f2c,
                      FC begin_face_centroids,
-                     const double* permeability,
-                     bool checkCellExistence=true);
+                     const double* permeability);
 
         WellsManager(const Opm::EclipseStateConstPtr eclipseState,
                      const size_t timeStep,
                      const UnstructuredGrid& grid,
-                     const double* permeability,
-                     bool checkCellExistence=true);
-
+                     const double* permeability);
         /// Destructor.
         ~WellsManager();
 
@@ -148,25 +140,26 @@ namespace Opm
 
 
     private:
-        template<class CC, class F2C, class FC>
-        void init(const Opm::EclipseGridParser& deck,
+        template<class CC, class C2F, class FC>
+        void init(const Opm::EclipseStateConstPtr eclipseState,
+                  const size_t timeStep,
                   int num_cells,
                   const int* global_cell,
                   const int* cart_dims,
                   int dimensions,
                   CC begin_cell_centroids,
-                  const F2C& f2c,
+                  const C2F& cell_to_faces,
                   FC begin_face_centroids,
                   const double* permeability);
         // Disable copying and assignment.
-        WellsManager(const WellsManager& other, bool checkCellExistence=true);
+        WellsManager(const WellsManager& other);
         WellsManager& operator=(const WellsManager& other);
         static void setupCompressedToCartesian(const int* global_cell, int number_of_cells, std::map<int,int>& cartesian_to_compressed );
         void setupWellControls(std::vector<WellConstPtr>& wells, size_t timeStep,
                                std::vector<std::string>& well_names, const PhaseUsage& phaseUsage);
         template<class C2F, class CC, class FC>
         void createWellsFromSpecs( std::vector<WellConstPtr>& wells, size_t timeStep,
-                                   const C2F& c2f, 
+                                   const C2F& cell_to_faces, 
                                    const int* cart_dims,
                                    FC begin_face_centroids, 
                                    CC begin_cell_centroids,
@@ -185,7 +178,6 @@ namespace Opm
         // Data
         Wells* w_;
         WellCollection well_collection_;
-        bool checkCellExistence_;
     };
 
 } // namespace Opm

opm/core/wells/WellsManager_impl.hpp

@@ -127,10 +127,8 @@ void WellsManager::createWellsFromSpecs(std::vector<WellConstPtr>& wells, size_t
                 int cart_grid_indx = i + cpgdim[0]*(j + cpgdim[1]*k);
                 std::map<int, int>::const_iterator cgit = cartesian_to_compressed.find(cart_grid_indx);
                 if (cgit == cartesian_to_compressed.end()) {
-                    if (checkCellExistence_)
                         OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << i << ' ' << j << ' '
                                   << k << " not found in grid (well = " << well->name() << ')');
-                    continue;
                 }
                 int cell = cgit->second;
                 PerfData pd;
@@ -206,522 +204,90 @@ void WellsManager::createWellsFromSpecs(std::vector<WellConstPtr>& wells, size_t
 }
 
 template<class CC, class C2F, class FC>
-WellsManager::WellsManager(const Opm::EclipseGridParser& deck,
+WellsManager::WellsManager(const Opm::EclipseStateConstPtr eclipseState,
+                           const size_t timeStep,
                            int number_of_cells,
                            const int* global_cell,
                            const int* cart_dims,
                            int dimensions,
                            CC begin_cell_centroids,
-                           const C2F& c2f,
+                           const C2F& cell_to_faces,
                            FC begin_face_centroids,
-                           const double* permeability,
-                           bool checkCellExistence)
-    : w_(0), checkCellExistence_(checkCellExistence)
+                           const double* permeability)
+    : w_(0)
 {
-    init(deck, number_of_cells, global_cell, cart_dims, dimensions,
-         begin_cell_centroids, c2f, begin_face_centroids, permeability);
+    init(eclipseState, timeStep, number_of_cells, global_cell, cart_dims, dimensions,
+         begin_cell_centroids, cell_to_faces, begin_face_centroids, permeability);
 }
 
-    /// Construct wells from deck.
+/// Construct wells from deck.
 template<class CC, class C2F, class FC>
-void WellsManager::init(const Opm::EclipseGridParser& deck,
+void WellsManager::init(const Opm::EclipseStateConstPtr eclipseState,
+                        const size_t timeStep,
                         int number_of_cells,
                         const int* global_cell,
                         const int* cart_dims,
                         int dimensions,
                         CC begin_cell_centroids,
-                        const C2F& c2f,
+                        const C2F& cell_to_faces,
                         FC begin_face_centroids,
                         const double* permeability)
 {
     if (dimensions != 3) {
         OPM_THROW(std::runtime_error, "We cannot initialize wells from a deck unless the corresponding grid is 3-dimensional.");
     }
-    // NOTE: Implementation copied and modified from dune-porsol's class BlackoilWells.
-    std::vector<std::string> keywords;
-    keywords.push_back("WELSPECS");
-    keywords.push_back("COMPDAT");
-    //      keywords.push_back("WELTARG");
-    if (!deck.hasFields(keywords)) {
-        OPM_MESSAGE("Missing well keywords in deck, initializing no wells.");
+
+    if (eclipseState->getSchedule()->numWells() == 0) {
+        OPM_MESSAGE("No wells specified in Schedule section, initializing no wells");
         return;
     }
-    if (!(deck.hasField("WCONINJE") || deck.hasField("WCONPROD")) ) {
-        OPM_THROW(std::runtime_error, "Needed field is missing in file");
-    }
+
+    std::map<int,int> cartesian_to_compressed;
+    setupCompressedToCartesian(global_cell, 
+                               number_of_cells, cartesian_to_compressed);
 
     // Obtain phase usage data.
-    PhaseUsage pu = phaseUsageFromDeck(deck);
+    PhaseUsage pu = phaseUsageFromDeck(eclipseState);
 
     // These data structures will be filled in this constructor,
     // then used to initialize the Wells struct.
     std::vector<std::string> well_names;
     std::vector<WellData> well_data;
-    std::vector<std::vector<PerfData> > wellperf_data;
+
 
     // For easy lookup:
     std::map<std::string, int> well_names_to_index;
-    typedef std::map<std::string, int>::const_iterator WNameIt;
 
-    // Get WELSPECS data.
-    // It is allowed to have multiple lines corresponding to
-    // the same well, in which case the last one encountered
-    // is the one used.
-    const WELSPECS& welspecs = deck.getWELSPECS();
-    const int num_welspecs = welspecs.welspecs.size();
-    well_names.reserve(num_welspecs);
-    well_data.reserve(num_welspecs);
-    for (int w = 0; w < num_welspecs; ++w) {
-        // First check if this well has already been encountered.
-        // If so, we modify it's data instead of appending a new well
-        // to the well_data and well_names vectors.
-        const std::string& name = welspecs.welspecs[w].name_;
-        const double refdepth = welspecs.welspecs[w].datum_depth_BHP_;
-        WNameIt wit = well_names_to_index.find(name);
-        if (wit == well_names_to_index.end()) {
-            // New well, append data.
-            well_names_to_index[welspecs.welspecs[w].name_] = well_data.size();
-            well_names.push_back(name);
-            WellData wd;
-            // If negative (defaulted), set refdepth to a marker
-            // value, will be changed after getting perforation
-            // data to the centroid of the cell of the top well
-            // perforation.
-            wd.reference_bhp_depth = (refdepth < 0.0) ? -1e100 : refdepth;
-            wd.welspecsline = w;
-            well_data.push_back(wd);
-        } else {
-            // Existing well, change data.
-            const int wix = wit->second;
-            well_data[wix].reference_bhp_depth = (refdepth < 0.0) ? -1e100 : refdepth;
-            well_data[wix].welspecsline = w;
-        }
-    }
-    const int num_wells = well_data.size();
-    wellperf_data.resize(num_wells);
+    ScheduleConstPtr schedule = eclipseState->getSchedule();
+    std::vector<WellConstPtr> wells = schedule->getWells(timeStep);
 
+    well_names.reserve(wells.size());
+    well_data.reserve(wells.size());
 
-    // global_cell is a map from compressed cells to Cartesian grid cells.
-    // We must make the inverse lookup.
-    const int* cpgdim = cart_dims;
-    std::map<int,int> cartesian_to_compressed;
+    createWellsFromSpecs(wells, timeStep, cell_to_faces,
+                         cart_dims,
+                         begin_face_centroids,
+                         begin_cell_centroids,
+                         dimensions,
+                         well_names, well_data, well_names_to_index, pu, cartesian_to_compressed, permeability);
+    setupWellControls(wells, timeStep, well_names, pu);
 
-    if (global_cell) {
-        for (int i = 0; i < number_of_cells; ++i) {
-            cartesian_to_compressed.insert(std::make_pair(global_cell[i], i));
-        }
-    }
-    else {
-        for (int i = 0; i < number_of_cells; ++i) {
-            cartesian_to_compressed.insert(std::make_pair(i, i));
-        }
+    {
+        GroupTreeNodeConstPtr fieldNode = eclipseState->getSchedule()->getGroupTree(timeStep)->getNode("FIELD");
+        GroupConstPtr fieldGroup = eclipseState->getSchedule()->getGroup(fieldNode->name());
+        well_collection_.addField(fieldGroup, timeStep, pu);
+        addChildGroups(fieldNode, eclipseState->getSchedule(), timeStep, pu);
     }
 
-    // Get COMPDAT data
-    // It is *not* allowed to have multiple lines corresponding to
-    // the same perforation!
-    const COMPDAT& compdat = deck.getCOMPDAT();
-    const int num_compdat  = compdat.compdat.size();
-    for (int kw = 0; kw < num_compdat; ++kw) {
-        // Extract well name, or the part of the well name that
-        // comes before the '*'.
-        std::string name = compdat.compdat[kw].well_;
-        std::string::size_type len = name.find('*');
-        if (len != std::string::npos) {
-            name = name.substr(0, len);
-        }
-        // Look for well with matching name.
-        bool found = false;
-        for (int wix = 0; wix < num_wells; ++wix) {
-            if (well_names[wix].compare(0,len, name) == 0) { // equal
-                // Extract corresponding WELSPECS defintion for
-                // purpose of default location specification.
-                const WelspecsLine& wspec = welspecs.welspecs[well_data[wix].welspecsline];
-
-                // We have a matching name.
-                int ix = compdat.compdat[kw].grid_ind_[0] - 1;
-                int jy = compdat.compdat[kw].grid_ind_[1] - 1;
-                int kz1 = compdat.compdat[kw].grid_ind_[2] - 1;
-                int kz2 = compdat.compdat[kw].grid_ind_[3] - 1;
-
-                if (ix < 0) {
-                    // Defaulted I location.  Extract from WELSPECS.
-                    ix = wspec.I_ - 1;
-                }
-                if (jy < 0) {
-                    // Defaulted J location.  Extract from WELSPECS.
-                    jy = wspec.J_ - 1;
-                }
-                if (kz1 < 0) {
-                    // Defaulted KZ1.  Use top layer.
-                    kz1 = 0;
-                }
-                if (kz2 < 0) {
-                    // Defaulted KZ2.  Use bottom layer.
-                    kz2 = cpgdim[2] - 1;
-                }
-
-                for (int kz = kz1; kz <= kz2; ++kz) {
-                    int cart_grid_indx = ix + cpgdim[0]*(jy + cpgdim[1]*kz);
-                    std::map<int, int>::const_iterator cgit =
-                        cartesian_to_compressed.find(cart_grid_indx);
-                    if (cgit == cartesian_to_compressed.end()) {
-                        if(checkCellExistence_)
-                            OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << ix << ' ' << jy << ' '
-                                      << kz << " not found in grid (well = " << name << ')');
-                        continue;
-                    }
-                    int cell = cgit->second;
-                    PerfData pd;
-                    pd.cell = cell;
-                    if (compdat.compdat[kw].connect_trans_fac_ > 0.0) {
-                        pd.well_index = compdat.compdat[kw].connect_trans_fac_;
-                    } else {
-                        double radius = 0.5*compdat.compdat[kw].diameter_;
-                        if (radius <= 0.0) {
-                            radius = 0.5*unit::feet;
-                            OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius);
-                        }
-                        std::array<double, 3> cubical = WellsManagerDetail::getCubeDim(c2f,
-                                                                                       begin_face_centroids,
-                                                                                       dimensions,
-                                                                                       cell);
-                        const double* cell_perm = &permeability[dimensions*dimensions*cell];
-                        pd.well_index = WellsManagerDetail::computeWellIndex(radius, cubical, cell_perm,
-                                                                             compdat.compdat[kw].skin_factor_);
-                    }
-                    wellperf_data[wix].push_back(pd);
-                }
-                found = true;
-                break;
-            }
-        }
-        if (!found) {
-            OPM_THROW(std::runtime_error, "Undefined well name: " << compdat.compdat[kw].well_
-                      << " in COMPDAT");
-        }
+    for (auto wellIter = wells.begin(); wellIter != wells.end(); ++wellIter ) {
+        well_collection_.addWell((*wellIter), timeStep, pu);
     }
 
-    // Set up reference depths that were defaulted. Count perfs.
-    int num_perfs = 0;
-    assert(dimensions == 3);
-    for (int w = 0; w < num_wells; ++w) {
-        num_perfs += wellperf_data[w].size();
-        if (well_data[w].reference_bhp_depth < 0.0) {
-            // It was defaulted. Set reference depth to minimum perforation depth.
-            double min_depth = 1e100;
-            int num_wperfs = wellperf_data[w].size();
-            for (int perf = 0; perf < num_wperfs; ++perf) {
-                double depth = UgGridHelpers::
-                    getCoordinate(UgGridHelpers::increment(begin_cell_centroids, 
-                                                           wellperf_data[w][perf].cell,
-                                                           dimensions), 2);
-                min_depth = std::min(min_depth, depth);
-            }
-            well_data[w].reference_bhp_depth = min_depth;
-        }
-    }
+    well_collection_.setWellsPointer(w_);
+    well_collection_.applyGroupControls();
 
-    // Create the well data structures.
-    w_ = create_wells(pu.num_phases, num_wells, num_perfs);
-    if (!w_) {
-        OPM_THROW(std::runtime_error, "Failed creating Wells struct.");
-    }
 
-    // Classify wells
-    if (deck.hasField("WCONINJE")) {
-        const std::vector<WconinjeLine>& lines = deck.getWCONINJE().wconinje;
-        for (size_t i = 0 ; i < lines.size(); ++i) {
-            const std::map<std::string, int>::const_iterator it = well_names_to_index.find(lines[i].well_);
-            if (it != well_names_to_index.end()) {
-                const int well_index = it->second;
-                well_data[well_index].type = INJECTOR;
-            } else {
-                OPM_THROW(std::runtime_error, "Unseen well name: " << lines[i].well_ << " first seen in WCONINJE");
-            }
-        }
-    }
-    if (deck.hasField("WCONPROD")) {
-        const std::vector<WconprodLine>& lines = deck.getWCONPROD().wconprod;
-        for (size_t i = 0; i < lines.size(); ++i) {
-            const std::map<std::string, int>::const_iterator it = well_names_to_index.find(lines[i].well_);
-            if (it != well_names_to_index.end()) {
-                const int well_index = it->second;
-                well_data[well_index].type = PRODUCER;
-            } else {
-                OPM_THROW(std::runtime_error, "Unseen well name: " << lines[i].well_ << " first seen in WCONPROD");
-            }
-
-        }
-    }
-
-    // Add wells.
-    for (int w = 0; w < num_wells; ++w) {
-        const int w_num_perf = wellperf_data[w].size();
-        std::vector<int> perf_cells(w_num_perf);
-        std::vector<double> perf_prodind(w_num_perf);
-        for (int perf = 0; perf < w_num_perf; ++perf) {
-            perf_cells[perf] = wellperf_data[w][perf].cell;
-            perf_prodind[perf] = wellperf_data[w][perf].well_index;
-        }
-        const double* comp_frac = NULL;
-        // We initialize all wells with a null component fraction,
-        // and must (for injection wells) overwrite it later.
-        int ok = add_well(well_data[w].type, well_data[w].reference_bhp_depth, w_num_perf,
-                          comp_frac, &perf_cells[0], &perf_prodind[0], well_names[w].c_str(), w_);
-        if (!ok) {
-            OPM_THROW(std::runtime_error, "Failed adding well " << well_names[w] << " to Wells data structure.");
-        }
-    }
-
-    // Get WCONINJE data, add injection controls to wells.
-    // It is allowed to have multiple lines corresponding to
-    // the same well, in which case the last one encountered
-    // is the one used.
-    if (deck.hasField("WCONINJE")) {
-        const WCONINJE& wconinjes = deck.getWCONINJE();
-        const int num_wconinjes = wconinjes.wconinje.size();
-        for (int kw = 0; kw < num_wconinjes; ++kw) {
-            const WconinjeLine& wci_line = wconinjes.wconinje[kw];
-            // Extract well name, or the part of the well name that
-            // comes before the '*'.
-            std::string name = wci_line.well_;
-            std::string::size_type len = name.find('*');
-            if (len != std::string::npos) {
-                name = name.substr(0, len);
-            }
-            bool well_found = false;
-            for (int wix = 0; wix < num_wells; ++wix) {
-                if (well_names[wix].compare(0,len, name) == 0) { //equal
-                    well_found = true;
-                    assert(well_data[wix].type == w_->type[wix]);
-                    if (well_data[wix].type != INJECTOR) {
-                        OPM_THROW(std::runtime_error, "Found WCONINJE entry for a non-injector well: " << well_names[wix]);
-                    }
-
-                    // Add all controls that are present in well.
-                    // First we must clear existing controls, in case the
-                    // current WCONINJE line is modifying earlier controls.
-                    clear_well_controls(wix, w_);
-                    int ok = 1;
-                    int control_pos[5] = { -1, -1, -1, -1, -1 };
-                    if (ok && wci_line.surface_flow_max_rate_ >= 0.0) {
-                        control_pos[WellsManagerDetail::InjectionControl::RATE] = well_controls_get_num(w_->ctrls[wix]);
-                        double distr[3] = { 0.0, 0.0, 0.0 };
-                        if (wci_line.injector_type_ == "WATER") {
-                            distr[pu.phase_pos[BlackoilPhases::Aqua]] = 1.0;
-                        } else if (wci_line.injector_type_ == "OIL") {
-                            distr[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
-                        } else if (wci_line.injector_type_ == "GAS") {
-                            distr[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
-                        } else {
-                            OPM_THROW(std::runtime_error, "Injector type " << wci_line.injector_type_ << " not supported."
-                                      "WellsManager only supports WATER, OIL and GAS injector types.");
-                        }
-                        ok = append_well_controls(SURFACE_RATE, wci_line.surface_flow_max_rate_,
-                                                  distr, wix, w_);
-                    }
-                    if (ok && wci_line.reservoir_flow_max_rate_ >= 0.0) {
-                        control_pos[WellsManagerDetail::InjectionControl::RESV] = well_controls_get_num(w_->ctrls[wix]);
-                        double distr[3] = { 0.0, 0.0, 0.0 };
-                        if (wci_line.injector_type_ == "WATER") {
-                            distr[pu.phase_pos[BlackoilPhases::Aqua]] = 1.0;
-                        } else if (wci_line.injector_type_ == "OIL") {
-                            distr[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
-                        } else if (wci_line.injector_type_ == "GAS") {
-                            distr[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
-                        } else {
-                            OPM_THROW(std::runtime_error, "Injector type " << wci_line.injector_type_ << " not supported."
-                                      "WellsManager only supports WATER, OIL and GAS injector types.");
-                        }
-                        ok = append_well_controls(RESERVOIR_RATE, wci_line.reservoir_flow_max_rate_,
-                                                  distr, wix, w_);
-                    }
-                    if (ok && wci_line.BHP_limit_ > 0.0) {
-                        control_pos[WellsManagerDetail::InjectionControl::BHP] = well_controls_get_num(w_->ctrls[wix]);
-                        ok = append_well_controls(BHP, wci_line.BHP_limit_,
-                                                  NULL, wix, w_);
-                    }
-                    if (ok && wci_line.THP_limit_ > 0.0) {
-                        OPM_THROW(std::runtime_error, "We cannot handle THP limit for well " << well_names[wix]);
-                    }
-                    if (!ok) {
-                        OPM_THROW(std::runtime_error, "Failure occured appending controls for well " << well_names[wix]);
-                    }
-                    WellsManagerDetail::InjectionControl::Mode mode = WellsManagerDetail::InjectionControl::mode(wci_line.control_mode_);
-                    int cpos = control_pos[mode];
-                    if (cpos == -1 && mode != WellsManagerDetail::InjectionControl::GRUP) {
-                        OPM_THROW(std::runtime_error, "Control for " << wci_line.control_mode_ << " not specified in well " << well_names[wix]);
-                    }
-                    // We need to check if the well is shut or not
-                    if (wci_line.open_shut_flag_ == "SHUT") {
-                        cpos = ~cpos;
-                    }
-                    set_current_control(wix, cpos, w_);
-
-                    // Set well component fraction.
-                    double cf[3] = { 0.0, 0.0, 0.0 };
-                    if (wci_line.injector_type_[0] == 'W') {
-                        if (!pu.phase_used[BlackoilPhases::Aqua]) {
-                            OPM_THROW(std::runtime_error, "Water phase not used, yet found water-injecting well.");
-                        }
-                        cf[pu.phase_pos[BlackoilPhases::Aqua]] = 1.0;
-                    } else if (wci_line.injector_type_[0] == 'O') {
-                        if (!pu.phase_used[BlackoilPhases::Liquid]) {
-                            OPM_THROW(std::runtime_error, "Oil phase not used, yet found oil-injecting well.");
-                        }
-                        cf[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
-                    } else if (wci_line.injector_type_[0] == 'G') {
-                        if (!pu.phase_used[BlackoilPhases::Vapour]) {
-                            OPM_THROW(std::runtime_error, "Gas phase not used, yet found gas-injecting well.");
-                        }
-                        cf[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
-                    }
-                    std::copy(cf, cf + pu.num_phases, w_->comp_frac + wix*pu.num_phases);
-                }
-            }
-            if (!well_found) {
-                OPM_THROW(std::runtime_error, "Undefined well name: " << wci_line.well_
-                          << " in WCONINJE");
-            }
-        }
-    }
-
-    // Get WCONPROD data
-    // It is allowed to have multiple lines corresponding to
-    // the same well, in which case the last one encountered
-    // is the one used.
-    if (deck.hasField("WCONPROD")) {
-        const WCONPROD& wconprods = deck.getWCONPROD();
-        const int num_wconprods   = wconprods.wconprod.size();
-        for (int kw = 0; kw < num_wconprods; ++kw) {
-            const WconprodLine& wcp_line = wconprods.wconprod[kw];
-            std::string name = wcp_line.well_;
-            std::string::size_type len = name.find('*');
-            if (len != std::string::npos) {
-                name = name.substr(0, len);
-            }
-            bool well_found = false;
-            for (int wix = 0; wix < num_wells; ++wix) {
-                if (well_names[wix].compare(0,len, name) == 0) { //equal
-                    well_found = true;
-                    assert(well_data[wix].type == w_->type[wix]);
-                    if (well_data[wix].type != PRODUCER) {
-                        OPM_THROW(std::runtime_error, "Found WCONPROD entry for a non-producer well: " << well_names[wix]);
-                    }
-                    // Add all controls that are present in well.
-                    // First we must clear existing controls, in case the
-                    // current WCONPROD line is modifying earlier controls.
-                    clear_well_controls(wix, w_);
-                    int control_pos[9] = { -1, -1, -1, -1, -1, -1, -1, -1, -1 };
-                    int ok = 1;
-                    if (ok && wcp_line.oil_max_rate_ >= 0.0) {
-                        if (!pu.phase_used[BlackoilPhases::Liquid]) {
-                            OPM_THROW(std::runtime_error, "Oil phase not active and ORAT control specified.");
-                        }
-                        control_pos[WellsManagerDetail::ProductionControl::ORAT] = well_controls_get_num(w_->ctrls[wix]);
-                        double distr[3] = { 0.0, 0.0, 0.0 };
-                        distr[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
-                        ok = append_well_controls(SURFACE_RATE, -wcp_line.oil_max_rate_,
-                                                  distr, wix, w_);
-                    }
-                    if (ok && wcp_line.water_max_rate_ >= 0.0) {
-                        if (!pu.phase_used[BlackoilPhases::Aqua]) {
-                            OPM_THROW(std::runtime_error, "Water phase not active and WRAT control specified.");
-                        }
-                        control_pos[WellsManagerDetail::ProductionControl::WRAT] = well_controls_get_num(w_->ctrls[wix]);
-                        double distr[3] = { 0.0, 0.0, 0.0 };
-                        distr[pu.phase_pos[BlackoilPhases::Aqua]] = 1.0;
-                        ok = append_well_controls(SURFACE_RATE, -wcp_line.water_max_rate_,
-                                                  distr, wix, w_);
-                    }
-                    if (ok && wcp_line.gas_max_rate_ >= 0.0) {
-                        if (!pu.phase_used[BlackoilPhases::Vapour]) {
-                            OPM_THROW(std::runtime_error, "Gas phase not active and GRAT control specified.");
-                        }
-                        control_pos[WellsManagerDetail::ProductionControl::GRAT] = well_controls_get_num(w_->ctrls[wix]);
-                        double distr[3] = { 0.0, 0.0, 0.0 };
-                        distr[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
-                        ok = append_well_controls(SURFACE_RATE, -wcp_line.gas_max_rate_,
-                                                  distr, wix, w_);
-                    }
-                    if (ok && wcp_line.liquid_max_rate_ >= 0.0) {
-                        if (!pu.phase_used[BlackoilPhases::Aqua]) {
-                            OPM_THROW(std::runtime_error, "Water phase not active and LRAT control specified.");
-                        }
-                        if (!pu.phase_used[BlackoilPhases::Liquid]) {
-                            OPM_THROW(std::runtime_error, "Oil phase not active and LRAT control specified.");
-                        }
-                        control_pos[WellsManagerDetail::ProductionControl::LRAT] = well_controls_get_num(w_->ctrls[wix]);
-                        double distr[3] = { 0.0, 0.0, 0.0 };
-                        distr[pu.phase_pos[BlackoilPhases::Aqua]] = 1.0;
-                        distr[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
-                        ok = append_well_controls(SURFACE_RATE, -wcp_line.liquid_max_rate_,
-                                                  distr, wix, w_);
-                    }
-                    if (ok && wcp_line.reservoir_flow_max_rate_ >= 0.0) {
-                        control_pos[WellsManagerDetail::ProductionControl::RESV] = well_controls_get_num(w_->ctrls[wix]);
-                        double distr[3] = { 1.0, 1.0, 1.0 };
-                        ok = append_well_controls(RESERVOIR_RATE, -wcp_line.reservoir_flow_max_rate_,
-                                                  distr, wix, w_);
-                    }
-                    if (ok && wcp_line.BHP_limit_ > 0.0) {
-                        control_pos[WellsManagerDetail::ProductionControl::BHP] = well_controls_get_num(w_->ctrls[wix]);
-                        ok = append_well_controls(BHP, wcp_line.BHP_limit_,
-                                                  NULL, wix, w_);
-                    }
-                    if (ok && wcp_line.THP_limit_ > 0.0) {
-                        OPM_THROW(std::runtime_error, "We cannot handle THP limit for well " << well_names[wix]);
-                    }
-                    if (!ok) {
-                        OPM_THROW(std::runtime_error, "Failure occured appending controls for well " << well_names[wix]);
-                    }
-                    WellsManagerDetail::ProductionControl::Mode mode = WellsManagerDetail::ProductionControl::mode(wcp_line.control_mode_);
-                    int cpos = control_pos[mode];
-                    if (cpos == -1 && mode != WellsManagerDetail::ProductionControl::GRUP) {
-                        OPM_THROW(std::runtime_error, "Control mode type " << mode << " not present in well " << well_names[wix]);
-                    }
-                    // If it's shut, we complement the cpos
-                    if (wcp_line.open_shut_flag_ == "SHUT") {
-                        cpos = ~cpos; // So we can easily retrieve the cpos later
-                    }
-                    set_current_control(wix, cpos, w_);
-                }
-            }
-            if (!well_found) {
-                OPM_THROW(std::runtime_error, "Undefined well name: " << wcp_line.well_
-                          << " in WCONPROD");
-            }
-        }
-    }
-
-    // Get WELTARG data
-    if (deck.hasField("WELTARG")) {
-        OPM_THROW(std::runtime_error, "We currently do not handle WELTARG.");
-        /*
-          const WELTARG& weltargs = deck.getWELTARG();
-          const int num_weltargs  = weltargs.weltarg.size();
-          for (int kw = 0; kw < num_weltargs; ++kw) {
-          std::string name = weltargs.weltarg[kw].well_;
-          std::string::size_type len = name.find('*');
-          if (len != std::string::npos) {
-          name = name.substr(0, len);
-          }
-          bool well_found = false;
-          for (int wix = 0; wix < num_wells; ++wix) {
-          if (well_names[wix].compare(0,len, name) == 0) { //equal
-          well_found = true;
-          well_data[wix].target = weltargs.weltarg[kw].new_value_;
-          break;
-          }
-          }
-          if (!well_found) {
-          OPM_THROW(std::runtime_error, "Undefined well name: " << weltargs.weltarg[kw].well_
-          << " in WELTARG");
-          }
-          }
-        */
-    }
+    setupGuideRates(wells, timeStep, well_data, well_names_to_index);
 
     // Debug output.
 #define EXTRA_OUTPUT
@@ -743,77 +309,6 @@ void WellsManager::init(const Opm::EclipseGridParser& deck,
       }
     */
 #endif
-
-    if (deck.hasField("WELOPEN")) {
-        const WELOPEN& welopen = deck.getWELOPEN();
-        for (size_t i = 0; i < welopen.welopen.size(); ++i) {
-            WelopenLine line = welopen.welopen[i];
-            std::string wellname = line.well_;
-            std::map<std::string, int>::const_iterator it = well_names_to_index.find(wellname);
-            if (it == well_names_to_index.end()) {
-                OPM_THROW(std::runtime_error, "Trying to open/shut well with name: \"" << wellname<<"\" but it's not registered under WELSPECS.");
-            }
-            const int index = it->second;
-            if (line.openshutflag_ == "SHUT") {
-                well_controls_shut_well( w_->ctrls[index] );
-            } else if (line.openshutflag_ == "OPEN") {
-                well_controls_open_well( w_->ctrls[index] );
-            } else {
-                OPM_THROW(std::runtime_error, "Unknown Open/close keyword: \"" << line.openshutflag_<< "\". Allowed values: OPEN, SHUT.");
-            }
-        }
-    }
-
-    // Build the well_collection_ well group hierarchy.
-    if (deck.hasField("GRUPTREE")) {
-        std::cout << "Found gruptree" << std::endl;
-        const GRUPTREE& gruptree = deck.getGRUPTREE();
-        std::map<std::string, std::string>::const_iterator it = gruptree.tree.begin();
-        for( ; it != gruptree.tree.end(); ++it) {
-            well_collection_.addChild(it->first, it->second, deck);
-        }
-    }
-    for (size_t i = 0; i < welspecs.welspecs.size(); ++i) {
-        WelspecsLine line = welspecs.welspecs[i];
-        well_collection_.addChild(line.name_, line.group_, deck);
-    }
-
-
-
-    // Set the guide rates:
-    if (deck.hasField("WGRUPCON")) {
-        std::cout << "Found Wgrupcon" << std::endl;
-        WGRUPCON wgrupcon = deck.getWGRUPCON();
-        const std::vector<WgrupconLine>& lines = wgrupcon.wgrupcon;
-        std::cout << well_collection_.getLeafNodes().size() << std::endl;
-        for (size_t i = 0; i < lines.size(); i++) {
-            std::string name = lines[i].well_;
-            const int wix = well_names_to_index[name];
-            WellNode& wellnode = *well_collection_.getLeafNodes()[wix];
-            assert(wellnode.name() == name);
-            if (well_data[wix].type == PRODUCER) {
-                wellnode.prodSpec().guide_rate_ = lines[i].guide_rate_;
-                if (lines[i].phase_ == "OIL") {
-                    wellnode.prodSpec().guide_rate_type_ = ProductionSpecification::OIL;
-                } else {
-                    OPM_THROW(std::runtime_error, "Guide rate type " << lines[i].phase_ << " specified for producer "
-                              << name << " in WGRUPCON, cannot handle.");
-                }
-            } else if (well_data[wix].type == INJECTOR) {
-                wellnode.injSpec().guide_rate_ = lines[i].guide_rate_;
-                if (lines[i].phase_ == "RAT") {
-                    wellnode.injSpec().guide_rate_type_ = InjectionSpecification::RAT;
-                } else {
-                    OPM_THROW(std::runtime_error, "Guide rate type " << lines[i].phase_ << " specified for injector "
-                              << name << " in WGRUPCON, cannot handle.");
-                }
-            } else {
-                OPM_THROW(std::runtime_error, "Unknown well type " << well_data[wix].type << " for well " << name);
-            }
-        }
-    }
-    well_collection_.setWellsPointer(w_);
-    well_collection_.applyGroupControls();
 }
 
 } // end namespace Opm

tests/test_linearsolver.cpp

@@ -0,0 +1,186 @@
+/*
+  Copyright 2014 Dr. Markus Blatt - HPC-Simulation-Software & Services
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <config.h>
+
+#if HAVE_DYNAMIC_BOOST_TEST
+#define BOOST_TEST_DYN_LINK
+#endif
+#define NVERBOSE // to suppress our messages when throwing
+
+#define BOOST_TEST_MODULE OPM-IterativeSolverTest
+#include <boost/test/unit_test.hpp>
+
+#include <opm/core/linalg/LinearSolverFactory.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+
+#include <dune/common/version.hh>
+#include <memory>
+#include <cstdlib>
+#include <string>
+
+struct MyMatrix
+{
+    MyMatrix(int rows, int nnz)
+        : data(nnz, 0.0), rowStart(rows+1, -1),
+          colIndex(nnz, -1)
+    {}
+    MyMatrix()
+        : data(), rowStart(), colIndex()
+    {}
+
+    std::vector<double> data;
+    std::vector<int> rowStart;
+    std::vector<int> colIndex;
+};
+
+std::shared_ptr<MyMatrix> createLaplacian(int N)
+{
+    MyMatrix* mm=new MyMatrix(N*N, N*N*5);
+    int nnz=0;
+    mm->rowStart[0]=0;
+    for(int row=0; row<N*N; row++)
+    {
+
+        int x=row%N;
+        int y=row/N;
+        if(y>0)
+        {
+            mm->colIndex[nnz]=row-N;
+            mm->data[nnz++]=-1;
+        }
+        if(x>0)
+        {
+            mm->colIndex[nnz]=row-1;
+            mm->data[nnz++]=-1;
+        }
+        mm->colIndex[nnz]=row;
+        mm->data[nnz++]=4;
+        if(x<N-1)
+        {
+            mm->colIndex[nnz]=row+1;
+            mm->data[nnz++]=-1;
+        }
+        if(y<N-1)
+        {
+            mm->colIndex[nnz]=row+N;
+            mm->data[nnz++]=-1;
+        }
+        mm->rowStart[row+1]=nnz;
+    }
+    mm->data.resize(nnz);
+    mm->colIndex.resize(nnz);
+    return std::shared_ptr<MyMatrix>(mm);
+}
+
+void createRandomVectors(int NN, std::vector<double>& x, std::vector<double>& b,
+                         const MyMatrix& mat)
+{
+    x.resize(NN);
+    for(auto entry=x.begin(), end =x.end(); entry!=end; ++entry)
+        *entry=((double) (rand()%100))/10.0;
+
+    b.resize(NN);
+    std::fill(b.begin(), b.end(), 0.0);
+
+    // Construct the right hand side as b=A*x
+    for(std::size_t row=0; row<mat.rowStart.size()-1; ++row)
+    {
+        for(int i=mat.rowStart[row], end=mat.rowStart[row+1]; i!=end; ++i)
+        {
+            b[row]+= mat.data[i]*x[mat.colIndex[i]];
+        }
+    }
+}
+
+void run_test(const Opm::parameter::ParameterGroup& param)
+{
+    int N=4;
+    auto mat = createLaplacian(N);
+    std::vector<double> x(N*N), b(N*N);
+    createRandomVectors(100*100, x, b, *mat);
+    std::vector<double> exact(x);
+    std::fill(x.begin(), x.end(), 0.0);
+    Opm::LinearSolverFactory ls(param);
+    ls.solve(N*N, mat->data.size(), &(mat->rowStart[0]),
+             &(mat->colIndex[0]), &(mat->data[0]), &(b[0]),
+             &(x[0]));
+}
+
+
+BOOST_AUTO_TEST_CASE(DefaultTest)
+{
+    Opm::parameter::ParameterGroup param;
+    param.insertParameter(std::string("linsolver_max_iterations"), std::string("200"));
+    param.insertParameter(std::string("linsolver_verbosity"), std::string("2"));
+    run_test(param);
+}
+
+#ifdef HAVE_DUNE_ISTL
+BOOST_AUTO_TEST_CASE(CGAMGTest)
+{
+    Opm::parameter::ParameterGroup param;
+    param.insertParameter(std::string("linsolver"), std::string("istl"));
+    param.insertParameter(std::string("linsolver_type"), std::string("1"));
+    param.insertParameter(std::string("linsolver_max_iterations"), std::string("200"));
+    param.insertParameter(std::string("linsolver_verbosity"), std::string("2"));
+    run_test(param);
+}
+
+BOOST_AUTO_TEST_CASE(CGILUTest)
+{
+    Opm::parameter::ParameterGroup param;
+    param.insertParameter(std::string("linsolver"), std::string("istl"));
+    param.insertParameter(std::string("linsolver_type"), std::string("0"));
+    param.insertParameter(std::string("linsolver_max_iterations"), std::string("200"));
+    param.insertParameter(std::string("linsolver_verbosity"), std::string("2"));
+    run_test(param);
+}
+
+BOOST_AUTO_TEST_CASE(BiCGILUTest)
+{
+    Opm::parameter::ParameterGroup param;
+    param.insertParameter(std::string("linsolver"), std::string("istl"));
+    param.insertParameter(std::string("linsolver_type"), std::string("2"));
+    param.insertParameter(std::string("linsolver_max_iterations"), std::string("200"));
+    param.insertParameter(std::string("linsolver_verbosity"), std::string("2"));
+    run_test(param);
+}
+
+#if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
+BOOST_AUTO_TEST_CASE(FastAMGTest)
+{
+    Opm::parameter::ParameterGroup param;
+    param.insertParameter(std::string("linsolver"), std::string("istl"));
+    param.insertParameter(std::string("linsolver_type"), std::string("3"));
+    param.insertParameter(std::string("linsolver_max_iterations"), std::string("200"));
+    param.insertParameter(std::string("linsolver_verbosity"), std::string("2"));
+    run_test(param);
+}
+
+BOOST_AUTO_TEST_CASE(KAMGTest)
+{
+    Opm::parameter::ParameterGroup param;
+    param.insertParameter(std::string("linsolver"), std::string("istl"));
+    param.insertParameter(std::string("linsolver_type"), std::string("4"));
+    param.insertParameter(std::string("linsolver_max_iterations"), std::string("200"));
+    run_test(param);
+}
+#endif
+#endif

tests/test_wellsgroup.cpp

@@ -59,16 +59,18 @@ BOOST_AUTO_TEST_CASE(ConstructGroupFromWell) {
         std::shared_ptr<WellsGroupInterface> wellsGroup = createWellWellsGroup(well, 2, pu);
         BOOST_CHECK_EQUAL(well->name(), wellsGroup->name());
         if (well->isInjector(2)) {
-            BOOST_CHECK_EQUAL(well->getSurfaceInjectionRate(2), wellsGroup->injSpec().surface_flow_max_rate_);
-            BOOST_CHECK_EQUAL(well->getBHPLimit(2), wellsGroup->injSpec().BHP_limit_);
-            BOOST_CHECK_EQUAL(well->getReservoirInjectionRate(2), wellsGroup->injSpec().reservoir_flow_max_rate_);
+            const WellInjectionProperties& properties = well->getInjectionProperties(2);
+            BOOST_CHECK_EQUAL(properties.surfaceInjectionRate, wellsGroup->injSpec().surface_flow_max_rate_);
+            BOOST_CHECK_EQUAL(properties.BHPLimit, wellsGroup->injSpec().BHP_limit_);
+            BOOST_CHECK_EQUAL(properties.reservoirInjectionRate, wellsGroup->injSpec().reservoir_flow_max_rate_);
             BOOST_CHECK_EQUAL(0.0, wellsGroup->prodSpec().guide_rate_);
         }
         if (well->isProducer(2)) {
-            BOOST_CHECK_EQUAL(well->getResVRate(2), wellsGroup->prodSpec().reservoir_flow_max_rate_);
-            BOOST_CHECK_EQUAL(well->getBHPLimit(2), wellsGroup->prodSpec().BHP_limit_);
-            BOOST_CHECK_EQUAL(well->getOilRate(2), wellsGroup->prodSpec().oil_max_rate_);
-            BOOST_CHECK_EQUAL(well->getWaterRate(2), wellsGroup->prodSpec().water_max_rate_);
+            const WellProductionProperties& properties = well->getProductionProperties(2);
+            BOOST_CHECK_EQUAL(properties.ResVRate, wellsGroup->prodSpec().reservoir_flow_max_rate_);
+            BOOST_CHECK_EQUAL(properties.BHPLimit, wellsGroup->prodSpec().BHP_limit_);
+            BOOST_CHECK_EQUAL(properties.OilRate, wellsGroup->prodSpec().oil_max_rate_);
+            BOOST_CHECK_EQUAL(properties.WaterRate, wellsGroup->prodSpec().water_max_rate_);
             BOOST_CHECK_EQUAL(0.0, wellsGroup->injSpec().guide_rate_);
         }
     }

tests/test_wellsmanager.cpp

@@ -173,30 +173,6 @@ void check_controls_epoch1( struct WellControls ** ctrls) {
     }
 }
 
-
-BOOST_AUTO_TEST_CASE(Constructor_Works) {
-    Opm::EclipseGridParser Deck("wells_manager_data.data");
-    Opm::GridManager gridManager(Deck);
-
-    Deck.setCurrentEpoch(0);
-    {
-        Opm::WellsManager wellsManager(Deck, *gridManager.c_grid(), NULL);
-        const Wells* wells = wellsManager.c_wells();
-        wells_static_check( wells );
-        check_controls_epoch0( wells->ctrls );
-    }
-
-
-    Deck.setCurrentEpoch(1);
-    {
-        Opm::WellsManager wellsManager(Deck, *gridManager.c_grid(), NULL);
-        const Wells* wells = wellsManager.c_wells();
-        
-        wells_static_check( wells );
-        check_controls_epoch1( wells->ctrls );
-    }
-}
-
 BOOST_AUTO_TEST_CASE(New_Constructor_Works) {
 
     Opm::ParserPtr parser(new Opm::Parser());
@@ -205,70 +181,16 @@ BOOST_AUTO_TEST_CASE(New_Constructor_Works) {
     Opm::EclipseGridParser Deck("wells_manager_data.data");
     Opm::GridManager gridManager(Deck);
 
-    Deck.setCurrentEpoch(0);
     {
         Opm::WellsManager wellsManager(eclipseState, 0, *gridManager.c_grid(), NULL);
-        Opm::WellsManager oldWellsManager(Deck, *gridManager.c_grid(), NULL);
-
-        std::cout << "Checking new well structure, epoch 0" << std::endl;
         wells_static_check( wellsManager.c_wells() );
-
-        std::cout << "Checking old well structure, epoch 0" << std::endl;
-        wells_static_check( oldWellsManager.c_wells() );
-
         check_controls_epoch0( wellsManager.c_wells()->ctrls );
-
-        BOOST_CHECK(wells_equal(wellsManager.c_wells(), oldWellsManager.c_wells() , false));
     }
 
-    Deck.setCurrentEpoch(1);
     {
         Opm::WellsManager wellsManager(eclipseState, 1, *gridManager.c_grid(), NULL);
-        Opm::WellsManager oldWellsManager(Deck, *gridManager.c_grid(), NULL);
-
-        std::cout << "Checking new well structure, epoch 1" << std::endl;
         wells_static_check( wellsManager.c_wells() );
-
-        std::cout << "Checking old well structure, epoch 1" << std::endl;
-        wells_static_check( oldWellsManager.c_wells() );
-
         check_controls_epoch1( wellsManager.c_wells()->ctrls );
-
-        BOOST_CHECK(wells_equal( wellsManager.c_wells(), oldWellsManager.c_wells(),false));
-    }
-}
-
-
-BOOST_AUTO_TEST_CASE(New_Constructor_Works_ExpandedData) {
-
-    Opm::ParserPtr parser(new Opm::Parser());
-    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(parser->parseFile("wells_manager_data_expanded.data")));
-
-    Opm::EclipseGridParser Deck("wells_manager_data_expanded.data");
-    Opm::GridManager gridManager(Deck);
-
-    Deck.setCurrentEpoch(0);
-    {
-        Opm::WellsManager wellsManager(eclipseState, 0, *gridManager.c_grid(), NULL);
-        Opm::WellsManager oldWellsManager(Deck, *gridManager.c_grid(), NULL);
-
-        BOOST_CHECK(wells_equal(wellsManager.c_wells(), oldWellsManager.c_wells(),false));
-    }
-
-    Deck.setCurrentEpoch(1);
-    {
-        Opm::WellsManager wellsManager(eclipseState, 1, *gridManager.c_grid(), NULL);
-        Opm::WellsManager oldWellsManager(Deck, *gridManager.c_grid(), NULL);
-
-        BOOST_CHECK(wells_equal( wellsManager.c_wells(), oldWellsManager.c_wells(), true));
-    }
-
-    Deck.setCurrentEpoch(2);
-    {
-        Opm::WellsManager wellsManager(eclipseState, 2, *gridManager.c_grid(), NULL);
-        Opm::WellsManager oldWellsManager(Deck, *gridManager.c_grid(), NULL);
-
-        BOOST_CHECK(wells_equal( wellsManager.c_wells(), oldWellsManager.c_wells(), true));
     }
 }
 
@@ -277,12 +199,11 @@ BOOST_AUTO_TEST_CASE(New_Constructor_Works_ExpandedData) {
 BOOST_AUTO_TEST_CASE(WellsEqual) {
     Opm::EclipseGridParser Deck("wells_manager_data.data");
     Opm::GridManager gridManager(Deck);
+    Opm::ParserPtr parser(new Opm::Parser());
+    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(parser->parseFile("wells_manager_data.data")));
 
-    Deck.setCurrentEpoch(0);
-    Opm::WellsManager wellsManager0(Deck, *gridManager.c_grid(), NULL);
-
-    Deck.setCurrentEpoch(1);
-    Opm::WellsManager wellsManager1(Deck, *gridManager.c_grid(), NULL);
+    Opm::WellsManager wellsManager0(eclipseState , 0 , *gridManager.c_grid(), NULL);
+    Opm::WellsManager wellsManager1(eclipseState , 1 , *gridManager.c_grid(), NULL);
 
     BOOST_CHECK(  wells_equal( wellsManager0.c_wells() , wellsManager0.c_wells(),false) );
     BOOST_CHECK( !wells_equal( wellsManager0.c_wells() , wellsManager1.c_wells(),false) );
@@ -293,11 +214,11 @@ BOOST_AUTO_TEST_CASE(ControlsEqual) {
     Opm::EclipseGridParser Deck("wells_manager_data.data");
     Opm::GridManager gridManager(Deck);
 
-    Deck.setCurrentEpoch(0);
-    Opm::WellsManager wellsManager0(Deck, *gridManager.c_grid(), NULL);
+    Opm::ParserPtr parser(new Opm::Parser());
+    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(parser->parseFile("wells_manager_data.data")));
 
-    Deck.setCurrentEpoch(1);
-    Opm::WellsManager wellsManager1(Deck, *gridManager.c_grid(), NULL);
+    Opm::WellsManager wellsManager0(eclipseState , 0 , *gridManager.c_grid(), NULL);
+    Opm::WellsManager wellsManager1(eclipseState , 1 , *gridManager.c_grid(), NULL);
 
     BOOST_CHECK(  well_controls_equal( wellsManager0.c_wells()->ctrls[0] , wellsManager0.c_wells()->ctrls[0] , false));
     BOOST_CHECK(  well_controls_equal( wellsManager0.c_wells()->ctrls[1] , wellsManager0.c_wells()->ctrls[1] , false));