AdaptiveSimulatorTimer: -improvement in time step adjustment near end of time interval

-max time step parameter

PIDTimeStepControl --> TimeStepControl:
- added simple iteration count time step control
- bug fix in PIDAndIterationCountTimeStepControl

AdaptiveTimeStepping: apply the above changes.

CMakeLists_files.cmake

@@ -103,7 +103,7 @@ list (APPEND MAIN_SOURCE_FILES
 	opm/core/props/satfunc/SaturationPropsFromDeck.cpp
 	opm/core/simulator/AdaptiveSimulatorTimer.cpp
 	opm/core/simulator/BlackoilState.cpp
-	opm/core/simulator/PIDTimeStepControl.cpp
+	opm/core/simulator/TimeStepControl.cpp
 	opm/core/simulator/SimulatorCompressibleTwophase.cpp
 	opm/core/simulator/SimulatorIncompTwophase.cpp
 	opm/core/simulator/SimulatorOutput.cpp
@@ -368,7 +368,7 @@ list (APPEND PUBLIC_HEADER_FILES
 	opm/core/simulator/AdaptiveTimeStepping_impl.hpp
 	opm/core/simulator/BlackoilState.hpp
 	opm/core/simulator/EquilibrationHelpers.hpp
-	opm/core/simulator/PIDTimeStepControl.hpp
+	opm/core/simulator/TimeStepControl.hpp
 	opm/core/simulator/SimulatorCompressibleTwophase.hpp
 	opm/core/simulator/SimulatorIncompTwophase.hpp
 	opm/core/simulator/SimulatorOutput.hpp

opm/core/simulator/AdaptiveSimulatorTimer.cpp

@@ -30,17 +30,18 @@
 namespace Opm
 {
     AdaptiveSimulatorTimer::
-    AdaptiveSimulatorTimer( const SimulatorTimerInterface& timer, const double lastStepTaken )
+    AdaptiveSimulatorTimer( const SimulatorTimerInterface& timer,
+                            const double lastStepTaken,
+                            const double maxTimeStep )
         : start_date_time_( timer.startDateTime() )
         , start_time_( timer.simulationTimeElapsed() )
         , total_time_( start_time_ + timer.currentStepLength() )
         , report_step_( timer.reportStepNum() )
+        , max_time_step_( maxTimeStep )
         , current_time_( start_time_ )
         , dt_( 0.0 )
         , current_step_( 0 )
         , steps_()
-        , suggestedMax_( 0.0 )
-        , suggestedAverage_( 0.0 )
     {
         // reserve memory for sub steps
         steps_.reserve( 10 );
@@ -61,31 +62,30 @@ namespace Opm
     void AdaptiveSimulatorTimer::
     provideTimeStepEstimate( const double dt_estimate )
     {
-        // store some information about the time steps suggested
-        suggestedMax_      = std::max( dt_estimate, suggestedMax_ );
-        suggestedAverage_ += dt_estimate;
-
         double remaining = (total_time_ - current_time_);
+        // apply max time step if it was set
+        dt_ = std::min( dt_estimate, max_time_step_ );
 
         if( remaining > 0 ) {
 
             // set new time step (depending on remaining time)
-            if( 1.5 * dt_estimate > remaining ) {
+            if( 1.05 * dt_ > remaining ) {
                 dt_ = remaining;
-                return ;
+                // check max time step again and use half remaining if to large
+                if( dt_ > max_time_step_ ) {
+                    dt_ = 0.5 * remaining;
+                }
+                return;
             }
 
             // check for half interval step to avoid very small step at the end
             // remaining *= 0.5;
 
-            if( 2.25 * dt_estimate > remaining ) {
+            if( 1.5 * dt_ > remaining ) {
                 dt_ = 0.5 * remaining;
-                return ;
+                return;
             }
         }
-
-        // otherwise set dt_estimate as is
-        dt_ = dt_estimate;
     }
 
     int AdaptiveSimulatorTimer::
@@ -137,16 +137,6 @@ namespace Opm
         return *(std::min_element( steps_.begin(), steps_.end() ));
     }
 
-    /// \brief return max suggested step length
-    double AdaptiveSimulatorTimer::suggestedMax () const { return suggestedMax_; }
-
-    /// \brief return average suggested step length
-    double AdaptiveSimulatorTimer::suggestedAverage () const
-    {
-        const int size = steps_.size();
-        return (size > 0 ) ? (suggestedAverage_ / double(size)) : suggestedAverage_;
-    }
-
     /// \brief report start and end time as well as used steps so far
     void AdaptiveSimulatorTimer::
     report(std::ostream& os) const

opm/core/simulator/AdaptiveSimulatorTimer.hpp

@@ -40,8 +40,12 @@ namespace Opm
     {
     public:
         /// \brief constructor taking a simulator timer to determine start and end time
-        ///  \param timer   in case of sub stepping this is the outer timer
+        ///  \param timer          in case of sub stepping this is the outer timer
-        explicit AdaptiveSimulatorTimer( const SimulatorTimerInterface& timer, const double lastStepTaken );
+        ///  \param lastStepTaken  last suggested time step
+        ///  \param maxTimeStep    maximum time step allowed
+        AdaptiveSimulatorTimer( const SimulatorTimerInterface& timer,
+                                const double lastStepTaken,
+                                const double maxTimeStep = std::numeric_limits<double>::max() );
 
         /// \brief advance time by currentStepLength
         AdaptiveSimulatorTimer& operator++ ();
@@ -79,12 +83,6 @@ namespace Opm
         /// \brief return min step length used so far
         double minStepLength () const;
 
-        /// \brief return max suggested step length
-        double suggestedMax () const;
-
-        /// \brief return average suggested step length
-        double suggestedAverage () const;
-
         /// \brief Previous step length. This is the length of the step that
         ///        was taken to arrive at this time.
         double stepLengthTaken () const;
@@ -100,14 +98,13 @@ namespace Opm
         const double start_time_;
         const double total_time_;
         const int report_step_;
+        const double max_time_step_;
 
         double current_time_;
         double dt_;
         int current_step_;
 
         std::vector< double > steps_;
-        double suggestedMax_;
-        double suggestedAverage_;
     };
 
 } // namespace Opm

opm/core/simulator/AdaptiveTimeStepping.hpp

@@ -77,6 +77,7 @@ namespace Opm {
         const double initial_fraction_;       //!< fraction to take as a guess for initial time interval
         const double restart_factor_;         //!< factor to multiply time step with when solver fails to converge
         const double growth_factor_;          //!< factor to multiply time step when solver recovered from failed convergence
+        const double max_time_step_;          //!< maximal allowed time step size
         const int solver_restart_max_;        //!< how many restart of solver are allowed
         const bool solver_verbose_;           //!< solver verbosity
         const bool timestep_verbose_;         //!< timestep verbosity

opm/core/simulator/AdaptiveTimeStepping_impl.hpp

@@ -25,7 +25,7 @@
 
 #include <opm/core/simulator/SimulatorTimer.hpp>
 #include <opm/core/simulator/AdaptiveSimulatorTimer.hpp>
-#include <opm/core/simulator/PIDTimeStepControl.hpp>
+#include <opm/core/simulator/TimeStepControl.hpp>
 
 namespace Opm {
 
@@ -37,13 +37,17 @@ namespace Opm {
         , initial_fraction_( param.getDefault("solver.initialfraction", double(0.25) ) )
         , restart_factor_( param.getDefault("solver.restartfactor", double(0.1) ) )
         , growth_factor_( param.getDefault("solver.growthfactor", double(1.25) ) )
+          // default is 1 year, convert to seconds
+        , max_time_step_( unit::convert::from(param.getDefault("timestep.max_timestep_in_days", 365.0 ), unit::day) )
         , solver_restart_max_( param.getDefault("solver.restart", int(3) ) )
         , solver_verbose_( param.getDefault("solver.verbose", bool(false) ) )
         , timestep_verbose_( param.getDefault("timestep.verbose", bool(false) ) )
         , last_timestep_( -1.0 )
     {
         // valid are "pid" and "pid+iteration"
-        std::string control = param.getDefault("timestep.control", std::string("pid") );
+        std::string control = param.getDefault("timestep.control", std::string("pid+iteration") );
+        // iterations is the accumulation of all linear iterations over all newton steops per time step
+        const int defaultTargetIterations = 30;
 
         const double tol = param.getDefault("timestep.control.tol", double(1e-3) );
         if( control == "pid" ) {
@@ -51,10 +55,17 @@ namespace Opm {
         }
         else if ( control == "pid+iteration" )
         {
-            const int iterations   = param.getDefault("timestep.control.targetiteration", int(25) );
+            const int iterations   = param.getDefault("timestep.control.targetiteration", defaultTargetIterations );
             const double maxgrowth = param.getDefault("timestep.control.maxgrowth", double(3.0) );
             timeStepControl_ = TimeStepControlType( new PIDAndIterationCountTimeStepControl( iterations, tol, maxgrowth ) );
         }
+        else if ( control == "iterationcount" )
+        {
+            const int iterations    = param.getDefault("timestep.control.targetiteration", defaultTargetIterations );
+            const double decayrate  = param.getDefault("timestep.control.decayrate",  double(0.75) );
+            const double growthrate = param.getDefault("timestep.control.growthrate", double(1.25) );
+            timeStepControl_ = TimeStepControlType( new SimpleIterationCountTimeStepControl( iterations, decayrate, growthrate ) );
+        }
         else
             OPM_THROW(std::runtime_error,"Unsupported time step control selected "<< control );
 
@@ -93,7 +104,7 @@ namespace Opm {
         }
 
         // create adaptive step timer with previously used sub step size
-        AdaptiveSimulatorTimer substepTimer( simulatorTimer, last_timestep_ );
+        AdaptiveSimulatorTimer substepTimer( simulatorTimer, last_timestep_, max_time_step_ );
 
         // copy states in case solver has to be restarted (to be revised)
         State  last_state( state );

opm/core/simulator/TimeStepControl.cpp

@@ -20,12 +20,66 @@
 #include <cassert>
 #include <cmath>
 #include <iostream>
+#include <stdexcept>
 
+#include <opm/core/utility/ErrorMacros.hpp>
 #include <opm/core/utility/Units.hpp>
-#include <opm/core/simulator/PIDTimeStepControl.hpp>
+#include <opm/core/simulator/TimeStepControl.hpp>
 
 namespace Opm
 {
+    ////////////////////////////////////////////////////////
+    //
+    //  InterationCountTimeStepControl Implementation
+    //
+    ////////////////////////////////////////////////////////
+
+    SimpleIterationCountTimeStepControl::
+    SimpleIterationCountTimeStepControl( const int target_iterations,
+                                         const double decayrate,
+                                         const double growthrate,
+                                         const bool verbose)
+        : target_iterations_( target_iterations )
+        , decayrate_( decayrate )
+        , growthrate_( growthrate )
+        , verbose_( verbose )
+    {
+        if( decayrate_  > 1.0 ) {
+            OPM_THROW(std::runtime_error,"SimpleIterationCountTimeStepControl: decay should be <= 1 " << decayrate_ );
+        }
+        if( growthrate_ < 1.0 ) {
+            OPM_THROW(std::runtime_error,"SimpleIterationCountTimeStepControl: growth should be >= 1 " << growthrate_ );
+        }
+    }
+
+    double SimpleIterationCountTimeStepControl::
+    computeTimeStepSize( const double dt, const int iterations, const SimulatorState& state ) const
+    {
+        double dtEstimate = dt ;
+
+        // reduce the time step size if we exceed the number of target iterations
+        if( iterations > target_iterations_ )
+        {
+            // scale dtEstimate down with a given rate
+            dtEstimate *= decayrate_;
+        }
+        // increase the time step size if we are below the number of target iterations
+        else if ( iterations < target_iterations_-1 )
+        {
+            dtEstimate *= growthrate_;
+        }
+
+        return dtEstimate;
+    }
+
+
+
+    ////////////////////////////////////////////////////////
+    //
+    //  PIDTimeStepControl Implementation
+    //
+    ////////////////////////////////////////////////////////
+
     PIDTimeStepControl::PIDTimeStepControl( const double tol, const bool verbose )
         : p0_()
         , sat0_()
@@ -99,6 +153,13 @@ namespace Opm
     }
 
 
+
+    ////////////////////////////////////////////////////////////
+    //
+    //  PIDAndIterationCountTimeStepControl  Implementation
+    //
+    ////////////////////////////////////////////////////////////
+
     PIDAndIterationCountTimeStepControl::
     PIDAndIterationCountTimeStepControl( const int target_iterations,
                                          const double tol,
@@ -122,7 +183,7 @@ namespace Opm
         }
 
         // limit the growth of the timestep size by the growth factor
-        dtEstimate = std::max( dtEstimate, double(maxgrowth_ * dt) );
+        dtEstimate = std::min( dtEstimate, double(maxgrowth_ * dt) );
 
         return dtEstimate;
     }

opm/core/simulator/TimeStepControl.hpp

@@ -16,8 +16,8 @@
   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_PIDTIMESTEPCONTROL_HEADER_INCLUDED
+#ifndef OPM_TIMESTEPCONTROL_HEADER_INCLUDED
-#define OPM_PIDTIMESTEPCONTROL_HEADER_INCLUDED
+#define OPM_TIMESTEPCONTROL_HEADER_INCLUDED
 
 #include <vector>
 
@@ -25,6 +25,34 @@
 
 namespace Opm
 {
+    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    ///
+    ///  A simple iteration count based adaptive time step control.
+    //
+    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    class SimpleIterationCountTimeStepControl : public TimeStepControlInterface
+    {
+    public:
+        /// \brief constructor
+        /// \param target_iterations  number of desired iterations (e.g. Newton iterations) per time step in one time step
+        //  \param decayrate          decayrate of time step when target iterations are not met (should be <= 1)
+        //  \param growthrate         growthrate of time step when target iterations are not met (should be >= 1)
+        /// \param verbose            if true get some output (default = false)
+        SimpleIterationCountTimeStepControl( const int target_iterations,
+                                             const double decayrate,
+                                             const double growthrate,
+                                             const bool verbose = false);
+
+        /// \brief \copydoc TimeStepControlInterface::computeTimeStepSize
+        double computeTimeStepSize( const double dt, const int iterations, const SimulatorState& state ) const;
+
+    protected:
+        const int     target_iterations_;
+        const double  decayrate_;
+        const double  growthrate_;
+        const bool    verbose_;
+    };
+
     ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
     ///
     ///  PID controller based adaptive time step control as suggested in:
@@ -90,7 +118,7 @@ namespace Opm
         /// \param target_iterations  number of desired iterations per time step
         /// \param tol        tolerance for the relative changes of the numerical solution to be accepted
         ///                   in one time step (default is 1e-3)
-        //  \param maxgrodth  max growth factor for new time step in relation of old time step (default = 3.0)
+        //  \param maxgrowth  max growth factor for new time step in relation of old time step (default = 3.0)
         /// \param verbose    if true get some output (default = false)
         PIDAndIterationCountTimeStepControl( const int target_iterations = 20,
                                              const double tol = 1e-3,