Merge pull request #381 from atgeirr/refactor-solver

Start refactoring the FullyImplicitBlackoilSolver class

CMakeLists_files.cmake

@@ -94,6 +94,8 @@ list (APPEND PUBLIC_HEADER_FILES
 	opm/autodiff/AutoDiffHelpers.hpp
 	opm/autodiff/AutoDiff.hpp
 	opm/autodiff/BackupRestore.hpp
+	opm/autodiff/BlackoilModel.hpp
+	opm/autodiff/BlackoilModel_impl.hpp
 	opm/autodiff/BlackoilPropsAdFromDeck.hpp
 	opm/autodiff/BlackoilPropsAdInterface.hpp
 	opm/autodiff/CPRPreconditioner.hpp
@@ -108,6 +110,8 @@ list (APPEND PUBLIC_HEADER_FILES
 	opm/autodiff/NewtonIterationBlackoilCPR.hpp
 	opm/autodiff/NewtonIterationBlackoilInterface.hpp
 	opm/autodiff/NewtonIterationBlackoilSimple.hpp
+	opm/autodiff/NewtonSolver.hpp
+	opm/autodiff/NewtonSolver_impl.hpp
 	opm/autodiff/LinearisedBlackoilResidual.hpp
 	opm/autodiff/RateConverter.hpp
 	opm/autodiff/RedistributeDataHandles.hpp

opm/autodiff/BlackoilModel.hpp

@@ -0,0 +1,441 @@
+/*
+  Copyright 2013, 2015 SINTEF ICT, Applied Mathematics.
+  Copyright 2014, 2015 Statoil ASA.
+  Copyright 2014, 2015 Dr. Markus Blatt - HPC-Simulation-Software & Services
+  Copyright 2015 NTNU
+
+  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_BLACKOILMODEL_HEADER_INCLUDED
+#define OPM_BLACKOILMODEL_HEADER_INCLUDED
+
+#include <cassert>
+
+#include <opm/autodiff/AutoDiffBlock.hpp>
+#include <opm/autodiff/AutoDiffHelpers.hpp>
+#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
+#include <opm/autodiff/LinearisedBlackoilResidual.hpp>
+#include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
+
+#include <array>
+
+struct UnstructuredGrid;
+struct Wells;
+
+namespace Opm {
+
+    namespace parameter { class ParameterGroup; }
+    class DerivedGeology;
+    class RockCompressibility;
+    class NewtonIterationBlackoilInterface;
+    class BlackoilState;
+    class WellStateFullyImplicitBlackoil;
+
+
+    /// A model implementation for three-phase black oil.
+    ///
+    /// The simulator is capable of handling three-phase problems
+    /// where gas can be dissolved in oil and vice versa. It
+    /// uses an industry-standard TPFA discretization with per-phase
+    /// upwind weighting of mobilities.
+    ///
+    /// It uses automatic differentiation via the class AutoDiffBlock
+    /// to simplify assembly of the jacobian matrix.
+    template<class Grid>
+    class BlackoilModel
+    {
+    public:
+        // ---------  Types and enums  ---------
+        typedef AutoDiffBlock<double> ADB;
+        typedef ADB::V V;
+        typedef ADB::M M;
+        typedef BlackoilState ReservoirState;
+        typedef WellStateFullyImplicitBlackoil WellState;
+
+        /// Model-specific solver parameters.
+        struct ModelParameters
+        {
+            double dp_max_rel_;
+            double ds_max_;
+            double dr_max_rel_;
+            double max_residual_allowed_;
+            double tolerance_mb_;
+            double tolerance_cnv_;
+            double tolerance_wells_;
+
+            explicit ModelParameters( const parameter::ParameterGroup& param );
+            ModelParameters();
+
+            void reset();
+        };
+
+        // ---------  Public methods  ---------
+
+        /// Construct the model. 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] 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
+        /// \param[in] has_disgas       turn on dissolved gas
+        /// \param[in] has_vapoil       turn on vaporized oil feature
+        /// \param[in] terminal_output  request output to cout/cerr
+        BlackoilModel(const ModelParameters&          param,
+                      const Grid&                     grid ,
+                      const BlackoilPropsAdInterface& fluid,
+                      const DerivedGeology&           geo  ,
+                      const RockCompressibility*      rock_comp_props,
+                      const Wells*                    wells,
+                      const NewtonIterationBlackoilInterface& linsolver,
+                      const bool has_disgas,
+                      const bool has_vapoil,
+                      const bool terminal_output);
+
+        /// \brief Set threshold pressures that prevent or reduce flow.
+        /// This prevents flow across faces if the potential
+        /// difference is less than the threshold. If the potential
+        /// difference is greater, the threshold value is subtracted
+        /// before calculating flow. This is treated symmetrically, so
+        /// flow is prevented or reduced in both directions equally.
+        /// \param[in]  threshold_pressures_by_face   array of size equal to the number of faces
+        ///                                   of the grid passed in the constructor.
+        void setThresholdPressures(const std::vector<double>& threshold_pressures_by_face);
+
+        /// Called once before each time step.
+        /// \param[in] dt                     time step size
+        /// \param[in, out] reservoir_state   reservoir state variables
+        /// \param[in, out] well_state        well state variables
+        void prepareStep(const double dt,
+                         ReservoirState& reservoir_state,
+                         WellState& well_state);
+
+        /// Called once after each time step.
+        /// In this class, this function does nothing.
+        /// \param[in] dt                     time step size
+        /// \param[in, out] reservoir_state   reservoir state variables
+        /// \param[in, out] well_state        well state variables
+        void afterStep(const double dt,
+                       ReservoirState& reservoir_state,
+                       WellState& well_state);
+
+        /// Assemble the residual and Jacobian of the nonlinear system.
+        /// \param[in]      reservoir_state   reservoir state variables
+        /// \param[in, out] well_state        well state variables
+        /// \param[in]      initial_assembly  pass true if this is the first call to assemble() in this timestep
+        void assemble(const BlackoilState& reservoir_state,
+                      WellStateFullyImplicitBlackoil& well_state,
+                      const bool initial_assembly);
+
+        /// \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;
+
+        /// The size (number of unknowns) of the nonlinear system of equations.
+        int sizeNonLinear() const;
+
+        /// Number of linear iterations used in last call to solveJacobianSystem().
+        int linearIterationsLastSolve() const;
+
+        /// Solve the Jacobian system Jx = r where J is the Jacobian and
+        /// r is the residual.
+        V solveJacobianSystem() const;
+
+        /// Apply an update to the primary variables, chopped if appropriate.
+        /// \param[in]      dx                updates to apply to primary variables
+        /// \param[in, out] reservoir_state   reservoir state variables
+        /// \param[in, out] well_state        well state variables
+        void updateState(const V& dx,
+                         BlackoilState& reservoir_state,
+                         WellStateFullyImplicitBlackoil& well_state);
+
+        /// Return true if output to cout is wanted.
+        bool terminalOutputEnabled() const;
+
+        /// Compute convergence based on total mass balance (tol_mb) and maximum
+        /// residual mass balance (tol_cnv).
+        /// \param[in]   dt          timestep length
+        /// \param[in]   iteration   current iteration number
+        bool getConvergence(const double dt, const int iteration);
+
+        /// The number of active phases in the model.
+        int numPhases() const;
+
+    private:
+
+        // ---------  Types and enums  ---------
+
+        typedef Eigen::Array<double,
+                             Eigen::Dynamic,
+                             Eigen::Dynamic,
+                             Eigen::RowMajor> DataBlock;
+
+        struct ReservoirResidualQuant {
+            ReservoirResidualQuant();
+            std::vector<ADB> accum; // Accumulations
+            ADB              mflux; // Mass flux (surface conditions)
+            ADB              b;     // Reciprocal FVF
+            ADB              head;  // Pressure drop across int. interfaces
+            ADB              mob;   // Phase mobility (per cell)
+        };
+
+        struct SolutionState {
+            SolutionState(const int np);
+            ADB              pressure;
+            ADB              temperature;
+            std::vector<ADB> saturation;
+            ADB              rs;
+            ADB              rv;
+            ADB              qs;
+            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);
+            M w2p;              // well -> perf (scatter)
+            M p2w;              // perf -> well (gather)
+        };
+
+        enum { Water        = BlackoilPropsAdInterface::Water,
+               Oil          = BlackoilPropsAdInterface::Oil  ,
+               Gas          = BlackoilPropsAdInterface::Gas  ,
+               MaxNumPhases = BlackoilPropsAdInterface::MaxNumPhases
+         };
+
+        enum PrimalVariables { Sg = 0, RS = 1, RV = 2 };
+
+        // ---------  Data members  ---------
+
+        const Grid&         grid_;
+        const BlackoilPropsAdInterface& fluid_;
+        const DerivedGeology&           geo_;
+        const RockCompressibility*      rock_comp_props_;
+        const Wells*                    wells_;
+        const NewtonIterationBlackoilInterface&    linsolver_;
+        // For each canonical phase -> true if active
+        const std::vector<bool>         active_;
+        // Size = # active phases. Maps active -> canonical phase indices.
+        const std::vector<int>          canph_;
+        const std::vector<int>          cells_;  // All grid cells
+        HelperOps                       ops_;
+        const WellOps                   wops_;
+        const bool has_disgas_;
+        const bool has_vapoil_;
+
+        ModelParameters                 param_;
+        bool use_threshold_pressure_;
+        V threshold_pressures_by_interior_face_;
+
+        std::vector<ReservoirResidualQuant> rq_;
+        std::vector<PhasePresence> phaseCondition_;
+        V well_perforation_pressure_diffs_; // Diff to bhp for each well perforation.
+
+        LinearisedBlackoilResidual residual_;
+
+        /// \brief Whether we print something to std::cout
+        bool terminal_output_;
+
+        std::vector<int>         primalVariable_;
+        V pvdt_;
+
+        // ---------  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 BlackoilState& x,
+                      const WellStateFullyImplicitBlackoil& xw) const;
+
+        void
+        makeConstantState(SolutionState& state) const;
+
+        SolutionState
+        variableState(const BlackoilState& x,
+                      const WellStateFullyImplicitBlackoil& xw) const;
+
+        void
+        computeAccum(const SolutionState& state,
+                     const int            aix  );
+
+        void computeWellConnectionPressures(const SolutionState& state,
+                                            const WellStateFullyImplicitBlackoil& xw);
+
+        void
+        addWellControlEq(const SolutionState& state,
+                         const WellStateFullyImplicitBlackoil& xw,
+                         const V& aliveWells);
+
+        void
+        addWellEq(const SolutionState& state,
+                  WellStateFullyImplicitBlackoil& xw,
+                  V& aliveWells);
+
+        void updateWellControls(WellStateFullyImplicitBlackoil& xw) const;
+
+        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;
+
+        void
+        computeMassFlux(const int               actph ,
+                        const V&                transi,
+                        const ADB&              kr    ,
+                        const ADB&              p     ,
+                        const SolutionState&    state );
+
+        void applyThresholdPressures(ADB& dp);
+
+        ADB
+        fluidViscosity(const int               phase,
+                       const ADB&              p    ,
+                       const ADB&              temp ,
+                       const ADB&              rs   ,
+                       const ADB&              rv   ,
+                       const std::vector<PhasePresence>& cond,
+                       const std::vector<int>& cells) const;
+
+        ADB
+        fluidReciprocFVF(const int               phase,
+                         const ADB&              p    ,
+                         const ADB&              temp ,
+                         const ADB&              rs   ,
+                         const ADB&              rv   ,
+                         const std::vector<PhasePresence>& cond,
+                         const std::vector<int>& cells) const;
+
+        ADB
+        fluidDensity(const int               phase,
+                     const ADB&              p    ,
+                     const ADB&              temp ,
+                     const ADB&              rs   ,
+                     const ADB&              rv   ,
+                     const std::vector<PhasePresence>& cond,
+                     const std::vector<int>& cells) const;
+
+        V
+        fluidRsSat(const V&                p,
+                   const V&                so,
+                   const std::vector<int>& cells) const;
+
+        ADB
+        fluidRsSat(const ADB&              p,
+                   const ADB&              so,
+                   const std::vector<int>& cells) const;
+
+        V
+        fluidRvSat(const V&                p,
+                   const V&                so,
+                   const std::vector<int>& cells) const;
+
+        ADB
+        fluidRvSat(const ADB&              p,
+                   const ADB&              so,
+                   const std::vector<int>& cells) const;
+
+        ADB
+        poroMult(const ADB& p) const;
+
+        ADB
+        transMult(const ADB& p) const;
+
+        void
+        classifyCondition(const SolutionState&        state,
+                          std::vector<PhasePresence>& cond ) const;
+
+        const std::vector<PhasePresence>
+        phaseCondition() const {return phaseCondition_;}
+
+        void
+        classifyCondition(const BlackoilState&        state);
+
+
+        /// update the primal variable for Sg, Rv or Rs. The Gas phase must
+        /// be active to call this method.
+        void
+        updatePrimalVariableFromState(const BlackoilState&        state);
+
+        /// Update the phaseCondition_ member based on the primalVariable_ member.
+        void
+        updatePhaseCondFromPrimalVariable();
+
+        /// \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[out] maxNormWell The maximum of the well equations for each phase.
+        /// \param[in]  nc    The number of cells of the local grid.
+        /// \param[in]  nw    The number of wells on the local grid.
+        /// \return The total pore volume over all cells.
+        double
+        convergenceReduction(const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases>& B,
+                             const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases>& tempV,
+                             const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases>& R,
+                             std::array<double,MaxNumPhases>& R_sum,
+                             std::array<double,MaxNumPhases>& maxCoeff,
+                             std::array<double,MaxNumPhases>& B_avg,
+                             std::vector<double>& maxNormWell,
+                             int nc,
+                             int nw) const;
+
+        double dpMaxRel() const { return param_.dp_max_rel_; }
+        double dsMax() const { return param_.ds_max_; }
+        double drMaxRel() const { return param_.dr_max_rel_; }
+        double maxResidualAllowed() const { return param_.max_residual_allowed_; }
+
+    };
+} // namespace Opm
+
+#include "BlackoilModel_impl.hpp"
+
+#endif // OPM_BLACKOILMODEL_HEADER_INCLUDED

opm/autodiff/NewtonSolver.hpp

@@ -0,0 +1,119 @@
+/*
+  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_NEWTONSOLVER_HEADER_INCLUDED
+#define OPM_NEWTONSOLVER_HEADER_INCLUDED
+
+#include <opm/autodiff/AutoDiffBlock.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+
+namespace Opm {
+
+
+    /// A Newton solver class suitable for general fully-implicit models.
+    template <class PhysicalModel>
+    class NewtonSolver
+    {
+    public:
+        // ---------  Types and enums  ---------
+        typedef AutoDiffBlock<double> ADB;
+        typedef ADB::V V;
+        typedef ADB::M M;
+
+        // The Newton relaxation scheme type
+        enum RelaxType { DAMPEN, SOR };
+
+        // Solver parameters controlling nonlinear Newton process.
+        struct SolverParameters
+        {
+            enum RelaxType relax_type_;
+            double         relax_max_;
+            double         relax_increment_;
+            double         relax_rel_tol_;
+            int            max_iter_; // max newton iterations
+            int            min_iter_; // min newton iterations
+
+            explicit SolverParameters( const parameter::ParameterGroup& param );
+            SolverParameters();
+
+            void reset();
+        };
+
+        // Forwarding types from PhysicalModel.
+        typedef typename PhysicalModel::ReservoirState ReservoirState;
+        typedef typename PhysicalModel::WellState WellState;
+
+        // ---------  Public methods  ---------
+
+        /// Construct solver for a given model.
+        /// \param[in]      param   parameters controlling nonlinear Newton process
+        /// \param[in, out] model   physical simulation model
+        explicit NewtonSolver(const SolverParameters& param,
+                              PhysicalModel& model);
+
+        /// Take a single forward step, after which the states will be modified
+        /// according to the physical model.
+        /// \param[in] dt                time step size
+        /// \param[in] reservoir_state   reservoir state variables
+        /// \param[in] well_state        well state variables
+        /// \return                      number of linear iterations used
+        int
+        step(const double dt,
+             ReservoirState& reservoir_state,
+             WellState& well_state);
+
+        /// Number of Newton iterations used in all calls to step().
+        unsigned int newtonIterations() const;
+
+        /// Number of linear solver iterations used in all calls to step().
+        unsigned int linearIterations() const;
+
+        /// Number of linear solver iterations used in the last call to step().
+        unsigned int newtonIterationsLastStep() const;
+
+        /// Number of linear solver iterations used in the last call to step().
+        unsigned int linearIterationsLastStep() const;
+
+    private:
+        // ---------  Data members  ---------
+        SolverParameters param_;
+        PhysicalModel& model_;
+        unsigned int newtonIterations_;
+        unsigned int linearIterations_;
+        unsigned int newtonIterationsLast_;
+        unsigned int linearIterationsLast_;
+
+        // ---------  Private methods  ---------
+        enum RelaxType relaxType() const { return param_.relax_type_; }
+        double relaxMax() const          { return param_.relax_max_; }
+        double relaxIncrement() const    { return param_.relax_increment_; }
+        double relaxRelTol() const       { return param_.relax_rel_tol_; }
+        double maxIter() const           { return param_.max_iter_; }
+        double minIter() const           { return param_.min_iter_; }
+        void detectNewtonOscillations(const std::vector<std::vector<double>>& residual_history,
+                                      const int it, const double relaxRelTol,
+                                      bool& oscillate, bool& stagnate) const;
+        void stabilizeNewton(V& dx, V& dxOld, const double omega, const RelaxType relax_type) const;
+    };
+} // namespace Opm
+
+#include "NewtonSolver_impl.hpp"
+
+#endif // OPM_NEWTONSOLVER_HEADER_INCLUDED

opm/autodiff/NewtonSolver_impl.hpp

@@ -0,0 +1,250 @@
+/*
+  Copyright 2013, 2015 SINTEF ICT, Applied Mathematics.
+  Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
+  Copyright 2015 NTNU
+  Copyright 2015 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_NEWTONSOLVER_IMPL_HEADER_INCLUDED
+#define OPM_NEWTONSOLVER_IMPL_HEADER_INCLUDED
+
+#include <opm/autodiff/NewtonSolver.hpp>
+
+namespace Opm
+{
+    template <class PhysicalModel>
+    NewtonSolver<PhysicalModel>::NewtonSolver(const SolverParameters& param,
+                                              PhysicalModel& model)
+        : param_(param),
+          model_(model),
+          newtonIterations_(0),
+          linearIterations_(0)
+    {
+    }
+
+    template <class PhysicalModel>
+    unsigned int NewtonSolver<PhysicalModel>::newtonIterations () const
+    {
+        return newtonIterations_;
+    }
+
+    template <class PhysicalModel>
+    unsigned int NewtonSolver<PhysicalModel>::linearIterations () const
+    {
+        return linearIterations_;
+    }
+
+
+    template <class PhysicalModel>
+    int
+    NewtonSolver<PhysicalModel>::
+    step(const double dt,
+         ReservoirState& reservoir_state,
+         WellState& well_state)
+    {
+        // Do model-specific once-per-step calculations.
+        model_.prepareStep(dt, reservoir_state, well_state);
+
+        // For each iteration we store in a vector the norms of the residual of
+        // the mass balance for each active phase, the well flux and the well equations.
+        std::vector<std::vector<double>> residual_norms_history;
+
+        // Assemble residual and Jacobian, store residual norms.
+        model_.assemble(reservoir_state, well_state, true);
+        residual_norms_history.push_back(model_.computeResidualNorms());
+
+        // Set up for main Newton loop.
+        double omega = 1.0;
+        int iteration = 0;
+        bool converged = model_.getConvergence(dt, iteration);
+        const int sizeNonLinear = model_.sizeNonLinear();
+        V dxOld = V::Zero(sizeNonLinear);
+        bool isOscillate = false;
+        bool isStagnate = false;
+        const enum RelaxType relaxtype = relaxType();
+        int linearIterations = 0;
+
+        // ----------  Main Newton loop  ----------
+        while ( (!converged && (iteration < maxIter())) || (minIter() > iteration)) {
+            // Compute the Newton update to the primary variables.
+            V dx = model_.solveJacobianSystem();
+
+            // Store number of linear iterations used.
+            linearIterations += model_.linearIterationsLastSolve();
+
+            // Stabilize the Newton update.
+            detectNewtonOscillations(residual_norms_history, iteration, relaxRelTol(), isOscillate, isStagnate);
+            if (isOscillate) {
+                omega -= relaxIncrement();
+                omega = std::max(omega, relaxMax());
+                if (model_.terminalOutputEnabled()) {
+                    std::cout << " Oscillating behavior detected: Relaxation set to " << omega << std::endl;
+                }
+            }
+            stabilizeNewton(dx, dxOld, omega, relaxtype);
+
+            // Apply the update, the model may apply model-dependent
+            // limitations and chopping of the update.
+            model_.updateState(dx, reservoir_state, well_state);
+
+            // Assemble residual and Jacobian, store residual norms.
+            model_.assemble(reservoir_state, well_state, false);
+            residual_norms_history.push_back(model_.computeResidualNorms());
+
+            // increase iteration counter
+            ++iteration;
+
+            converged = model_.getConvergence(dt, iteration);
+        }
+
+        if (!converged) {
+            if (model_.terminalOutputEnabled()) {
+                std::cerr << "WARNING: Failed to compute converged solution in " << iteration << " iterations." << std::endl;
+            }
+            return -1; // -1 indicates that the solver has to be restarted
+        }
+
+        linearIterations_ += linearIterations;
+        newtonIterations_ += iteration;
+        linearIterationsLast_ = linearIterations;
+        newtonIterationsLast_ = iteration;
+
+        // Do model-specific post-step actions.
+        model_.afterStep(dt, reservoir_state, well_state);
+
+        return linearIterations;
+    }
+
+
+
+    template <class PhysicalModel>
+    void NewtonSolver<PhysicalModel>::SolverParameters::
+    reset()
+    {
+        // default values for the solver parameters
+        relax_type_      = DAMPEN;
+        relax_max_       = 0.5;
+        relax_increment_ = 0.1;
+        relax_rel_tol_   = 0.2;
+        max_iter_        = 15;
+        min_iter_        = 1;
+    }
+
+    template <class PhysicalModel>
+    NewtonSolver<PhysicalModel>::SolverParameters::
+    SolverParameters()
+    {
+        // set default values
+        reset();
+    }
+
+    template <class PhysicalModel>
+    NewtonSolver<PhysicalModel>::SolverParameters::
+    SolverParameters( const parameter::ParameterGroup& param )
+    {
+        // set default values
+        reset();
+
+        // overload with given parameters
+        relax_max_   = param.getDefault("relax_max", relax_max_);
+        max_iter_    = param.getDefault("max_iter", max_iter_);
+        min_iter_    = param.getDefault("min_iter", min_iter_);
+
+        std::string relaxation_type = param.getDefault("relax_type", std::string("dampen"));
+        if (relaxation_type == "dampen") {
+            relax_type_ = DAMPEN;
+        } else if (relaxation_type == "sor") {
+            relax_type_ = SOR;
+        } else {
+            OPM_THROW(std::runtime_error, "Unknown Relaxtion Type " << relaxation_type);
+        }
+    }
+
+    template <class PhysicalModel>
+    void
+    NewtonSolver<PhysicalModel>::detectNewtonOscillations(const std::vector<std::vector<double>>& residual_history,
+                                                          const int it, const double relaxRelTol,
+                                                          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 = residual_history[it];
+        const std::vector<double>& F1 = residual_history[it - 1];
+        const std::vector<double>& F2 = residual_history[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);
+    }
+
+
+    template <class PhysicalModel>
+    void
+    NewtonSolver<PhysicalModel>::stabilizeNewton(V& dx, V& dxOld, const double omega,
+                                                 const RelaxType relax_type) const
+    {
+        // The dxOld is updated with dx.
+        // If omega is equal to 1., no relaxtion will be appiled.
+
+        const V tempDxOld = dxOld;
+        dxOld = dx;
+
+        switch (relax_type) {
+            case DAMPEN:
+                if (omega == 1.) {
+                    return;
+                }
+                dx = dx*omega;
+                return;
+            case SOR:
+                if (omega == 1.) {
+                    return;
+                }
+                dx = dx*omega + (1.-omega)*tempDxOld;
+                return;
+            default:
+                OPM_THROW(std::runtime_error, "Can only handle DAMPEN and SOR relaxation type.");
+        }
+
+        return;
+    }
+
+
+} // namespace Opm
+
+
+#endif // OPM_FULLYIMPLICITSOLVER_IMPL_HEADER_INCLUDED

opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp

@@ -24,7 +24,8 @@
 #include <opm/core/utility/ErrorMacros.hpp>
 
 #include <opm/autodiff/GeoProps.hpp>
-#include <opm/autodiff/FullyImplicitBlackoilSolver.hpp>
+#include <opm/autodiff/NewtonSolver.hpp>
+#include <opm/autodiff/BlackoilModel.hpp>
 #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 #include <opm/autodiff/RateConverter.hpp>
@@ -231,7 +232,13 @@ namespace Opm
         std::string tstep_filename = output_writer_.outputDirectory() + "/step_timing.txt";
         std::ofstream tstep_os(tstep_filename.c_str());
 
-        typename FullyImplicitBlackoilSolver<T>::SolverParameter solverParam( param_ );
+        typedef T Grid;
+        typedef BlackoilModel<Grid> Model;
+        typedef typename Model::ModelParameters ModelParams;
+        ModelParams modelParams( param_ );
+        typedef NewtonSolver<Model> Solver;
+        typedef typename Solver::SolverParameters SolverParams;
+        SolverParams solverParams( param_ );
 
         // adaptive time stepping
         std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping;
@@ -291,10 +298,11 @@ namespace Opm
             // Run a multiple steps of the solver depending on the time step control.
             solver_timer.start();
 
-            FullyImplicitBlackoilSolver<T> solver(solverParam, grid_, props_, geo_, rock_comp_props_, wells, solver_, has_disgas_, has_vapoil_, terminal_output_);
+            Model model(modelParams, grid_, props_, geo_, rock_comp_props_, wells, solver_, has_disgas_, has_vapoil_, terminal_output_);
             if (!threshold_pressures_by_face_.empty()) {
-                solver.setThresholdPressures(threshold_pressures_by_face_);
+                model.setThresholdPressures(threshold_pressures_by_face_);
             }
+            Solver solver(solverParams, model);
 
             // If sub stepping is enabled allow the solver to sub cycle
             // in case the report steps are to large for the solver to converge