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

Manually resolved conficts in the following files
	examples/sim_fibo_ad.cpp
	opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
	opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp

In additions examples/sim_fibo_ad_cp.cpp was adapted to compile again.

CMakeLists_files.cmake

@@ -109,4 +109,5 @@ list (APPEND PUBLIC_HEADER_FILES
 	opm/autodiff/SimulatorIncompTwophaseAd.hpp
 	opm/autodiff/TransportSolverTwophaseAd.hpp
 	opm/autodiff/WellDensitySegmented.hpp
+	opm/autodiff/WellStateFullyImplicitBlackoil.hpp
 	)

examples/sim_fibo_ad.cpp

@@ -39,7 +39,7 @@
 #include <opm/core/linalg/LinearSolverFactory.hpp>
 
 #include <opm/core/simulator/BlackoilState.hpp>
-#include <opm/core/simulator/WellState.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 
 #include <opm/autodiff/SimulatorFullyImplicitBlackoil.hpp>
 #include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
@@ -103,10 +103,6 @@ try
 #if USE_NEW_PARSER
     Opm::ParserPtr newParser(new Opm::Parser() );
     Opm::DeckConstPtr newParserDeck = newParser->parseFile( deck_filename );
-
-#warning "HACK: required until the WellsManager and the EclipseWriter don't require the old parser anymore"
-    std::shared_ptr<EclipseGridParser> deck;
-    deck.reset(new EclipseGridParser(deck_filename));
 #else
     std::shared_ptr<EclipseGridParser> deck;
     deck.reset(new EclipseGridParser(deck_filename));
@@ -119,8 +115,7 @@ try
     grid.reset(new GridManager(*deck));
 #endif
 
-#warning "HACK: required until the WellsManager and the EclipseWriter don't require the old parser anymore"
+#if USE_NEW_PARSER
-#if 0 // USE_NEW_PARSER
     Opm::EclipseWriter outputWriter(param, newParserDeck, share_obj(*grid->c_grid()));
 #else
     Opm::EclipseWriter outputWriter(param, deck, share_obj(*grid->c_grid()));
@@ -200,21 +195,45 @@ try
         param.writeParam(output_dir + "/simulation.param");
     }
 
-#warning "TODO: convert the well handling code to the new parser"
 #if USE_NEW_PARSER
     std::cout << "\n\n================    Starting main simulation loop     ===============\n"
               << std::flush;
 
-    WellState well_state;
+    WellStateFullyImplicitBlackoil well_state;
     Opm::TimeMapPtr timeMap(new Opm::TimeMap(newParserDeck));
     SimulatorTimer simtimer;
+    std::shared_ptr<EclipseState> eclipseState(new EclipseState(newParserDeck));
 
-    // Create new wells, well_state
+    // initialize variables
-    WellsManager wells(*deck, *grid->c_grid(), props->permeability());
+    simtimer.init(timeMap, /*beginReportStepIdx=*/0, /*endReportStepIdx=*/0);
+
+    SimulatorReport fullReport;
+    for (size_t reportStepIdx = 0; reportStepIdx < timeMap->numTimesteps(); ++reportStepIdx) {
+        // Report on start of a report step.
+        std::cout << "\n"
+                  << "---------------------------------------------------------------\n"
+                  << "--------------    Starting report step " << reportStepIdx << "    --------------\n"
+                  << "---------------------------------------------------------------\n"
+                  << "\n";
+
+        WellsManager wells(eclipseState,
+                           reportStepIdx,
+                           *grid->c_grid(),
+                           props->permeability());
+
+        if (reportStepIdx == 0) {
             // @@@ 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.
             well_state.init(wells.c_wells(), state);
+        }
+
+        simtimer.init(timeMap,
+                      /*beginReportStepIdx=*/reportStepIdx,
+                      /*endReportStepIdx=*/reportStepIdx + 1);
+
+        if (reportStepIdx == 0)
+            outputWriter.writeInit(simtimer, state, well_state.basicWellState());
 
         // Create and run simulator.
     SimulatorFullyImplicitBlackoil<UnstructuredGrid> simulator(param,
@@ -225,21 +244,24 @@ try
                                                                linsolver,
                                                                grav,
                                                                outputWriter);
-    simtimer.init(timeMap);
+        SimulatorReport episodeReport = simulator.run(simtimer, state, well_state);
-    SimulatorReport report = simulator.run(simtimer, state, well_state);
+
+        outputWriter.writeTimeStep(simtimer, state, well_state.basicWellState());
+        fullReport += episodeReport;
 
         if (output) {
-        report.reportParam(outStream);
+            episodeReport.reportParam(outStream);
-        warnIfUnusedParams(param);
+        }
     }
 
     std::cout << "\n\n================    End of simulation     ===============\n\n";
-    report.report(std::cout);
+    fullReport.report(std::cout);
 
     if (output) {
         std::string filename = output_dir + "/walltime.param";
         std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out);
-        report.reportParam(tot_os);
+        fullReport.reportParam(tot_os);
+        warnIfUnusedParams(param);
     }
 #else
     std::cout << "\n\n================    Starting main simulation loop     ===============\n"
@@ -248,12 +270,13 @@ try
 
     SimulatorReport rep;
     // With a deck, we may have more epochs etc.
-    WellState well_state;
+    WellStateFullyImplicitBlackoil well_state;
     int step = 0;
     SimulatorTimer simtimer;
     // Use timer for last epoch to obtain total time.
     deck->setCurrentEpoch(deck->numberOfEpochs() - 1);
     simtimer.init(*deck);
+
     const double total_time = simtimer.totalTime();
     for (int epoch = 0; epoch < deck->numberOfEpochs(); ++epoch) {
         // Set epoch index.
@@ -285,6 +308,9 @@ try
             well_state.init(wells.c_wells(), state);
         }
 
+        if (epoch == 0)
+            outputWriter.writeInit(simtimer, state, well_state.basicWellState());
+
         // Create and run simulator.
         SimulatorFullyImplicitBlackoil<UnstructuredGrid> simulator(param,
                                                  *grid->c_grid(),
@@ -294,6 +320,8 @@ try
                                                  linsolver,
                                                  grav,
                                                  outputWriter);
+        outputWriter.writeTimeStep(simtimer, state, well_state.basicWellState());
+
         if (epoch == 0) {
             warnIfUnusedParams(param);
         }

examples/sim_fibo_ad_cp.cpp

@@ -48,17 +48,20 @@
 #include <opm/core/linalg/LinearSolverFactory.hpp>
 
 #include <opm/core/simulator/BlackoilState.hpp>
-#include <opm/core/simulator/WellState.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 
 #include <opm/autodiff/SimulatorFullyImplicitBlackoil.hpp>
 #include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
 #include <opm/autodiff/GridHelpers.hpp>
 #include <opm/core/utility/share_obj.hpp>
 
-#include <boost/scoped_ptr.hpp>
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+#include <opm/parser/eclipse/Parser/Parser.hpp>
+
 #include <boost/filesystem.hpp>
 #include <boost/algorithm/string.hpp>
 
+#include <memory>
 #include <algorithm>
 #include <iostream>
 #include <vector>
@@ -97,23 +100,24 @@ try
         OPM_THROW(std::runtime_error, "This program must be run with an input deck. "
                   "Specify the deck with deck_filename=deckname.data (for example).");
     }
-    boost::scoped_ptr<EclipseGridParser> deck;
+    std::shared_ptr<Dune::CpGrid> grid;
-    boost::scoped_ptr<Dune::CpGrid> grid;
+    std::shared_ptr<BlackoilPropertiesInterface> props;
-    boost::scoped_ptr<BlackoilPropertiesInterface> props;
+    std::shared_ptr<BlackoilPropsAdInterface> new_props;
-    boost::scoped_ptr<BlackoilPropsAdInterface> new_props;
+    std::shared_ptr<RockCompressibility> rock_comp;
-    boost::scoped_ptr<RockCompressibility> rock_comp;
     BlackoilState state;
     // bool check_well_controls = false;
-    // int max_well_control_iterations = 0
+    // int max_well_control_iterations = 0;
     double gravity[3] = { 0.0 };
     std::string deck_filename = param.get<std::string>("deck_filename");
+    std::shared_ptr<EclipseGridParser> deck;
     deck.reset(new EclipseGridParser(deck_filename));
+    std::shared_ptr<const EclipseGridParser> cdeck(deck);
     // Grid init
     grid.reset(new Dune::CpGrid());
     
     grid->processEclipseFormat(*deck, 2e-12, false);
 
-    Opm::EclipseWriter outputWriter(param, share_obj(*deck),
+    Opm::EclipseWriter outputWriter(param, cdeck,
                                     Opm::UgGridHelpers::numCells(*grid),
                                     Opm::UgGridHelpers::globalCell(*grid),
                                     Opm::UgGridHelpers::cartDims(*grid),
@@ -194,12 +198,13 @@ try
 
     SimulatorReport rep;
     // With a deck, we may have more epochs etc.
-    WellState well_state;
+    WellStateFullyImplicitBlackoil well_state;
     int step = 0;
     SimulatorTimer simtimer;
     // Use timer for last epoch to obtain total time.
     deck->setCurrentEpoch(deck->numberOfEpochs() - 1);
     simtimer.init(*deck);
+
     const double total_time = simtimer.totalTime();
     for (int epoch = 0; epoch < deck->numberOfEpochs(); ++epoch) {
         // Set epoch index.

examples/test_implicit_ad.cpp

@@ -23,6 +23,7 @@
 #include <opm/autodiff/FullyImplicitBlackoilSolver.hpp>
 #include <opm/autodiff/GeoProps.hpp>
 #include <opm/autodiff/BlackoilPropsAd.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 
 #include <opm/core/grid.h>
 #include <opm/core/wells.h>
@@ -37,7 +38,6 @@
 #include <opm/core/utility/Units.hpp>
 
 #include <opm/core/simulator/BlackoilState.hpp>
-#include <opm/core/simulator/WellState.hpp>
 #include <opm/core/simulator/initState.hpp>
 
 #include <boost/shared_ptr.hpp>
@@ -110,7 +110,7 @@ try
     initStateBasic(*g, props0, param, 0.0, state);
     initBlackoilSurfvol(*g, props0, state);
 
-    Opm::WellState well_state;
+    Opm::WellStateFullyImplicitBlackoil well_state;
     well_state.init(wells.get(), state);
 
     solver.step(1.0, state, well_state);

opm/autodiff/FullyImplicitBlackoilSolver.hpp

@@ -33,7 +33,7 @@ namespace Opm {
     class RockCompressibility;
     class LinearSolverInterface;
     class BlackoilState;
-    class WellState;
+    class WellStateFullyImplicitBlackoil;
 
 
     /// A fully implicit solver for the black-oil problem.
@@ -79,7 +79,7 @@ namespace Opm {
         void
         step(const double   dt    ,
              BlackoilState& state ,
-             WellState&     wstate);
+             WellStateFullyImplicitBlackoil&     wstate);
 
     private:
         // Types and enums
@@ -151,11 +151,11 @@ namespace Opm {
         // Private methods.
         SolutionState
         constantState(const BlackoilState& x,
-                      const WellState&     xw);
+                      const WellStateFullyImplicitBlackoil& xw);
 
         SolutionState
         variableState(const BlackoilState& x,
-                      const WellState&     xw);
+                      const WellStateFullyImplicitBlackoil& xw);
 
         void
         computeAccum(const SolutionState& state,
@@ -164,13 +164,13 @@ namespace Opm {
         void
         assemble(const V&             dtpv,
                  const BlackoilState& x,
-                 const WellState&     xw  );
+                 const WellStateFullyImplicitBlackoil& xw);
 
         V solveJacobianSystem() const;
 
         void updateState(const V& dx,
                          BlackoilState& state,
-                         WellState& well_state);
+                         WellStateFullyImplicitBlackoil& well_state);
 
         std::vector<ADB>
         computePressures(const SolutionState& state) const;

opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp

@@ -25,12 +25,12 @@
 #include <opm/autodiff/GridHelpers.hpp>
 #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
 #include <opm/autodiff/GeoProps.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 
 #include <opm/core/grid.h>
 #include <opm/core/linalg/LinearSolverInterface.hpp>
 #include <opm/core/props/rock/RockCompressibility.hpp>
 #include <opm/core/simulator/BlackoilState.hpp>
-#include <opm/core/simulator/WellState.hpp>
 #include <opm/core/utility/ErrorMacros.hpp>
 #include <opm/core/well_controls.h>
 
@@ -231,7 +231,7 @@ namespace {
     FullyImplicitBlackoilSolver<T>::
     step(const double   dt,
          BlackoilState& x ,
-         WellState&     xw)
+         WellStateFullyImplicitBlackoil& xw)
     {
         const V pvdt = geo_.poreVolume() / dt;
 
@@ -343,7 +343,7 @@ namespace {
     template<class T>
     typename FullyImplicitBlackoilSolver<T>::SolutionState
     FullyImplicitBlackoilSolver<T>::constantState(const BlackoilState& x,
-                                               const WellState&     xw)
+                                                  const WellStateFullyImplicitBlackoil&     xw)
     {
         using namespace Opm::AutoDiffGrid;
         const int nc = numCells(grid_);
@@ -441,7 +441,7 @@ namespace {
     template<class T>
     typename FullyImplicitBlackoilSolver<T>::SolutionState
     FullyImplicitBlackoilSolver<T>::variableState(const BlackoilState& x,
-                                               const WellState&     xw)
+                                                  const WellStateFullyImplicitBlackoil&     xw)
     {
         using namespace Opm::AutoDiffGrid;
         const int nc = numCells(grid_);
@@ -633,7 +633,7 @@ namespace {
     FullyImplicitBlackoilSolver<T>::
     assemble(const V&             pvdt,
              const BlackoilState& x   ,
-             const WellState&     xw  )
+             const WellStateFullyImplicitBlackoil& xw  )
     {
         using namespace Opm::AutoDiffGrid;
         // Create the primary variables.
@@ -898,7 +898,7 @@ namespace {
     template<class T>
     void FullyImplicitBlackoilSolver<T>::updateState(const V& dx,
                                                   BlackoilState& state,
-                                                  WellState& well_state)
+                                                  WellStateFullyImplicitBlackoil& well_state)
     {
         using namespace Opm::AutoDiffGrid;
         const int np = fluid_.numPhases();

opm/autodiff/SimulatorCompressibleAd.cpp

@@ -496,13 +496,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/autodiff/SimulatorFullyImplicitBlackoil.hpp

@@ -37,7 +37,7 @@ namespace Opm
     class LinearSolverInterface;
     class SimulatorTimer;
     class BlackoilState;
-    class WellState;
+    class WellStateFullyImplicitBlackoil;
     struct SimulatorReport;
 
     /// Class collecting all necessary components for a two-phase simulation.
@@ -87,7 +87,7 @@ namespace Opm
         /// \return                    simulation report, with timing data
         SimulatorReport run(SimulatorTimer& timer,
                             BlackoilState& state,
-                            WellState& well_state);
+                            WellStateFullyImplicitBlackoil& well_state);
 
     private:
         class Impl;

opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp

@@ -17,6 +17,7 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 
+#include  <opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp>
 #include <opm/autodiff/SimulatorFullyImplicitBlackoil.hpp>
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
 #include <opm/core/utility/ErrorMacros.hpp>
@@ -24,7 +25,7 @@
 #include <opm/autodiff/GeoProps.hpp>
 #include <opm/autodiff/FullyImplicitBlackoilSolver.hpp>
 #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
-#include <opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 #include <opm/core/grid.h>
 #include <opm/core/wells.h>
 #include <opm/core/pressure/flow_bc.h>
@@ -33,6 +34,7 @@
 #include <opm/core/simulator/SimulatorTimer.hpp>
 #include <opm/core/utility/StopWatch.hpp>
 #include <opm/core/io/eclipse/EclipseWriter.hpp>
+#include <opm/core/io/vtk/writeVtkData.hpp>
 #include <opm/core/utility/miscUtilities.hpp>
 #include <opm/core/utility/miscUtilitiesBlackoil.hpp>
 
@@ -42,7 +44,6 @@
 
 #include <opm/core/grid/ColumnExtract.hpp>
 #include <opm/core/simulator/BlackoilState.hpp>
-#include <opm/core/simulator/WellState.hpp>
 #include <opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp>
 
 #include <boost/filesystem.hpp>
@@ -71,7 +72,7 @@ namespace Opm
 
         SimulatorReport run(SimulatorTimer& timer,
                             BlackoilState& state,
-                            WellState& well_state);
+                            WellStateFullyImplicitBlackoil& well_state);
 
     private:
         // Data.
@@ -123,13 +124,69 @@ namespace Opm
     template<class T>
     SimulatorReport SimulatorFullyImplicitBlackoil<T>::run(SimulatorTimer& timer,
                                                         BlackoilState& state,
-                                                        WellState& well_state)
+                                                        WellStateFullyImplicitBlackoil& well_state)
     {
         return pimpl_->run(timer, state, well_state);
     }
 
 
 
+    static void outputWellStateMatlab(const Opm::WellStateFullyImplicitBlackoil& well_state,
+                                  const int step,
+                                  const std::string& output_dir)
+    {
+        Opm::DataMap dm;
+        dm["bhp"] = &well_state.bhp();
+        dm["wellrates"] = &well_state.wellRates();
+
+        // Write data (not grid) in Matlab format
+        for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) {
+            std::ostringstream fname;
+            fname << output_dir << "/" << it->first;
+            boost::filesystem::path fpath = fname.str();
+            try {
+                create_directories(fpath);
+            }
+            catch (...) {
+                OPM_THROW(std::runtime_error,"Creating directories failed: " << fpath);
+            }
+            fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt";
+            std::ofstream file(fname.str().c_str());
+            if (!file) {
+                OPM_THROW(std::runtime_error,"Failed to open " << fname.str());
+            }
+            file.precision(15);
+            const std::vector<double>& d = *(it->second);
+            std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
+        }
+    }
+
+#if 0
+    static void outputWaterCut(const Opm::Watercut& watercut,
+                               const std::string& output_dir)
+    {
+        // Write water cut curve.
+        std::string fname = output_dir  + "/watercut.txt";
+        std::ofstream os(fname.c_str());
+        if (!os) {
+            OPM_THROW(std::runtime_error, "Failed to open " << fname);
+        }
+        watercut.write(os);
+    }
+
+    static void outputWellReport(const Opm::WellReport& wellreport,
+                                 const std::string& output_dir)
+    {
+        // Write well report.
+        std::string fname = output_dir  + "/wellreport.txt";
+        std::ofstream os(fname.c_str());
+        if (!os) {
+            OPM_THROW(std::runtime_error, "Failed to open " << fname);
+        }
+        wellreport.write(os);
+    }
+#endif
+
 
     // \TODO: Treat bcs.
     template<class T>
@@ -187,11 +244,8 @@ namespace Opm
     template<class T>
     SimulatorReport SimulatorFullyImplicitBlackoil<T>::Impl::run(SimulatorTimer& timer,
                                                               BlackoilState& state,
-                                                              WellState& well_state)
+                                                              WellStateFullyImplicitBlackoil& well_state)
     {
-        eclipseWriter_.writeInit(timer, state, well_state);
-        eclipseWriter_.writeTimeStep(timer, state, well_state);
-
         // Initialisation.
         std::vector<double> porevol;
         if (rock_comp_props_ && rock_comp_props_->isActive()) {
@@ -259,13 +313,9 @@ namespace Opm
 
                 // Optionally, check if well controls are satisfied.
                 if (check_well_controls_) {
-                    Opm::computePhaseFlowRatesPerWell(*wells_,
-                                                      well_state.perfRates(),
-                                                      fractional_flows,
-                                                      well_resflows_phase);
                     std::cout << "Checking well conditions." << std::endl;
                     // For testing we set surface := reservoir
-                    well_control_passed = wells_manager_.conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase);
+                    well_control_passed = wells_manager_.conditionsMet(well_state.bhp(), well_state.wellRates(), well_state.wellRates());
                     ++well_control_iteration;
                     if (!well_control_passed && well_control_iteration > max_well_control_iterations_) {
                         OPM_THROW(std::runtime_error, "Could not satisfy well conditions in " << max_well_control_iterations_ << " tries.");
@@ -299,13 +349,13 @@ namespace Opm
                 tstep_os.close();
             }
 
-            // write an output file for later inspection
-            if (output_) {
-                eclipseWriter_.writeTimeStep(timer, state, well_state);
-            }
-
             // advance to next timestep before reporting at this location
             ++timer;
+
+            // write an output file for later inspection
+            if (output_) {
+                eclipseWriter_.writeTimeStep(timer, state, well_state.basicWellState());
+            }
         }
 
         total_timer.stop();

opm/autodiff/SimulatorIncompTwophaseAd.cpp

@@ -582,12 +582,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/autodiff/WellStateFullyImplicitBlackoil.hpp

@@ -0,0 +1,101 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_WELLSTATEFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
+#define OPM_WELLSTATEFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
+
+
+#include <opm/core/wells.h>
+#include <opm/core/well_controls.h>
+#include <opm/core/simulator/WellState.hpp>
+#include <vector>
+#include <cassert>
+
+namespace Opm
+{
+
+    /// The state of a set of wells, tailored for use by the fully
+    /// implicit blackoil simulator.
+    class WellStateFullyImplicitBlackoil
+    {
+    public:
+        /// Allocate and initialize if wells is non-null.  Also tries
+        /// to give useful initial values to the bhp(), wellRates()
+        /// and perfPhaseRates() fields, depending on controls
+        template <class State>
+        void init(const Wells* wells, const State& state)
+        {
+            if (wells == 0) {
+                return;
+            }
+
+            // We use the WellState::init() function to do bhp and well rates init.
+            // The alternative would be to copy that function wholesale.
+            basic_well_state_.init(wells, state);
+
+            // Initialize perfphaserates_, which must be done here.
+            const int nw = wells->number_of_wells;
+            const int np = wells->number_of_phases;
+            const int nperf = wells->well_connpos[nw];
+            perfphaserates_.resize(nperf * 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];
+                if (well_controls_well_is_shut(ctrl)) {
+                    // Shut well: perfphaserates_ are all zero.
+                } else {
+                    // Open well: Initialize perfphaserates_ to well
+                    // rates divided by the number of perforations.
+                    const int num_perf_this_well = wells->well_connpos[w + 1] - wells->well_connpos[w];
+                    for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) {
+                        for (int p = 0; p < np; ++p) {
+                            perfphaserates_[np*perf + p] = wellRates()[np*w + p] / double(num_perf_this_well);
+                        }
+                    }
+                }
+            }
+        }
+
+        /// One bhp pressure per well.
+        std::vector<double>& bhp() { return basic_well_state_.bhp(); }
+        const std::vector<double>& bhp() const { return basic_well_state_.bhp(); }
+
+        /// One rate per well and phase.
+        std::vector<double>& wellRates() { return basic_well_state_.wellRates(); }
+        const std::vector<double>& wellRates() const { return basic_well_state_.wellRates(); }
+
+        /// One rate per phase and well connection.
+        std::vector<double>& perfPhaseRates() { return perfphaserates_; }
+        const std::vector<double>& perfPhaseRates() const { return perfphaserates_; }
+
+        /// For interfacing with functions that take a WellState.
+        const WellState& basicWellState() const
+        {
+            return basic_well_state_;
+        }
+
+    private:
+        WellState basic_well_state_;
+        std::vector<double> perfphaserates_;
+    };
+
+} // namespace Opm
+
+
+#endif // OPM_WELLSTATEFULLYIMPLICITBLACKOIL_HEADER_INCLUDED