Merge pull request #1491 from andlaus/remove_more_state

Remove more state

CMakeLists_files.cmake

@@ -322,7 +322,6 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/autodiff/WellDensitySegmented.hpp
   opm/autodiff/WellStateFullyImplicitBlackoil.hpp
   opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp
-  opm/autodiff/BlackoilOutputEbos.hpp
   opm/autodiff/ThreadHandle.hpp
   opm/autodiff/VFPProperties.hpp
   opm/autodiff/VFPHelpers.hpp

opm/autodiff/BlackoilModelEbos.hpp

@@ -146,8 +146,7 @@ namespace Opm {
                           BlackoilWellModel<TypeTag>& well_model,
                           BlackoilAquiferModel<TypeTag>& aquifer_model,
                           const NewtonIterationBlackoilInterface& linsolver,
-                          const bool terminal_output
+                          const bool terminal_output)
-                          )
         : ebosSimulator_(ebosSimulator)
         , grid_(ebosSimulator_.vanguard().grid())
         , istlSolver_( dynamic_cast< const ISTLSolverType* > (&linsolver) )
@@ -180,11 +179,7 @@ namespace Opm {
 
         /// Called once before each time step.
         /// \param[in] timer                  simulation timer
-        /// \param[in, out] reservoir_state   reservoir state variables
+        void prepareStep(const SimulatorTimerInterface& timer)
-        /// \param[in, out] well_state        well state variables
-        void prepareStep(const SimulatorTimerInterface& timer,
-                         const ReservoirState& /*reservoir_state*/,
-                         const WellState& /* well_state */)
         {
 
             // update the solution variables in ebos
@@ -228,9 +223,7 @@ namespace Opm {
         template <class NonlinearSolverType>
         SimulatorReport nonlinearIteration(const int iteration,
                                            const SimulatorTimerInterface& timer,
-                                           NonlinearSolverType& nonlinear_solver,
+                                           NonlinearSolverType& nonlinear_solver)
-                                           ReservoirState& /*reservoir_state*/,
-                                           WellState& /*well_state*/)
         {
             SimulatorReport report;
             failureReport_ = SimulatorReport();
@@ -341,16 +334,8 @@ namespace Opm {
         /// Called once after each time step.
         /// In this class, this function does nothing.
         /// \param[in] timer                  simulation timer
-        /// \param[in, out] reservoir_state   reservoir state variables
+        void afterStep(const SimulatorTimerInterface& OPM_UNUSED timer)
-        /// \param[in, out] well_state        well state variables
-        void afterStep(const SimulatorTimerInterface& timer,
-                       const ReservoirState& reservoir_state,
-                       WellState& well_state)
         {
-            DUNE_UNUSED_PARAMETER(timer);
-            DUNE_UNUSED_PARAMETER(reservoir_state);
-            DUNE_UNUSED_PARAMETER(well_state);
-
             wellModel().timeStepSucceeded();
             aquiferModel().timeStepSucceeded(timer);
             ebosSimulator_.problem().endTimeStep();
@@ -411,8 +396,7 @@ namespace Opm {
         }
 
         // compute the "relative" change of the solution between time steps
-        template <class Dummy>
+        double relativeChange() const
-        double relativeChange(const Dummy&, const Dummy&) const
         {
             Scalar resultDelta = 0.0;
             Scalar resultDenom = 0.0;
@@ -484,7 +468,7 @@ namespace Opm {
             resultDenom = gridView.comm().sum(resultDenom);
 
             if (resultDenom > 0.0)
-              return resultDelta/resultDenom;
+                return resultDelta/resultDenom;
             return 0.0;
         }
 
@@ -1069,6 +1053,9 @@ namespace Opm {
         const Simulator& ebosSimulator() const
         { return ebosSimulator_; }
 
+        Simulator& ebosSimulator()
+        { return ebosSimulator_; }
+
         /// return the statistics if the nonlinearIteration() method failed
         const SimulatorReport& failureReport() const
         { return failureReport_; }
@@ -1145,4 +1132,4 @@ namespace Opm {
     };
 } // namespace Opm
 
-#endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED
+#endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED

opm/autodiff/BlackoilOutputEbos.hpp

@@ -1,198 +0,0 @@
-/*
-  Copyright (c) 2017 IRIS AS
-
-  This file is part of the Open Porous Media project (OPM).
-
-  OPM is free software: you can redistribute it and/or modify
-  it under the terms of the GNU General Public License as published by
-  the Free Software Foundation, either version 3 of the License, or
-  (at your option) any later version.
-
-  OPM is distributed in the hope that it will be useful,
-  but WITHOUT ANY WARRANTY; without even the implied warranty of
-  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-  GNU General Public License for more details.
-
-  You should have received a copy of the GNU General Public License
-  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
-*/
-#ifndef OPM_BLACKOILOUTPUTEBOS_HEADER_INCLUDED
-#define OPM_BLACKOILOUTPUTEBOS_HEADER_INCLUDED
-
-
-#include <ebos/eclproblem.hh>
-#include <ewoms/common/start.hh>
-
-#include <opm/grid/UnstructuredGrid.h>
-#include <opm/simulators/timestepping/SimulatorTimerInterface.hpp>
-#include <opm/core/utility/DataMap.hpp>
-#include <opm/common/ErrorMacros.hpp>
-#include <opm/common/OpmLog/OpmLog.hpp>
-#include <opm/core/utility/miscUtilities.hpp>
-#include <opm/common/utility/parameters/ParameterGroup.hpp>
-#include <opm/core/wells/DynamicListEconLimited.hpp>
-#include <opm/core/simulator/SimulatorReport.hpp>
-#include <opm/core/wells/WellsManager.hpp>
-#include <opm/output/data/Cells.hpp>
-#include <opm/output/data/Solution.hpp>
-
-#include <opm/autodiff/Compat.hpp>
-
-#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
-
-#include <opm/parser/eclipse/Units/UnitSystem.hpp>
-#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
-#include <opm/parser/eclipse/EclipseState/InitConfig/InitConfig.hpp>
-
-#include <string>
-#include <sstream>
-#include <iomanip>
-#include <fstream>
-#include <thread>
-#include <map>
-
-#include <boost/filesystem.hpp>
-
-#ifdef HAVE_OPM_GRID
-#include <opm/grid/CpGrid.hpp>
-#endif
-namespace Opm
-{
-
-
-    /// Extra data to read/write for OPM restarting
-    struct ExtraData
-    {
-        double suggested_step = -1.0;
-    };
-
-
-    /** \brief Wrapper ECL output. */
-    template<class TypeTag>
-    class BlackoilOutputEbos
-    {
-    public:
-
-        typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
-        typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
-        typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
-        typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
-        // constructor creating different sub writers
-        BlackoilOutputEbos(Simulator& ebosSimulator,
-                           const ParameterGroup& param)
-            : output_( [ &param ] () -> bool {
-                    // If output parameter is true or all, then we do output
-                    const std::string outputString = param.getDefault("output", std::string("all"));
-                    return ( outputString == "all" ||  outputString == "true" );
-                }()
-                ),
-            ebosSimulator_(ebosSimulator),
-            phaseUsage_(phaseUsageFromDeck(eclState()))
-        {}
-
-        /*!
-         * \brief Write a blackoil reservoir state to disk for later inspection with
-         *        visualization tools like ResInsight. This function will extract the
-         *        requested output cell properties specified by the RPTRST keyword
-         *        and write these to file.
-         */
-        template<class SimulationDataContainer, class Model>
-        void writeTimeStep(const SimulatorTimerInterface& timer,
-                           const SimulationDataContainer& /*reservoirStateDummy*/,
-                           const Opm::WellStateFullyImplicitBlackoil& /*wellStateDummy*/,
-                           const Model& physicalModel,
-                           const bool substep = false,
-                           const double nextstep = -1.0,
-                           const SimulatorReport& simulatorReport = SimulatorReport())
-        {
-            if( output_ )
-            {
-                // Add TCPU if simulatorReport is not defaulted.
-                const double totalSolverTime = simulatorReport.solver_time;
-
-                const Opm::WellStateFullyImplicitBlackoil& localWellState = physicalModel.wellModel().wellState();
-
-                // The writeOutput expects a local data::solution vector and a local data::well vector.
-                auto localWellData = localWellState.report(phaseUsage_, Opm::UgGridHelpers::globalCell(grid()) );
-                ebosSimulator_.problem().writeOutput(localWellData, timer.simulationTimeElapsed(), substep, totalSolverTime, nextstep);
-            }
-        }
-
-        template <class SimulationDataContainer, class WellState>
-        void initFromRestartFile(const PhaseUsage& /*phaseUsage*/,
-                                 const Grid& /*grid */,
-                                 SimulationDataContainer& simulatorstate,
-                                 WellState& wellstate,
-                                 ExtraData& extra)   {
-
-            std::vector<RestartKey> extra_keys = {
-              {"OPMEXTRA" , Opm::UnitSystem::measure::identity, false}
-            };
-
-            // gives a dummy dynamic_list_econ_limited
-            DynamicListEconLimited dummy_list_econ_limited;
-            const auto& defunct_well_names = ebosSimulator_.vanguard().defunctWellNames();
-            WellsManager wellsmanager(eclState(),
-                                      schedule(),
-                                      // The restart step value is used to identify wells present at the given
-                                      // time step. Wells that are added at the same time step as RESTART is initiated
-                                      // will not be present in a restart file. Use the previous time step to retrieve
-                                      // wells that have information written to the restart file.
-                                      std::max(eclState().getInitConfig().getRestartStep() - 1, 0),
-                                      Opm::UgGridHelpers::numCells(grid()),
-                                      Opm::UgGridHelpers::globalCell(grid()),
-                                      Opm::UgGridHelpers::cartDims(grid()),
-                                      Opm::UgGridHelpers::dimensions(grid()),
-                                      Opm::UgGridHelpers::cell2Faces(grid()),
-                                      Opm::UgGridHelpers::beginFaceCentroids(grid()),
-                                      dummy_list_econ_limited,
-                                      grid().comm().size() > 1,
-                                      defunct_well_names);
-
-            const Wells* wells = wellsmanager.c_wells();
-
-            std::vector<RestartKey> solution_keys = {};
-            auto restart_values = ebosSimulator_.problem().eclIO().loadRestart(solution_keys, extra_keys);
-
-            const int nw = wells->number_of_wells;
-            if (nw > 0) {
-                wellstate.resize(wells, simulatorstate, phaseUsage_ ); //Resize for restart step
-                wellsToState( restart_values.wells, phaseUsage_, wellstate );
-            }
-
-            if (restart_values.hasExtra("OPMEXTRA")) {
-                std::vector<double> opmextra = restart_values.getExtra("OPMEXTRA");
-                assert(opmextra.size() == 1);
-                extra.suggested_step = opmextra[0];
-            } else {
-                OpmLog::warning("Restart data is missing OPMEXTRA field, restart run may deviate from original run.");
-                extra.suggested_step = -1.0;
-            }
-        }
-
-        const Grid& grid()
-        { return ebosSimulator_.vanguard().grid(); }
-
-        const Schedule& schedule() const
-        { return ebosSimulator_.vanguard().schedule(); }
-
-        const EclipseState& eclState() const
-        { return ebosSimulator_.vanguard().eclState(); }
-
-        bool isRestart() const {
-            const auto& initconfig = eclState().getInitConfig();
-            return initconfig.restartRequested();
-        }   
-
-    protected:
-        const bool output_;
-        Simulator& ebosSimulator_;
-        Opm::PhaseUsage phaseUsage_;
-    };
-
-
-
-
-
-}
-#endif

opm/autodiff/BlackoilWellModel.hpp

@@ -49,6 +49,7 @@
 #include <opm/autodiff/WellInterface.hpp>
 #include <opm/autodiff/StandardWell.hpp>
 #include <opm/autodiff/MultisegmentWell.hpp>
+#include <opm/autodiff/Compat.hpp>
 #include<opm/autodiff/SimFIBODetails.hpp>
 #include<dune/common/fmatrix.hh>
 #include<dune/istl/bcrsmatrix.hh>
@@ -103,6 +104,40 @@ namespace Opm {
                               const ModelParameters& param,
                               const bool terminal_output);
 
+            void initFromRestartFile(const RestartValue& restartValues)
+            {
+                // gives a dummy dynamic_list_econ_limited
+                DynamicListEconLimited dummyListEconLimited;
+                const auto& defunctWellNames = ebosSimulator_.vanguard().defunctWellNames();
+                WellsManager wellsmanager(eclState(),
+                                          schedule(),
+                                          // The restart step value is used to identify wells present at the given
+                                          // time step. Wells that are added at the same time step as RESTART is initiated
+                                          // will not be present in a restart file. Use the previous time step to retrieve
+                                          // wells that have information written to the restart file.
+                                          std::max(eclState().getInitConfig().getRestartStep() - 1, 0),
+                                          Opm::UgGridHelpers::numCells(grid()),
+                                          Opm::UgGridHelpers::globalCell(grid()),
+                                          Opm::UgGridHelpers::cartDims(grid()),
+                                          Opm::UgGridHelpers::dimensions(grid()),
+                                          Opm::UgGridHelpers::cell2Faces(grid()),
+                                          Opm::UgGridHelpers::beginFaceCentroids(grid()),
+                                          dummyListEconLimited,
+                                          grid().comm().size() > 1,
+                                          defunctWellNames);
+
+                const Wells* wells = wellsmanager.c_wells();
+
+                const int nw = wells->number_of_wells;
+                if (nw > 0) {
+                    auto phaseUsage = phaseUsageFromDeck(eclState());
+                    size_t numCells = Opm::UgGridHelpers::numCells(grid());
+                    well_state_.resize(wells, numCells, phaseUsage); //Resize for restart step
+                    wellsToState(restartValues.wells, phaseUsage, well_state_);
+                    previous_well_state_ = well_state_;
+                }
+            }
+
             // compute the well fluxes and assemble them in to the reservoir equations as source terms
             // and in the well equations.
             void assemble(const int iterationIdx,
@@ -161,6 +196,32 @@ namespace Opm {
             }
 
         protected:
+            void extractLegacyPressure_(std::vector<double>& cellPressure) const
+            {
+                size_t nc = number_of_cells_;
+                std::vector<double> cellPressures(nc, 0.0);
+                ElementContext elemCtx(ebosSimulator_);
+                const auto& gridView = ebosSimulator_.vanguard().gridView();
+                const auto& elemEndIt = gridView.template end</*codim=*/0>();
+                for (auto elemIt = gridView.template begin</*codim=*/0>();
+                     elemIt != elemEndIt;
+                     ++elemIt)
+                {
+                    const auto& elem = *elemIt;
+                    if (elem.partitionType() != Dune::InteriorEntity) {
+                        continue;
+                    }
+                    elemCtx.updatePrimaryStencil(elem);
+                    elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
+
+                    const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
+                    const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
+                    const auto& fs = intQuants.fluidState();
+
+                    const double p = fs.pressure(FluidSystem::oilPhaseIdx).value();
+                    cellPressures[cellIdx] = p;
+                }
+            }
 
             Simulator& ebosSimulator_;
             std::unique_ptr<WellsManager> wells_manager_;
@@ -204,6 +265,12 @@ namespace Opm {
 
             const Wells* wells() const { return wells_manager_->c_wells(); }
 
+            const Grid& grid() const
+            { return ebosSimulator_.vanguard().grid(); }
+
+            const EclipseState& eclState() const
+            { return ebosSimulator_.vanguard().eclState(); }
+
             const Schedule& schedule() const
             { return ebosSimulator_.vanguard().schedule(); }
 

opm/autodiff/BlackoilWellModel_impl.hpp

@@ -88,7 +88,7 @@ namespace Opm {
             const double p = fs.pressure(FluidSystem::oilPhaseIdx).value();
             cellPressures[cellIdx] = p;
         }
-        well_state_.init(wells(), cellPressures, previous_well_state_, phase_usage_);
+        well_state_.init(wells(), cellPressures, &previous_well_state_, phase_usage_);
 
         // handling MS well related
         if (param_.use_multisegment_well_) { // if we use MultisegmentWell model

opm/autodiff/FlowMainEbos.hpp

@@ -85,7 +85,6 @@ namespace Opm
 
         typedef Opm::SimulatorFullyImplicitBlackoilEbos<TypeTag> Simulator;
         typedef typename Simulator::ReservoirState ReservoirState;
-        typedef typename Simulator::OutputWriter OutputWriter;
 
         /// This is the main function of Flow.
         /// It runs a complete simulation, with the given grid and
@@ -107,7 +106,6 @@ namespace Opm
                 setupLogging();
                 printPRTHeader();
                 runDiagnostics();
-                setupOutputWriter();
                 setupLinearSolver();
                 createSimulator();
 
@@ -480,20 +478,6 @@ namespace Opm
             diagnostic.diagnosis(eclState(), deck(), this->grid());
         }
 
-        // Setup output writer.
-        // Writes to:
-        //   output_writer_
-        void setupOutputWriter()
-        {
-            // create output writer after grid is distributed, otherwise the parallel output
-            // won't work correctly since we need to create a mapping from the distributed to
-            // the global view
-
-            output_writer_.reset(new OutputWriter(*ebosSimulator_,
-                                                   param_));
-
-        }
-
         // Run the simulator.
         // Returns EXIT_SUCCESS if it does not throw.
         int runSimulator()
@@ -567,10 +551,7 @@ namespace Opm
             // Create the simulator instance.
             simulator_.reset(new Simulator(*ebosSimulator_,
                                            param_,
-                                           *fis_solver_,
+                                           *fis_solver_));
-                                           FluidSystem::enableDissolvedGas(),
-                                           FluidSystem::enableVaporizedOil(),
-                                           *output_writer_));
         }
 
     private:
@@ -619,7 +600,6 @@ namespace Opm
         bool must_distribute_ = false;
         ParameterGroup param_;
         bool output_to_files_ = false;
-        std::unique_ptr<OutputWriter> output_writer_;
         boost::any parallel_information_;
         std::unique_ptr<NewtonIterationBlackoilInterface> fis_solver_;
         std::unique_ptr<Simulator> simulator_;

opm/autodiff/NonlinearSolverEbos.hpp

@@ -0,0 +1,350 @@
+/*
+  Copyright 2015 SINTEF ICT, Applied Mathematics.
+  Copyright 2015 Statoil ASA.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_NON_LINEAR_SOLVER_EBOS_HPP
+#define OPM_NON_LINEAR_SOLVER_EBOS_HPP
+
+#include <opm/core/simulator/SimulatorReport.hpp>
+#include <opm/common/utility/parameters/ParameterGroup.hpp>
+#include <opm/common/ErrorMacros.hpp>
+#include <opm/common/Exceptions.hpp>
+#include <opm/simulators/timestepping/SimulatorTimerInterface.hpp>
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <memory>
+
+namespace Opm {
+
+
+    /// A nonlinear solver class suitable for general fully-implicit models,
+    /// as well as pressure, transport and sequential models.
+    template <class PhysicalModel>
+    class NonlinearSolverEbos
+    {
+    public:
+        // Available relaxation scheme types.
+        enum RelaxType {
+            Dampen,
+            SOR
+        };
+
+        // Solver parameters controlling nonlinear process.
+        struct SolverParametersEbos
+        {
+            RelaxType relaxType_;
+            double relaxMax_;
+            double relaxIncrement_;
+            double relaxRelTol_;
+            int maxIter_; // max nonlinear iterations
+            int minIter_; // min nonlinear iterations
+
+            explicit SolverParametersEbos(const ParameterGroup& param)
+            {
+                // set default values
+                reset();
+
+                // overload with given parameters
+                relaxMax_ = param.getDefault("relax_max", relaxMax_);
+                maxIter_ = param.getDefault("max_iter", maxIter_);
+                minIter_ = param.getDefault("min_iter", minIter_);
+
+                std::string relaxationType = param.getDefault("relax_type", std::string("dampen"));
+                if (relaxationType == "dampen") {
+                    relaxType_ = Dampen;
+                } else if (relaxationType == "sor") {
+                    relaxType_ = SOR;
+                } else {
+                    OPM_THROW(std::runtime_error, "Unknown Relaxtion Type " << relaxationType);
+                }
+            }
+
+            SolverParametersEbos()
+            { reset(); }
+
+            void reset()
+            {
+                // default values for the solver parameters
+                relaxType_ = Dampen;
+                relaxMax_ = 0.5;
+                relaxIncrement_ = 0.1;
+                relaxRelTol_ = 0.2;
+                maxIter_ = 10;
+                minIter_ = 1;
+            }
+
+        };
+
+        // Forwarding types from PhysicalModel.
+        typedef typename PhysicalModel::ReservoirState ReservoirState;
+        typedef typename PhysicalModel::WellState WellState;
+
+        // ---------  Public methods  ---------
+
+        /// Construct solver for a given model.
+        ///
+        /// The model is a std::unique_ptr because the object to which model points to is
+        /// not allowed to be deleted as long as the NonlinearSolver object exists.
+        ///
+        /// \param[in]      param   parameters controlling nonlinear process
+        /// \param[in, out] model   physical simulation model.
+        NonlinearSolverEbos(const SolverParametersEbos& param,
+                            std::unique_ptr<PhysicalModel> model)
+            : param_(param)
+            , model_(std::move(model))
+            , linearizations_(0)
+            , nonlinearIterations_(0)
+            , linearIterations_(0)
+            , wellIterations_(0)
+            , nonlinearIterationsLast_(0)
+            , linearIterationsLast_(0)
+            , wellIterationsLast_(0)
+        {
+            if (!model_) {
+                OPM_THROW(std::logic_error, "Must provide a non-null model argument for NonlinearSolver.");
+            }
+        }
+
+
+        SimulatorReport step(const SimulatorTimerInterface& timer)
+        {
+            SimulatorReport iterReport;
+            SimulatorReport report;
+            failureReport_ = SimulatorReport();
+
+            // Do model-specific once-per-step calculations.
+            model_->prepareStep(timer);
+
+            int iteration = 0;
+
+            // Let the model do one nonlinear iteration.
+
+            // Set up for main solver loop.
+            bool converged = false;
+
+            // ----------  Main nonlinear solver loop  ----------
+            do {
+                try {
+                    // Do the nonlinear step. If we are in a converged state, the
+                    // model will usually do an early return without an expensive
+                    // solve, unless the minIter() count has not been reached yet.
+                    iterReport = model_->nonlinearIteration(iteration, timer, *this);
+
+                    report += iterReport;
+                    report.converged = iterReport.converged;
+
+                    converged = report.converged;
+                    iteration += 1;
+                }
+                catch (...) {
+                    // if an iteration fails during a time step, all previous iterations
+                    // count as a failure as well
+                    failureReport_ += report;
+                    failureReport_ += model_->failureReport();
+                    throw;
+                }
+            }
+            while ( (!converged && (iteration <= maxIter())) || (iteration <= minIter()));
+
+            if (!converged) {
+                failureReport_ += report;
+
+                std::string msg = "Solver convergence failure - Failed to complete a time step within " + std::to_string(maxIter()) + " iterations.";
+                OPM_THROW_NOLOG(Opm::TooManyIterations, msg);
+            }
+
+            // Do model-specific post-step actions.
+            model_->afterStep(timer);
+            report.converged = true;
+
+            return report;
+        }
+
+        /// return the statistics if the step() method failed
+        const SimulatorReport& failureReport() const
+        { return failureReport_; }
+
+        /// Number of linearizations used in all calls to step().
+        int linearizations() const
+        { return linearizations_; }
+
+        /// Number of full nonlinear solver iterations used in all calls to step().
+        int nonlinearIterations() const
+        { return nonlinearIterations_; }
+
+        /// Number of linear solver iterations used in all calls to step().
+        int linearIterations() const
+        { return linearIterations_; }
+
+        /// Number of well iterations used in all calls to step().
+        int wellIterations() const
+        { return wellIterations_; }
+
+        /// Number of nonlinear solver iterations used in the last call to step().
+        int nonlinearIterationsLastStep() const
+        { return nonlinearIterationsLast_; }
+
+        /// Number of linear solver iterations used in the last call to step().
+        int linearIterationsLastStep() const
+        { return linearIterationsLast_; }
+
+        /// Number of well iterations used in all calls to step().
+        int wellIterationsLastStep() const
+        { return wellIterationsLast_; }
+
+        /// Compute fluid in place.
+        /// \param[in]    ReservoirState
+        /// \param[in]    FIPNUM for active cells not global cells.
+        /// \return fluid in place, number of fip regions, each region contains 5 values which are liquid, vapour, water, free gas and dissolved gas.
+        std::vector<std::vector<double> >
+        computeFluidInPlace(const ReservoirState& x, const std::vector<int>& fipnum) const
+        { return model_->computeFluidInPlace(x, fipnum); }
+
+        std::vector<std::vector<double> >
+        computeFluidInPlace(const std::vector<int>& fipnum) const
+        { return model_->computeFluidInPlace(fipnum); }
+
+        /// Reference to physical model.
+        const PhysicalModel& model() const
+        { return *model_; }
+
+        /// Mutable reference to physical model.
+        PhysicalModel& model()
+        { return *model_; }
+
+        /// Detect oscillation or stagnation in a given residual history.
+        void detectOscillations(const std::vector<std::vector<double>>& residualHistory,
+                                const int it, bool& oscillate, bool& stagnate) const
+        {
+            // The detection of oscillation in two primary variable results in the report of the detection
+            // of oscillation for the solver.
+            // Only the saturations are used for oscillation detection for the black oil model.
+            // Stagnate is not used for any treatment here.
+
+            if ( it < 2 ) {
+                oscillate = false;
+                stagnate = false;
+                return;
+            }
+
+            stagnate = true;
+            int oscillatePhase = 0;
+            const std::vector<double>& F0 = residualHistory[it];
+            const std::vector<double>& F1 = residualHistory[it - 1];
+            const std::vector<double>& F2 = residualHistory[it - 2];
+            for (int p= 0; p < model_->numPhases(); ++p){
+                const double d1 = std::abs((F0[p] - F2[p]) / F0[p]);
+                const double d2 = std::abs((F0[p] - F1[p]) / F0[p]);
+
+                oscillatePhase += (d1 < relaxRelTol()) && (relaxRelTol() < d2);
+
+                // Process is 'stagnate' unless at least one phase
+                // exhibits significant residual change.
+                stagnate = (stagnate && !(std::abs((F1[p] - F2[p]) / F2[p]) > 1.0e-3));
+            }
+
+            oscillate = (oscillatePhase > 1);
+        }
+
+
+        /// Apply a stabilization to dx, depending on dxOld and relaxation parameters.
+        /// Implemention for Dune block vectors.
+        template <class BVector>
+        void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld, const double omega) const
+        {
+            // The dxOld is updated with dx.
+            // If omega is equal to 1., no relaxtion will be appiled.
+
+            BVector tempDxOld = dxOld;
+            dxOld = dx;
+
+            switch (relaxType()) {
+            case Dampen: {
+                if (omega == 1.) {
+                    return;
+                }
+                auto i = dx.size();
+                for (i = 0; i < dx.size(); ++i) {
+                    dx[i] *= omega;
+                }
+                return;
+            }
+            case SOR: {
+                if (omega == 1.) {
+                    return;
+                }
+                auto i = dx.size();
+                for (i = 0; i < dx.size(); ++i) {
+                    dx[i] *= omega;
+                    tempDxOld[i] *= (1.-omega);
+                    dx[i] += tempDxOld[i];
+                }
+                return;
+            }
+            default:
+                OPM_THROW(std::runtime_error, "Can only handle Dampen and SOR relaxation type.");
+            }
+
+            return;
+        }
+
+        /// The greatest relaxation factor (i.e. smallest factor) allowed.
+        double relaxMax() const
+        { return param_.relaxMax_; }
+
+        /// The step-change size for the relaxation factor.
+        double relaxIncrement() const
+        { return param_.relaxIncrement_; }
+
+        /// The relaxation type (Dampen or SOR).
+        enum RelaxType relaxType() const
+        { return param_.relaxType_; }
+
+        /// The relaxation relative tolerance.
+        double relaxRelTol() const
+        { return param_.relaxRelTol_; }
+
+        /// The maximum number of nonlinear iterations allowed.
+        int maxIter() const
+        { return param_.maxIter_; }
+
+        /// The minimum number of nonlinear iterations allowed.
+        int minIter() const
+        { return param_.minIter_; }
+
+        /// Set parameters to override those given at construction time.
+        void setParameters(const SolverParametersEbos& param)
+        { param_ = param; }
+
+    private:
+        // ---------  Data members  ---------
+        SimulatorReport failureReport_;
+        SolverParametersEbos param_;
+        std::unique_ptr<PhysicalModel> model_;
+        int linearizations_;
+        int nonlinearIterations_;
+        int linearIterations_;
+        int wellIterations_;
+        int nonlinearIterationsLast_;
+        int linearIterationsLast_;
+        int wellIterationsLast_;
+    };
+} // namespace Opm
+
+#endif // OPM_NON_LINEAR_SOLVER_EBOS_HPP

opm/autodiff/ParallelDebugOutput.hpp

@@ -596,8 +596,9 @@ namespace Opm
         };
 
         // gather solution to rank 0 for EclipseWriter
+        template <class WellState>
         bool collectToIORank( const SimulationDataContainer& /*localReservoirState*/,
-                              const WellStateFullyImplicitBlackoil& localWellState,
+                              const WellState& localWellState,
                               const data::Solution& localCellData,
                               const int wellStateStepNumber )
         {
@@ -627,7 +628,7 @@ namespace Opm
                                            std::unordered_set<std::string>());
 
                 const Wells* wells = wells_manager.c_wells();
-                globalWellState_.init(wells, *globalReservoirState_, globalWellState_, phaseUsage_ );
+                globalWellState_.initLegacy(wells, *globalReservoirState_, globalWellState_, phaseUsage_ );
                 globalCellData_->clear();
             }
 

opm/autodiff/SimulatorBase_impl.hpp

@@ -174,7 +174,7 @@ namespace Opm
                                        defunct_well_names_);
             const Wells* wells = wells_manager.c_wells();
             WellState well_state;
-            well_state.init(wells, state, prev_well_state, props_.phaseUsage());
+            well_state.initLegacy(wells, state, prev_well_state, props_.phaseUsage());
 
             // give the polymer and surfactant simulators the chance to do their stuff
             asImpl().handleAdditionalWellInflow(timer, wells_manager, well_state, wells);

opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp

@@ -22,16 +22,15 @@
 #ifndef OPM_SIMULATORFULLYIMPLICITBLACKOILEBOS_HEADER_INCLUDED
 #define OPM_SIMULATORFULLYIMPLICITBLACKOILEBOS_HEADER_INCLUDED
 
-#include <opm/autodiff/BlackoilOutputEbos.hpp>
 #include <opm/autodiff/IterationReport.hpp>
-#include <opm/autodiff/NonlinearSolver.hpp>
+#include <opm/autodiff/NonlinearSolverEbos.hpp>
 #include <opm/autodiff/BlackoilModelEbos.hpp>
 #include <opm/autodiff/BlackoilModelParameters.hpp>
 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 #include <opm/autodiff/BlackoilWellModel.hpp>
 #include <opm/autodiff/BlackoilAquiferModel.hpp>
 #include <opm/autodiff/moduleVersion.hpp>
-#include <opm/simulators/timestepping/AdaptiveTimeStepping.hpp>
+#include <opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp>
 #include <opm/grid/utility/StopWatch.hpp>
 
 #include <opm/common/Exceptions.hpp>
@@ -61,10 +60,9 @@ public:
 
     typedef WellStateFullyImplicitBlackoil WellState;
     typedef BlackoilState ReservoirState;
-    typedef BlackoilOutputEbos<TypeTag> OutputWriter;
     typedef BlackoilModelEbos<TypeTag> Model;
     typedef BlackoilModelParameters ModelParameters;
-    typedef NonlinearSolver<Model> Solver;
+    typedef NonlinearSolverEbos<Model> Solver;
     typedef BlackoilWellModel<TypeTag> WellModel;
     typedef BlackoilAquiferModel<TypeTag> AquiferModel;
 
@@ -94,30 +92,21 @@ public:
     /// \param[in] threshold_pressures_by_face   if nonempty, threshold pressures that inhibit flow
     SimulatorFullyImplicitBlackoilEbos(Simulator& ebosSimulator,
                                        const ParameterGroup& param,
-                                       NewtonIterationBlackoilInterface& linsolver,
+                                       NewtonIterationBlackoilInterface& linsolver)
-                                       const bool has_disgas,
+        : ebosSimulator_(ebosSimulator)
-                                       const bool has_vapoil,
+        , param_(param)
-                                       OutputWriter& output_writer)
+        , modelParam_(param)
-        : ebosSimulator_(ebosSimulator),
+        , solverParam_(param)
-          param_(param),
+        , solver_(linsolver)
-          model_param_(param),
+        , phaseUsage_(phaseUsageFromDeck(eclState()))
-          solver_param_(param),
+        , terminalOutput_(param.getDefault("output_terminal", true))
-          solver_(linsolver),
-          phaseUsage_(phaseUsageFromDeck(eclState())),
-          has_disgas_(has_disgas),
-          has_vapoil_(has_vapoil),
-          terminal_output_(param.getDefault("output_terminal", true)),
-          output_writer_(output_writer),
-          is_parallel_run_( false )
     {
 #if HAVE_MPI
-        if ( solver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
+        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_ = terminal_output_ && ( info.communicator().rank() == 0 );
+            terminalOutput_ = terminalOutput_ && (info.communicator().rank() == 0);
-            is_parallel_run_ = ( info.communicator().size() > 1 );
         }
 #endif
     }
@@ -130,43 +119,51 @@ public:
     /// \return                    simulation report, with timing data
     SimulatorReport run(SimulatorTimer& timer)
     {
-
-        ReservoirState dummy_state(0,0,0);
-
-        WellState prev_well_state;
-
-        ExtraData extra;
-
         failureReport_ = SimulatorReport();
 
-        if (output_writer_.isRestart()) {
+        // handle restarts
-            // This is a restart, populate WellState
+        std::unique_ptr<RestartValue> restartValues;
-            ReservoirState stateInit(Opm::UgGridHelpers::numCells(grid()),
+        if (isRestart()) {
-                                     Opm::UgGridHelpers::numFaces(grid()),
+            std::vector<RestartKey> extraKeys = {
-                                     phaseUsage_.num_phases);
+                {"OPMEXTRA" , Opm::UnitSystem::measure::identity, false}
-            output_writer_.initFromRestartFile(phaseUsage_, grid(), stateInit, prev_well_state, extra);
+            };
+
+            std::vector<RestartKey> solutionKeys = {};
+            restartValues.reset(new RestartValue(ebosSimulator_.problem().eclIO().loadRestart(solutionKeys, extraKeys)));
         }
 
         // Create timers and file for writing timing info.
-        Opm::time::StopWatch solver_timer;
+        Opm::time::StopWatch solverTimer;
-        Opm::time::StopWatch total_timer;
+        Opm::time::StopWatch totalTimer;
-        total_timer.start();
+        totalTimer.start();
 
         // adaptive time stepping
         const auto& events = schedule().getEvents();
-        std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping;
+        std::unique_ptr< AdaptiveTimeSteppingEbos > adaptiveTimeStepping;
         const bool useTUNING = param_.getDefault("use_TUNING", false);
-        if( param_.getDefault("timestep.adaptive", true ) )
+        if (param_.getDefault("timestep.adaptive", true)) {
-        {
             if (useTUNING) {
-                adaptiveTimeStepping.reset( new AdaptiveTimeStepping( schedule().getTuning(), timer.currentStepNum(), param_, terminal_output_ ) );
+                adaptiveTimeStepping.reset(new AdaptiveTimeSteppingEbos(schedule().getTuning(), timer.currentStepNum(), param_, terminalOutput_));
-            } else {
+            }
-                adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_, terminal_output_ ) );
+            else {
+                adaptiveTimeStepping.reset(new AdaptiveTimeSteppingEbos(param_, terminalOutput_));
             }
 
-            if (output_writer_.isRestart()) {
+            double suggestedStepSize = -1.0;
-                if (extra.suggested_step > 0.0) {
+            if (isRestart()) {
-                    adaptiveTimeStepping->setSuggestedNextStep(extra.suggested_step);
+                // This is a restart, determine the time step size from the restart data
+                if (restartValues->hasExtra("OPMEXTRA")) {
+                    std::vector<double> opmextra = restartValues->getExtra("OPMEXTRA");
+                    assert(opmextra.size() == 1);
+                    suggestedStepSize = opmextra[0];
+                }
+                else {
+                    OpmLog::warning("Restart data is missing OPMEXTRA field, restart run may deviate from original run.");
+                    suggestedStepSize = -1.0;
+                }
+
+                if (suggestedStepSize > 0.0) {
+                    adaptiveTimeStepping->setSuggestedNextStep(suggestedStepSize);
                 }
             }
         }
@@ -174,15 +171,13 @@ public:
         SimulatorReport report;
         SimulatorReport stepReport;
 
-        WellModel well_model(ebosSimulator_, model_param_, terminal_output_);
+        WellModel wellModel(ebosSimulator_, modelParam_, terminalOutput_);
-        if (output_writer_.isRestart()) {
+        if (isRestart()) {
-            well_model.setRestartWellState(prev_well_state); // Neccessary for perfect restarts
+            wellModel.initFromRestartFile(*restartValues);
         }
 
-        WellState wellStateDummy; //not used. Only passed to make the old interfaces happy
+        if (modelParam_.matrix_add_well_contributions_ ||
-
+             modelParam_.preconditioner_add_well_contributions_)
-        if ( model_param_.matrix_add_well_contributions_ ||
-             model_param_.preconditioner_add_well_contributions_ )
         {
             ebosSimulator_.model().clearAuxiliaryModules();
             auto auxMod = std::make_shared<WellConnectionAuxiliaryModule<TypeTag> >(schedule(), grid());
@@ -194,19 +189,18 @@ public:
         // Main simulation loop.
         while (!timer.done()) {
             // Report timestep.
-            if ( terminal_output_ )
+            if (terminalOutput_) {
-            {
                 std::ostringstream ss;
                 timer.report(ss);
                 OpmLog::debug(ss.str());
             }
 
             // Run a multiple steps of the solver depending on the time step control.
-            solver_timer.start();
+            solverTimer.start();
 
-            well_model.beginReportStep(timer.currentStepNum());
+            wellModel.beginReportStep(timer.currentStepNum());
 
-            auto solver = createSolver(well_model, aquifer_model);
+            auto solver = createSolver(wellModel, aquifer_model);
 
             // write the inital state at the report stage
             if (timer.initialStep()) {
@@ -215,22 +209,26 @@ public:
 
                 // No per cell data is written for initial step, but will be
                 // for subsequent steps, when we have started simulating
-                output_writer_.writeTimeStep( timer, dummy_state, well_model.wellState(), solver->model() );
+                auto localWellData = wellModel.wellState().report(phaseUsage_, Opm::UgGridHelpers::globalCell(grid()));
+                ebosSimulator_.problem().writeOutput(localWellData,
+                                                     timer.simulationTimeElapsed(),
+                                                     /*isSubstep=*/false,
+                                                     totalTimer.secsSinceStart(),
+                                                     /*nextStepSize=*/-1.0);
 
                 report.output_write_time += perfTimer.stop();
             }
 
-            if( terminal_output_ )
+            if (terminalOutput_) {
-            {
+                std::ostringstream stepMsg;
-                std::ostringstream step_msg;
                 boost::posix_time::time_facet* facet = new boost::posix_time::time_facet("%d-%b-%Y");
-                step_msg.imbue(std::locale(std::locale::classic(), facet));
+                stepMsg.imbue(std::locale(std::locale::classic(), facet));
-                step_msg << "\nReport step " << std::setw(2) <<timer.currentStepNum()
+                stepMsg << "\nReport step " << std::setw(2) <<timer.currentStepNum()
                          << "/" << timer.numSteps()
                          << " at day " << (double)unit::convert::to(timer.simulationTimeElapsed(), unit::day)
                          << "/" << (double)unit::convert::to(timer.totalTime(), unit::day)
                          << ", date = " << timer.currentDateTime();
-                OpmLog::info(step_msg.str());
+                OpmLog::info(stepMsg.str());
             }
 
             solver->model().beginReportStep();
@@ -240,9 +238,9 @@ public:
             //
             // \Note: The report steps are met in any case
             // \Note: The sub stepping will require a copy of the state variables
-            if( adaptiveTimeStepping ) {
+            if (adaptiveTimeStepping) {
                 if (useTUNING) {
-                    if(events.hasEvent(ScheduleEvents::TUNING_CHANGE,timer.currentStepNum())) {
+                    if (events.hasEvent(ScheduleEvents::TUNING_CHANGE,timer.currentStepNum())) {
                         adaptiveTimeStepping->updateTUNING(schedule().getTuning(), timer.currentStepNum());
                     }
                 }
@@ -251,18 +249,17 @@ public:
                         events.hasEvent(ScheduleEvents::PRODUCTION_UPDATE, timer.currentStepNum()) ||
                         events.hasEvent(ScheduleEvents::INJECTION_UPDATE, timer.currentStepNum()) ||
                         events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE, timer.currentStepNum());
-                stepReport = adaptiveTimeStepping->step( timer, *solver, dummy_state, wellStateDummy, event, output_writer_, nullptr );
+                stepReport = adaptiveTimeStepping->step(timer, *solver, event, nullptr);
                 report += stepReport;
                 failureReport_ += adaptiveTimeStepping->failureReport();
             }
             else {
                 // solve for complete report step
-                stepReport = solver->step(timer, dummy_state, wellStateDummy);
+                stepReport = solver->step(timer);
                 report += stepReport;
                 failureReport_ += solver->failureReport();
 
-                if( terminal_output_ )
+                if (terminalOutput_) {
-                {
                     std::ostringstream ss;
                     stepReport.reportStep(ss);
                     OpmLog::info(ss.str());
@@ -270,20 +267,19 @@ public:
             }
 
             solver->model().endReportStep();
-            well_model.endReportStep();
+            wellModel.endReportStep();
 
             // take time that was used to solve system for this reportStep
-            solver_timer.stop();
+            solverTimer.stop();
 
             // update timing.
-            report.solver_time += solver_timer.secsSinceStart();
+            report.solver_time += solverTimer.secsSinceStart();
 
             // Increment timer, remember well state.
             ++timer;
 
 
-            if (terminal_output_ )
+            if (terminalOutput_) {
-            {
                 if (!timer.initialStep()) {
                     const std::string version = moduleVersionName();
                     outputTimestampFIP(timer, version);
@@ -295,13 +291,17 @@ public:
             perfTimer.start();
             const double nextstep = adaptiveTimeStepping ? adaptiveTimeStepping->suggestedNextStep() : -1.0;
 
-            output_writer_.writeTimeStep( timer, dummy_state, well_model.wellState(), solver->model(), false, nextstep, report);
+            auto localWellData = wellModel.wellState().report(phaseUsage_, Opm::UgGridHelpers::globalCell(grid()));
+            ebosSimulator_.problem().writeOutput(localWellData,
+                                                 timer.simulationTimeElapsed(),
+                                                 /*isSubstep=*/false,
+                                                 totalTimer.secsSinceStart(),
+                                                 nextstep);
             report.output_write_time += perfTimer.stop();
 
-            if (terminal_output_ )
+            if (terminalOutput_) {
-            {
                 std::string msg =
-                    "Time step took " + std::to_string(solver_timer.secsSinceStart()) + " seconds; "
+                    "Time step took " + std::to_string(solverTimer.secsSinceStart()) + " seconds; "
                     "total solver time " + std::to_string(report.solver_time) + " seconds.";
                 OpmLog::debug(msg);
             }
@@ -309,8 +309,8 @@ public:
         }
 
         // Stop timer and create timing report
-        total_timer.stop();
+        totalTimer.stop();
-        report.total_time = total_timer.secsSinceStart();
+        report.total_time = totalTimer.secsSinceStart();
         report.converged = true;
 
         return report;
@@ -318,23 +318,24 @@ public:
 
     /** \brief Returns the simulator report for the failed substeps of the simulation.
      */
-    const SimulatorReport& failureReport() const { return failureReport_; };
+    const SimulatorReport& failureReport() const
+    { return failureReport_; };
 
     const Grid& grid() const
     { return ebosSimulator_.vanguard().grid(); }
 
 protected:
 
-    std::unique_ptr<Solver> createSolver(WellModel& well_model, AquiferModel& aquifer_model)
+    std::unique_ptr<Solver> createSolver(WellModel& wellModel, AquiferModel& aquifer_model)
     {
         auto model = std::unique_ptr<Model>(new Model(ebosSimulator_,
-                                                      model_param_,
+                                                      modelParam_,
-                                                      well_model,
+                                                      wellModel,
                                                       aquifer_model,
                                                       solver_,
-                                                      terminal_output_));
+                                                      terminalOutput_));
 
-        return std::unique_ptr<Solver>(new Solver(solver_param_, std::move(model)));
+        return std::unique_ptr<Solver>(new Solver(solverParam_, std::move(model)));
     }
 
     void outputTimestampFIP(const SimulatorTimer& timer, const std::string version)
@@ -358,31 +359,28 @@ protected:
     const Schedule& schedule() const
     { return ebosSimulator_.vanguard().schedule(); }
 
+    bool isRestart() const
+    {
+        const auto& initconfig = eclState().getInitConfig();
+        return initconfig.restartRequested();
+    }
 
     // Data.
     Simulator& ebosSimulator_;
 
-    typedef typename Solver::SolverParameters SolverParameters;
+    typedef typename Solver::SolverParametersEbos SolverParametersEbos;
 
     SimulatorReport failureReport_;
 
     const ParameterGroup param_;
-    ModelParameters model_param_;
+    ModelParameters modelParam_;
-    SolverParameters solver_param_;
+    SolverParametersEbos solverParam_;
 
     // Observed objects.
     NewtonIterationBlackoilInterface& solver_;
     PhaseUsage phaseUsage_;
     // Misc. data
-    const bool has_disgas_;
+    bool       terminalOutput_;
-    const bool has_vapoil_;
-    bool       terminal_output_;
-    // output_writer
-    OutputWriter& output_writer_;
-
-    // Whether this a parallel simulation or not
-    bool is_parallel_run_;
-
 };
 
 } // namespace Opm

opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp

@@ -459,7 +459,8 @@ namespace Opm
                                   std::unordered_set<std::string>());
 
         const Wells* wells = wellsmanager.c_wells();
-        wellstate.resize(wells, simulatorstate, phaseUsage ); //Resize for restart step
+        size_t numCells = Opm::UgGridHelpers::numCells(grid);
+        wellstate.resize(wells, numCells, phaseUsage ); //Resize for restart step
         auto restart_values = eclIO_->loadRestart(solution_keys, extra_keys);
 
         solutionToSim( restart_values, phaseUsage, simulatorstate );

opm/autodiff/StandardWell_impl.hpp

@@ -658,9 +658,8 @@ namespace Opm
             {
                 const double target_rate = well_controls_get_current_target(well_controls_); // surface rate target
                 if (well_type_ == INJECTOR) {
-                    const WellInjectionProperties& injection = well_ecl_->getInjectionProperties(current_step_);
                     // only handles single phase injection now
-                    assert(injection.injectorType != WellInjector::MULTI);
+                    assert(well_ecl_->getInjectionProperties(current_step_).injectorType != WellInjector::MULTI);
                     control_eq = getWQTotal() - target_rate;
                 } else if (well_type_ == PRODUCER) {
                     if (target_rate != 0.) {
@@ -704,9 +703,8 @@ namespace Opm
             {
                 const double target_rate = well_controls_get_current_target(well_controls_); // reservoir rate target
                 if (well_type_ == INJECTOR) {
-                    const WellInjectionProperties& injection = well_ecl_->getInjectionProperties(current_step_);
                     // only handles single phase injection now
-                    assert(injection.injectorType != WellInjector::MULTI);
+                    assert(well_ecl_->getInjectionProperties(current_step_).injectorType != WellInjector::MULTI);
                     const double* distr = well_controls_get_current_distr(well_controls_);
                     for (int phase = 0; phase < number_of_phases_; ++phase) {
                         if (distr[phase] > 0.0) {

opm/autodiff/WellStateFullyImplicitBlackoil.hpp

@@ -59,17 +59,10 @@ namespace Opm
         using BaseType :: numWells;
         using BaseType :: numPhases;
 
-        template <class State, class PrevWellState>
-        void init(const Wells* wells, const State& state, const PrevWellState& prevState, const PhaseUsage& pu)
-        {
-            init(wells, state.pressure(), prevState, pu);
-        }
-
         /// Allocate and initialize if wells is non-null.  Also tries
         /// to give useful initial values to the bhp(), wellRates()
         /// and perfPhaseRates() fields, depending on controls
-        template <class PrevWellState>
+        void init(const Wells* wells, const std::vector<double>& cellPressures, const WellStateFullyImplicitBlackoil* prevState, const PhaseUsage& pu)
-        void init(const Wells* wells, const std::vector<double>& cellPressures , const PrevWellState& prevState, const PhaseUsage& pu)
         {
             // call init on base class
             BaseType :: init(wells, cellPressures);
@@ -92,10 +85,10 @@ namespace Opm
             well_vaporized_oil_rates_.resize(nw, 0.0);
 
             is_new_well_.resize(nw, true);
-            if ( !prevState.wellMap().empty() ) {
+            if (prevState && !prevState->wellMap().empty()) {
-                const auto& end = prevState.wellMap().end();
+                const auto& end = prevState->wellMap().end();
                 for (int w = 0; w < nw; ++w) {
-                    const auto& it = prevState.wellMap().find( wells->name[w]);
+                    const auto& it = prevState->wellMap().find( wells->name[w]);
                     if (it != end) {
                         is_new_well_[w] = false;
                     }
@@ -135,6 +128,175 @@ namespace Opm
             perfRateSolvent_.clear();
             perfRateSolvent_.resize(nperf, 0.0);
 
+            // intialize wells that have been there before
+            // order may change so the mapping is based on the well name
+            if(prevState && !prevState->wellMap().empty()) {
+                typedef typename WellMapType :: const_iterator const_iterator;
+                const_iterator end = prevState->wellMap().end();
+                for (int w = 0; w < nw; ++w) {
+                    std::string name( wells->name[ w ] );
+                    const_iterator it = prevState->wellMap().find( name );
+                    if( it != end )
+                    {
+                        const int oldIndex = (*it).second[ 0 ];
+                        const int newIndex = w;
+
+                        // bhp
+                        bhp()[ newIndex ] = prevState->bhp()[ oldIndex ];
+
+                        // thp
+                        thp()[ newIndex ] = prevState->thp()[ oldIndex ];
+
+                        // wellrates
+                        for( int i=0, idx=newIndex*np, oldidx=oldIndex*np; i<np; ++i, ++idx, ++oldidx )
+                        {
+                            wellRates()[ idx ] = prevState->wellRates()[ oldidx ];
+                        }
+
+                        // perfPhaseRates
+                        const int oldPerf_idx_beg = (*it).second[ 1 ];
+                        const int num_perf_old_well = (*it).second[ 2 ];
+                        const int num_perf_this_well = wells->well_connpos[newIndex + 1] - wells->well_connpos[newIndex];
+                        // copy perforation rates when the number of perforations is equal,
+                        // otherwise initialize perfphaserates to well rates divided by the number of perforations.
+                        if( num_perf_old_well == num_perf_this_well )
+                        {
+                            int old_perf_phase_idx = oldPerf_idx_beg *np;
+                            for (int perf_phase_idx = wells->well_connpos[ newIndex ]*np;
+                                 perf_phase_idx < wells->well_connpos[ newIndex + 1]*np; ++perf_phase_idx, ++old_perf_phase_idx )
+                            {
+                                perfPhaseRates()[ perf_phase_idx ] = prevState->perfPhaseRates()[ old_perf_phase_idx ];
+                            }
+                        } else {
+                            for (int perf = wells->well_connpos[newIndex]; perf < wells->well_connpos[newIndex + 1]; ++perf) {
+                                for (int p = 0; p < np; ++p) {
+                                    perfPhaseRates()[np*perf + p] = wellRates()[np*newIndex + p] / double(num_perf_this_well);
+                                }
+                            }
+                        }
+                        // perfPressures
+                        if( num_perf_old_well == num_perf_this_well )
+                        {
+                            int oldPerf_idx = oldPerf_idx_beg;
+                            for (int perf = wells->well_connpos[ newIndex ];
+                                 perf < wells->well_connpos[ newIndex + 1]; ++perf, ++oldPerf_idx )
+                            {
+                                perfPress()[ perf ] = prevState->perfPress()[ oldPerf_idx ];
+                            }
+                        }
+                        // perfSolventRates
+                        if (pu.has_solvent) {
+                            if( num_perf_old_well == num_perf_this_well )
+                            {
+                                int oldPerf_idx = oldPerf_idx_beg;
+                                for (int perf = wells->well_connpos[ newIndex ];
+                                     perf < wells->well_connpos[ newIndex + 1]; ++perf, ++oldPerf_idx )
+                                {
+                                    perfRateSolvent()[ perf ] = prevState->perfRateSolvent()[ oldPerf_idx ];
+                                }
+                            }
+                        }
+                    }
+
+
+                    // If in the new step, there is no THP related target/limit anymore, its thp value should be
+                    // set to zero.
+                    const WellControls* ctrl = wells->ctrls[w];
+                    const int nwc = well_controls_get_num(ctrl);
+                    int ctrl_index = 0;
+                    for (; ctrl_index < nwc; ++ctrl_index) {
+                        if (well_controls_iget_type(ctrl, ctrl_index) == THP) {
+                            break;
+                        }
+                    }
+                    // not finding any thp related control/limits
+                    if (ctrl_index == nwc) {
+                        thp()[w] = 0.;
+                    }
+                }
+            }
+
+            {
+                // we need to create a trival segment related values to avoid there will be some
+                // multi-segment wells added later.
+                top_segment_index_.reserve(nw);
+                for (int w = 0; w < nw; ++w) {
+                    top_segment_index_.push_back(w);
+                }
+                segpress_ = bhp();
+                segrates_ = wellRates();
+            }
+        }
+
+        void resize(const Wells* wells, size_t numCells, const PhaseUsage& pu)
+        {
+            std::vector<double> tmp(numCells, 0.0); // <- UGLY HACK to pass the size
+            init(wells, tmp, nullptr, pu);
+        }
+
+        /// Allocate and initialize if wells is non-null.  Also tries
+        /// to give useful initial values to the bhp(), wellRates()
+        /// and perfPhaseRates() fields, depending on controls.
+        ///
+        /// this method is only for flow_legacy!
+        template <class PrevWellState>
+        void initLegacy(const Wells* wells, const std::vector<double>& cellPressures , const PrevWellState& prevState, const PhaseUsage& pu)
+        {
+            // call init on base class
+            BaseType :: init(wells, cellPressures);
+
+            // if there are no well, do nothing in init
+            if (wells == 0) {
+                return;
+            }
+
+            const int nw = wells->number_of_wells;
+
+            if( nw == 0 ) return ;
+
+            // Initialize perfphaserates_, which must be done here.
+            const int np = wells->number_of_phases;
+            const int nperf = wells->well_connpos[nw];
+
+            well_reservoir_rates_.resize(nw * np, 0.0);
+            well_dissolved_gas_rates_.resize(nw, 0.0);
+            well_vaporized_oil_rates_.resize(nw, 0.0);
+
+            // Ensure that we start out with zero rates by default.
+            perfphaserates_.clear();
+            perfphaserates_.resize(nperf * np, 0.0);
+            for (int w = 0; w < nw; ++w) {
+                assert((wells->type[w] == INJECTOR) || (wells->type[w] == PRODUCER));
+                const WellControls* ctrl = wells->ctrls[w];
+
+                if (well_controls_well_is_stopped(ctrl)) {
+                    // Shut well: perfphaserates_ are all zero.
+                } else {
+                    const int num_perf_this_well = wells->well_connpos[w + 1] - wells->well_connpos[w];
+                    // Open well: Initialize perfphaserates_ to well
+                    // rates divided by the number of perforations.
+                    for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) {
+                        for (int p = 0; p < np; ++p) {
+                            perfphaserates_[np*perf + p] = wellRates()[np*w + p] / double(num_perf_this_well);
+                        }
+                        perfPress()[perf] = cellPressures[wells->well_cells[perf]];
+                    }
+                }
+            }
+
+            // Initialize current_controls_.
+            // The controls set in the Wells object are treated as defaults,
+            // and also used for initial values.
+            current_controls_.resize(nw);
+            for (int w = 0; w < nw; ++w) {
+                current_controls_[w] = well_controls_get_current(wells->ctrls[w]);
+            }
+
+            is_new_well_.resize(nw, true);
+
+            perfRateSolvent_.clear();
+            perfRateSolvent_.resize(nperf, 0.0);
+
             // intialize wells that have been there before
             // order may change so the mapping is based on the well name
             if( ! prevState.wellMap().empty() )
@@ -146,6 +308,9 @@ namespace Opm
                     const_iterator it = prevState.wellMap().find( name );
                     if( it != end )
                     {
+                        // this is not a new added well
+                        is_new_well_[w] = false;
+
                         const int oldIndex = (*it).second[ 0 ];
                         const int newIndex = w;
 
@@ -236,10 +401,19 @@ namespace Opm
             }
         }
 
+        // this method is only for flow_legacy!
+        template <class State, class PrevWellState>
+        void initLegacy(const Wells* wells, const State& state, const PrevWellState& prevState, const PhaseUsage& pu)
+        {
+            initLegacy(wells, state.pressure(), prevState, pu);
+        }
+
+        // this method is only for flow_legacy!
         template <class State>
-        void resize(const Wells* wells, const State& state, const PhaseUsage& pu) {
+        void resizeLegacy(const Wells* wells, const State& state, const PhaseUsage& pu)
+        {
             const WellStateFullyImplicitBlackoil dummy_state{}; // Init with an empty previous state only resizes
-            init(wells, state, dummy_state, pu) ;
+            initLegacy(wells, state, dummy_state, pu) ;
         }
 
         /// One rate per phase and well connection.
@@ -250,7 +424,8 @@ namespace Opm
         std::vector<int>& currentControls() { return current_controls_; }
         const std::vector<int>& currentControls() const { return current_controls_; }
 
-        data::Wells report(const PhaseUsage &pu, const int* globalCellIdxMap) const override {
+        data::Wells report(const PhaseUsage &pu, const int* globalCellIdxMap) const override
+        {
             data::Wells res = WellState::report(pu, globalCellIdxMap);
 
             const int nw = this->numWells();

opm/simulators/WellSwitchingLogger.hpp

@@ -69,7 +69,7 @@ public:
         if( cc_.size() > 1 )
         {
             using Pair = typename SwitchMap::value_type;
-            switchMap_.insert(Pair(name, {char(from), char(to)}));
+            switchMap_.insert(Pair(name, {{char(from), char(to)}}));
         }
         else
         {

opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp

@@ -0,0 +1,383 @@
+/*
+*/
+#ifndef OPM_ADAPTIVE_TIME_STEPPING_EBOS_HPP
+#define OPM_ADAPTIVE_TIME_STEPPING_EBOS_HPP
+
+#include <iostream>
+#include <utility>
+
+#include <opm/grid/utility/StopWatch.hpp>
+#include <opm/common/OpmLog/OpmLog.hpp>
+#include <opm/common/utility/parameters/ParameterGroup.hpp>
+#include <opm/common/ErrorMacros.hpp>
+#include <opm/simulators/timestepping/SimulatorTimer.hpp>
+#include <opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp>
+#include <opm/simulators/timestepping/TimeStepControlInterface.hpp>
+#include <opm/simulators/timestepping/TimeStepControl.hpp>
+
+namespace Opm {
+    // AdaptiveTimeStepping
+    //---------------------
+
+    class AdaptiveTimeSteppingEbos
+    {
+        template <class Solver>
+        class SolutionTimeErrorSolverWrapperEbos : public RelativeChangeInterface
+        {
+            const Solver& solver_;
+        public:
+            SolutionTimeErrorSolverWrapperEbos(const Solver& solver)
+              : solver_(solver)
+            {}
+
+            /// return || u^n+1 - u^n || / || u^n+1 ||
+            double relativeChange() const
+            { return solver_.model().relativeChange(); }
+        };
+
+        template<class E>
+        void logException_(const E& exception, bool verbose)
+        {
+            if (verbose) {
+                std::string message;
+                message = "Caught Exception: ";
+                message += exception.what();
+                OpmLog::debug(message);
+            }
+        }
+
+    public:
+        //! \brief contructor taking parameter object
+        AdaptiveTimeSteppingEbos(const ParameterGroup& param,
+                                 const bool terminalOutput = true)
+            : timeStepControl_()
+            , restartFactor_( param.getDefault("solver.restartfactor", double(0.33)))
+            , growthFactor_( param.getDefault("solver.growthfactor", double(2)))
+            , maxGrowth_( param.getDefault("timestep.control.maxgrowth", double(3.0)))
+              // default is 1 year, convert to seconds
+            , maxTimeStep_( unit::convert::from(param.getDefault("timestep.max_timestep_in_days", 365.0), unit::day))
+            , solverRestartMax_( param.getDefault("solver.restart", int(10)))
+            , solverVerbose_( param.getDefault("solver.verbose", bool(true)) && terminalOutput)
+            , timestepVerbose_( param.getDefault("timestep.verbose", bool(true)) && terminalOutput)
+            , suggestedNextTimestep_( unit::convert::from(param.getDefault("timestep.initial_timestep_in_days", 1.0), unit::day))
+            , fullTimestepInitially_( param.getDefault("full_timestep_initially", bool(false)))
+            , timestepAfterEvent_( unit::convert::from(param.getDefault("timestep.timestep_in_days_after_event", -1.0), unit::day))
+            , useNewtonIteration_(false)
+        {
+            init_(param);
+        }
+
+
+
+        //! \brief contructor taking parameter object
+        //! \param tuning Pointer to ecl TUNING keyword
+        //! \param timeStep current report step
+        //! \param param The parameter object
+        AdaptiveTimeSteppingEbos(const Tuning& tuning,
+                                 size_t timeStep,
+                                 const ParameterGroup& param,
+                                 const bool terminalOutput = true)
+            : timeStepControl_()
+            , restartFactor_(tuning.getTSFCNV(timeStep))
+            , growthFactor_(tuning.getTFDIFF(timeStep))
+            , maxGrowth_(tuning.getTSFMAX(timeStep))
+              // default is 1 year, convert to seconds
+            , maxTimeStep_(tuning.getTSMAXZ(timeStep))
+            , solverRestartMax_(param.getDefault("solver.restart", int(10)))
+            , solverVerbose_(param.getDefault("solver.verbose", bool(true)) && terminalOutput)
+            , timestepVerbose_(param.getDefault("timestep.verbose", bool(true)) && terminalOutput)
+            , suggestedNextTimestep_(tuning.getTSINIT(timeStep))
+            , fullTimestepInitially_(param.getDefault("full_timestep_initially", bool(false)))
+            , timestepAfterEvent_(tuning.getTMAXWC(timeStep))
+            , useNewtonIteration_(false)
+        {
+            init_(param);
+        }
+
+        /** \brief  step method that acts like the solver::step method
+                    in a sub cycle of time steps
+        */
+        template <class Solver>
+        SimulatorReport step(const SimulatorTimer& simulatorTimer,
+                             Solver& solver,
+                             const bool isEvent,
+                             const std::vector<int>* fipnum = nullptr)
+        {
+            SimulatorReport report;
+            const double timestep = simulatorTimer.currentStepLength();
+
+            // init last time step as a fraction of the given time step
+            if (suggestedNextTimestep_ < 0) {
+                suggestedNextTimestep_ = restartFactor_ * timestep;
+            }
+
+            if (fullTimestepInitially_) {
+                suggestedNextTimestep_ = timestep;
+            }
+
+            // use seperate time step after event
+            if (isEvent && timestepAfterEvent_ > 0) {
+                suggestedNextTimestep_ = timestepAfterEvent_;
+            }
+
+            auto& ebosSimulator = solver.model().ebosSimulator();
+            auto& ebosProblem = ebosSimulator.problem();
+            auto phaseUsage = Opm::phaseUsageFromDeck(ebosSimulator.vanguard().eclState());
+
+            // create adaptive step timer with previously used sub step size
+            AdaptiveSimulatorTimer substepTimer(simulatorTimer, suggestedNextTimestep_, maxTimeStep_);
+
+            // reset the statistics for the failed substeps
+            failureReport_ = SimulatorReport();
+
+            // counter for solver restarts
+            int restarts = 0;
+
+            // sub step time loop
+            while (!substepTimer.done()) {
+                // get current delta t
+                const double dt = substepTimer.currentStepLength() ;
+                if (timestepVerbose_) {
+                    std::ostringstream ss;
+                    ss <<"\nTime step " << substepTimer.currentStepNum() << ", stepsize "
+                       << unit::convert::to(substepTimer.currentStepLength(), unit::day) << " days.";
+                    OpmLog::info(ss.str());
+                }
+
+                SimulatorReport substepReport;
+                std::string causeOfFailure = "";
+                try {
+                    substepReport = solver.step(substepTimer);
+                    report += substepReport;
+
+                    if (solverVerbose_) {
+                        // report number of linear iterations
+                        OpmLog::debug("Overall linear iterations used: " + std::to_string(substepReport.total_linear_iterations));
+                    }
+                }
+                catch (const Opm::TooManyIterations& e) {
+                    substepReport += solver.failureReport();
+                    causeOfFailure = "Solver convergence failure - Iteration limit reached";
+
+                    logException_(e, solverVerbose_);
+                    // since linearIterations is < 0 this will restart the solver
+                }
+                catch (const Opm::LinearSolverProblem& e) {
+                    substepReport += solver.failureReport();
+                    causeOfFailure = "Linear solver convergence failure";
+
+                    logException_(e, solverVerbose_);
+                    // since linearIterations is < 0 this will restart the solver
+                }
+                catch (const Opm::NumericalIssue& e) {
+                    substepReport += solver.failureReport();
+                    causeOfFailure = "Solver convergence failure - Numerical problem encountered";
+
+                    logException_(e, solverVerbose_);
+                    // since linearIterations is < 0 this will restart the solver
+                }
+                catch (const std::runtime_error& e) {
+                    substepReport += solver.failureReport();
+
+                    logException_(e, solverVerbose_);
+                    // also catch linear solver not converged
+                }
+                catch (const Dune::ISTLError& e) {
+                    substepReport += solver.failureReport();
+
+                    logException_(e, solverVerbose_);
+                    // also catch errors in ISTL AMG that occur when time step is too large
+                }
+                catch (const Dune::MatrixBlockError& e) {
+                    substepReport += solver.failureReport();
+
+                    logException_(e, solverVerbose_);
+                    // this can be thrown by ISTL's ILU0 in block mode, yet is not an ISTLError
+                }
+
+                if (substepReport.converged) {
+                    // advance by current dt
+                    ++substepTimer;
+
+                    // create object to compute the time error, simply forwards the call to the model
+                    SolutionTimeErrorSolverWrapperEbos<Solver> relativeChange(solver);
+
+                    // compute new time step estimate
+                    const int iterations = useNewtonIteration_ ? substepReport.total_newton_iterations
+                        : substepReport.total_linear_iterations;
+                    double dtEstimate = timeStepControl_->computeTimeStepSize(dt, iterations, relativeChange,
+                                                                               substepTimer.simulationTimeElapsed());
+
+                    // limit the growth of the timestep size by the growth factor
+                    dtEstimate = std::min(dtEstimate, double(maxGrowth_ * dt));
+
+                    // further restrict time step size growth after convergence problems
+                    if (restarts > 0) {
+                        dtEstimate = std::min(growthFactor_ * dt, dtEstimate);
+                        // solver converged, reset restarts counter
+                        restarts = 0;
+                    }
+
+                    if (timestepVerbose_) {
+                        std::ostringstream ss;
+                        substepReport.reportStep(ss);
+                        OpmLog::info(ss.str());
+                    }
+
+                    // write data if outputWriter was provided
+                    // if the time step is done we do not need
+                    // to write it as this will be done by the simulator
+                    // anyway.
+                    if (!substepTimer.done()) {
+                        if (fipnum) {
+                            solver.computeFluidInPlace(*fipnum);
+                        }
+                        Opm::time::StopWatch perfTimer;
+                        perfTimer.start();
+
+                        // The writeOutput expects a local data::solution vector and a local data::well vector.
+                        auto localWellData = solver.model().wellModel().wellState().report(phaseUsage, Opm::UgGridHelpers::globalCell(ebosSimulator.vanguard().grid()));
+                        ebosProblem.writeOutput(localWellData,
+                                                substepTimer.simulationTimeElapsed(),
+                                                /*isSubstep=*/true,
+                                                substepReport.total_time,
+                                                /*nextStepSize=*/-1.0);
+
+                        report.output_write_time += perfTimer.secsSinceStart();
+                    }
+
+                    // set new time step length
+                    substepTimer.provideTimeStepEstimate(dtEstimate);
+
+                    report.converged = substepTimer.done();
+                    substepTimer.setLastStepFailed(false);
+
+                }
+                else { // in case of no convergence (linearIterations < 0)
+                    substepTimer.setLastStepFailed(true);
+
+                    failureReport_ += substepReport;
+
+                    // increase restart counter
+                    if (restarts >= solverRestartMax_) {
+                        const auto msg = std::string("Solver failed to converge after cutting timestep ")
+                            + std::to_string(restarts) + " times.";
+                        if (solverVerbose_) {
+                            OpmLog::error(msg);
+                        }
+                        OPM_THROW_NOLOG(Opm::NumericalIssue, msg);
+                    }
+
+                    const double newTimeStep = restartFactor_ * dt;
+                    // we need to revise this
+                    substepTimer.provideTimeStepEstimate(newTimeStep);
+                    if (solverVerbose_) {
+                        std::string msg;
+                        msg = causeOfFailure + "\nTimestep chopped to "
+                            + std::to_string(unit::convert::to(substepTimer.currentStepLength(), unit::day)) + " days\n";
+                        OpmLog::problem(msg);
+                    }
+
+                    ++restarts;
+                }
+            }
+
+
+            // store estimated time step for next reportStep
+            suggestedNextTimestep_ = substepTimer.currentStepLength();
+            if (timestepVerbose_) {
+                std::ostringstream ss;
+                substepTimer.report(ss);
+                ss << "Suggested next step size = " << unit::convert::to(suggestedNextTimestep_, unit::day) << " (days)" << std::endl;
+                OpmLog::debug(ss.str());
+            }
+
+            if (! std::isfinite(suggestedNextTimestep_)) { // check for NaN
+                suggestedNextTimestep_ = timestep;
+            }
+            return report;
+        }
+
+        /** \brief Returns the simulator report for the failed substeps of the last
+         *         report step.
+         */
+        const SimulatorReport& failureReport() const
+        { return failureReport_; };
+
+        double suggestedNextStep() const
+        { return suggestedNextTimestep_; }
+
+        void setSuggestedNextStep(const double x)
+        { suggestedNextTimestep_ = x; }
+
+        void updateTUNING(const Tuning& tuning, size_t timeStep)
+        {
+            restartFactor_ = tuning.getTSFCNV(timeStep);
+            growthFactor_ = tuning.getTFDIFF(timeStep);
+            maxGrowth_ = tuning.getTSFMAX(timeStep);
+            maxTimeStep_ = tuning.getTSMAXZ(timeStep);
+            suggestedNextTimestep_ = tuning.getTSINIT(timeStep);
+            timestepAfterEvent_ = tuning.getTMAXWC(timeStep);
+        }
+
+
+    protected:
+        void init_(const ParameterGroup& param)
+        {
+            // valid are "pid" and "pid+iteration"
+            std::string control = param.getDefault("timestep.control", std::string("pid"));
+            // iterations is the accumulation of all linear iterations over all newton steops per time step
+            const int defaultTargetIterations = 30;
+            const int defaultTargetNewtonIterations = 8;
+
+            const double tol = param.getDefault("timestep.control.tol", double(1e-1));
+            if (control == "pid") {
+                timeStepControl_ = TimeStepControlType(new PIDTimeStepControl(tol));
+            }
+            else if (control == "pid+iteration") {
+                const int iterations = param.getDefault("timestep.control.targetiteration", defaultTargetIterations);
+                timeStepControl_ = TimeStepControlType(new PIDAndIterationCountTimeStepControl(iterations, tol));
+            }
+            else if (control == "pid+newtoniteration") {
+                const int iterations = param.getDefault("timestep.control.targetiteration", defaultTargetNewtonIterations);
+                timeStepControl_ = TimeStepControlType(new PIDAndIterationCountTimeStepControl(iterations, tol));
+                useNewtonIteration_ = true;
+            }
+            else if (control == "iterationcount") {
+                const int iterations = param.getDefault("timestep.control.targetiteration", defaultTargetIterations);
+                const double decayrate = param.getDefault("timestep.control.decayrate",  double(0.75));
+                const double growthrate = param.getDefault("timestep.control.growthrate", double(1.25));
+                timeStepControl_ = TimeStepControlType(new SimpleIterationCountTimeStepControl(iterations, decayrate, growthrate));
+            }
+            else if (control == "hardcoded") {
+                const std::string filename = param.getDefault("timestep.control.filename", std::string("timesteps"));
+                timeStepControl_ = TimeStepControlType(new HardcodedTimeStepControl(filename));
+
+            }
+            else
+                OPM_THROW(std::runtime_error,"Unsupported time step control selected "<< control);
+
+            // make sure growth factor is something reasonable
+            assert(growthFactor_ >= 1.0);
+        }
+
+
+        typedef std::unique_ptr<TimeStepControlInterface> TimeStepControlType;
+
+        SimulatorReport failureReport_;       //!< statistics for the failed substeps of the last timestep
+        TimeStepControlType timeStepControl_; //!< time step control object
+        double restartFactor_;               //!< factor to multiply time step with when solver fails to converge
+        double growthFactor_;                //!< factor to multiply time step when solver recovered from failed convergence
+        double maxGrowth_;                   //!< factor that limits the maximum growth of a time step
+        double maxTimeStep_;                //!< maximal allowed time step size
+        const int solverRestartMax_;        //!< how many restart of solver are allowed
+        const bool solverVerbose_;           //!< solver verbosity
+        const bool timestepVerbose_;         //!< timestep verbosity
+        double suggestedNextTimestep_;      //!< suggested size of next timestep
+        bool fullTimestepInitially_;        //!< beginning with the size of the time step from data file
+        double timestepAfterEvent_;         //!< suggested size of timestep after an event
+        bool useNewtonIteration_;           //!< use newton iteration count for adaptive time step control
+    };
+}
+
+#endif // OPM_ADAPTIVE_TIME_STEPPING_EBOS_HPP