some revision, time step control is now completly in the Simulator run method.

The solver simply returns a number of iterations.

opm/autodiff/FullyImplicitBlackoilSolver.hpp

@@ -25,7 +25,6 @@
 #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
 #include <opm/autodiff/LinearisedBlackoilResidual.hpp>
 #include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
-#include <opm/autodiff/TimeStepControl.hpp>
 
 struct UnstructuredGrid;
 struct Wells;
@@ -115,8 +114,8 @@ namespace Opm {
         /// \param[in] dt        time step size
         /// \param[in] state     reservoir state
         /// \param[in] wstate    well state
-        /// \return suggested time step for next step call
+        /// \return              number of iterations used
-        double
+        int
         step(const double   dt    ,
              BlackoilState& state ,
              WellStateFullyImplicitBlackoil&     wstate);
@@ -191,9 +190,6 @@ namespace Opm {
 
         std::vector<int>         primalVariable_;
 
-        //IterationCountTimeStepControl timeStepControl_;
-        PIDTimeStepControl timeStepControl_;
-
         // Private methods.
         SolutionState
         constantState(const BlackoilState& x,

opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp

@@ -235,7 +235,6 @@ namespace {
         , residual_ ( { std::vector<ADB>(fluid.numPhases(), ADB::null()),
                         ADB::null(),
                         ADB::null() } )
-        , timeStepControl_()
     {
     }
 
@@ -263,7 +262,7 @@ namespace {
 
 
     template<class T>
-    double
+    int
     FullyImplicitBlackoilSolver<T>::
     step(const double   dt,
          BlackoilState& x ,
@@ -279,9 +278,6 @@ namespace {
             computeWellConnectionPressures(state, xw);
         }
 
-        // initialize time step control 
-        timeStepControl_.initialize( x );
-
         std::vector<std::vector<double>> residual_history;
 
         assemble(pvdt, x, xw);
@@ -335,18 +331,19 @@ namespace {
 
             converged = getConvergence(dt);
 
-            it += 1;
+            // increase iteration counter
+            ++it;
             std::cout << std::setw(9) << it << std::setprecision(9)
                       << std::setw(18) << r << std::endl;
         }
 
         if (!converged) {
-            std::cerr << "Failed to compute converged solution in " << it << " iterations. Ignoring!\n";
+            std::cerr << "ERROR: Failed to compute converged solution in " << it << " iterations." << std::endl;
             // OPM_THROW(std::runtime_error, "Failed to compute converged solution in " << it << " iterations.");
+            return -1;
         }
 
-        std::cout << "Linear iterations: " << linearIterations << std::endl;
+        return linearIterations;
-        return timeStepControl_.computeTimeStepSize( dt, linearIterations, x );
     }
 
 
@@ -1704,6 +1701,7 @@ namespace {
         for (; quantityIt != endQuantityIt; ++quantityIt) {
             const double quantityResid = (*quantityIt).value().matrix().norm();
             if (!std::isfinite(quantityResid)) {
+                //std::cout << quantityResid << " quantity" << std::endl;
                 const int trouble_phase = quantityIt - residual_.material_balance_eq.begin();
                 OPM_THROW(Opm::NumericalProblem,
                           "Encountered a non-finite residual in material balance equation "

opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp

@@ -28,6 +28,7 @@
 #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 #include <opm/autodiff/RateConverter.hpp>
+#include <opm/autodiff/TimeStepControl.hpp>
 
 #include <opm/core/grid.h>
 #include <opm/core/wells.h>
@@ -37,6 +38,7 @@
 #include <opm/core/io/eclipse/EclipseWriter.hpp>
 #include <opm/core/simulator/SimulatorReport.hpp>
 #include <opm/core/simulator/SimulatorTimer.hpp>
+#include <opm/core/simulator/AdaptiveSimulatorTimer.hpp>
 #include <opm/core/utility/StopWatch.hpp>
 #include <opm/core/io/vtk/writeVtkData.hpp>
 #include <opm/core/utility/miscUtilities.hpp>
@@ -298,6 +300,9 @@ namespace Opm
 
         const bool subStepping = param_.getDefault("substepping", bool(false) );
 
+        std::unique_ptr< TimeStepControlInterface > 
+            timeStepControl( new PIDAndIterationCountTimeStepControl( 1e-3, 50 ) );
+
         // Main simulation loop.
         while (!timer.done()) {
             // Report timestep.
@@ -353,31 +358,67 @@ namespace Opm
                 solver.setThresholdPressures(threshold_pressures_by_face_);
             }
 
+            // If sub stepping is enabled allow the solver to sub cycle
+            // in case the report steps are to large for the solver to converge
+            // \Note: The report steps are met in any case
             if( subStepping )
             {
                 // create sub step simulator timer with previously used sub step size
-                SubStepSimulatorTimer subStepper( timer, lastSubStep );
+                const double start_time = timer.simulationTimeElapsed();
+                const double end_time   = start_time + timer.currentStepLength();
+                AdaptiveSimulatorTimer subStepper( start_time, end_time, lastSubStep );
 
+                // copy states in case solver has to be restarted
+                //BlackoilState last_state( state );
+                //WellStateFullyImplicitBlackoil last_well_state( well_state );
+
+                // sub step time loop
                 while( ! subStepper.done() )
                 {
-                    const double dt_new = solver.step(subStepper.currentStepLength(), state, well_state);
+                    // initialize time step control in case current state is needed later
-                    subStepper.next( dt_new );
+                    timeStepControl->initialize( state );
+
+                    // keep last state for solver restart and time step control
+                    // (linearIterations < 0 means on convergence in solver)
+                    const int linearIterations = solver.step(subStepper.currentStepLength(), state, well_state);
+
+                    // (linearIterations < 0 means on convergence in solver)
+                    if( linearIterations >= 0 )
+                    {
+                        // advance by current dt
+                        subStepper.advance();
+
+                        // compute new time step estimate
+                        const double dtEstimate = 
+                            timeStepControl->computeTimeStepSize( subStepper.currentStepLength(), linearIterations, state );
+                        // set new time step length
+                        subStepper.provideTimeStepEstimate( dtEstimate );
+                    }
+                    else // in case of no convergence
+                    {
+                        // we need to revise this
+                        abort();
+                        subStepper.halfTimeStep();
+                        std::cerr << "Solver convergence failed, restarting solver with half time step ("<< subStepper.currentStepLength()<<" days)." << std::endl;
+                        // reset states 
+                        // state      = last_state;
+                        // well_state = last_well_state;
+                    }
                 }
 
                 subStepper.report( std::cout );
 
                 // store last small time step for next reportStep
-                //lastSubStep = subStepper.maxStepLength();
-                //lastSubStep = subStepper.averageStepLength();
-                //lastSubStep = subStepper.suggestedMax();
                 lastSubStep = subStepper.suggestedAverage();
+                std::cout << "Last suggested step size = " << lastSubStep/86400.0 << " (days)" << std::endl;
 
-                std::cout << "Last suggested step size = " << lastSubStep << std::endl;
+                if( ! std::isfinite( lastSubStep ) ) // check for NaN
-                if( lastSubStep != lastSubStep )
                     lastSubStep = timer.currentStepLength();
             }
-            else
+            else {
+                // solve for complete report step
                 solver.step(timer.currentStepLength(), state, well_state);
+            }
 
             // take time that was used to solve system for this reportStep
             solver_timer.stop();

opm/autodiff/TimeStepControl.hpp

@@ -22,117 +22,130 @@
 namespace Opm
 {
 
+    ///////////////////////////////////////////////////////////////////
+    ///
+    ///  TimeStepControlInterface 
+    /// 
+    ///////////////////////////////////////////////////////////////////
     class TimeStepControlInterface 
     {
     protected:    
         TimeStepControlInterface() {}
     public:
+        /// \param state simulation state before computing update in the solver (default is empty)
         virtual void initialize( const SimulatorState& state ) {}
+
+        /// compute new time step size suggestions based on the PID controller
+        /// \param dt          time step size used in the current step
+        /// \param iterations  number of iterations used (linear/nonlinear) 
+        /// \param state       new solution state
+        ///
+        /// \return suggested time step size for the next step
         virtual double computeTimeStepSize( const double dt, const int iterations, const SimulatorState& ) const = 0;
+
+        /// virtual destructor (empty)
         virtual ~TimeStepControlInterface () {}
     };
 
-    class IterationCountTimeStepControl : public TimeStepControlInterface
+    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
-    {
+    ///
-    protected:    
+    ///  PID controller based adaptive time step control as suggested in: 
-        mutable double prevDt_;
+    ///  Turek and Kuzmin. Algebraic Flux Correction III. Incompressible Flow Problems. Uni Dortmund. 
-        mutable int prevIterations_;
+    /// 
-        const int targetIterationCount_;
+    ///  See also: 
-        const double adjustmentFactor_;
+    ///  D. Kuzmin and S.Turek. Numerical simulation of turbulent bubbly flows. Techreport Uni Dortmund. 2004
-
+    ///  
-        const int upperTargetIterationCount_;
+    ///  and the original article: 
-        const int lowerTargetIterationCount_;
+    ///  Valli, Coutinho, and Carey. Adaptive Control for Time Step Selection in Finite Element
-
+    ///  Simulation of Coupled Viscous Flow and Heat Transfer. Proc of the 10th 
-    public:
+    ///  International Conference on Numerical Methods in Fluids. 1998.
-        IterationCountTimeStepControl() 
+    ///
-            : prevDt_( 0.0 ), prevIterations_( 0 ), 
+    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
-              targetIterationCount_( 100 ), adjustmentFactor_( 1.25 ),
-              upperTargetIterationCount_( 200 ), lowerTargetIterationCount_( 30 )
-        {}
-
-        double computeTimeStepSize( const double dt, const int iterations, const SimulatorState& ) 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;
-        }
-    };
-
     class PIDTimeStepControl : public TimeStepControlInterface
     {
     protected:   
         mutable std::vector<double> p0_;
         mutable std::vector<double> sat0_;
 
-        mutable double prevDt_;
-        mutable int prevIterations_;
-        const int targetIterationCount_;
-        const double adjustmentFactor_;
-
-        const int upperTargetIterationCount_;
-        const int lowerTargetIterationCount_;
-
         const double tol_;
-
         mutable std::vector< double > errors_;
 
+        const bool verbose_;
     public:
-        PIDTimeStepControl( const double tol = 8e-4 ) 
+        /// \brief constructor 
-            : p0_(), sat0_(), prevDt_( 0.0 ), prevIterations_( 0 ), 
+        /// \param tol  tolerance for the relative changes of the numerical solution to be accepted 
-              targetIterationCount_( 100 ), adjustmentFactor_( 1.25 ),
+        ///             in one time step (default is 1e-3)
-              upperTargetIterationCount_( 200 ), lowerTargetIterationCount_( 30 ),
+        PIDTimeStepControl( const double tol = 1e-3, const bool verbose = false ) 
-              tol_( tol ),
+            : p0_()
-              errors_( 3, tol_ )
+            , sat0_() 
+            , tol_( tol )
+            , errors_( 3, tol_ )
+            , verbose_( verbose )
         {}
 
+        /// \brief \copydoc TimeStepControlInterface::initialize
         void initialize( const SimulatorState& state ) 
         {
-            // store current state
+            // store current state for later time step computation
             p0_   = state.pressure();
             sat0_ = state.saturation();
         }
 
+        /// \brief \copydoc TimeStepControlInterface::computeTimeStepSize
+        double computeTimeStepSize( const double dt, const int /* iterations */, const SimulatorState& state ) const
+        {
+            const size_t size = p0_.size();
+            assert( state.pressure().size() == size );
+            assert( state.saturation().size() == size );
+            assert( sat0_.size() == size );
+
+            // compute u^n - u^n+1 
+            for( size_t i=0; i<size; ++i ) 
+            {
+                p0_[ i ]   -= state.pressure()[ i ];
+                sat0_[ i ] -= state.saturation()[ i ];
+            }
+
+            // compute || u^n - u^n+1 || 
+            const double stateOld  = inner_product( p0_.begin(),   p0_.end() ) +
+                                     inner_product( sat0_.begin(), sat0_.end() );
+
+            // compute || u^n+1 || 
+            const double stateNew  = inner_product( state.pressure().begin(), state.pressure().end() ) +
+                                     inner_product( state.saturation().begin(), state.saturation().end() );
+
+            // shift errors
+            for( int i=0; i<2; ++i )
+              errors_[ i ] = errors_[i+1];
+
+            // store new error
+            const double error = stateOld / stateNew;
+            errors_[ 2 ] =  error ;
+
+            if( error > tol_ )
+            {
+                // adjust dt by given tolerance 
+                if( verbose_ )
+                    std::cout << "Computed step size (tol): " << (dt * tol_ / error )/86400.0 << " (days)" << std::endl;
+                return (dt * tol_ / error );
+            }
+            else
+            {
+                // values taking from turek time stepping paper
+                const double kP = 0.075 ;
+                const double kI = 0.175 ;
+                const double kD = 0.01 ;
+                double newDt = (dt * std::pow( errors_[ 1 ] / errors_[ 2 ], kP ) *
+                             std::pow( tol_         / errors_[ 2 ], kI ) *
+                             std::pow( (errors_[0]*errors_[0]/errors_[ 1 ]*errors_[ 2 ]), kD ));
+                if( verbose_ )
+                    std::cout << "Computed step size (pow): " << newDt/86400.0 << " (days)" << std::endl;
+                return newDt;
+            }
+        }
+
+    protected:    
+        // return inner product for given container, here std::vector
         template <class Iterator>
         double inner_product( Iterator it, const Iterator end ) const 
         {
@@ -142,55 +155,37 @@ namespace Opm
             return product;
         }
 
+    };
+
+    class PIDAndIterationCountTimeStepControl : public PIDTimeStepControl
+    {
+        typedef PIDTimeStepControl BaseType;
+    protected:   
+        const int targetIterationCount_;
+
+    public:
+        PIDAndIterationCountTimeStepControl( const int target_iterations = 20,
+                                             const double tol = 1e-3, 
+                                             const bool verbose = false) 
+            : BaseType( tol, verbose )
+            , targetIterationCount_( target_iterations )
+        {}
+
         double computeTimeStepSize( const double dt, const int iterations, const SimulatorState& state ) const
         {
-            const size_t size = p0_.size();
+            double dtEstimate = BaseType :: computeTimeStepSize( dt, iterations, state );
-            assert( state.pressure().size() == size );
+
-            // compute u^n - u^n+1 
+            // further reduce step size if to many iterations were used
-            for( size_t i=0; i<size; ++i ) 
+            if( iterations > targetIterationCount_ ) 
             {
-                p0_[ i ]   -= state.pressure()[ i ];
+                // if iterations was the same or dts were the same, do some magic
-                sat0_[ i ] -= state.saturation()[ i ];
+                dtEstimate *= double( targetIterationCount_ ) / double(iterations);
             }
 
-            // compute || u^n - u^n+1 || 
+            return dtEstimate;
-            double stateN0  = inner_product( p0_.begin(),   p0_.end() ) +
-                              inner_product( sat0_.begin(), sat0_.end() );
-
-            // compute || u^n+1 || 
-            double stateN   = inner_product( state.pressure().begin(), state.pressure().end() ) +
-                              inner_product( state.saturation().begin(), state.saturation().end() );
-
-
-            for( int i=0; i<2; ++i )
-              errors_[ i ] = errors_[i+1];
-
-            double error = stateN0 / stateN ;
-            errors_[ 2 ] =  error ;
-
-            prevDt_ = dt;
-            prevIterations_ = iterations;
-
-            if( error > tol_ )
-            {
-                // adjust dt by given tolerance 
-                std::cout << "Computed step size (tol): " << (dt * tol_ / error ) << std::endl;
-                return (dt * tol_ / error );
-            }
-            else
-            {
-                const double kP = 0.075 ;
-                const double kI = 0.175 ;
-                const double kD = 0.01 ;
-                double newDt = (dt * std::pow( errors_[ 1 ] / errors_[ 2 ], kP ) *
-                             std::pow( tol_         / errors_[ 2 ], kI ) *
-                             std::pow( (errors_[0]*errors_[0]/errors_[ 1 ]*errors_[ 2 ]), kD ));
-                std::cout << "Computed step size (pow): " << newDt << std::endl;
-                return newDt;
-            }
         }
     };
 
-} // end namespace OPM
 
+} // end namespace OPM
 #endif