Make BlackoilPolymerModel usable with NewtonSolver from opm-autodiff.

opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp

@@ -43,15 +43,15 @@ namespace Opm {
     class RockCompressibility;
     class NewtonIterationBlackoilInterface;
     class PolymerBlackoilState;
-    class WellStateFullyImplicitBlackoil;
+    class WellStateFullyImplicitBlackoilPolymer;
 
 
-    /// A fully implicit solver for the black-oil-polymer problem.
+    /// A model implementation for three-phase black oil with polymer.
     ///
     /// The simulator is capable of handling three-phase problems
-    /// where gas can be dissolved in oil (but not vice versa). It
+    /// where gas can be dissolved in oil and vice versa, with polymer
-    /// uses an industry-standard TPFA discretization with per-phase
+    /// in the water phase. It uses an industry-standard TPFA
-    /// upwind weighting of mobilities.
+    /// discretization with per-phase upwind weighting of mobilities.
     ///
     /// It uses automatic differentiation via the class AutoDiffBlock
     /// to simplify assembly of the jacobian matrix.
@@ -59,33 +59,35 @@ namespace Opm {
     class BlackoilPolymerModel
     {
     public:
-        // the Newton relaxation type
-        enum RelaxType { DAMPEN, SOR };
 
-        // class holding the solver parameters
+        // ---------  Types and enums  ---------
-        struct SolverParameter
+
+        typedef AutoDiffBlock<double> ADB;
+        typedef ADB::V V;
+        typedef ADB::M M;
+        typedef PolymerBlackoilState ReservoirState;
+        typedef WellStateFullyImplicitBlackoilPolymer WellState;
+
+        /// Model-specific solver parameters.
+        struct ModelParameter
         {
             double dp_max_rel_;
             double ds_max_;
             double dr_max_rel_;
-            enum RelaxType                  relax_type_;
-            double                          relax_max_;
-            double                          relax_increment_;
-            double                          relax_rel_tol_;
             double max_residual_allowed_;
             double tolerance_mb_;
             double tolerance_cnv_;
             double tolerance_wells_;
-            int                             max_iter_; // max newton iterations
-            int                             min_iter_; // min newton iterations
 
-            SolverParameter( const parameter::ParameterGroup& param );
+            ModelParameter( const parameter::ParameterGroup& param );
-            SolverParameter();
+            ModelParameter();
 
             void reset();
         };
 
-        /// Construct a solver. It will retain references to the
+        // ---------  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
@@ -95,7 +97,11 @@ namespace Opm {
         /// \param[in] rock_comp_props  if non-null, rock compressibility properties
         /// \param[in] wells            well structure
         /// \param[in] linsolver        linear solver
-        BlackoilPolymerModel(const SolverParameter&          param,
+        /// \param[in] has_disgas       turn on dissolved gas
+        /// \param[in] has_vapoil       turn on vaporized oil feature
+        /// \param[in] has_polymer      turn on polymer feature
+        /// \param[in] terminal_output  request output to cout/cerr
+        BlackoilPolymerModel(const ModelParameter&          param,
                              const Grid&                     grid ,
                              const BlackoilPropsAdInterface& fluid,
                              const DerivedGeology&           geo  ,
@@ -118,30 +124,70 @@ namespace Opm {
         ///                                   of the grid passed in the constructor.
         void setThresholdPressures(const std::vector<double>& threshold_pressures_by_face);
 
-        /// Take a single forward step, modifiying
+        /// Called once before each time step.
-        ///   state.pressure()
-        ///   state.faceflux()
-        ///   state.saturation()
-        ///   state.gasoilratio()
-        ///   wstate.bhp()
         /// \param[in] dt                time step size
-        /// \param[in] state     reservoir state
+        /// \param[in] reservoir_state   reservoir state variables
-        /// \param[in] wstate    well state
+        /// \param[in] well_state        well state variables
-        /// \return              number of linear iterations used
+        void prepareStep(const double dt,
-        int
+                         ReservoirState& reservoir_state,
-        step(const double   dt    ,
+                         WellState& well_state);
-             PolymerBlackoilState& state ,
-             WellStateFullyImplicitBlackoil& wstate,
-             const std::vector<double>& polymer_inflow);
 
-        unsigned int newtonIterations () const { return newtonIterations_; }
+        /// Called once after each time step.
-        unsigned int linearIterations () const { return linearIterations_; }
+        /// \param[in] dt                time step size
+        /// \param[in] reservoir_state   reservoir state variables
+        /// \param[in] 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 ReservoirState& reservoir_state,
+                      WellState& 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,
+                         ReservoirState& reservoir_state,
+                         WellState& well_state);
+
+        /// Return true if output to cout is wanted.
+        bool terminalOutput() 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 AutoDiffBlock<double> ADB;
+        // ---------  Types and enums  ---------
-        typedef ADB::V V;
+
-        typedef ADB::M M;
         typedef Eigen::Array<double,
                              Eigen::Dynamic,
                              Eigen::Dynamic,
@@ -184,7 +230,8 @@ namespace Opm {
 
         enum PrimalVariables { Sg = 0, RS = 1, RV = 2 };
 
-        // Member data
+        // ---------  Data members  ---------
+
         const Grid&         grid_;
         const BlackoilPropsAdInterface& fluid_;
         const DerivedGeology&           geo_;
@@ -205,7 +252,7 @@ namespace Opm {
         const bool has_polymer_;
         const int  poly_pos_;
 
-        SolverParameter                 param_;
+        ModelParameter                 param_;
         bool use_threshold_pressure_;
         V threshold_pressures_by_interior_face_;
 
@@ -221,8 +268,9 @@ namespace Opm {
         unsigned int linearIterations_;
 
         std::vector<int>         primalVariable_;
+        V pvdt_;
 
-        // Private methods.
+        // ---------  Private methods  ---------
 
         // return true if wells are available
         bool wellsActive() const { return wells_ ? wells_->number_of_wells > 0 : false ; }
@@ -231,47 +279,33 @@ namespace Opm {
 
         SolutionState
         constantState(const PolymerBlackoilState& x,
-                      const WellStateFullyImplicitBlackoil& xw) const;
+                      const WellStateFullyImplicitBlackoilPolymer& xw) const;
 
         void
         makeConstantState(SolutionState& state) const;
 
         SolutionState
         variableState(const PolymerBlackoilState& x,
-                      const WellStateFullyImplicitBlackoil& xw) const;
+                      const WellStateFullyImplicitBlackoilPolymer& xw) const;
 
         void
         computeAccum(const SolutionState& state,
                      const int            aix  );
 
         void computeWellConnectionPressures(const SolutionState& state,
-                                            const WellStateFullyImplicitBlackoil& xw);
+                                            const WellStateFullyImplicitBlackoilPolymer& xw);
 
         void
         addWellControlEq(const SolutionState& state,
-                         const WellStateFullyImplicitBlackoil& xw,
+                         const WellStateFullyImplicitBlackoilPolymer& xw,
                          const V& aliveWells);
 
         void
         addWellEq(const SolutionState& state,
-                  WellStateFullyImplicitBlackoil& xw,
+                  WellStateFullyImplicitBlackoilPolymer& xw,
-                  V& aliveWells,
+                  V& aliveWells);
-                  const std::vector<double>& polymer_inflow);
 
-        void updateWellControls(WellStateFullyImplicitBlackoil& xw) const;
+        void updateWellControls(WellStateFullyImplicitBlackoilPolymer& xw) const;
-
-        void
-        assemble(const V&             dtpv,
-                 const PolymerBlackoilState& x,
-                 const bool initial_assembly,
-                 WellStateFullyImplicitBlackoil& xw,
-                 const std::vector<double>& polymer_inflow);
-
-        V solveJacobianSystem() const;
-
-        void updateState(const V& dx,
-                         PolymerBlackoilState& state,
-                         WellStateFullyImplicitBlackoil& well_state);
 
         std::vector<ADB>
         computePressures(const SolutionState& state) const;
@@ -306,15 +340,6 @@ namespace Opm {
 
         void applyThresholdPressures(ADB& dp);
 
-        double
-        residualNorm() const;
-
-        /// \brief Compute the residual norms of the mass balance for each phase,
-        /// the well flux, and the well equation.
-        /// \return a vector that contains for each phase the norm of the mass balance
-        /// and afterwards the norm of the residual of the well flux and the well equation.
-        std::vector<double> computeResidualNorms() const;
-
         ADB
         fluidViscosity(const int               phase,
                        const ADB&              p    ,
@@ -388,10 +413,6 @@ namespace Opm {
         void
         updatePhaseCondFromPrimalVariable();
 
-        /// Compute convergence based on total mass balance (tol_mb) and maximum
-        /// residual mass balance (tol_cnv).
-        bool getConvergence(const double dt, const int iteration);
-
         /// \brief Compute the reduction within the convergence check.
         /// \param[in] B     A matrix with MaxNumPhases columns and the same number rows
         ///                  as the number of cells of the grid. B.col(i) contains the values
@@ -420,21 +441,9 @@ namespace Opm {
                              std::array<double,MaxNumPhases+1>& B_avg,
                              int nc) const;
 
-        void detectNewtonOscillations(const std::vector<std::vector<double>>& residual_history,
-                                      const int it, const double relaxRelTol,
-                                      bool& oscillate, bool& stagnate) const;
-
-        void stablizeNewton(V& dx, V& dxOld, const double omega, const RelaxType relax_type) const;
-
         double dpMaxRel() const { return param_.dp_max_rel_; }
         double dsMax() const { return param_.ds_max_; }
         double drMaxRel() const { return param_.dr_max_rel_; }
-        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_; }
         double maxResidualAllowed() const { return param_.max_residual_allowed_; }
 
     };

opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp

@@ -25,6 +25,7 @@
 
 #include <opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp>
 #include <opm/polymer/PolymerBlackoilState.hpp>
+#include <opm/polymer/fullyimplicit/WellStateFullyImplicitBlackoilPolymer.hpp>
 
 #include <opm/autodiff/AutoDiffBlock.hpp>
 #include <opm/autodiff/AutoDiffHelpers.hpp>
@@ -32,7 +33,6 @@
 #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
 #include <opm/autodiff/GeoProps.hpp>
 #include <opm/autodiff/WellDensitySegmented.hpp>
-#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 
 #include <opm/core/grid.h>
 #include <opm/core/linalg/LinearSolverInterface.hpp>
@@ -151,19 +151,13 @@ namespace detail {
 } // namespace detail
 
     template <class Grid>
-    void BlackoilPolymerModel<Grid>::SolverParameter::
+    void BlackoilPolymerModel<Grid>::ModelParameter::
     reset()
     {
         // default values for the solver parameters
         dp_max_rel_      = 1.0e9;
         ds_max_          = 0.2;
         dr_max_rel_      = 1.0e9;
-        relax_type_      = DAMPEN;
-        relax_max_       = 0.5;
-        relax_increment_ = 0.1;
-        relax_rel_tol_   = 0.2;
-        max_iter_        = 15; // not more then 15 its by default
-        min_iter_        = 1;  // Default to always do at least one nonlinear iteration.
         max_residual_allowed_ = 1e7;
         tolerance_mb_    = 1.0e-5;
         tolerance_cnv_   = 1.0e-2;
@@ -172,16 +166,16 @@ namespace detail {
     }
 
     template <class Grid>
-    BlackoilPolymerModel<Grid>::SolverParameter::
+    BlackoilPolymerModel<Grid>::ModelParameter::
-    SolverParameter()
+    ModelParameter()
     {
         // set default values
         reset();
     }
 
     template <class Grid>
-    BlackoilPolymerModel<Grid>::SolverParameter::
+    BlackoilPolymerModel<Grid>::ModelParameter::
-    SolverParameter( const parameter::ParameterGroup& param )
+    ModelParameter( const parameter::ParameterGroup& param )
     {
         // set default values
         reset();
@@ -190,29 +184,16 @@ namespace detail {
         dp_max_rel_  = param.getDefault("dp_max_rel", dp_max_rel_);
         ds_max_      = param.getDefault("ds_max", ds_max_);
         dr_max_rel_  = param.getDefault("dr_max_rel", dr_max_rel_);
-        relax_max_   = param.getDefault("relax_max", relax_max_);
-        max_iter_    = param.getDefault("max_iter", max_iter_);
-        min_iter_    = param.getDefault("min_iter", min_iter_);
         max_residual_allowed_ = param.getDefault("max_residual_allowed", max_residual_allowed_);
-
         tolerance_mb_    = param.getDefault("tolerance_mb", tolerance_mb_);
         tolerance_cnv_   = param.getDefault("tolerance_cnv", tolerance_cnv_);
         tolerance_wells_ = param.getDefault("tolerance_wells", tolerance_wells_ );
-
-        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 Grid>
     BlackoilPolymerModel<Grid>::
-    BlackoilPolymerModel(const SolverParameter&          param,
+    BlackoilPolymerModel(const ModelParameter&          param,
                          const Grid&                     grid ,
                          const BlackoilPropsAdInterface& fluid,
                          const DerivedGeology&           geo  ,
@@ -276,6 +257,81 @@ namespace detail {
 
 
 
+    template <class Grid>
+    void
+    BlackoilPolymerModel<Grid>::
+    prepareStep(const double dt,
+                ReservoirState& reservoir_state,
+                WellState& /* well_state */)
+    {
+        pvdt_ = geo_.poreVolume() / dt;
+        if (active_[Gas]) {
+            updatePrimalVariableFromState(reservoir_state);
+        }
+        // Initial max concentration of this time step from PolymerBlackoilState.
+        cmax_ = Eigen::Map<const V>(reservoir_state.maxconcentration().data(), Opm::AutoDiffGrid::numCells(grid_));
+    }
+
+
+
+
+    template <class Grid>
+    void
+    BlackoilPolymerModel<Grid>::
+    afterStep(const double /* dt */,
+              ReservoirState& reservoir_state,
+              WellState& /* well_state */)
+    {
+        computeCmax(reservoir_state);
+    }
+
+
+
+
+    template <class Grid>
+    int
+    BlackoilPolymerModel<Grid>::
+    sizeNonLinear() const
+    {
+        return residual_.sizeNonLinear();
+    }
+
+
+
+
+    template <class Grid>
+    int
+    BlackoilPolymerModel<Grid>::
+    linearIterationsLastSolve() const
+    {
+        return linsolver_.iterations();
+    }
+
+
+
+
+    template <class Grid>
+    bool
+    BlackoilPolymerModel<Grid>::
+    terminalOutput() const
+    {
+        return terminal_output_;
+    }
+
+
+
+
+    template <class Grid>
+    int
+    BlackoilPolymerModel<Grid>::
+    numPhases() const
+    {
+        return fluid_.numPhases();
+    }
+
+
+
+
     template <class Grid>
     void
     BlackoilPolymerModel<Grid>::
@@ -297,89 +353,6 @@ namespace detail {
 
 
 
-    template <class Grid>
-    int
-    BlackoilPolymerModel<Grid>::
-    step(const double   dt,
-         PolymerBlackoilState& x ,
-         WellStateFullyImplicitBlackoil& xw,
-         const std::vector<double>& polymer_inflow)
-    {
-        const V pvdt = geo_.poreVolume() / dt;
-
-        // Initial max concentration of this time step from PolymerBlackoilState.
-        cmax_ = Eigen::Map<V>(&x.maxconcentration()[0], Opm::AutoDiffGrid::numCells(grid_));
-        if (active_[Gas]) { updatePrimalVariableFromState(x); }
-
-        // 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(pvdt, x, true, xw, polymer_inflow);
-
-
-        bool converged = false;
-        double omega = 1.;
-
-        residual_norms_history.push_back(computeResidualNorms());
-
-        int          it  = 0;
-        converged = getConvergence(dt,it);
-        const int sizeNonLinear = residual_.sizeNonLinear();
-
-        V dxOld = V::Zero(sizeNonLinear);
-
-        bool isOscillate = false;
-        bool isStagnate = false;
-        const enum RelaxType relaxtype = relaxType();
-        int linearIterations = 0;
-
-        while ((!converged && (it < maxIter())) || (minIter() > it)) {
-            V dx = solveJacobianSystem();
-
-            // store number of linear iterations used
-            linearIterations += linsolver_.iterations();
-
-            detectNewtonOscillations(residual_norms_history, it, relaxRelTol(), isOscillate, isStagnate);
-
-            if (isOscillate) {
-                omega -= relaxIncrement();
-                omega = std::max(omega, relaxMax());
-                if (terminal_output_)
-                {
-                    std::cout << " Oscillating behavior detected: Relaxation set to " << omega << std::endl;
-                }
-            }
-
-            stablizeNewton(dx, dxOld, omega, relaxtype);
-
-            updateState(dx, x, xw);
-
-            assemble(pvdt, x, false, xw, polymer_inflow);
-
-            residual_norms_history.push_back(computeResidualNorms());
-
-            // increase iteration counter
-            ++it;
-
-            converged = getConvergence(dt,it);
-        }
-
-        if (!converged) {
-            std::cerr << "WARNING: Failed to compute converged solution in " << it << " iterations." << std::endl;
-            return -1; // -1 indicates that the solver has to be restarted
-        }
-
-        linearIterations_ += linearIterations;
-        newtonIterations_ += it;
-        // Update max concentration.
-        computeCmax(x);
-        return linearIterations;
-    }
-
-
-
-
 
     template <class Grid>
     BlackoilPolymerModel<Grid>::ReservoirResidualQuant::ReservoirResidualQuant()
@@ -452,7 +425,7 @@ namespace detail {
     template <class Grid>
     typename BlackoilPolymerModel<Grid>::SolutionState
     BlackoilPolymerModel<Grid>::constantState(const PolymerBlackoilState& x,
-                                                         const WellStateFullyImplicitBlackoil&     xw) const
+                                              const WellStateFullyImplicitBlackoilPolymer& xw) const
     {
         auto state = variableState(x, xw);
         makeConstantState(state);
@@ -496,7 +469,7 @@ namespace detail {
     template <class Grid>
     typename BlackoilPolymerModel<Grid>::SolutionState
     BlackoilPolymerModel<Grid>::variableState(const PolymerBlackoilState& x,
-                                                         const WellStateFullyImplicitBlackoil&     xw) const
+                                              const WellStateFullyImplicitBlackoilPolymer& xw) const
     {
         using namespace Opm::AutoDiffGrid;
         const int nc = numCells(grid_);
@@ -738,7 +711,7 @@ namespace detail {
 
     template <class Grid>
     void BlackoilPolymerModel<Grid>::computeWellConnectionPressures(const SolutionState& state,
-                                                                               const WellStateFullyImplicitBlackoil& xw)
+                                                                    const WellStateFullyImplicitBlackoilPolymer& xw)
     {
         if( ! wellsActive() ) return ;
 
@@ -832,11 +805,9 @@ namespace detail {
     template <class Grid>
     void
     BlackoilPolymerModel<Grid>::
-    assemble(const V&             pvdt,
+    assemble(const PolymerBlackoilState& x,
-             const PolymerBlackoilState& x   ,
+             WellStateFullyImplicitBlackoilPolymer& xw,
-             const bool initial_assembly,
+             const bool initial_assembly)
-             WellStateFullyImplicitBlackoil& xw,
-             const std::vector<double>& polymer_inflow)
     {
         using namespace Opm::AutoDiffGrid;
 
@@ -883,7 +854,7 @@ namespace detail {
         for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
             computeMassFlux(phaseIdx, transi, kr[canph_[phaseIdx]], state.canonical_phase_pressures[canph_[phaseIdx]], state);
             residual_.material_balance_eq[ phaseIdx ] =
-                pvdt*(rq_[phaseIdx].accum[1] - rq_[phaseIdx].accum[0])
+                pvdt_ * (rq_[phaseIdx].accum[1] - rq_[phaseIdx].accum[0])
                 + ops_.div*rq_[phaseIdx].mflux;
         }
 
@@ -913,13 +884,13 @@ namespace detail {
 
         // Add polymer equation.
         if (has_polymer_) {
-            residual_.material_balance_eq[poly_pos_] = pvdt*(rq_[poly_pos_].accum[1] - rq_[poly_pos_].accum[0])
+            residual_.material_balance_eq[poly_pos_] = pvdt_ * (rq_[poly_pos_].accum[1] - rq_[poly_pos_].accum[0])
                                                + ops_.div*rq_[poly_pos_].mflux;
         }
 
         // Add contribution from wells and set up the well equations.
         V aliveWells;
-        addWellEq(state, xw, aliveWells, polymer_inflow);
+        addWellEq(state, xw, aliveWells);
         addWellControlEq(state, xw, aliveWells);
     }
 
@@ -929,9 +900,8 @@ namespace detail {
 
     template <class Grid>
     void BlackoilPolymerModel<Grid>::addWellEq(const SolutionState& state,
-                                                          WellStateFullyImplicitBlackoil& xw,
+                                               WellStateFullyImplicitBlackoilPolymer& xw,
-                                                          V& aliveWells,
+                                               V& aliveWells)
-                                                          const std::vector<double>& polymer_inflow)
     {
         if( ! wellsActive() ) return ;
 
@@ -1065,7 +1035,7 @@ namespace detail {
         // Add well contributions to polymer mass balance equation
         if (has_polymer_) {
             const ADB mc = computeMc(state);
-            const V polyin = Eigen::Map<const V>(&polymer_inflow[0], nc);
+            const V polyin = Eigen::Map<const V>(xw.polymerInflow().data(), nc);
             const V poly_in_perf = subset(polyin, well_cells);
             const V poly_mc_perf = subset(mc, well_cells).value();
             residual_.material_balance_eq[poly_pos_] -= superset(cq_ps[pu.phase_pos[Water]] * poly_mc_perf
@@ -1183,7 +1153,7 @@ namespace detail {
 
 
     template <class Grid>
-    void BlackoilPolymerModel<Grid>::updateWellControls(WellStateFullyImplicitBlackoil& xw) const
+    void BlackoilPolymerModel<Grid>::updateWellControls(WellStateFullyImplicitBlackoilPolymer& xw) const
     {
         if( ! wellsActive() ) return ;
 
@@ -1259,7 +1229,7 @@ namespace detail {
 
     template <class Grid>
     void BlackoilPolymerModel<Grid>::addWellControlEq(const SolutionState& state,
-                                                                 const WellStateFullyImplicitBlackoil& xw,
+                                                      const WellStateFullyImplicitBlackoilPolymer& xw,
                                                       const V& aliveWells)
     {
         if( ! wellsActive() ) return;
@@ -1360,7 +1330,7 @@ namespace detail {
     template <class Grid>
     void BlackoilPolymerModel<Grid>::updateState(const V& dx,
                                                  PolymerBlackoilState& state,
-                                                            WellStateFullyImplicitBlackoil& well_state)
+                                                 WellStateFullyImplicitBlackoilPolymer& well_state)
     {
         using namespace Opm::AutoDiffGrid;
         const int np = fluid_.numPhases();
@@ -1835,30 +1805,6 @@ namespace detail {
 
 
 
-    template <class Grid>
-    double
-    BlackoilPolymerModel<Grid>::residualNorm() const
-    {
-        double globalNorm = 0;
-        std::vector<ADB>::const_iterator quantityIt = residual_.material_balance_eq.begin();
-        const std::vector<ADB>::const_iterator endQuantityIt = residual_.material_balance_eq.end();
-        for (; quantityIt != endQuantityIt; ++quantityIt) {
-            const double quantityResid = (*quantityIt).value().matrix().norm();
-            if (!std::isfinite(quantityResid)) {
-                const int trouble_phase = quantityIt - residual_.material_balance_eq.begin();
-                OPM_THROW(Opm::NumericalProblem,
-                          "Encountered a non-finite residual in material balance equation "
-                          << trouble_phase);
-            }
-            globalNorm = std::max(globalNorm, quantityResid);
-        }
-        globalNorm = std::max(globalNorm, residual_.well_flux_eq.value().matrix().norm());
-        globalNorm = std::max(globalNorm, residual_.well_eq.value().matrix().norm());
-
-        return globalNorm;
-    }
-
-
     template <class Grid>
     std::vector<double>
     BlackoilPolymerModel<Grid>::computeResidualNorms() const
@@ -1895,73 +1841,8 @@ namespace detail {
         return residualNorms;
     }
 
-    template <class Grid>
-    void
-    BlackoilPolymerModel<Grid>::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 < fluid_.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 Grid>
-    void
-    BlackoilPolymerModel<Grid>::stablizeNewton(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;
-    }
 
     template <class Grid>
     double

opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp

@@ -22,6 +22,7 @@
 #include <opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp>
 #include <opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp>
 #include <opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp>
+#include <opm/polymer/fullyimplicit/WellStateFullyImplicitBlackoilPolymer.hpp>
 #include <opm/polymer/PolymerBlackoilState.hpp>
 #include <opm/polymer/PolymerInflow.hpp>
 
@@ -30,7 +31,6 @@
 
 #include <opm/autodiff/GeoProps.hpp>
 #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
-#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 #include <opm/autodiff/RateConverter.hpp>
 #include <opm/autodiff/NewtonSolver.hpp>
 
@@ -135,7 +135,7 @@ namespace Opm
         computeRESV(const std::size_t               step,
                     const Wells*                    wells,
                     const BlackoilState&     x,
-                    WellStateFullyImplicitBlackoil& xw);
+                    WellStateFullyImplicitBlackoilPolymer& xw);
     };
 
 
@@ -236,7 +236,7 @@ namespace Opm
     SimulatorReport SimulatorFullyImplicitBlackoilPolymer<T>::Impl::run(SimulatorTimer& timer,
                                                                         PolymerBlackoilState& state)
     {
-        WellStateFullyImplicitBlackoil prev_well_state;
+        WellStateFullyImplicitBlackoilPolymer prev_well_state;
 
         // Create timers and file for writing timing info.
         Opm::time::StopWatch solver_timer;
@@ -249,11 +249,10 @@ namespace Opm
 
         typedef T Grid;
         typedef BlackoilPolymerModel<Grid> Model;
-        // typedef typename Model::ModelParameter ModelParam;
+        typedef typename Model::ModelParameter ModelParam;
-        // ModelParam modelParam( param_ );
+        ModelParam modelParam( param_ );
-        // typedef NewtonSolver<Model> Solver;
+        typedef NewtonSolver<Model> Solver;
-        // typedef typename Solver::SolverParameter SolverParam;
+        typedef typename Solver::SolverParameter SolverParam;
-        typedef typename Model::SolverParameter SolverParam;
         SolverParam solverParam( param_ );
 
         //adaptive time stepping
@@ -297,7 +296,7 @@ namespace Opm
                                        Opm::UgGridHelpers::beginFaceCentroids(grid_),
                                        props_.permeability());
             const Wells* wells = wells_manager.c_wells();
-            WellStateFullyImplicitBlackoil well_state;
+            WellStateFullyImplicitBlackoilPolymer well_state;
             well_state.init(wells, state.blackoilState(), prev_well_state);
 
             // compute polymer inflow
@@ -316,6 +315,7 @@ namespace Opm
             polymer_inflow_ptr->getInflowValues(timer.simulationTimeElapsed(), 
                                                 timer.simulationTimeElapsed() + timer.currentStepLength(),
                                                 polymer_inflow_c);
+            well_state.polymerInflow() = polymer_inflow_c;
 
             // write simulation state at the report stage
             output_writer_.writeTimeStep( timer, state.blackoilState(), well_state );
@@ -330,10 +330,11 @@ namespace Opm
             // Run a multiple steps of the solver depending on the time step control.
             solver_timer.start();
 
-            Model solver(solverParam, grid_, props_, geo_, rock_comp_props_, polymer_props_, wells, solver_, has_disgas_, has_vapoil_, has_polymer_, terminal_output_);
+            Model model(modelParam, grid_, props_, geo_, rock_comp_props_, polymer_props_, wells, solver_, has_disgas_, has_vapoil_, has_polymer_, terminal_output_);
             if (!threshold_pressures_by_face_.empty()) {
-                solver.setThresholdPressures(threshold_pressures_by_face_);
+                model.setThresholdPressures(threshold_pressures_by_face_);
             }
+            Solver solver(solverParam, 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
@@ -344,7 +345,7 @@ namespace Opm
             //                adaptiveTimeStepping->step( timer, solver, state, well_state,  output_writer_ );
             //            } else {
                 // solve for complete report step
-            solver.step(timer.currentStepLength(), state, well_state, polymer_inflow_c);
+            solver.step(timer.currentStepLength(), state, well_state);
                 //            }
 
             // take time that was used to solve system for this reportStep
@@ -511,7 +512,7 @@ namespace Opm
     Impl::computeRESV(const std::size_t               step,
                       const Wells*                    wells,
                       const BlackoilState&            x,
-                      WellStateFullyImplicitBlackoil& xw)
+                      WellStateFullyImplicitBlackoilPolymer& xw)
     {
         typedef SimFIBODetails::WellMap WellMap;
 

opm/polymer/fullyimplicit/WellStateFullyImplicitBlackoilPolymer.hpp

@@ -0,0 +1,39 @@
+/*
+  Copyright 2015 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_WELLSTATEFULLYIMPLICITBLACKOILPOLYMER_HEADER_INCLUDED
+#define OPM_WELLSTATEFULLYIMPLICITBLACKOILPOLYMER_HEADER_INCLUDED
+
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
+
+namespace Opm
+{
+
+    class WellStateFullyImplicitBlackoilPolymer : public WellStateFullyImplicitBlackoil
+    {
+    public:
+        std::vector<double>& polymerInflow() { return polymer_inflow_; }
+        const std::vector<double>& polymerInflow() const { return polymer_inflow_; }
+    private:
+        std::vector<double> polymer_inflow_;
+    };
+
+} // namespace Opm
+
+#endif // OPM_WELLSTATEFULLYIMPLICITBLACKOILPOLYMER_HEADER_INCLUDED