From fbaa790e70bb270f236bb3597010795186e47445 Mon Sep 17 00:00:00 2001 From: Robert K Date: Thu, 5 Feb 2015 14:43:43 +0100 Subject: [PATCH] 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 | 4 +- opm/core/simulator/AdaptiveSimulatorTimer.cpp | 38 ++++------- opm/core/simulator/AdaptiveSimulatorTimer.hpp | 17 ++--- opm/core/simulator/AdaptiveTimeStepping.hpp | 1 + .../simulator/AdaptiveTimeStepping_impl.hpp | 19 ++++-- ...imeStepControl.cpp => TimeStepControl.cpp} | 65 ++++++++++++++++++- ...imeStepControl.hpp => TimeStepControl.hpp} | 34 +++++++++- 7 files changed, 133 insertions(+), 45 deletions(-) rename opm/core/simulator/{PIDTimeStepControl.cpp => TimeStepControl.cpp} (67%) rename opm/core/simulator/{PIDTimeStepControl.hpp => TimeStepControl.hpp} (73%) diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index edd8c8b5..abca78ac 100644 --- a/CMakeLists_files.cmake +++ b/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 diff --git a/opm/core/simulator/AdaptiveSimulatorTimer.cpp b/opm/core/simulator/AdaptiveSimulatorTimer.cpp index 7e5c200e..fd906eec 100644 --- a/opm/core/simulator/AdaptiveSimulatorTimer.cpp +++ b/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 diff --git a/opm/core/simulator/AdaptiveSimulatorTimer.hpp b/opm/core/simulator/AdaptiveSimulatorTimer.hpp index 870fffd8..ac4bb6ed 100644 --- a/opm/core/simulator/AdaptiveSimulatorTimer.hpp +++ b/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 - explicit AdaptiveSimulatorTimer( const SimulatorTimerInterface& timer, const double lastStepTaken ); + /// \param timer in case of sub stepping this is the outer timer + /// \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::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 diff --git a/opm/core/simulator/AdaptiveTimeStepping.hpp b/opm/core/simulator/AdaptiveTimeStepping.hpp index fe5f860a..1cd38c53 100644 --- a/opm/core/simulator/AdaptiveTimeStepping.hpp +++ b/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 diff --git a/opm/core/simulator/AdaptiveTimeStepping_impl.hpp b/opm/core/simulator/AdaptiveTimeStepping_impl.hpp index c77ec1d8..7f467c98 100644 --- a/opm/core/simulator/AdaptiveTimeStepping_impl.hpp +++ b/opm/core/simulator/AdaptiveTimeStepping_impl.hpp @@ -25,7 +25,7 @@ #include #include -#include +#include 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 ); diff --git a/opm/core/simulator/PIDTimeStepControl.cpp b/opm/core/simulator/TimeStepControl.cpp similarity index 67% rename from opm/core/simulator/PIDTimeStepControl.cpp rename to opm/core/simulator/TimeStepControl.cpp index f4d51073..98114eb2 100644 --- a/opm/core/simulator/PIDTimeStepControl.cpp +++ b/opm/core/simulator/TimeStepControl.cpp @@ -20,12 +20,66 @@ #include #include #include +#include +#include #include -#include +#include 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; } diff --git a/opm/core/simulator/PIDTimeStepControl.hpp b/opm/core/simulator/TimeStepControl.hpp similarity index 73% rename from opm/core/simulator/PIDTimeStepControl.hpp rename to opm/core/simulator/TimeStepControl.hpp index 19ab28f2..ed8f1369 100644 --- a/opm/core/simulator/PIDTimeStepControl.hpp +++ b/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 . */ -#ifndef OPM_PIDTIMESTEPCONTROL_HEADER_INCLUDED -#define OPM_PIDTIMESTEPCONTROL_HEADER_INCLUDED +#ifndef OPM_TIMESTEPCONTROL_HEADER_INCLUDED +#define OPM_TIMESTEPCONTROL_HEADER_INCLUDED #include @@ -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,