First version of time step control. Not really fully functioning yet.

opm/autodiff/FullyImplicitBlackoilSolver.hpp

@@ -25,6 +25,7 @@
 #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
 #include <opm/autodiff/LinearisedBlackoilResidual.hpp>
 #include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
+#include <opm/autodiff/TimeStepControl.hpp>
 
 struct UnstructuredGrid;
 struct Wells;
@@ -93,7 +94,8 @@ namespace Opm {
         /// \param[in] dt        time step size
         /// \param[in] state     reservoir state
         /// \param[in] wstate    well state
-        void
+        /// \return new suggested time step
+        double 
         step(const double   dt    ,
              BlackoilState& state ,
              WellStateFullyImplicitBlackoil&     wstate);
@@ -176,6 +178,8 @@ namespace Opm {
 
         std::vector<int>         primalVariable_;
 
+        IterationCountTimeStepControl timeStepControl_;
+
         // Private methods.
         SolutionState
         constantState(const BlackoilState& x,

opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp

@@ -196,12 +196,13 @@ namespace {
         , residual_ ( { std::vector<ADB>(fluid.numPhases(), ADB::null()),
                         ADB::null(),
                         ADB::null() } )
+        , timeStepControl_() 
     {
-        dp_max_rel_ = param.getDefault("dp_max_rel", dp_max_rel_);
-        ds_max_ = param.getDefault("ds_max", ds_max_);
+        dp_max_rel_  = param.getDefault("dp_max_rel", dp_max_rel_);
+        ds_max_      = param.getDefault("ds_max", ds_max_);
         drs_max_rel_ = param.getDefault("drs_max_rel", drs_max_rel_);
-        relax_max_ = param.getDefault("relax_max", relax_max_);
-        max_iter_ = param.getDefault("max_iter", max_iter_);
+        relax_max_   = param.getDefault("relax_max", relax_max_);
+        max_iter_    = param.getDefault("max_iter", max_iter_);
 
         std::string relaxation_type = param.getDefault("relax_type", std::string("dampen"));
         if (relaxation_type == "dampen") {
@@ -234,10 +235,8 @@ namespace {
     }
 
 
-
-
     template<class T>
-    void
+    double
     FullyImplicitBlackoilSolver<T>::
     step(const double   dt,
          BlackoilState& x ,
@@ -279,10 +278,14 @@ namespace {
         bool isOscillate = false;
         bool isStagnate = false;
         const enum RelaxType relaxtype = relaxType();
+        int linearIterations = 0 ;
 
         while ((!converged) && (it < maxIter())) {
             V dx = solveJacobianSystem();
 
+            // store number of linear iterations used
+            linearIterations += linsolver_.iterationCount();
+
             detectNewtonOscillations(residual_history, it, relaxRelTol(), isOscillate, isStagnate);
 
             if (isOscillate) {
@@ -312,6 +315,9 @@ namespace {
             std::cerr << "Failed to compute converged solution in " << it << " iterations. Ignoring!\n";
             // OPM_THROW(std::runtime_error, "Failed to compute converged solution in " << it << " iterations.");
         }
+
+        std::cout << "Iterations count " << linearIterations << std::endl;
+        return timeStepControl_.computeTimeStepSize( dt, linearIterations );
     }
 
 

opm/autodiff/NewtonIterationBlackoilCPR.cpp

@@ -108,6 +108,7 @@ namespace Opm
 
     /// Construct a system solver.
     NewtonIterationBlackoilCPR::NewtonIterationBlackoilCPR(const parameter::ParameterGroup& param)
+        : iteration_count_( 0 )
     {
         use_amg_ = param.getDefault("cpr_use_amg", false);
         use_bicgstab_ = param.getDefault("cpr_use_bicgstab", true);
@@ -191,7 +192,7 @@ namespace Opm
         // Construct linear solver.
         const double tolerance = 1e-3;
         const int maxit = 5000;
-        const int verbosity = 1;
+        const int verbosity = 0;
         const int restart = 40;
         Dune::RestartedGMResSolver<Vector> linsolve(opA, sp, precond, tolerance, restart, maxit, verbosity);
 
@@ -199,6 +200,9 @@ namespace Opm
         Dune::InverseOperatorResult result;
         linsolve.apply(x, istlb, result);
 
+        // store number of iterations used
+        iteration_count_ = result.iterations;
+
         // Check for failure of linear solver.
         if (!result.converged) {
             OPM_THROW(std::runtime_error, "Convergence failure for linear solver.");

opm/autodiff/NewtonIterationBlackoilCPR.hpp

@@ -54,7 +54,9 @@ namespace Opm
         /// \return               the solution x
         virtual SolutionVector computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const;
 
+        virtual int iterationCount () const { return iteration_count_; }
     private:
+        mutable int iteration_count_;
         bool use_amg_;
         bool use_bicgstab_;
     };

opm/autodiff/NewtonIterationBlackoilInterface.hpp

@@ -38,6 +38,8 @@ namespace Opm
         /// \param[in] residual   residual object containing A and b.
         /// \return               the solution x
         virtual SolutionVector computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const = 0;
+
+        virtual int iterationCount () const { return 0; }
     };
 
 } // namespace Opm

opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp

@@ -291,6 +291,8 @@ namespace Opm
         std::string tstep_filename = output_dir_ + "/step_timing.txt";
         std::ofstream tstep_os(tstep_filename.c_str());
 
+        double lastSubStep = timer.currentStepLength();
+
         // Main simulation loop.
         while (!timer.done()) {
             // Report timestep.
@@ -338,13 +340,40 @@ namespace Opm
             // Compute reservoir volumes for RESV controls.
             computeRESV(timer.currentStepNum(), wells, state, well_state);
 
-            // Run a single step of the solver.
+            // Run a multiple steps of the solver depending on the time step control.
             solver_timer.start();
+
             FullyImplicitBlackoilSolver<T> solver(param_, grid_, props_, geo_, rock_comp_props_, *wells, solver_, has_disgas_, has_vapoil_);
             if (!threshold_pressures_by_face_.empty()) {
                 solver.setThresholdPressures(threshold_pressures_by_face_);
             }
-            solver.step(timer.currentStepLength(), state, well_state);
+
+            const bool subStepping = true;
+            if( subStepping )
+            {
+
+                // create sub step simulator timer with previously used sub step size
+                SubStepSimulatorTimer subStepper( timer, lastSubStep );
+
+                while( ! subStepper.done() )
+                {
+                    const double dt_new = solver.step(subStepper.currentStepLength(), state, well_state);
+                    subStepper.next( dt_new );
+                }
+
+                subStepper.report( std::cout );
+
+                // store last small time step for next reportStep
+                lastSubStep = subStepper.currentStepLength();
+
+                std::cout << "Last suggested step size = " << lastSubStep << std::endl;
+                if( lastSubStep != lastSubStep )
+                    lastSubStep = timer.currentStepLength();
+            }
+            else
+                solver.step(timer.currentStepLength(), state, well_state);
+
+            // take time that was used to solve system for this reportStep
             solver_timer.stop();
 
             // Report timing.

opm/autodiff/TimeStepControl.hpp

@@ -0,0 +1,103 @@
+/*
+  Copyright 2014 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_TIMESTEPCONTROL_HEADER_INCLUDED
+#define OPM_TIMESTEPCONTROL_HEADER_INCLUDED
+
+namespace Opm
+{
+
+    class TimeStepControlInterface 
+    {
+    protected:    
+        TimeStepControlInterface() {}
+    public:
+        virtual double computeTimeStepSize( const double dt, const int iterations ) const = 0;
+        virtual ~TimeStepControlInterface () {}
+    };
+
+    class IterationCountTimeStepControl : public TimeStepControlInterface
+    {
+    protected:    
+        mutable double prevDt_;
+        mutable int prevIterations_;
+        const int targetIterationCount_;
+        const double adjustmentFactor_;
+
+        const int upperTargetIterationCount_;
+        const int lowerTargetIterationCount_;
+
+    public:
+        IterationCountTimeStepControl() 
+            : prevDt_( 0.0 ), prevIterations_( 0 ), 
+              targetIterationCount_( 100 ), adjustmentFactor_( 1.25 ),
+              upperTargetIterationCount_( 200 ), lowerTargetIterationCount_( 30 )
+        {}
+
+        double computeTimeStepSize( const double dt, const int iterations ) const
+        {
+            // make sure dt is somewhat reliable
+            assert( dt > 0 && dt == dt );
+            double newDt = dt;
+            double derivation = double(std::abs( iterations - targetIterationCount_ )) / double(targetIterationCount_);
+
+            if( derivation > 0.1 )
+            {
+                if( iterations < targetIterationCount_ ) 
+                {
+                    newDt = dt * adjustmentFactor_;
+                }
+                else 
+                    newDt = dt / adjustmentFactor_;
+            }
+                           
+            /*
+            if( prevDt_ > 0 && std::abs( dt - prevDt_ ) > 1e-12 ) {
+                const double dFdt  = double(iterations - prevIterations_) / ( dt - prevDt_ );
+                if( std::abs( dFdt ) > 1e-12 )
+                    newDt = dt + (targetIterationCount_ - iterations) / dFdt;
+                else
+                    // if iterations was the same or dts were the same, do some magic
+                    newDt = dt * double( targetIterationCount_ ) / double(targetIterationCount_ - iterations);
+            }
+
+            if( newDt < 0 || ! (prevDt_ > 0) || ( iterations == prevIterations_) ) 
+            {
+                if( iterations > upperTargetIterationCount_ )
+                    newDt = dt / adjustmentFactor_;
+                else if( iterations < lowerTargetIterationCount_ )
+                    newDt = dt * adjustmentFactor_;
+                else
+                    newDt = dt;
+            }
+            */
+
+            assert( newDt == newDt );
+
+            //std::cout << "dt = " << dt << "  " << prevDt_ << std::endl;
+
+            prevDt_ = dt;
+            prevIterations_ = iterations;
+
+            return newDt;
+        }
+    };
+
+} // end namespace OPM
+
+#endif