improvement of the time step adjustment.

opm/simulators/timestepping/SimulatorTimer.hpp

@@ -1,5 +1,6 @@
 /*
   Copyright 2012 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 IRIS AS
 
   This file is part of the Open Porous Media project (OPM).
 
@@ -24,6 +25,9 @@
 
 #include <iosfwd>
 #include <vector>
+#include <algorithm>
+#include <numeric>
+
 #include <boost/date_time/gregorian/gregorian.hpp>
 #include <boost/date_time/posix_time/posix_time_types.hpp>
 
@@ -120,14 +124,32 @@ namespace Opm
         int sub_step_;
 
         std::vector< double > steps_;
+        double suggestedMax_;
+        double suggestedAverage_;
+
+        double computeInitialTimeStep( const double lastDt ) const 
+        {
+            const double maxTimeStep = simulator_timer_.currentStepLength();
+            const double fraction = (lastDt / maxTimeStep);
+            // when lastDt and maxTimeStep are close together, choose the max time step
+            if( fraction > 0.95 ) return       maxTimeStep;
+            // if lastDt is still pretty large, choose half step size since we have to
+            // do two steps anyway
+            // if( fraction > 0.85 ) return 0.5 * maxTimeStep;
+
+            // otherwise choose lastDt
+            return std::min( lastDt, maxTimeStep );
+        }
     public:
         SubStepSimulatorTimer( const SimulatorTimer& st, const double lastDt )
             : simulator_timer_( st )
-            , dt_( lastDt )
+            , dt_( computeInitialTimeStep( lastDt ) )
             , current_time_( simulator_timer_.simulationTimeElapsed() )
             , total_time_( current_time_ + simulator_timer_.currentStepLength() )
             , sub_step_( 0 )
             , steps_()
+            , suggestedMax_( 0.0 )
+            , suggestedAverage_( 0.0 )
         {
             steps_.reserve( 10 );
         }
@@ -139,9 +161,31 @@ namespace Opm
             // store used time step sizes
             steps_.push_back( dt_ );
 
-            double remaining = total_time_ - current_time_;
+            // store some information about the time steps suggested 
-            // set new time step (depending on remaining time)
+            suggestedMax_      = std::max( new_dt, suggestedMax_ );
-            dt_ = ( new_dt > remaining && remaining > 0 ) ? remaining : new_dt ;
+            suggestedAverage_ += new_dt;
+
+            double remaining = (total_time_ - current_time_);
+
+            if( remaining > 0 ) {
+
+                // set new time step (depending on remaining time)
+                if( 1.5 * new_dt > remaining ) {
+                    dt_ = remaining; 
+                    return ;
+                }
+
+                // check for half interval step to avoid very small step at the end
+                // remaining *= 0.5;
+
+                if( 2.25 * new_dt > remaining ) {
+                    dt_ = 0.5 * remaining ;
+                    return ;
+                }
+            }
+
+            // otherwise set new_dt as is
+            dt_ = new_dt;
         }
 
         int currentStepNum () const { return sub_step_; }
@@ -158,14 +202,43 @@ namespace Opm
 
         bool done () const { return (current_time_ >= total_time_) ; }
 
+        double averageStepLength() const
+        {
+            const int size = steps_.size();
+            if( size == 0 ) return 0.0;
+
+            const double sum = std::accumulate(steps_.begin(), steps_.end(), 0.0);
+            return sum / double(size);
+        }
+
+        double maxStepLength () const
+        {
+            if( steps_.size() == 0 ) return 0.0;
+            return *(std::max_element( steps_.begin(), steps_.end() ));
+        }
+
+        double minStepLength () const
+        {
+            if( steps_.size() == 0 ) return 0.0;
+            return *(std::min_element( steps_.begin(), steps_.end() ));
+        }
+
+        double suggestedMax () const { return suggestedMax_; }
+        double suggestedAverage () const 
+        { 
+            const int size = steps_.size();
+            return (size > 0 ) ? (suggestedAverage_ / double(size)) : suggestedAverage_; 
+        }
+
         void report(std::ostream& os) const
         {
-            os << "Sub steps started at time = " << simulator_timer_.simulationTimeElapsed() << std::endl;
+            const double factor = 24.0 * 3600.0;
+            os << "Sub steps started at time = " << simulator_timer_.simulationTimeElapsed()/factor << " (days)" << std::endl;
             for( size_t i=0; i<steps_.size(); ++i )
             {
-                os << " step[ " << i << " ] = " << steps_[ i ] << std::endl;
+                os << " step[ " << i << " ] = " << steps_[ i ]/factor << " (days)" << std::endl;
             }
-            std::cout << "sub steps end time = " << simulationTimeElapsed() << std::endl;
+            std::cout << "sub steps end time = " << simulationTimeElapsed()/factor << " (days)" << std::endl;
         }
     };