Merge branch 'master' into support-resv

Conflicts:
	examples/sim_fibo_ad.cpp
	examples/sim_fibo_ad_cp.cpp
	opm/autodiff/SimulatorFullyImplicitBlackoil.hpp
	opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp

examples/sim_fibo_ad.cpp

@@ -163,9 +163,9 @@ try
 
     // Write parameters used for later reference.
     bool output = param.getDefault("output", true);
-    std::ofstream outStream;
     std::string output_dir;
     if (output) {
+        // Create output directory if needed.
         output_dir =
             param.getDefault("output_dir", std::string("output"));
         boost::filesystem::path fpath(output_dir);
@@ -175,20 +175,11 @@ try
         catch (...) {
             OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
         }
-        std::string filename = output_dir + "/timing.param";
+        // Write simulation parameters.
-        outStream.open(filename.c_str(), std::fstream::trunc | std::fstream::out);
-        // open file to clean it. The file is appended to in SimulatorTwophase
-        filename = output_dir + "/step_timing.param";
-        std::fstream step_os(filename.c_str(), std::fstream::trunc | std::fstream::out);
-        step_os.close();
         param.writeParam(output_dir + "/simulation.param");
     }
 
-    std::cout << "\n\n================    Starting main simulation loop     ===============\n"
+    Opm::TimeMapConstPtr timeMap(eclipseState->getSchedule()->getTimeMap());
-              << std::flush;
-
-    WellStateFullyImplicitBlackoil well_state;
-    Opm::TimeMapPtr timeMap(new Opm::TimeMap(deck));
     SimulatorTimer simtimer;
 
     // initialize variables
@@ -196,60 +187,28 @@ try
 
     Opm::DerivedGeology geology(*grid->c_grid(), *new_props, eclipseState, grav);
 
-    SimulatorReport fullReport;
+    SimulatorFullyImplicitBlackoil<UnstructuredGrid> simulator(param,
-    for (size_t reportStepIdx = 0; reportStepIdx < timeMap->numTimesteps(); ++reportStepIdx) {
+                                             *grid->c_grid(),
-        // Report on start of a report step.
+                                             geology,
-        std::cout << "\n"
+                                             *new_props,
-                  << "---------------------------------------------------------------\n"
+                                             rock_comp->isActive() ? rock_comp.get() : 0,
-                  << "--------------    Starting report step " << reportStepIdx << "    --------------\n"
+                                             *fis_solver,
-                  << "---------------------------------------------------------------\n"
+                                             grav,
-                  << "\n";
+                                             deck->hasKeyword("DISGAS"),
+                                             deck->hasKeyword("VAPOIL"),
+                                             eclipseState,
+                                             outputWriter);
 
-        WellsManager wells(eclipseState,
+    std::cout << "\n\n================ Starting main simulation loop ===============\n"
-                           reportStepIdx,
+              << std::flush;
-                           *grid->c_grid(),
-                           props->permeability());
 
-        if (reportStepIdx == 0) {
+    SimulatorReport fullReport = simulator.run(simtimer, state);
-            // @@@ 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.setCurrentStepNum(reportStepIdx);
-
-        if (reportStepIdx == 0) {
-            outputWriter.writeInit(simtimer);
-            outputWriter.writeTimeStep(simtimer, state, well_state.basicWellState());
-        }
-
-        // Create and run simulator.
-        SimulatorFullyImplicitBlackoil<UnstructuredGrid> simulator(param,
-                                                                   eclipseState->getSchedule(),
-                                                 *grid->c_grid(),
-                                                 geology,
-                                                 *new_props,
-                                                 rock_comp->isActive() ? rock_comp.get() : 0,
-                                                 wells.c_wells(),
-                                                 *fis_solver,
-                                                 grav,
-                                                 deck->hasKeyword("DISGAS"),
-                                                 deck->hasKeyword("VAPOIL") );
-
-        SimulatorReport episodeReport = simulator.run(simtimer, state, well_state);
-
-        ++simtimer;
-
-        outputWriter.writeTimeStep(simtimer, state, well_state.basicWellState());
-        fullReport += episodeReport;
-    }
 
     std::cout << "\n\n================    End of simulation     ===============\n\n";
     fullReport.report(std::cout);
 
     if (output) {
-        std::string filename = output_dir + "/walltime.param";
+        std::string filename = output_dir + "/walltime.txt";
         std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out);
         fullReport.reportParam(tot_os);
         warnIfUnusedParams(param);

examples/sim_fibo_ad_cp.cpp

@@ -48,6 +48,7 @@
 #include <opm/core/wells/WellsManager.hpp>
 #include <opm/core/utility/ErrorMacros.hpp>
 #include <opm/core/simulator/initState.hpp>
+#include <opm/core/simulator/initStateEquil.hpp>
 #include <opm/core/simulator/SimulatorReport.hpp>
 #include <opm/core/simulator/SimulatorTimer.hpp>
 #include <opm/core/utility/miscUtilities.hpp>
@@ -130,8 +131,9 @@ try
     std::string deck_filename = param.get<std::string>("deck_filename");
 
     Opm::ParserPtr newParser(new Opm::Parser() );
-    Opm::DeckConstPtr deck = newParser->parseFile( deck_filename );
+    bool strict_parsing = param.getDefault("strict_parsing", true);
-    Opm::EclipseStateConstPtr eclipseState(new EclipseState(deck));
+    Opm::DeckConstPtr deck = newParser->parseFile(deck_filename, strict_parsing);
+    std::shared_ptr<EclipseState> eclipseState(new EclipseState(deck));
 
     // Grid init
     grid.reset(new Dune::CpGrid());
@@ -184,6 +186,8 @@ try
                     / state.surfacevol()[c*np + pu.phase_pos[Oil]];
             }
         }
+    } else if (deck->hasKeyword("EQUIL") && props->numPhases() == 3) {
+        OPM_THROW(std::logic_error, "sim_fibo_ad_cp does not support EQUIL initialization.");
     } else {
         initBlackoilStateFromDeck(grid->numCells(), &(grid->globalCell())[0],
                                   grid->numFaces(), AutoDiffGrid::faceCells(*grid),
@@ -205,9 +209,9 @@ try
 
     // Write parameters used for later reference.
     bool output = param.getDefault("output", true);
-    std::ofstream outStream;
     std::string output_dir;
     if (output) {
+        // Create output directory if needed.
         output_dir =
             param.getDefault("output_dir", std::string("output"));
         boost::filesystem::path fpath(output_dir);
@@ -217,19 +221,10 @@ try
         catch (...) {
             OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
         }
-        std::string filename = output_dir + "/timing.param";
+        // Write simulation parameters.
-        outStream.open(filename.c_str(), std::fstream::trunc | std::fstream::out);
-        // open file to clean it. The file is appended to in SimulatorTwophase
-        filename = output_dir + "/step_timing.param";
-        std::fstream step_os(filename.c_str(), std::fstream::trunc | std::fstream::out);
-        step_os.close();
         param.writeParam(output_dir + "/simulation.param");
     }
 
-    std::cout << "\n\n================    Starting main simulation loop     ===============\n"
-              << std::flush;
-
-    WellStateFullyImplicitBlackoil well_state;
     Opm::TimeMapConstPtr timeMap(eclipseState->getSchedule()->getTimeMap());
     SimulatorTimer simtimer;
 
@@ -238,64 +233,28 @@ try
 
     Opm::DerivedGeology geology(*grid, *new_props, eclipseState, grav);
 
-    SimulatorReport fullReport;
+    SimulatorFullyImplicitBlackoil<Dune::CpGrid> simulator(param,
-    for (size_t reportStepIdx = 0; reportStepIdx < timeMap->numTimesteps(); ++reportStepIdx) {
+                                                           *grid,
-        // Report on start of a report step.
+                                                           geology,
-        std::cout << "\n"
+                                                           *new_props,
-                  << "---------------------------------------------------------------\n"
+                                                           rock_comp->isActive() ? rock_comp.get() : 0,
-                  << "--------------    Starting report step " << reportStepIdx << "    --------------\n"
+                                                           *fis_solver,
-                  << "---------------------------------------------------------------\n"
+                                                           grav,
-                  << "\n";
+                                                           deck->hasKeyword("DISGAS"),
+                                                           deck->hasKeyword("VAPOIL"),
+                                                           eclipseState,
+                                                           outputWriter);
 
-        // Create new wells, well_state
+    std::cout << "\n\n================ Starting main simulation loop ===============\n"
-        WellsManager wells(eclipseState, reportStepIdx, Opm::UgGridHelpers::numCells(*grid), 
+              << std::flush;
-                           Opm::UgGridHelpers::globalCell(*grid),
-                           Opm::UgGridHelpers::cartDims(*grid),
-                           Opm::UgGridHelpers::dimensions(*grid),
-                           Opm::UgGridHelpers::beginCellCentroids(*grid),
-                           Opm::UgGridHelpers::cell2Faces(*grid),
-                           Opm::UgGridHelpers::beginFaceCentroids(*grid),
-                           props->permeability());
 
-        if (reportStepIdx == 0) {
+    SimulatorReport fullReport = simulator.run(simtimer, state);
-            // @@@ 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.setCurrentStepNum(reportStepIdx);
-
-        if (reportStepIdx == 0) {
-            outputWriter.writeInit(simtimer);
-            outputWriter.writeTimeStep(simtimer, state, well_state.basicWellState());
-        }
-
-        SimulatorFullyImplicitBlackoil<Dune::CpGrid> simulator(param,
-                                                               eclipseState->getSchedule(),
-                                                               *grid,
-                                                               geology,
-                                                               *new_props,
-                                                               rock_comp->isActive() ? rock_comp.get() : 0,
-                                                               wells.c_wells(),
-                                                               *fis_solver,
-                                                               grav,
-                                                               deck->hasKeyword("DISGAS"),
-                                                               deck->hasKeyword("VAPOIL") );
-
-        SimulatorReport episodeReport = simulator.run(simtimer, state, well_state);
-
-        ++simtimer;
-
-        outputWriter.writeTimeStep(simtimer, state, well_state.basicWellState());
-        fullReport += episodeReport;
-    }
 
     std::cout << "\n\n================    End of simulation     ===============\n\n";
     fullReport.report(std::cout);
 
     if (output) {
-        std::string filename = output_dir + "/walltime.param";
+        std::string filename = output_dir + "/walltime.txt";
         std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out);
         fullReport.reportParam(tot_os);
         warnIfUnusedParams(param);

opm/autodiff/BlackoilPropsAd.cpp

@@ -609,6 +609,23 @@ namespace Opm
         OPM_THROW(std::runtime_error, "Method rsMax() not implemented.");
     }
 
+    /// Bubble point curve for Rs as function of oil pressure.
+    /// \param[in]  po     Array of n oil pressure values.
+    /// \param[in]  so     Array of n oil saturation values.
+    /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+    /// \return            Array of n bubble point values for Rs.
+    V BlackoilPropsAd::rsSat(const V& po,
+                             const V& so,
+                             const Cells& cells) const
+    {
+        // Suppress warning about "unused parameters".
+        static_cast<void>(po);
+        static_cast<void>(so);
+        static_cast<void>(cells);
+
+        OPM_THROW(std::runtime_error, "Method rsMax() not implemented.");
+    }
+
     /// Bubble point curve for Rs as function of oil pressure.
     /// \param[in]  po     Array of n oil pressure values.
     /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
@@ -623,6 +640,23 @@ namespace Opm
         OPM_THROW(std::runtime_error, "Method rsMax() not implemented.");
     }
 
+    /// Bubble point curve for Rs as function of oil pressure.
+    /// \param[in]  po     Array of n oil pressure values.
+    /// \param[in]  so     Array of n oil saturation values.
+    /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+    /// \return            Array of n bubble point values for Rs.
+    ADB BlackoilPropsAd::rsSat(const ADB& po,
+                               const ADB& so,
+                               const Cells& cells) const
+    {
+        // Suppress warning about "unused parameters".
+        static_cast<void>(po);
+        static_cast<void>(so);
+        static_cast<void>(cells);
+
+        OPM_THROW(std::runtime_error, "Method rsMax() not implemented.");
+    }
+
     // ------ Rs bubble point curve ------
 
     /// Bubble point curve for Rs as function of oil pressure.
@@ -639,6 +673,23 @@ namespace Opm
         OPM_THROW(std::runtime_error, "Method rsMax() not implemented.");
     }
 
+    /// Bubble point curve for Rs as function of oil pressure.
+    /// \param[in]  po     Array of n oil pressure values.
+    /// \param[in]  so     Array of n oil saturation values.
+    /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+    /// \return            Array of n bubble point values for Rs.
+    V BlackoilPropsAd::rvSat(const V& po,
+                             const V& so,
+                             const Cells& cells) const
+    {
+        // Suppress warning about "unused parameters".
+        static_cast<void>(po);
+        static_cast<void>(so);
+        static_cast<void>(cells);
+
+        OPM_THROW(std::runtime_error, "Method rsMax() not implemented.");
+    }
+
     /// Bubble point curve for Rs as function of oil pressure.
     /// \param[in]  po     Array of n oil pressure values.
     /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
@@ -653,6 +704,23 @@ namespace Opm
         OPM_THROW(std::runtime_error, "Method rsMax() not implemented.");
     }
 
+    /// Bubble point curve for Rs as function of oil pressure.
+    /// \param[in]  po     Array of n oil pressure values.
+    /// \param[in]  so     Array of n oil saturation values.
+    /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+    /// \return            Array of n bubble point values for Rs.
+    ADB BlackoilPropsAd::rvSat(const ADB& po,
+                               const ADB& so,
+                               const Cells& cells) const
+    {
+        // Suppress warning about "unused parameters".
+        static_cast<void>(po);
+        static_cast<void>(so);
+        static_cast<void>(cells);
+
+        OPM_THROW(std::runtime_error, "Method rsMax() not implemented.");
+    }
+
     // ------ Relative permeability ------
 
     /// Relative permeabilities for all phases.
@@ -829,6 +897,11 @@ namespace Opm
         OPM_THROW(std::logic_error, "BlackoilPropsAd class does not support hysteresis.");
     }
     
+    /// Update for max oil saturation.
+    void BlackoilPropsAd::updateSatOilMax(const std::vector<double>& saturation)
+    {
+        OPM_THROW(std::logic_error, "BlackoilPropsAd class does not support this functionality.");
+    }
 
 } // namespace Opm
 

opm/autodiff/BlackoilPropsAd.hpp

@@ -254,6 +254,15 @@ namespace Opm
         V rsSat(const V& po,
                 const Cells& cells) const;
 
+        /// Solution gas/oil ratio and its derivatives at saturated condition as a function of p.
+        /// \param[in]  po     Array of n oil pressure values.
+        /// \param[in]  so     Array of n oil saturation values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n bubble point values for Rs.
+        V rsSat(const V& po,
+                const V& so,
+                const Cells& cells) const;
+
         /// Solution gas/oil ratio and its derivatives at saturated condition as a function of p.
         /// \param[in]  po     Array of n oil pressure values.
         /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
@@ -261,6 +270,15 @@ namespace Opm
         ADB rsSat(const ADB& po,
                   const Cells& cells) const;
 
+        /// Solution gas/oil ratio and its derivatives at saturated condition as a function of p.
+        /// \param[in]  po     Array of n oil pressure values.
+        /// \param[in]  so     Array of n oil saturation values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n bubble point values for Rs.
+        ADB rsSat(const ADB& po,
+                  const ADB& so,
+                  const Cells& cells) const;
+
         // ------ Rv condensation curve ------
 
         /// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
@@ -270,6 +288,15 @@ namespace Opm
         V rvSat(const V& po,
                 const Cells& cells) const;
 
+        /// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
+        /// \param[in]  po     Array of n oil pressure values.
+        /// \param[in]  so     Array of n oil saturation values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n bubble point values for Rs.
+        V rvSat(const V& po,
+                const V& so,
+                const Cells& cells) const;
+
         /// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
         /// \param[in]  po     Array of n oil pressure values.
         /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
@@ -277,6 +304,15 @@ namespace Opm
         ADB rvSat(const ADB& po,
                   const Cells& cells) const;
 
+        /// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
+        /// \param[in]  po     Array of n oil pressure values.
+        /// \param[in]  so     Array of n oil saturation values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n bubble point values for Rs.
+        ADB rvSat(const ADB& po,
+                  const ADB& so,
+                  const Cells& cells) const;
+
         // ------ Relative permeability ------
 
         /// Relative permeabilities for all phases.
@@ -320,6 +356,11 @@ namespace Opm
         /// \param[in]  cells       Array of n cell indices to be associated with the saturation values.
         void updateSatHyst(const std::vector<double>& saturation,
                            const std::vector<int>& cells);
+
+
+        /// Update for max oil saturation.
+        void updateSatOilMax(const std::vector<double>& saturation);
+
     private:
         const BlackoilPropertiesInterface& props_;
         PhaseUsage pu_;

opm/autodiff/BlackoilPropsAdFromDeck.cpp

@@ -188,6 +188,15 @@ namespace Opm
                 OPM_THROW(std::runtime_error, "Input is missing PVDG or PVTG\n");
             }
         }
+        // Oil vaporization controls (kw VAPPARS)
+        vap1_ = vap2_ = 0.0;
+        if (deck->hasKeyword("VAPPARS") && deck->hasKeyword("VAPOIL") && deck->hasKeyword("DISGAS")) {
+            vap1_ = deck->getKeyword("VAPPARS")->getRecord(0)->getItem(0)->getRawDouble(0);
+            vap2_ = deck->getKeyword("VAPPARS")->getRecord(0)->getItem(1)->getRawDouble(0);
+            satOilMax_.resize(number_of_cells, 0.0);
+        } else if (deck->hasKeyword("VAPPARS")) {
+            OPM_THROW(std::runtime_error, "Input has VAPPARS, but missing VAPOIL and/or DISGAS\n");
+        }
 
         SaturationPropsFromDeck<SatFuncGwsegNonuniform>* ptr
             = new SaturationPropsFromDeck<SatFuncGwsegNonuniform>();
@@ -199,6 +208,7 @@ namespace Opm
                       "Inconsistent number of phases in pvt data (" << phase_usage_.num_phases
                       << ") and saturation-dependent function data (" << satprops_->numPhases() << ").");
         }
+        vap_satmax_guard_ = 0.01;
     }
 
     ////////////////////////////
@@ -746,6 +756,20 @@ namespace Opm
         return rbub;
     }
 
+    /// Bubble point curve for Rs as function of oil pressure.
+    /// \param[in]  po     Array of n oil pressure values.
+    /// \param[in]  so     Array of n oil saturation values.
+    /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+    /// \return            Array of n bubble point values for Rs.
+    V BlackoilPropsAdFromDeck::rsSat(const V& po,
+                                     const V& so,
+                                     const Cells& cells) const
+    {
+        V rs = rsSat(po, cells);
+        applyVap(rs, so, cells, vap2_);
+        return rs;
+    }
+
     /// Bubble point curve for Rs as function of oil pressure.
     /// \param[in]  po     Array of n oil pressure values.
     /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
@@ -770,6 +794,20 @@ namespace Opm
         return ADB::function(rbub, jacs);
     }
 
+    /// Bubble point curve for Rs as function of oil pressure.
+    /// \param[in]  po     Array of n oil pressure values.
+    /// \param[in]  so     Array of n oil saturation values.
+    /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+    /// \return            Array of n bubble point values for Rs.
+    ADB BlackoilPropsAdFromDeck::rsSat(const ADB& po,
+                                       const ADB& so,
+                                       const Cells& cells) const
+    {
+        ADB rs = rsSat(po, cells);
+        applyVap(rs, so, cells, vap2_);
+        return rs;
+    }
+
     // ------ Condensation curve ------
 
     /// Condensation curve for Rv as function of oil pressure.
@@ -790,6 +828,20 @@ namespace Opm
         return rv;
     }
 
+    /// Condensation curve for Rv as function of oil pressure.
+    /// \param[in]  po     Array of n oil pressure values.
+    /// \param[in]  so     Array of n oil saturation values.
+    /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+    /// \return            Array of n bubble point values for Rs.
+    V BlackoilPropsAdFromDeck::rvSat(const V& po,
+                                     const V& so,
+                                     const Cells& cells) const
+    {
+        V rv = rvSat(po, cells);
+        applyVap(rv, so, cells, vap1_);
+        return rv;
+    }
+
     /// Condensation curve for Rv as function of oil pressure.
     /// \param[in]  po     Array of n oil pressure values.
     /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
@@ -814,6 +866,20 @@ namespace Opm
         return ADB::function(rv, jacs);
     }
 
+    /// Condensation curve for Rv as function of oil pressure.
+    /// \param[in]  po     Array of n oil pressure values.
+    /// \param[in]  so     Array of n oil saturation values.
+    /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+    /// \return            Array of n bubble point values for Rs.
+    ADB BlackoilPropsAdFromDeck::rvSat(const ADB& po,
+                                       const ADB& so,
+                                       const Cells& cells) const
+    {
+        ADB rv = rvSat(po, cells);
+        applyVap(rv, so, cells, vap1_);
+        return rv;
+    }
+
     // ------ Relative permeability ------
 
     /// Relative permeabilities for all phases.
@@ -988,5 +1054,75 @@ namespace Opm
         satprops_->updateSatHyst(n, cells.data(), saturation.data());
     }
     
+    /// Update for max oil saturation.
+    void BlackoilPropsAdFromDeck::updateSatOilMax(const std::vector<double>& saturation)
+    {
+        if (!satOilMax_.empty()) {
+            const int n = satOilMax_.size();
+            const int np = phase_usage_.num_phases;
+            const int posOil = phase_usage_.phase_pos[Oil];
+            const double* s = saturation.data();
+            for (int i=0; i<n; ++i) {
+                if (satOilMax_[i] < s[np*i+posOil]) {
+                    satOilMax_[i] = s[np*i+posOil];
+                }
+            }
+        }
+    }
+    
+    /// Apply correction to rs/rv according to kw VAPPARS
+    /// \param[in/out] r     Array of n rs/rv values.
+    /// \param[in]     so    Array of n oil saturation values.
+    /// \param[in]     cells Array of n cell indices to be associated with the r and so values.
+    /// \param[in]     vap   Correction parameter.
+    void BlackoilPropsAdFromDeck::applyVap(V& r,
+                                           const V& so,
+                                           const std::vector<int>& cells,
+                                           const double vap) const
+    {
+        if (!satOilMax_.empty() && vap > 0.0) { 
+            const int n = cells.size();
+            V factor = V::Ones(n, 1);
+            for (int i=0; i<n; ++i) {
+                if (satOilMax_[cells[i]] > vap_satmax_guard_ && so[i] < satOilMax_[cells[i]]) {
+                    factor[i] = std::pow(so[i]/satOilMax_[cells[i]], vap);
+                } 
+            }
+            r = factor*r;
+        }
+    }
+    
+    /// Apply correction to rs/rv according to kw VAPPARS
+    /// \param[in/out] r     Array of n rs/rv values.
+    /// \param[in]     so    Array of n oil saturation values.
+    /// \param[in]     cells Array of n cell indices to be associated with the r and so values.
+    /// \param[in]     vap   Correction parameter.
+    void BlackoilPropsAdFromDeck::applyVap(ADB& r,
+                                           const ADB& so,
+                                           const std::vector<int>& cells,
+                                           const double vap) const
+    {
+        if (!satOilMax_.empty() && vap > 0.0) { 
+            const int n = cells.size();
+            V factor = V::Ones(n, 1);
+            //V dfactor_dso = V::Zero(n, 1);  TODO: Consider effect of complete jacobian (including so-derivatives)
+            for (int i=0; i<n; ++i) {
+                if (satOilMax_[cells[i]] > vap_satmax_guard_ && so.value()[i] < satOilMax_[cells[i]]) {
+                    factor[i] = std::pow(so.value()[i]/satOilMax_[cells[i]], vap);
+                    //dfactor_dso[i] = vap*std::pow(so.value()[i]/satOilMax_[cells[i]], vap-1.0)/satOilMax_[cells[i]];
+                }
+            }
+            //ADB::M dfactor_dso_diag = spdiag(dfactor_dso);
+            //const int num_blocks = so.numBlocks();
+            //std::vector<ADB::M> jacs(num_blocks);
+            //for (int block = 0; block < num_blocks; ++block) {
+            //    jacs[block] = dfactor_dso_diag * so.derivative()[block];
+            //}
+            //r = ADB::function(factor, jacs)*r;
+            
+            r = factor*r;
+        }
+    }
+
 } // namespace Opm
 

opm/autodiff/BlackoilPropsAdFromDeck.hpp

@@ -277,6 +277,15 @@ namespace Opm
         V rsSat(const V& po,
                 const Cells& cells) const;
 
+        /// Bubble point curve for Rs as function of oil pressure.
+        /// \param[in]  po     Array of n oil pressure values.
+        /// \param[in]  so     Array of n oil saturation values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n bubble point values for Rs.
+        V rsSat(const V& po,
+                const V& so,
+                const Cells& cells) const;
+
         /// Bubble point curve for Rs as function of oil pressure.
         /// \param[in]  po     Array of n oil pressure values.
         /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
@@ -284,6 +293,15 @@ namespace Opm
         ADB rsSat(const ADB& po,
                   const Cells& cells) const;
 
+        /// Bubble point curve for Rs as function of oil pressure.
+        /// \param[in]  po     Array of n oil pressure values.
+        /// \param[in]  so     Array of n oil saturation values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n bubble point values for Rs.
+        ADB rsSat(const ADB& po,
+                  const ADB& so,
+                  const Cells& cells) const;
+
         // ------ Rv condensation curve ------
 
         /// Condensation curve for Rv as function of oil pressure.
@@ -293,6 +311,15 @@ namespace Opm
         V rvSat(const V& po,
                 const Cells& cells) const;
 
+        /// Condensation curve for Rv as function of oil pressure.
+        /// \param[in]  po     Array of n oil pressure values.
+        /// \param[in]  so     Array of n oil saturation values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n bubble point values for Rs.
+        V rvSat(const V& po,
+                const V& so,
+                const Cells& cells) const;
+
         /// Condensation curve for Rv as function of oil pressure.
         /// \param[in]  po     Array of n oil pressure values.
         /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
@@ -300,6 +327,15 @@ namespace Opm
         ADB rvSat(const ADB& po,
                   const Cells& cells) const;
 
+        /// Condensation curve for Rv as function of oil pressure.
+        /// \param[in]  po     Array of n oil pressure values.
+        /// \param[in]  so     Array of n oil saturation values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n bubble point values for Rs.
+        ADB rvSat(const ADB& po,
+                  const ADB& so,
+                  const Cells& cells) const;
+
         // ------ Relative permeability ------
 
         /// Relative permeabilities for all phases.
@@ -344,6 +380,9 @@ namespace Opm
         void updateSatHyst(const std::vector<double>& saturation,
                            const std::vector<int>& cells);
 
+        /// Update for max oil saturation.                  
+        void updateSatOilMax(const std::vector<double>& saturation);
+
     private:
         /// Initializes the properties.
         template <class CentroidIterator>
@@ -356,6 +395,17 @@ namespace Opm
                   int dimension,
                   const bool init_rock);
 
+        /// Correction to rs/rv according to kw VAPPARS
+        void applyVap(V& r,
+                      const V& so,
+                      const std::vector<int>& cells,
+                      const double vap) const;
+
+        void applyVap(ADB& r,
+                      const ADB& so,
+                      const std::vector<int>& cells,
+                      const double vap) const;
+
         RockFromDeck rock_;
         std::unique_ptr<SaturationPropsInterface> satprops_;
 
@@ -380,6 +430,13 @@ namespace Opm
         std::vector<int> pvtTableIdx_;
 
         std::vector<std::array<double, BlackoilPhases::MaxNumPhases> > densities_;
+        
+        // VAPPARS
+        double vap1_;
+        double vap2_;
+        std::vector<double> satOilMax_;
+        double vap_satmax_guard_;  //Threshold value to promote stability
+
     };
 } // namespace Opm
 

opm/autodiff/BlackoilPropsAdInterface.hpp

@@ -247,6 +247,16 @@ namespace Opm
         V rsSat(const V& po,
                 const Cells& cells) const = 0;
 
+        /// Bubble point curve for Rs as function of oil pressure.
+        /// \param[in]  po     Array of n oil pressure values.
+        /// \param[in]  so     Array of n oil saturation values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n bubble point values for Rs.
+        virtual
+        V rsSat(const V& po,
+                const V& so,
+                const Cells& cells) const = 0;
+
         /// Bubble point curve for Rs as function of oil pressure.
         /// \param[in]  po     Array of n oil pressure values.
         /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
@@ -255,6 +265,16 @@ namespace Opm
         ADB rsSat(const ADB& po,
                   const Cells& cells) const = 0;
 
+        /// Bubble point curve for Rs as function of oil pressure.
+        /// \param[in]  po     Array of n oil pressure values.
+        /// \param[in]  so     Array of n oil saturation values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n bubble point values for Rs.
+        virtual
+        ADB rsSat(const ADB& po,
+                  const ADB& so,
+                  const Cells& cells) const = 0;
+
         // ------ Rs bubble point curve ------
 
         /// Bubble point curve for Rs as function of oil pressure.
@@ -265,6 +285,16 @@ namespace Opm
         V rvSat(const V& po,
                 const Cells& cells) const = 0;
 
+        /// Bubble point curve for Rs as function of oil pressure.
+        /// \param[in]  po     Array of n oil pressure values.
+        /// \param[in]  so     Array of n oil saturation values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n bubble point values for Rs.
+        virtual
+        V rvSat(const V& po,
+                const V& so,
+                const Cells& cells) const = 0;
+
         /// Bubble point curve for Rs as function of oil pressure.
         /// \param[in]  po     Array of n oil pressure values.
         /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
@@ -273,6 +303,16 @@ namespace Opm
         ADB rvSat(const ADB& po,
                   const Cells& cells) const = 0;
 
+        /// Bubble point curve for Rs as function of oil pressure.
+        /// \param[in]  po     Array of n oil pressure values.
+        /// \param[in]  so     Array of n oil saturation values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n bubble point values for Rs.
+        virtual
+        ADB rvSat(const ADB& po,
+                  const ADB& so,
+                  const Cells& cells) const = 0;
+
         // ------ Relative permeability ------
 
         /// Relative permeabilities for all phases.
@@ -321,6 +361,10 @@ namespace Opm
         virtual
         void updateSatHyst(const std::vector<double>& saturation,
                            const std::vector<int>& cells) = 0;
+
+        /// Update for max oil saturation.
+        virtual                  
+        void updateSatOilMax(const std::vector<double>& saturation) = 0;
     };
 
 } // namespace Opm

opm/autodiff/CPRPreconditioner.hpp

@@ -37,6 +37,7 @@
 
 #include "reenable_warning_pragmas.h"
 
+#include <opm/core/utility/ErrorMacros.hpp>
 
 namespace Opm
 {
@@ -161,8 +162,8 @@ namespace Opm
             // Solve system.
             Dune::InverseOperatorResult result;
             linsolve.apply(x, de, result);
-            if (result.converged) {
+            if (!result.converged) {
-                // std::cout << "solveElliptic() successful!" << std::endl;
+                OPM_THROW(std::runtime_error, "CPRPreconditioner failed to solve elliptic subsystem.");
             }
             return x;
         }

opm/autodiff/FullyImplicitBlackoilSolver.hpp

@@ -254,18 +254,22 @@ namespace Opm {
 
         V
         fluidRsSat(const V&                p,
+                   const V&                so,
                    const std::vector<int>& cells) const;
 
         ADB
         fluidRsSat(const ADB&              p,
+                   const ADB&              so,
                    const std::vector<int>& cells) const;
 
         V
         fluidRvSat(const V&                p,
+                   const V&                so,
                    const std::vector<int>& cells) const;
 
         ADB
         fluidRvSat(const ADB&              p,
+                   const ADB&              so,
                    const std::vector<int>& cells) const;
 
         ADB

opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp

@@ -480,12 +480,12 @@ namespace {
 
             if (active_[ Gas]) {
                 // Define Sg Rs and Rv in terms of xvar.
-                const ADB& rsSat = fluidRsSat(state.pressure, cells_);
-                const ADB& rvSat = fluidRvSat(state.pressure, cells_);
                 const ADB& xvar = vars[ nextvar++ ];
                 const ADB& sg = isSg*xvar + isRv* so;
                 state.saturation[ pu.phase_pos[ Gas ] ] = sg;
                 so = so - sg;
+                const ADB rsSat = fluidRsSat(state.pressure, so, cells_);
+                const ADB rvSat = fluidRvSat(state.pressure, so, cells_);
 
                 if (has_disgas_) {
                     state.rs = (1-isRs) * rsSat + isRs*xvar;
@@ -587,18 +587,20 @@ namespace {
             const ADB bw = fluid_.bWat(perf_press, well_cells);
             b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw.value();
         }
+        assert(active_[Oil]);
+        const ADB perf_so =  subset(state.saturation[pu.phase_pos[Oil]], well_cells);
         if (pu.phase_used[BlackoilPhases::Liquid]) {
             const ADB perf_rs = subset(state.rs, well_cells);
             const ADB bo = fluid_.bOil(perf_press, perf_rs, perf_cond, well_cells);
             b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo.value();
-            const V rssat = fluidRsSat(perf_press.value(), well_cells);
+            const V rssat = fluidRsSat(perf_press.value(), perf_so.value(), well_cells);
             rssat_perf.assign(rssat.data(), rssat.data() + nperf);
         }
         if (pu.phase_used[BlackoilPhases::Vapour]) {
             const ADB perf_rv = subset(state.rv, well_cells);
             const ADB bg = fluid_.bGas(perf_press, perf_rv, perf_cond, well_cells);
             b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg.value();
-            const V rvsat = fluidRvSat(perf_press.value(), well_cells);
+            const V rvsat = fluidRvSat(perf_press.value(), perf_so.value(), well_cells);
             rvsat_perf.assign(rvsat.data(), rvsat.data() + nperf);
         }
         // b is row major, so can just copy data.
@@ -644,6 +646,15 @@ namespace {
         // Create the primary variables.
         SolutionState state = variableState(x, xw);
 
+        // DISKVAL(state.pressure);
+        // DISKVAL(state.saturation[0]);
+        // DISKVAL(state.saturation[1]);
+        // DISKVAL(state.saturation[2]);
+        // DISKVAL(state.rs);
+        // DISKVAL(state.rv);
+        // DISKVAL(state.qs);
+        // DISKVAL(state.bhp);
+
         // -------- Mass balance equations --------
 
         // Compute b_p and the accumulation term b_p*s_p for each phase,
@@ -1303,8 +1314,8 @@ namespace {
 
 
         // phase translation sg <-> rs
-        const V rsSat0 = fluidRsSat(p_old, cells_);
+        const V rsSat0 = fluidRsSat(p_old, s_old.col(pu.phase_pos[Oil]), cells_);
-        const V rsSat = fluidRsSat(p, cells_);
+        const V rsSat = fluidRsSat(p, so, cells_);
 
         std::fill(primalVariable_.begin(), primalVariable_.end(), PrimalVariables::Sg);
 
@@ -1328,8 +1339,8 @@ namespace {
         }
 
         // phase transitions so <-> rv
-        const V rvSat0 = fluidRvSat(p_old, cells_);
+        const V rvSat0 = fluidRvSat(p_old, s_old.col(pu.phase_pos[Oil]), cells_);
-        const V rvSat = fluidRvSat(p, cells_);
+        const V rvSat = fluidRvSat(p, so, cells_);
 
         if (has_vapoil_) {
             // The obvious case
@@ -1373,7 +1384,7 @@ namespace {
         for (int c = 0; c < nc; ++c) {
             if (ixw[c]) {
                 so[c] = so[c] / (1-sw[c]);
-                sg[c] = sg[c] / (1-so[c]);
+                sg[c] = sg[c] / (1-sw[c]);
                 sw[c] = 0;
             }
         }
@@ -1591,8 +1602,10 @@ namespace {
         for (; quantityIt != endQuantityIt; ++quantityIt) {
             const double quantityResid = (*quantityIt).value().matrix().norm();
             if (!std::isfinite(quantityResid)) {
+                const int trouble_phase = quantityIt - residual_.material_balance_eq.begin();
                 OPM_THROW(Opm::NumericalProblem,
-                          "Encountered a non-finite residual");
+                          "Encountered a non-finite residual in material balance equation "
+                          << trouble_phase);
             }
             globalNorm = std::max(globalNorm, quantityResid);
         }
@@ -1881,9 +1894,10 @@ namespace {
     template<class T>
     V
     FullyImplicitBlackoilSolver<T>::fluidRsSat(const V&                p,
+                                            const V&                satOil,
                                             const std::vector<int>& cells) const
     {
-        return fluid_.rsSat(p, cells);
+        return fluid_.rsSat(p, satOil, cells);
     }
 
 
@@ -1893,17 +1907,19 @@ namespace {
     template<class T>
     ADB
     FullyImplicitBlackoilSolver<T>::fluidRsSat(const ADB&              p,
+                                            const ADB&              satOil,
                                             const std::vector<int>& cells) const
     {
-        return fluid_.rsSat(p, cells);
+        return fluid_.rsSat(p, satOil, cells);
     }
 
     template<class T>
     V
     FullyImplicitBlackoilSolver<T>::fluidRvSat(const V&                p,
+                                            const V&              satOil,
                                             const std::vector<int>& cells) const
     {
-        return fluid_.rvSat(p, cells);
+        return fluid_.rvSat(p, satOil, cells);
     }
 
 
@@ -1913,9 +1929,10 @@ namespace {
     template<class T>
     ADB
     FullyImplicitBlackoilSolver<T>::fluidRvSat(const ADB&              p,
+                                            const ADB&              satOil,
                                             const std::vector<int>& cells) const
     {
-        return fluid_.rvSat(p, cells);
+        return fluid_.rvSat(p, satOil, cells);
     }
 
 

opm/autodiff/SimulatorFullyImplicitBlackoil.hpp

@@ -37,6 +37,8 @@ namespace Opm
     class SimulatorTimer;
     class BlackoilState;
     class WellStateFullyImplicitBlackoil;
+    class EclipseState;
+    class EclipseWriter;
     struct SimulatorReport;
 
     class Schedule;
@@ -72,17 +74,19 @@ namespace Opm
         /// \param[in] wells         Collection of wells.  Null if no wells.
         /// \param[in] linsolver     linear solver
         /// \param[in] gravity       if non-null, gravity vector
+        /// \param[in] eclipse_state
+        /// \param[in] output_writer
         SimulatorFullyImplicitBlackoil(const parameter::ParameterGroup& param,
-                                       ScheduleConstPtr schedule,
                                        const Grid& grid,
                                        const DerivedGeology& geo,
                                        BlackoilPropsAdInterface& props,
                                        const RockCompressibility* rock_comp_props,
-                                       const Wells* wells,
                                        NewtonIterationBlackoilInterface& linsolver,
                                        const double* gravity,
                                        const bool disgas,
-                                       const bool vapoil );
+                                       const bool vapoil,
+                                       std::shared_ptr<EclipseState> eclipse_state,
+                                       EclipseWriter& output_writer);
 
         /// Run the simulation.
         /// This will run succesive timesteps until timer.done() is true. It will
@@ -92,8 +96,7 @@ namespace Opm
         /// \param[in,out] well_state  state of wells: bhp, perforation rates
         /// \return                    simulation report, with timing data
         SimulatorReport run(SimulatorTimer& timer,
-                            BlackoilState& state,
+                            BlackoilState& state);
-                            WellStateFullyImplicitBlackoil& well_state);
 
     private:
         class Impl;

opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp

@@ -33,6 +33,7 @@
 #include <opm/core/well_controls.h>
 #include <opm/core/pressure/flow_bc.h>
 
+#include <opm/core/io/eclipse/EclipseWriter.hpp>
 #include <opm/core/simulator/SimulatorReport.hpp>
 #include <opm/core/simulator/SimulatorTimer.hpp>
 #include <opm/core/utility/StopWatch.hpp>
@@ -73,49 +74,54 @@ namespace Opm
     {
     public:
         Impl(const parameter::ParameterGroup& param,
-             ScheduleConstPtr schedule,
              const Grid& grid,
              const DerivedGeology& geo,
              BlackoilPropsAdInterface& props,
              const RockCompressibility* rock_comp_props,
-             const Wells* wells,
              NewtonIterationBlackoilInterface& linsolver,
              const double* gravity,
              bool has_disgas,
-             bool has_vapoil );
+             bool has_vapoil,
+             std::shared_ptr<EclipseState> eclipse_state,
+             EclipseWriter& output_writer);
 
         SimulatorReport run(SimulatorTimer& timer,
-                            BlackoilState& state,
+                            BlackoilState& state);
-                            WellStateFullyImplicitBlackoil& well_state);
 
     private:
         // Data.
-
         typedef RateConverter::
         SurfaceToReservoirVoidage< BlackoilPropsAdInterface,
                                    std::vector<int> > RateConverterType;
 
+        const parameter::ParameterGroup param_;
+
         // Parameters for output.
         bool output_;
         bool output_vtk_;
         std::string output_dir_;
         int output_interval_;
         // Observed objects.
-        ScheduleConstPtr schedule_;
         const Grid& grid_;
         BlackoilPropsAdInterface& props_;
         const RockCompressibility* rock_comp_props_;
-        std::shared_ptr<Wells> wells_;
         const double* gravity_;
         // Solvers
-        const DerivedGeology &geo_;
+        const DerivedGeology& geo_;
-        FullyImplicitBlackoilSolver<Grid> solver_;
+        NewtonIterationBlackoilInterface& solver_;
         // Misc. data
-        RateConverterType rateConverter_;
         std::vector<int> allcells_;
+        const bool has_disgas_;
+        const bool has_vapoil_;
+        // eclipse_state
+        std::shared_ptr<EclipseState> eclipse_state_;
+        // output_writer
+        EclipseWriter& output_writer_;
+        RateConverterType rateConverter_;
 
         void
         computeRESV(const std::size_t               step,
+                    const Wells*                    wells,
                     const BlackoilState&            x,
                     WellStateFullyImplicitBlackoil& xw);
     };
@@ -125,19 +131,20 @@ namespace Opm
 
     template<class T>
     SimulatorFullyImplicitBlackoil<T>::SimulatorFullyImplicitBlackoil(const parameter::ParameterGroup& param,
-                                                                      ScheduleConstPtr schedule,
                                                                    const Grid& grid,
                                                                    const DerivedGeology& geo,
                                                                    BlackoilPropsAdInterface& props,
                                                                    const RockCompressibility* rock_comp_props,
-                                                                   const Wells* wells,
                                                                    NewtonIterationBlackoilInterface& linsolver,
                                                                    const double* gravity,
                                                                    const bool has_disgas,
-                                                                   const bool has_vapoil )
+                                                                   const bool has_vapoil,
+                                                                   std::shared_ptr<EclipseState> eclipse_state,
+                                                                   EclipseWriter& output_writer)
 
     {
-        pimpl_.reset(new Impl(param, schedule, grid, geo, props, rock_comp_props, wells, linsolver, gravity, has_disgas, has_vapoil));
+        pimpl_.reset(new Impl(param, grid, geo, props, rock_comp_props, linsolver, gravity, has_disgas, has_vapoil,
+                              eclipse_state, output_writer));
     }
 
 
@@ -146,10 +153,9 @@ namespace Opm
 
     template<class T>
     SimulatorReport SimulatorFullyImplicitBlackoil<T>::run(SimulatorTimer& timer,
-                                                        BlackoilState& state,
+                                                        BlackoilState& state)
-                                                        WellStateFullyImplicitBlackoil& well_state)
     {
-        return pimpl_->run(timer, state, well_state);
+        return pimpl_->run(timer, state);
     }
 
 
@@ -214,25 +220,28 @@ namespace Opm
     // \TODO: Treat bcs.
     template<class T>
     SimulatorFullyImplicitBlackoil<T>::Impl::Impl(const parameter::ParameterGroup& param,
-                                                  ScheduleConstPtr schedule,
                                                const Grid& grid,
                                                const DerivedGeology& geo,
                                                BlackoilPropsAdInterface& props,
                                                const RockCompressibility* rock_comp_props,
-                                               const Wells* wells,
                                                NewtonIterationBlackoilInterface& linsolver,
                                                const double* gravity,
                                                const bool has_disgas,
-                                               const bool has_vapoil)
+                                               const bool has_vapoil,
-        : schedule_(schedule)
+                                               std::shared_ptr<EclipseState> eclipse_state,
-        , grid_(grid),
+                                               EclipseWriter& output_writer)
+        : param_(param),
+          grid_(grid),
           props_(props),
           rock_comp_props_(rock_comp_props),
-          wells_(clone_wells(wells), & destroy_wells),
           gravity_(gravity),
           geo_(geo),
-          solver_(param, grid_, props_, geo_, rock_comp_props, *wells_, linsolver, has_disgas, has_vapoil)
+          solver_(linsolver),
-        , rateConverter_(props_, std::vector<int>(AutoDiffGrid::numCells(grid_), 0))
+          has_disgas_(has_disgas),
+          has_vapoil_(has_vapoil),
+          eclipse_state_(eclipse_state),
+          output_writer_(output_writer),
+          rateConverter_(props_, std::vector<int>(AutoDiffGrid::numCells(grid_), 0))
     {
         // For output.
         output_ = param.getDefault("output", true);
@@ -260,93 +269,100 @@ namespace Opm
 
     template<class T>
     SimulatorReport SimulatorFullyImplicitBlackoil<T>::Impl::run(SimulatorTimer& timer,
-                                                              BlackoilState& state,
+                                                                 BlackoilState& state)
-                                                              WellStateFullyImplicitBlackoil& well_state)
     {
-        // Initialisation.
+        WellStateFullyImplicitBlackoil prev_well_state;
-        std::vector<double> porevol;
-        if (rock_comp_props_ && rock_comp_props_->isActive()) {
-            computePorevolume(AutoDiffGrid::numCells(grid_), AutoDiffGrid::beginCellVolumes(grid_), props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
-        } else {
-            computePorevolume(AutoDiffGrid::numCells(grid_), AutoDiffGrid::beginCellVolumes(grid_), props_.porosity(), porevol);
-        }
-        // const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
-        std::vector<double> initial_porevol = porevol;
 
-        // Main simulation loop.
+        // Create timers and file for writing timing info.
         Opm::time::StopWatch solver_timer;
         double stime = 0.0;
         Opm::time::StopWatch step_timer;
         Opm::time::StopWatch total_timer;
         total_timer.start();
-        std::vector<double> fractional_flows;
+        std::string tstep_filename = output_dir_ + "/step_timing.txt";
-        std::vector<double> well_resflows_phase;
+        std::ofstream tstep_os(tstep_filename.c_str());
-        if (wells_) {
+
-            well_resflows_phase.resize((wells_->number_of_phases)*(wells_->number_of_wells), 0.0);
+        // Main simulation loop.
-        }
-        std::fstream tstep_os;
-        if (output_) {
-            std::string filename = output_dir_ + "/step_timing.param";
-            tstep_os.open(filename.c_str(), std::fstream::out | std::fstream::app);
-        }
         while (!timer.done()) {
-            // Report timestep and (optionally) write state to disk.
+            // Report timestep.
             step_timer.start();
             timer.report(std::cout);
+
+            // Create wells and well state.
+            WellsManager wells_manager(eclipse_state_,
+                                       timer.currentStepNum(),
+                                       Opm::UgGridHelpers::numCells(grid_),
+                                       Opm::UgGridHelpers::globalCell(grid_),
+                                       Opm::UgGridHelpers::cartDims(grid_),
+                                       Opm::UgGridHelpers::dimensions(grid_),
+                                       Opm::UgGridHelpers::beginCellCentroids(grid_),
+                                       Opm::UgGridHelpers::cell2Faces(grid_),
+                                       Opm::UgGridHelpers::beginFaceCentroids(grid_),
+                                       props_.permeability());
+            const Wells* wells = wells_manager.c_wells();
+            WellStateFullyImplicitBlackoil well_state;
+            well_state.init(wells, state);
+            if (timer.currentStepNum() != 0) {
+                // Transfer previous well state to current.
+                well_state.partialCopy(prev_well_state, *wells, prev_well_state.numWells());
+            }
+
+            // Output state at start of time step.
             if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
                 if (output_vtk_) {
                     outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
                 }
                 outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
                 outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_);
-
+            }
+            if (output_) {
+                if (timer.currentStepNum() == 0) {
+                    output_writer_.writeInit(timer);
+                }
+                output_writer_.writeTimeStep(timer, state, well_state.basicWellState());
             }
 
-            SimulatorReport sreport;
+            // Max oil saturation (for VPPARS), hysteresis update.
-            {
+            props_.updateSatOilMax(state.saturation());
-                computeRESV(timer.currentStepNum(), state, well_state);
-
-                // Run solver.
-                solver_timer.start();
-                solver_.step(timer.currentStepLength(), state, well_state);
-
-                // Stop timer and report.
-                solver_timer.stop();
-                const double st = solver_timer.secsSinceStart();
-                std::cout << "Fully implicit solver took:  " << st << " seconds." << std::endl;
-
-                stime += st;
-                sreport.pressure_time = st;
-            }
-
-            // Update pore volumes if rock is compressible.
-            if (rock_comp_props_ && rock_comp_props_->isActive()) {
-                initial_porevol = porevol;
-                computePorevolume(AutoDiffGrid::numCells(grid_), AutoDiffGrid::beginCellVolumes(grid_), props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
-            }
-
-            // Hysteresis
             props_.updateSatHyst(state.saturation(), allcells_);
 
-            sreport.total_time =  step_timer.secsSinceStart();
+            // Compute reservoir volumes for RESV controls.
-            if (output_) {
+            computeRESV(timer.currentStepNum(), wells, state, well_state);
-                sreport.reportParam(tstep_os);
 
-                if (output_vtk_) {
+            // Run a single step of the solver.
-                    outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
+            solver_timer.start();
-                }
+            FullyImplicitBlackoilSolver<T> solver(param_, grid_, props_, geo_, rock_comp_props_, *wells, solver_, has_disgas_, has_vapoil_);
-                outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
+            solver.step(timer.currentStepLength(), state, well_state);
-                outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_);
+            solver_timer.stop();
-                tstep_os.close();
+
+            // Report timing.
+            const double st = solver_timer.secsSinceStart();
+            std::cout << "Fully implicit solver took: " << st << " seconds." << std::endl;
+            stime += st;
+            if (output_) {
+                SimulatorReport step_report;
+                step_report.pressure_time = st;
+                step_report.total_time =  step_timer.secsSinceStart();
+                step_report.reportParam(tstep_os);
             }
 
-            // advance to next timestep before reporting at this location
+            // Increment timer, remember well state.
-            // ++timer; // Commented out since this has temporarily moved to the main() function.
+            ++timer;
-            break; // this is a temporary measure
+            prev_well_state = well_state;
         }
 
-        total_timer.stop();
+        // Write final simulation state.
+        if (output_) {
+            if (output_vtk_) {
+                outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
+            }
+            outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
+            outputWellStateMatlab(prev_well_state, timer.currentStepNum(), output_dir_);
+            output_writer_.writeTimeStep(timer, state, prev_well_state.basicWellState());
+        }
 
+        // Stop timer and create timing report
+        total_timer.stop();
         SimulatorReport report;
         report.pressure_time = stime;
         report.transport_time = 0.0;
@@ -470,16 +486,17 @@ namespace Opm
     void
     SimulatorFullyImplicitBlackoil<T>::
     Impl::computeRESV(const std::size_t               step,
+                      const Wells*                    wells,
                       const BlackoilState&            x,
                       WellStateFullyImplicitBlackoil& xw)
     {
         typedef SimFIBODetails::WellMap WellMap;
 
-        const std::vector<WellConstPtr>& w_ecl = schedule_->getWells(step);
+        const std::vector<WellConstPtr>& w_ecl = eclipse_state_->getSchedule()->getWells(step);
         const WellMap& wmap = SimFIBODetails::mapWells(w_ecl);
 
         const std::vector<int>& resv_prod =
-            SimFIBODetails::resvProducers(*wells_, step, wmap);
+            SimFIBODetails::resvProducers(*wells, step, wmap);
 
         if (! resv_prod.empty()) {
             const PhaseUsage&                    pu = props_.phaseUsage();
@@ -495,7 +512,7 @@ namespace Opm
                      rp = resv_prod.begin(), e = resv_prod.end();
                  rp != e; ++rp)
             {
-                WellControls* ctrl = wells_->ctrls[*rp];
+                WellControls* ctrl = wells->ctrls[*rp];
 
                 // RESV control mode, all wells
                 {
@@ -519,8 +536,8 @@ namespace Opm
 
                 // RESV control, WCONHIST wells.  A bit of duplicate
                 // work, regrettably.
-                if (wells_->name[*rp] != 0) {
+                if (wells->name[*rp] != 0) {
-                    WellMap::const_iterator i = wmap.find(wells_->name[*rp]);
+                    WellMap::const_iterator i = wmap.find(wells->name[*rp]);
 
                     if (i != wmap.end()) {
                         WellConstPtr wp = i->second;

opm/autodiff/WellStateFullyImplicitBlackoil.hpp

@@ -24,6 +24,7 @@
 #include <opm/core/wells.h>
 #include <opm/core/well_controls.h>
 #include <opm/core/simulator/WellState.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
 #include <vector>
 #include <cassert>
 
@@ -102,6 +103,53 @@ namespace Opm
             return basic_well_state_;
         }
 
+        /// The number of wells present.
+        int numWells() const
+        {
+            return bhp().size();
+        }
+
+        /// The number of phases present.
+        int numPhases() const
+        {
+            return wellRates().size() / numWells();
+        }
+
+        /// Copy data for the first num_wells_to_copy from source,
+        /// overwriting any data in this object associated with those
+        /// wells. Assumes that the number of phases are the same,
+        /// that the number of perforations associated with the wells
+        /// is unchanging, and that both objects contain at least
+        /// num_wells_to_copy wells.
+        void partialCopy(const WellStateFullyImplicitBlackoil& source,
+                         const Wells& wells,
+                         const int num_wells_to_copy)
+        {
+            if (numPhases() != source.numPhases()) {
+                OPM_THROW(std::logic_error, "partialCopy(): source and destination have different number of phases.");
+            }
+            if (num_wells_to_copy > numWells() || num_wells_to_copy > source.numWells()) {
+                OPM_THROW(std::logic_error, "partialCopy(): trying to copy too many wells.");
+            }
+            // bhp
+            std::copy(source.bhp().begin(),
+                      source.bhp().begin() + num_wells_to_copy,
+                      bhp().begin());
+            // wellRates
+            std::copy(source.wellRates().begin(),
+                      source.wellRates().begin() + numPhases()*num_wells_to_copy,
+                      wellRates().begin());
+            // perfPhaseRates
+            const int num_perfs_to_copy = wells.well_connpos[num_wells_to_copy];
+            std::copy(source.perfPhaseRates().begin(),
+                      source.perfPhaseRates().begin() + numPhases()*num_perfs_to_copy,
+                      perfPhaseRates().begin());
+            // currentControls
+            std::copy(source.currentControls().begin(),
+                      source.currentControls().begin() + num_wells_to_copy,
+                      currentControls().begin());
+        }
+
     private:
         WellState basic_well_state_;
         std::vector<double> perfphaserates_;