Updates FIBOPOLYMER simulator based on flow simulator.

Rename sim_poly_fibo_ad to flow_polymer.

examples/sim_poly_fibo_ad.cpp

@@ -1,252 +0,0 @@
-/*
-  Copyright 2014 SINTEF ICT, Applied Mathematics.
-  Copyright 2014 STATOIL ASA.
-
-  This file is part of the Open Porous Media project (OPM).
-
-  OPM is free software: you can redistribute it and/or modify
-  it under the terms of the GNU General Public License as published by
-  the Free Software Foundation, either version 3 of the License, or
-  (at your option) any later version.
-
-  OPM is distributed in the hope that it will be useful,
-  but WITHOUT ANY WARRANTY; without even the implied warranty of
-  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-  GNU General Public License for more details.
-
-  You should have received a copy of the GNU General Public License
-  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
-*/
-#include "config.h"
-
-#include <opm/core/pressure/FlowBCManager.hpp>
-
-#include <opm/core/grid.h>
-#include <opm/core/grid/GridManager.hpp>
-#include <opm/core/wells.h>
-#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>
-#include <opm/core/utility/parameters/ParameterGroup.hpp>
-#include <opm/core/utility/thresholdPressures.hpp>
-
-#include <opm/core/io/eclipse/EclipseWriter.hpp>
-#include <opm/core/props/BlackoilPropertiesBasic.hpp>
-#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
-#include <opm/core/props/rock/RockCompressibility.hpp>
-
-#include <opm/core/linalg/LinearSolverFactory.hpp>
-#include <opm/autodiff/NewtonIterationBlackoilSimple.hpp>
-#include <opm/autodiff/NewtonIterationBlackoilCPR.hpp>
-
-#include <opm/polymer/PolymerBlackoilState.hpp>
-#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
-
-#include <opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp>
-#include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
-#include <opm/polymer/PolymerProperties.hpp>
-#include <opm/polymer/PolymerInflow.hpp>
-#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
-#include <opm/autodiff/BlackoilPropsAdInterface.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>
-#include <numeric>
-
-
-namespace
-{
-    void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param)
-    {
-        if (param.anyUnused()) {
-            std::cout << "--------------------   Unused parameters:   --------------------\n";
-            param.displayUsage();
-            std::cout << "----------------------------------------------------------------" << std::endl;
-        }
-    }
-} // anon namespace
-
-
-
-// ----------------- Main program -----------------
-int
-main(int argc, char** argv)
-try
-{
-    using namespace Opm;
-
-    std::cout << "\n================    Test program for fully implicit three-phase black-oil flow     ===============\n\n";
-    parameter::ParameterGroup param(argc, argv, false);
-    std::cout << "---------------    Reading parameters     ---------------" << std::endl;
-
-    // If we have a "deck_filename", grid and props will be read from that.
-    bool use_deck = param.has("deck_filename");
-    if (!use_deck) {
-        OPM_THROW(std::runtime_error, "This program must be run with an input deck. "
-                  "Specify the deck with deck_filename=deckname.data (for example).");
-    }
-    std::shared_ptr<GridManager> grid;
-    std::shared_ptr<BlackoilPropertiesInterface> props;
-    std::shared_ptr<BlackoilPropsAdInterface> new_props;
-    std::shared_ptr<RockCompressibility> rock_comp;
-    PolymerBlackoilState state;
-    // bool check_well_controls = false;
-    // int max_well_control_iterations = 0;
-    double gravity[3] = { 0.0 };
-    std::string deck_filename = param.get<std::string>("deck_filename");
-
-    Opm::ParserPtr newParser(new Opm::Parser() );
-    Opm::DeckConstPtr deck = newParser->parseFile(deck_filename);
-    std::shared_ptr<EclipseState> eclipseState(new EclipseState(deck));
-
-    // Grid init
-    std::vector<double> porv;
-    if (eclipseState->hasDoubleGridProperty("PORV")) {
-        porv = eclipseState->getDoubleGridProperty("PORV")->getData();
-    }
-    grid.reset(new GridManager(eclipseState->getEclipseGrid(), porv));
-    auto &cGrid = *grid->c_grid();
-    const PhaseUsage pu = Opm::phaseUsageFromDeck(deck);
-    Opm::EclipseWriter outputWriter(param,
-                                    eclipseState,
-                                    pu,
-                                    cGrid.number_of_cells,
-                                    cGrid.global_cell);
-
-    // Rock and fluid init
-    props.reset(new BlackoilPropertiesFromDeck(deck, eclipseState, *grid->c_grid(), param));
-    new_props.reset(new BlackoilPropsAdFromDeck(deck, eclipseState, *grid->c_grid()));
-    const bool polymer = deck->hasKeyword("POLYMER");
-    const bool use_wpolymer = deck->hasKeyword("WPOLYMER");
-    PolymerProperties polymer_props(deck, eclipseState);
-    PolymerPropsAd polymer_props_ad(polymer_props);
-    // check_well_controls = param.getDefault("check_well_controls", false);
-    // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
-    // Rock compressibility.
-    rock_comp.reset(new RockCompressibility(deck, eclipseState));
-
-    // Gravity.
-    gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
-
-    // Init state variables (saturation and pressure).
-    if (param.has("init_saturation")) {
-        initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
-        initBlackoilSurfvol(*grid->c_grid(), *props, state);
-        enum { Oil = BlackoilPhases::Liquid, Gas = BlackoilPhases::Vapour };
-        if (pu.phase_used[Oil] && pu.phase_used[Gas]) {
-            const int np = props->numPhases();
-            const int nc = grid->c_grid()->number_of_cells;
-            for (int c = 0; c < nc; ++c) {
-                state.gasoilratio()[c] = state.surfacevol()[c*np + pu.phase_pos[Gas]]
-                    / state.surfacevol()[c*np + pu.phase_pos[Oil]];
-            }
-        }
-    } else if (deck->hasKeyword("EQUIL") && props->numPhases() == 3) {
-        state.init(*grid->c_grid(), props->numPhases());
-        const double grav = param.getDefault("gravity", unit::gravity);
-        initStateEquil(*grid->c_grid(), *props, deck, eclipseState, grav, state.blackoilState());
-        state.faceflux().resize(grid->c_grid()->number_of_faces, 0.0);
-    } else {
-        state.init(*grid->c_grid(), props->numPhases());
-        initBlackoilStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state.blackoilState());
-    }
-
-    bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
-    const double *grav = use_gravity ? &gravity[0] : 0;
-
-    // Solver for Newton iterations.
-    std::unique_ptr<NewtonIterationBlackoilInterface> fis_solver;
-    if (param.getDefault("use_cpr", true)) {
-        fis_solver.reset(new NewtonIterationBlackoilCPR(param));
-    } else {
-        fis_solver.reset(new NewtonIterationBlackoilSimple(param));
-    }
-
-    // Write parameters used for later reference.
-    bool output = param.getDefault("output", true);
-    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);
-        try {
-            create_directories(fpath);
-        }
-        catch (...) {
-            OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
-        }
-        // Write simulation parameters.
-        param.writeParam(output_dir + "/simulation.param");
-    }
-
-    Opm::TimeMapConstPtr timeMap(eclipseState->getSchedule()->getTimeMap());
-    SimulatorTimer simtimer;
-
-    // initialize variables
-    simtimer.init(timeMap);
-    if (polymer){
-        if (!use_wpolymer) {
-            OPM_MESSAGE("Warning: simulate polymer injection without WPOLYMER.");
-        } else {
-            if (param.has("polymer_start_days")) {
-                OPM_MESSAGE("Warning: Using WPOLYMER to control injection since it was found in deck."
-                            "You seem to be trying to control it via parameter poly_start_days (etc.) as well.");
-            }
-        }
-    } else {
-        if (use_wpolymer) {
-            OPM_MESSAGE("Warning: use WPOLYMER in a non-polymer scenario.");
-        }
-    }
-    std::cout << "\n\n================ Starting main simulation loop ===============\n"
-              << std::flush;
-    SimulatorReport fullReport;
-    Opm::DerivedGeology geology(*grid->c_grid(), *new_props, eclipseState, grav);
-
-    std::vector<double> threshold_pressures = thresholdPressures(eclipseState, *grid->c_grid());
-    SimulatorFullyImplicitBlackoilPolymer<UnstructuredGrid> simulator(param,
-                                             *grid->c_grid(),
-                                             geology,
-                                             *new_props,
-                                             polymer_props_ad,
-                                             rock_comp->isActive() ? rock_comp.get() : 0,
-                                             *fis_solver,
-                                             grav,
-                                             deck->hasKeyword("DISGAS"),
-                                             deck->hasKeyword("VAPOIL"),
-                                             polymer,
-                                             eclipseState,
-                                             outputWriter,
-                                             deck,
-                                             threshold_pressures);
-
-
-    fullReport = simulator.run(simtimer, state);
-    std::cout << "\n\n================    End of simulation     ===============\n\n";
-    fullReport.report(std::cout);
-
-    if (output) {
-        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);
-    }
-}
-catch (const std::exception &e) {
-    std::cerr << "Program threw an exception: " << e.what() << "\n";
-    throw;
-}
-

opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver.hpp

@@ -1,7 +1,6 @@
 /*
   Copyright 2013 SINTEF ICT, Applied Mathematics.
   Copyright 2014 STATOIL ASA.
-  
 
   This file is part of the Open Porous Media project (OPM).
 
@@ -22,6 +21,8 @@
 #ifndef OPM_FULLYIMPLICITBLACKOILPOLYMERSOLVER_HEADER_INCLUDED
 #define OPM_FULLYIMPLICITBLACKOILPOLYMERSOLVER_HEADER_INCLUDED
 
+#include <cassert>
+
 #include <opm/autodiff/AutoDiffBlock.hpp>
 #include <opm/autodiff/AutoDiffHelpers.hpp>
 #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
@@ -30,6 +31,8 @@
 #include <opm/polymer/PolymerProperties.hpp>
 #include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
 
+#include <array>
+
 struct UnstructuredGrid;
 struct Wells;
 
@@ -56,29 +59,55 @@ namespace Opm {
     class FullyImplicitBlackoilPolymerSolver
     {
     public:
+        // the Newton relaxation type
+        enum RelaxType { DAMPEN, SOR };
+
+        // class holding the solver parameters
+        struct SolverParameter
+        {
+            double                          dp_max_rel_;
+            double                          ds_max_;
+            double                          dr_max_rel_;
+            enum RelaxType                  relax_type_;
+            double                          relax_max_;
+            double                          relax_increment_;
+            double                          relax_rel_tol_;
+            double                          max_residual_allowed_;
+            double                          tolerance_mb_;
+            double                          tolerance_cnv_;
+            double                          tolerance_wells_;
+            int                             max_iter_;
+
+            SolverParameter( const parameter::ParameterGroup& param );
+            SolverParameter();
+
+            void reset();
+        };
+
         /// \brief The type of the grid that we use.
         typedef T Grid;
         /// Construct a solver. It will retain references to the
         /// arguments of this functions, and they are expected to
         /// remain in scope for the lifetime of the solver.
-        /// \param[in] param       parameters
+        /// \param[in] param            parameters
         /// \param[in] grid             grid data structure
         /// \param[in] fluid            fluid properties
         /// \param[in] geo              rock properties
         /// \param[in] rock_comp_props  if non-null, rock compressibility properties
         /// \param[in] wells            well structure
         /// \param[in] linsolver        linear solver
-        FullyImplicitBlackoilPolymerSolver(const parameter::ParameterGroup& param,
+        FullyImplicitBlackoilPolymerSolver(const SolverParameter&          param,
                                            const Grid&                     grid ,
                                            const BlackoilPropsAdInterface& fluid,
                                            const DerivedGeology&           geo  ,
                                            const RockCompressibility*      rock_comp_props,
                                            const PolymerPropsAd&           polymer_props_ad,
-                                           const Wells&                    wells,
+                                           const Wells*                    wells,
                                            const NewtonIterationBlackoilInterface& linsolver,
                                            const bool has_disgas,
                                            const bool has_vapoil,
-                                           const bool has_polymer);
+                                           const bool has_polymer,
+                                           const bool terminal_output);
 
         /// \brief Set threshold pressures that prevent or reduce flow.
         /// This prevents flow across faces if the potential
@@ -99,12 +128,16 @@ namespace Opm {
         /// \param[in] dt        time step size
         /// \param[in] state     reservoir state
         /// \param[in] wstate    well state
-        void
+        /// \return              number of linear iterations used
+        int
         step(const double   dt    ,
              PolymerBlackoilState& state ,
              WellStateFullyImplicitBlackoil& wstate,
              const std::vector<double>& polymer_inflow);
 
+        unsigned int newtonIterations () const { return newtonIterations_; }
+        unsigned int linearIterations () const { return linearIterations_; }
+
     private:
         // Types and enums
         typedef AutoDiffBlock<double> ADB;
@@ -133,21 +166,23 @@ namespace Opm {
             ADB              rv;
             ADB              concentration;
             ADB              qs;
-            ADB              bhp;       
+            ADB              bhp;
+            // Below are quantities stored in the state for optimization purposes.
+            std::vector<ADB> canonical_phase_pressures; // Always has 3 elements, even if only 2 phases active.
         };
 
         struct WellOps {
-            WellOps(const Wells& wells);
+            WellOps(const Wells* wells);
             M w2p;              // well -> perf (scatter)
             M p2w;              // perf -> well (gather)
         };
 
-        enum { Water = BlackoilPropsAdInterface::Water,
+        enum { Water        = BlackoilPropsAdInterface::Water,
-               Oil   = BlackoilPropsAdInterface::Oil  ,
+               Oil          = BlackoilPropsAdInterface::Oil  ,
-               Gas   = BlackoilPropsAdInterface::Gas  };
+               Gas          = BlackoilPropsAdInterface::Gas  ,
+               MaxNumPhases = BlackoilPropsAdInterface::MaxNumPhases
+         };
 
-        // the Newton relaxation type
-        enum RelaxType { DAMPEN, SOR };
         enum PrimalVariables { Sg = 0, RS = 1, RV = 2 };
 
         // Member data
@@ -155,8 +190,8 @@ namespace Opm {
         const BlackoilPropsAdInterface& fluid_;
         const DerivedGeology&           geo_;
         const RockCompressibility*      rock_comp_props_;
-        const PolymerPropsAd&           polymer_props_ad_; 
+        const PolymerPropsAd&           polymer_props_ad_;
-        const Wells&                    wells_;
+        const Wells*                    wells_;
         const NewtonIterationBlackoilInterface&    linsolver_;
         // For each canonical phase -> true if active
         const std::vector<bool>         active_;
@@ -170,14 +205,8 @@ namespace Opm {
         const bool has_vapoil_;
         const bool has_polymer_;
         const int  poly_pos_;
-        double                          dp_max_rel_;
+
-        double                          ds_max_;
+        SolverParameter                 param_;
-        double                          drs_max_rel_;
-        enum RelaxType                  relax_type_;
-        double                          relax_max_;
-        double                          relax_increment_;
-        double                          relax_rel_tol_;
-        int                             max_iter_;
         bool use_threshold_pressure_;
         V threshold_pressures_by_interior_face_;
 
@@ -187,16 +216,30 @@ namespace Opm {
 
         LinearisedBlackoilResidual residual_;
 
+        /// \brief Whether we print something to std::cout
+        bool terminal_output_;
+        unsigned int newtonIterations_;
+        unsigned int linearIterations_;
+
         std::vector<int>         primalVariable_;
 
         // Private methods.
+
+        // return true if wells are available
+        bool wellsActive() const { return wells_ ? wells_->number_of_wells > 0 : false ; }
+        // return wells object
+        const Wells& wells () const { assert( bool(wells_ != 0) ); return *wells_; }
+
         SolutionState
         constantState(const PolymerBlackoilState& x,
-                      const WellStateFullyImplicitBlackoil& xw);
+                      const WellStateFullyImplicitBlackoil& xw) const;
+
+        void
+        makeConstantState(SolutionState& state) const;
 
         SolutionState
         variableState(const PolymerBlackoilState& x,
-                      const WellStateFullyImplicitBlackoil& xw);
+                      const WellStateFullyImplicitBlackoil& xw) const;
 
         void
         computeAccum(const SolutionState& state,
@@ -223,6 +266,7 @@ namespace Opm {
         void
         assemble(const V&             dtpv,
                  const PolymerBlackoilState& x,
+                 const bool initial_assembly,
                  WellStateFullyImplicitBlackoil& xw,
                  const std::vector<double>& polymer_inflow);
 
@@ -235,6 +279,18 @@ namespace Opm {
         std::vector<ADB>
         computePressures(const SolutionState& state) const;
 
+        std::vector<ADB>
+        computePressures(const ADB& po,
+                         const ADB& sw,
+                         const ADB& so,
+                         const ADB& sg) const;
+
+        V
+        computeGasPressure(const V& po,
+                           const V& sw,
+                           const V& so,
+                           const V& sg) const;
+
         std::vector<ADB>
         computeRelPerm(const SolutionState& state) const;
 
@@ -244,35 +300,26 @@ namespace Opm {
                             const std::vector<int>& well_cells) const;
 
         void
-        computeMassFlux(const int               actph ,
+        computeMassFlux(const V&                transi,
-                        const V&                transi,
+                        const std::vector<ADB>& kr    ,
-                        const ADB&              kr    ,
-                        const ADB&              p     ,
-                        const SolutionState&    state );
-
-        void
-        computeMassFlux(const V&                trans,
-                        const std::vector<ADB>& kr,
                         const std::vector<ADB>& phasePressure,
-                        const SolutionState&    state);
+                        const SolutionState&    state );
-
         void
         computeCmax(PolymerBlackoilState& state);
 
         ADB
         computeMc(const SolutionState& state) const;
 
-        ADB
-        rockPorosity(const ADB& p) const;
-        
-        ADB
-        rockPermeability(const ADB& p) const;
         void applyThresholdPressures(ADB& dp);
 
         double
         residualNorm() const;
 
-        std::vector<double> residuals() const;
+        /// \brief Compute the residual norms of the mass balance for each phase,
+        /// the well flux, and the well equation.
+        /// \return a vector that contains for each phase the norm of the mass balance
+        /// and afterwards the norm of the residual of the well flux and the well equation.
+        std::vector<double> computeResidualNorms() const;
 
         ADB
         fluidViscosity(const int               phase,
@@ -321,7 +368,6 @@ namespace Opm {
                    const ADB&              so,
                    const std::vector<int>& cells) const;
 
-
         ADB
         poroMult(const ADB& p) const;
 
@@ -350,7 +396,35 @@ namespace Opm {
 
         /// Compute convergence based on total mass balance (tol_mb) and maximum
         /// residual mass balance (tol_cnv).
-        bool getConvergence(const double dt, const int it);
+        bool getConvergence(const double dt, const int iteration);
+
+        /// \brief Compute the reduction within the convergence check.
+        /// \param[in] B     A matrix with MaxNumPhases columns and the same number rows
+        ///                  as the number of cells of the grid. B.col(i) contains the values
+        ///                  for phase i.
+        /// \param[in] tempV A matrix with MaxNumPhases columns and the same number rows
+        ///                  as the number of cells of the grid. tempV.col(i) contains the
+        ///                   values
+        ///                  for phase i.
+        /// \param[in] R     A matrix with MaxNumPhases columns and the same number rows
+        ///                  as the number of cells of the grid. B.col(i) contains the values
+        ///                  for phase i.
+        /// \param[out] R_sum An array of size MaxNumPhases where entry i contains the sum
+        ///                   of R for the phase i.
+        /// \param[out] maxCoeff An array of size MaxNumPhases where entry i contains the
+        ///                   maximum of tempV for the phase i.
+        /// \param[out] B_avg An array of size MaxNumPhases where entry i contains the average
+        ///                   of B for the phase i.
+        /// \param[in]  nc    The number of cells of the local grid.
+        /// \return The total pore volume over all cells.
+        double
+        convergenceReduction(const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases+1>& B,
+                             const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases+1>& tempV,
+                             const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases+1>& R,
+                             std::array<double,MaxNumPhases+1>& R_sum,
+                             std::array<double,MaxNumPhases+1>& maxCoeff,
+                             std::array<double,MaxNumPhases+1>& B_avg,
+                             int nc) const;
 
         void detectNewtonOscillations(const std::vector<std::vector<double>>& residual_history,
                                       const int it, const double relaxRelTol,
@@ -358,14 +432,15 @@ namespace Opm {
 
         void stablizeNewton(V& dx, V& dxOld, const double omega, const RelaxType relax_type) const;
 
-        double dpMaxRel() const { return dp_max_rel_; }
+        double dpMaxRel() const { return param_.dp_max_rel_; }
-        double dsMax() const { return ds_max_; }
+        double dsMax() const { return param_.ds_max_; }
-        double drsMaxRel() const { return drs_max_rel_; }
+        double drMaxRel() const { return param_.dr_max_rel_; }
-        enum RelaxType relaxType() const { return relax_type_; }
+        enum RelaxType relaxType() const { return param_.relax_type_; }
-        double relaxMax() const { return relax_max_; };
+        double relaxMax() const { return param_.relax_max_; };
-        double relaxIncrement() const { return relax_increment_; };
+        double relaxIncrement() const { return param_.relax_increment_; };
-        double relaxRelTol() const { return relax_rel_tol_; };
+        double relaxRelTol() const { return param_.relax_rel_tol_; };
-        double maxIter() const { return max_iter_; }
+        double maxIter() const     { return param_.max_iter_; }
+        double maxResidualAllowed() const { return param_.max_residual_allowed_; }
 
     };
 } // namespace Opm

opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp

@@ -1,5 +1,5 @@
 /*
-  Copyright 2014 SINTEF ICT, Applied Mathematics.
+  Copyright 2013 SINTEF ICT, Applied Mathematics.
   Copyright 2014 STATOIL ASA.
 
   This file is part of the Open Porous Media project (OPM).
@@ -26,6 +26,7 @@
 
 struct UnstructuredGrid;
 struct Wells;
+struct FlowBoundaryConditions;
 
 namespace Opm
 {
@@ -37,9 +38,8 @@ namespace Opm
     class SimulatorTimer;
     class PolymerBlackoilState;
     class WellStateFullyImplicitBlackoil;
-    class WellsManager;
     class EclipseState;
-    class EclipseWriter;
+    class BlackoilOutputWriter;
     class PolymerPropsAd;
     class PolymerInflowInterface;
     struct SimulatorReport;
@@ -70,7 +70,6 @@ namespace Opm
         /// \param[in] grid          grid data structure
         /// \param[in] geo           derived geological properties
         /// \param[in] props         fluid and rock properties
-        /// \param[in] polymer_props polymer properties
         /// \param[in] rock_comp_props     if non-null, rock compressibility properties
         /// \param[in] linsolver     linear solver
         /// \param[in] gravity       if non-null, gravity vector
@@ -91,7 +90,7 @@ namespace Opm
                                               const bool vapoil,
                                               const bool polymer,
                                               std::shared_ptr<EclipseState> eclipse_state,
-                                              EclipseWriter& output_writer,
+                                              BlackoilOutputWriter& output_writer,
                                               Opm::DeckConstPtr& deck,
                                               const std::vector<double>& threshold_pressures_by_face);
 

opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymerOutput.cpp

@@ -1,211 +0,0 @@
-#include "config.h"
-
-#include "SimulatorFullyImplicitBlackoilPolymerOutput.hpp"
-
-#include <opm/core/utility/DataMap.hpp>
-#include <opm/core/io/vtk/writeVtkData.hpp>
-#include <opm/core/utility/ErrorMacros.hpp>
-#include <opm/core/utility/miscUtilities.hpp>
-
-#include <opm/autodiff/GridHelpers.hpp>
-
-#include <sstream>
-#include <iomanip>
-#include <fstream>
-
-#include <boost/filesystem.hpp>
-
-#ifdef HAVE_DUNE_CORNERPOINT
-#include <opm/core/utility/platform_dependent/disable_warnings.h>
-#include <dune/common/version.hh>
-#include <dune/grid/io/file/vtk/vtkwriter.hh>
-#include <opm/core/utility/platform_dependent/reenable_warnings.h>
-#endif
-namespace Opm
-{
-
-
-    void outputStateVtk(const UnstructuredGrid& grid,
-                        const PolymerBlackoilState& state,
-                        const int step,
-                        const std::string& output_dir)
-    {
-        // Write data in VTK format.
-        std::ostringstream vtkfilename;
-        vtkfilename << output_dir << "/vtk_files";
-        boost::filesystem::path fpath(vtkfilename.str());
-        try {
-            create_directories(fpath);
-        }
-        catch (...) {
-            OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
-        }
-        vtkfilename << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu";
-        std::ofstream vtkfile(vtkfilename.str().c_str());
-        if (!vtkfile) {
-            OPM_THROW(std::runtime_error, "Failed to open " << vtkfilename.str());
-        }
-        Opm::DataMap dm;
-        dm["saturation"] = &state.saturation();
-        dm["pressure"] = &state.pressure();
-        dm["concentration"] = &state.concentration();
-        dm["maxconcentration"] = &state.maxconcentration();
-        std::vector<double> cell_velocity;
-        Opm::estimateCellVelocity(AutoDiffGrid::numCells(grid),
-                                  AutoDiffGrid::numFaces(grid),
-                                  AutoDiffGrid::beginFaceCentroids(grid),
-                                  AutoDiffGrid::faceCells(grid),
-                                  AutoDiffGrid::beginCellCentroids(grid),
-                                  AutoDiffGrid::beginCellVolumes(grid),
-                                  AutoDiffGrid::dimensions(grid),
-                                  state.faceflux(), cell_velocity);
-        dm["velocity"] = &cell_velocity;
-        Opm::writeVtkData(grid, dm, vtkfile);
-    }
-
-
-    void outputStateMatlab(const UnstructuredGrid& grid,
-                           const PolymerBlackoilState& state,
-                           const int step,
-                           const std::string& output_dir)
-    {
-        Opm::DataMap dm;
-        dm["saturation"] = &state.saturation();
-        dm["pressure"] = &state.pressure();
-        dm["surfvolume"] = &state.surfacevol();
-        dm["rs"] = &state.gasoilratio();
-        dm["rv"] = &state.rv();
-        dm["concentration"] = &state.concentration();
-        dm["maxconcentration"] = &state.maxconcentration();
-        std::vector<double> cell_velocity;
-        Opm::estimateCellVelocity(AutoDiffGrid::numCells(grid),
-                                  AutoDiffGrid::numFaces(grid),
-                                  AutoDiffGrid::beginFaceCentroids(grid),
-                                  UgGridHelpers::faceCells(grid),
-                                  AutoDiffGrid::beginCellCentroids(grid),
-                                  AutoDiffGrid::beginCellVolumes(grid),
-                                  AutoDiffGrid::dimensions(grid),
-                                  state.faceflux(), cell_velocity);
-        dm["velocity"] = &cell_velocity;
-
-        // 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"));
-        }
-    }
-    void outputWellStateMatlab(const Opm::WellState& 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
-    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);
-    }
-
-    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
-
-#ifdef HAVE_DUNE_CORNERPOINT
-    void outputStateVtk(const Dune::CpGrid& grid,
-                        const PolymerBlackoilState& state,
-                        const int step,
-                        const std::string& output_dir)
-    {
-        // Write data in VTK format.
-        std::ostringstream vtkfilename;
-        std::ostringstream vtkpath;
-        vtkpath << output_dir << "/vtk_files";
-        vtkpath << "/output-" << std::setw(3) << std::setfill('0') << step;
-        boost::filesystem::path fpath(vtkpath.str());
-        try {
-            create_directories(fpath);
-        }
-        catch (...) {
-            OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
-        }
-        vtkfilename << "output-" << std::setw(3) << std::setfill('0') << step;
-#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 3)
-        Dune::VTKWriter<Dune::CpGrid::LeafGridView> writer(grid.leafGridView(), Dune::VTK::nonconforming);
-#else
-        Dune::VTKWriter<Dune::CpGrid::LeafGridView> writer(grid.leafView(), Dune::VTK::nonconforming);
-#endif
-        writer.addCellData(state.saturation(), "saturation", state.numPhases());
-        writer.addCellData(state.pressure(), "pressure", 1);
-        writer.addCellData(state.concentration(), "concentration", 1);
-        writer.addCellData(state.maxconcentration(), "maxconcentration", 1);
-        
-        std::vector<double> cell_velocity;
-        Opm::estimateCellVelocity(AutoDiffGrid::numCells(grid),
-                                  AutoDiffGrid::numFaces(grid),
-                                  AutoDiffGrid::beginFaceCentroids(grid),
-                                  AutoDiffGrid::faceCells(grid),
-                                  AutoDiffGrid::beginCellCentroids(grid),
-                                  AutoDiffGrid::beginCellVolumes(grid),
-                                  AutoDiffGrid::dimensions(grid),
-                                  state.faceflux(), cell_velocity);
-        writer.addCellData(cell_velocity, "velocity", Dune::CpGrid::dimension);
-        writer.pwrite(vtkfilename.str(), vtkpath.str(), std::string("."), Dune::VTK::ascii);
-    }
-#endif
-
-}

opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymerOutput.hpp

@@ -1,96 +0,0 @@
-#ifndef OPM_SIMULATORFULLYIMPLICITBLACKOILPOLYMEROUTPUT_HEADER_INCLUDED
-#define OPM_SIMULATORFULLYIMPLICITBLACKOILPOLYMEROUTPUT_HEADER_INCLUDED
-#include <opm/core/grid.h>
-#include <opm/polymer/PolymerBlackoilState.hpp>
-#include <opm/core/simulator/WellState.hpp>
-#include <opm/core/utility/DataMap.hpp>
-#include <opm/core/utility/ErrorMacros.hpp>
-#include <opm/core/utility/miscUtilities.hpp>
-
-#include <opm/autodiff/GridHelpers.hpp>
-
-#include <string>
-#include <sstream>
-#include <iomanip>
-#include <fstream>
-
-#include <boost/filesystem.hpp>
-
-#ifdef HAVE_DUNE_CORNERPOINT
-#include <dune/grid/CpGrid.hpp>
-#endif
-namespace Opm
-{
-
-    void outputStateVtk(const UnstructuredGrid& grid,
-                        const PolymerBlackoilState& state,
-                        const int step,
-                        const std::string& output_dir);
-
-
-    void outputStateMatlab(const UnstructuredGrid& grid,
-                           const PolymerBlackoilState& state,
-                           const int step,
-                           const std::string& output_dir);
-
-    void outputWellStateMatlab(const Opm::WellState& well_state,
-                               const int step,
-                               const std::string& output_dir);
-#ifdef HAVE_DUNE_CORNERPOINT
-    void outputStateVtk(const Dune::CpGrid& grid,
-                        const PolymerBlackoilState& state,
-                        const int step,
-                        const std::string& output_dir);
-#endif
-
-    template<class Grid>
-    void outputStateMatlab(const Grid& grid,
-                           const PolymerBlackoilState& state,
-                           const int step,
-                           const std::string& output_dir)
-    {
-        Opm::DataMap dm;
-        dm["saturation"] = &state.saturation();
-        dm["pressure"] = &state.pressure();
-        dm["surfvolume"] = &state.surfacevol();
-        dm["rs"] = &state.gasoilratio();
-        dm["rv"] = &state.rv();
-        dm["concentration"] = &state.concentration();
-        dm["maxconcentration"] = &state.maxconcentration();
-
-        std::vector<double> cell_velocity;
-        Opm::estimateCellVelocity(AutoDiffGrid::numCells(grid),
-                                  AutoDiffGrid::numFaces(grid),
-                                  AutoDiffGrid::beginFaceCentroids(grid),
-                                  UgGridHelpers::faceCells(grid),
-                                  AutoDiffGrid::beginCellCentroids(grid),
-                                  AutoDiffGrid::beginCellVolumes(grid),
-                                  AutoDiffGrid::dimensions(grid),
-                                  state.faceflux(), cell_velocity);
-        dm["velocity"] = &cell_velocity;
-
-        // 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"));
-        }
-    }
-
-
-}
-#endif

opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp

@@ -1,5 +1,6 @@
 /*
-  Copyright 2014 SINTEF ICT, Applied Mathematics.
+  Copyright 2013 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 IRIS AS
   Copyright 2014 STATOIL ASA.
 
   This file is part of the Open Porous Media project (OPM).
@@ -18,11 +19,13 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-#include <opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymerOutput.hpp>
+//#include <opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilOutput.hpp>
+#include <opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp>
 #include <opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp>
 #include <opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver.hpp>
 #include <opm/polymer/PolymerBlackoilState.hpp>
 #include <opm/polymer/PolymerInflow.hpp>
+
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
 #include <opm/core/utility/ErrorMacros.hpp>
 
@@ -36,9 +39,9 @@
 #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/simulator/AdaptiveSimulatorTimer.hpp>
 #include <opm/core/utility/StopWatch.hpp>
 #include <opm/core/io/vtk/writeVtkData.hpp>
 #include <opm/core/utility/miscUtilities.hpp>
@@ -46,6 +49,7 @@
 
 #include <opm/core/props/rock/RockCompressibility.hpp>
 
+//#include <opm/core/simulator/AdaptiveTimeStepping.hpp>
 #include <opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp>
 
 #include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
@@ -80,7 +84,7 @@ namespace Opm
              const Grid& grid,
              const DerivedGeology& geo,
              BlackoilPropsAdInterface& props,
-             const PolymerPropsAd&  polymer_props,
+             const PolymerPropsAd& polymer_props,
              const RockCompressibility* rock_comp_props,
              NewtonIterationBlackoilInterface& linsolver,
              const double* gravity,
@@ -88,7 +92,7 @@ namespace Opm
              bool has_vapoil,
              bool has_polymer,
              std::shared_ptr<EclipseState> eclipse_state,
-             EclipseWriter& output_writer,
+             BlackoilOutputWriter& output_writer,
              Opm::DeckConstPtr& deck,
              const std::vector<double>& threshold_pressures_by_face);
 
@@ -103,11 +107,6 @@ namespace Opm
 
         const parameter::ParameterGroup param_;
 
-        // Parameters for output.
-        bool output_;
-        bool output_vtk_;
-        std::string output_dir_;
-        int output_interval_;
         // Observed objects.
         const Grid& grid_;
         BlackoilPropsAdInterface& props_;
@@ -122,10 +121,11 @@ namespace Opm
         const bool has_disgas_;
         const bool has_vapoil_;
         const bool has_polymer_;
+        bool       terminal_output_;
         // eclipse_state
         std::shared_ptr<EclipseState> eclipse_state_;
         // output_writer
-        EclipseWriter& output_writer_;
+        BlackoilOutputWriter& output_writer_;
         Opm::DeckConstPtr& deck_;
         RateConverterType rateConverter_;
         // Threshold pressures.
@@ -134,7 +134,7 @@ namespace Opm
         void
         computeRESV(const std::size_t               step,
                     const Wells*                    wells,
-                    const PolymerBlackoilState&     x,
+                    const BlackoilState&     x,
                     WellStateFullyImplicitBlackoil& xw);
     };
 
@@ -154,13 +154,13 @@ namespace Opm
                                                                                     const bool has_vapoil,
                                                                                     const bool has_polymer,
                                                                                     std::shared_ptr<EclipseState> eclipse_state,
-                                                                                    EclipseWriter& output_writer,
+                                                                                    BlackoilOutputWriter& output_writer,
                                                                                     Opm::DeckConstPtr& deck,
                                                                                     const std::vector<double>& threshold_pressures_by_face)
 
     {
-        pimpl_.reset(new Impl(param, grid, geo, props, polymer_props, rock_comp_props, linsolver, gravity, has_disgas, 
+        pimpl_.reset(new Impl(param, grid, geo, props, polymer_props, rock_comp_props, linsolver, gravity, has_disgas, has_vapoil, has_polymer,
-                              has_vapoil, has_polymer, eclipse_state, output_writer, deck, threshold_pressures_by_face));
+                              eclipse_state, output_writer, deck, threshold_pressures_by_face));
     }
 
 
@@ -175,64 +175,6 @@ namespace Opm
     }
 
 
-
-    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>
     SimulatorFullyImplicitBlackoilPolymer<T>::Impl::Impl(const parameter::ParameterGroup& param,
@@ -247,7 +189,7 @@ namespace Opm
                                                          const bool has_vapoil,
                                                          const bool has_polymer,
                                                          std::shared_ptr<EclipseState> eclipse_state,
-                                                         EclipseWriter& output_writer,
+                                                         BlackoilOutputWriter& output_writer,
                                                          Opm::DeckConstPtr& deck,
                                                          const std::vector<double>& threshold_pressures_by_face)
         : param_(param),
@@ -261,34 +203,30 @@ namespace Opm
           has_disgas_(has_disgas),
           has_vapoil_(has_vapoil),
           has_polymer_(has_polymer),
+          terminal_output_(param.getDefault("output_terminal", true)),
           eclipse_state_(eclipse_state),
           output_writer_(output_writer),
           deck_(deck),
           rateConverter_(props_, std::vector<int>(AutoDiffGrid::numCells(grid_), 0)),
           threshold_pressures_by_face_(threshold_pressures_by_face)
     {
-        // For output.
-        output_ = param.getDefault("output", true);
-        if (output_) {
-            output_vtk_ = param.getDefault("output_vtk", true);
-            output_dir_ = param.getDefault("output_dir", std::string("output"));
-            // Ensure that output dir exists
-            boost::filesystem::path fpath(output_dir_);
-            try {
-                create_directories(fpath);
-            }
-            catch (...) {
-                OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
-            }
-            output_interval_ = param.getDefault("output_interval", 1);
-        }
-
         // Misc init.
         const int num_cells = AutoDiffGrid::numCells(grid);
         allcells_.resize(num_cells);
         for (int cell = 0; cell < num_cells; ++cell) {
             allcells_[cell] = cell;
         }
+#if HAVE_MPI
+        if ( terminal_output_ ) {
+            if ( solver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
+            {
+                const ParallelISTLInformation& info =
+                    boost::any_cast<const ParallelISTLInformation&>(solver_.parallelInformation());
+                // Only rank 0 does print to std::cout
+                terminal_output_= (info.communicator().rank()==0);
+            }
+        }
+#endif
     }
 
 
@@ -306,14 +244,40 @@ namespace Opm
         Opm::time::StopWatch step_timer;
         Opm::time::StopWatch total_timer;
         total_timer.start();
-        std::string tstep_filename = output_dir_ + "/step_timing.txt";
+        std::string tstep_filename = output_writer_.outputDirectory() + "/step_timing.txt";
         std::ofstream tstep_os(tstep_filename.c_str());
 
+        typename FullyImplicitBlackoilPolymerSolver<T>::SolverParameter solverParam( param_ );
+
+        //adaptive time stepping
+        //        std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping;
+        //        if( param_.getDefault("timestep.adaptive", bool(false) ) )
+        //        {
+        //            adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_ ) );
+        //        }
+
+        // init output writer
+        output_writer_.writeInit( timer );
+
+        std::string restorefilename = param_.getDefault("restorefile", std::string("") );
+        if( ! restorefilename.empty() )
+        {
+            // -1 means that we'll take the last report step that was written
+            const int desiredRestoreStep = param_.getDefault("restorestep", int(-1) );
+            output_writer_.restore( timer, state.blackoilState(), prev_well_state, restorefilename, desiredRestoreStep );
+        }
+
+        unsigned int totalNewtonIterations = 0;
+        unsigned int totalLinearIterations = 0;
+
         // Main simulation loop.
         while (!timer.done()) {
             // Report timestep.
             step_timer.start();
-            timer.report(std::cout);
+            if ( terminal_output_ )
+            {
+                timer.report(std::cout);
+            }
 
             // Create wells and well state.
             WellsManager wells_manager(eclipse_state_,
@@ -328,6 +292,7 @@ namespace Opm
             const Wells* wells = wells_manager.c_wells();
             WellStateFullyImplicitBlackoil well_state;
             well_state.init(wells, state.blackoilState(), prev_well_state);
+
             // compute polymer inflow
             std::unique_ptr<PolymerInflowInterface> polymer_inflow_ptr;
             if (deck_->hasKeyword("WPOLYMER")) {
@@ -344,42 +309,54 @@ namespace Opm
             polymer_inflow_ptr->getInflowValues(timer.simulationTimeElapsed(), 
                                                 timer.simulationTimeElapsed() + timer.currentStepLength(),
                                                 polymer_inflow_c);
-            // Output state at start of time step.
+            
-            if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
+            // write simulation state at the report stage
-                if (output_vtk_) {
+            output_writer_.writeTimeStep( timer, state.blackoilState(), well_state );
-                    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.blackoilState(), well_state);
-            }
 
             // Max oil saturation (for VPPARS), hysteresis update.
             props_.updateSatOilMax(state.saturation());
             props_.updateSatHyst(state.saturation(), allcells_);
 
             // Compute reservoir volumes for RESV controls.
-            computeRESV(timer.currentStepNum(), wells, state, well_state);
+            computeRESV(timer.currentStepNum(), wells, state.blackoilState(), well_state);
 
-            // Run a single step of the solver.
+            // Run a multiple steps of the solver depending on the time step control.
             solver_timer.start();
-            FullyImplicitBlackoilPolymerSolver<T> solver(param_, grid_, props_, geo_, rock_comp_props_, polymer_props_, *wells, solver_, has_disgas_, has_vapoil_, has_polymer_);
+
+            FullyImplicitBlackoilPolymerSolver<T> solver(solverParam, grid_, props_, geo_, rock_comp_props_, polymer_props_, wells, solver_, has_disgas_, has_vapoil_, has_polymer_, terminal_output_);
             if (!threshold_pressures_by_face_.empty()) {
                 solver.setThresholdPressures(threshold_pressures_by_face_);
             }
+
+            // If sub stepping is enabled allow the solver to sub cycle
+            // in case the report steps are to large for the solver to converge
+            //
+            // \Note: The report steps are met in any case
+            // \Note: The sub stepping will require a copy of the state variables
+            //            if( adaptiveTimeStepping ) {
+            //                adaptiveTimeStepping->step( timer, solver, state, well_state,  output_writer_ );
+            //            } else {
+                // solve for complete report step
             solver.step(timer.currentStepLength(), state, well_state, polymer_inflow_c);
+                //            }
+
+            // take time that was used to solve system for this reportStep
             solver_timer.stop();
 
+            // accumulate the number of Newton and Linear Iterations
+            totalNewtonIterations += solver.newtonIterations();
+            totalLinearIterations += solver.linearIterations();
+
             // Report timing.
             const double st = solver_timer.secsSinceStart();
-            std::cout << "Fully implicit solver took: " << st << " seconds." << std::endl;
+
+            if ( terminal_output_ )
+            {
+                std::cout << "Fully implicit solver took: " << st << " seconds." << std::endl;
+            }
+
             stime += st;
-            if (output_) {
+            if ( output_writer_.output() ) {
                 SimulatorReport step_report;
                 step_report.pressure_time = st;
                 step_report.total_time =  step_timer.secsSinceStart();
@@ -392,14 +369,7 @@ namespace Opm
         }
 
         // Write final simulation state.
-        if (output_) {
+        output_writer_.writeTimeStep( timer, state.blackoilState(), prev_well_state );
-            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.blackoilState(), prev_well_state);
-        } 
 
         // Stop timer and create timing report
         total_timer.stop();
@@ -407,6 +377,8 @@ namespace Opm
         report.pressure_time = stime;
         report.transport_time = 0.0;
         report.total_time = total_timer.secsSinceStart();
+        report.total_newton_iterations = totalNewtonIterations;
+        report.total_linear_iterations = totalLinearIterations;
         return report;
     }
 
@@ -444,17 +416,16 @@ namespace Opm
         }
 
         inline bool
-        is_resv_prod(const Wells& wells,
+        is_resv(const Wells& wells,
-                     const int    w)
+                const int    w)
         {
-            return ((wells.type[w] == PRODUCER) &&
+            return (0 <= resv_control(wells.ctrls[w]));
-                    (0 <= resv_control(wells.ctrls[w])));
         }
 
         inline bool
-        is_resv_prod(const WellMap&     wmap,
+        is_resv(const WellMap&     wmap,
-                     const std::string& name,
+                const std::string& name,
-                     const std::size_t  step)
+                const std::size_t  step)
         {
             bool match = false;
 
@@ -465,29 +436,34 @@ namespace Opm
 
                 match = (wp->isProducer(step) &&
                          wp->getProductionProperties(step)
-                         .hasProductionControl(WellProducer::RESV));
+                         .hasProductionControl(WellProducer::RESV))
+                    ||  (wp->isInjector(step) &&
+                         wp->getInjectionProperties(step)
+                         .hasInjectionControl(WellInjector::RESV));
             }
 
             return match;
         }
 
         inline std::vector<int>
-        resvProducers(const Wells&      wells,
+        resvWells(const Wells*      wells,
-                      const std::size_t step,
+                  const std::size_t step,
-                      const WellMap&    wmap)
+                  const WellMap&    wmap)
         {
-            std::vector<int> resv_prod;
+            std::vector<int> resv_wells;
-
+            if( wells )
-            for (int w = 0, nw = wells.number_of_wells; w < nw; ++w) {
+            {
-                if (is_resv_prod(wells, w) ||
+                for (int w = 0, nw = wells->number_of_wells; w < nw; ++w) {
-                    ((wells.name[w] != 0) &&
+                    if (is_resv(*wells, w) ||
-                     is_resv_prod(wmap, wells.name[w], step)))
+                        ((wells->name[w] != 0) &&
-                {
+                         is_resv(wmap, wells->name[w], step)))
-                    resv_prod.push_back(w);
+                    {
+                        resv_wells.push_back(w);
+                    }
                 }
             }
 
-            return resv_prod;
+            return resv_wells;
         }
 
         inline void
@@ -527,7 +503,7 @@ namespace Opm
     SimulatorFullyImplicitBlackoilPolymer<T>::
     Impl::computeRESV(const std::size_t               step,
                       const Wells*                    wells,
-                      const PolymerBlackoilState&     x,
+                      const BlackoilState&            x,
                       WellStateFullyImplicitBlackoil& xw)
     {
         typedef SimFIBODetails::WellMap WellMap;
@@ -535,24 +511,24 @@ namespace Opm
         const std::vector<WellConstPtr>& w_ecl = eclipse_state_->getSchedule()->getWells(step);
         const WellMap& wmap = SimFIBODetails::mapWells(w_ecl);
 
-        const std::vector<int>& resv_prod =
+        const std::vector<int>& resv_wells = SimFIBODetails::resvWells(wells, step, wmap);
-            SimFIBODetails::resvProducers(*wells, step, wmap);
 
-        if (! resv_prod.empty()) {
+        if (! resv_wells.empty()) {
             const PhaseUsage&                    pu = props_.phaseUsage();
             const std::vector<double>::size_type np = props_.numPhases();
 
-            rateConverter_.defineState(x.blackoilState());
+            rateConverter_.defineState(x);
 
             std::vector<double> distr (np);
             std::vector<double> hrates(np);
             std::vector<double> prates(np);
 
             for (std::vector<int>::const_iterator
-                     rp = resv_prod.begin(), e = resv_prod.end();
+                     rp = resv_wells.begin(), e = resv_wells.end();
                  rp != e; ++rp)
             {
                 WellControls* ctrl = wells->ctrls[*rp];
+                const bool is_producer = wells->type[*rp] == PRODUCER;
 
                 // RESV control mode, all wells
                 {
@@ -561,11 +537,17 @@ namespace Opm
                     if (0 <= rctrl) {
                         const std::vector<double>::size_type off = (*rp) * np;
 
-                        // Convert to positive rates to avoid issues
+                        if (is_producer) {
-                        // in coefficient calculations.
+                            // Convert to positive rates to avoid issues
-                        std::transform(xw.wellRates().begin() + (off + 0*np),
+                            // in coefficient calculations.
-                                       xw.wellRates().begin() + (off + 1*np),
+                            std::transform(xw.wellRates().begin() + (off + 0*np),
-                                       prates.begin(), std::negate<double>());
+                                           xw.wellRates().begin() + (off + 1*np),
+                                           prates.begin(), std::negate<double>());
+                        } else {
+                            std::copy(xw.wellRates().begin() + (off + 0*np),
+                                      xw.wellRates().begin() + (off + 1*np),
+                                      prates.begin());
+                        }
 
                         const int fipreg = 0; // Hack.  Ignore FIP regions.
                         rateConverter_.calcCoeff(prates, fipreg, distr);
@@ -576,7 +558,7 @@ namespace Opm
 
                 // RESV control, WCONHIST wells.  A bit of duplicate
                 // work, regrettably.
-                if (wells->name[*rp] != 0) {
+                if (is_producer && wells->name[*rp] != 0) {
                     WellMap::const_iterator i = wmap.find(wells->name[*rp]);
 
                     if (i != wmap.end()) {
@@ -603,11 +585,18 @@ namespace Opm
                             well_controls_clear(ctrl);
                             well_controls_assert_number_of_phases(ctrl, int(np));
 
-                            const int ok =
+                            const int ok_resv =
                                 well_controls_add_new(RESERVOIR_RATE, target,
                                                       & distr[0], ctrl);
 
-                            if (ok != 0) {
+                            // For WCONHIST/RESV the BHP limit is set to 1 atm.
+                            // TODO: Make it possible to modify the BHP limit using
+                            // the WELTARG keyword
+                            const int ok_bhp =
+                                well_controls_add_new(BHP, unit::convert::from(1.0, unit::atm),
+                                                      NULL, ctrl);
+
+                            if (ok_resv != 0 && ok_bhp != 0) {
                                 xw.currentControls()[*rp] = 0;
                                 well_controls_set_current(ctrl, 0);
                             }