the adaptive time stepping utility classes.

This commit is contained in:
Robert K 2014-10-06 14:06:07 +02:00
parent 8a17e5e5d6
commit 294b899ee8
4 changed files with 231 additions and 1 deletions

View File

@ -57,12 +57,13 @@ namespace Opm
}
/// \brief advance time by currentStepLength
void advance()
AdaptiveSimulatorTimer& operator++ ()
{
++current_step_;
current_time_ += dt_;
// store used time step sizes
steps_.push_back( dt_ );
return *this;
}
/// \brief provide and estimate for new time step size

View File

@ -0,0 +1,48 @@
#ifndef OPM_SUBSTEPPING_HEADER_INCLUDED
#define OPM_SUBSTEPPING_HEADER_INCLUDED
#include <iostream>
#include <utility>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/simulator/TimeStepControlInterface.hpp>
namespace Opm {
// AdaptiveTimeStepping
//---------------------
class AdaptiveTimeStepping
{
public:
//! \brief contructor taking parameter object
AdaptiveTimeStepping( const parameter::ParameterGroup& param );
/** \brief step method that acts like the solver::step method
in a sub cycle of time steps
\param solver solver object that must implement a method step( dt, state, well_state )
\param state current state of the solution variables
\param well_state additional well state object
\param time current simulation time
\param timestep current time step length that is to be sub cycled
*/
template <class Solver, class State, class WellState>
void step( Solver& solver, State& state, WellState& well_state,
const double time, const double timestep );
protected:
typedef std::unique_ptr< TimeStepControlInterface > TimeStepControlType;
TimeStepControlType timeStepControl_; //!< time step control object
const double restart_factor_; //!< factor to multiply time step with when solver fails to converge
const int solver_restart_max_; //!< how many restart of solver are allowed
const bool solver_verbose_; //!< solver verbosity
const bool timestep_verbose_; //!< timestep verbosity
double last_timestep_; //!< size of last timestep
};
}
#include <opm/core/simulator/AdaptiveTimeStepping_impl.hpp>
#endif

View File

@ -0,0 +1,130 @@
#ifndef OPM_ADAPTIVETIMESTEPPING_IMPL_HEADER_INCLUDED
#define OPM_ADAPTIVETIMESTEPPING_IMPL_HEADER_INCLUDED
#include <iostream>
#include <string>
#include <utility>
#include <opm/core/simulator/AdaptiveSimulatorTimer.hpp>
#include <opm/core/simulator/PIDTimeStepControl.hpp>
namespace Opm {
// AdaptiveTimeStepping
//---------------------
AdaptiveTimeStepping::AdaptiveTimeStepping( const parameter::ParameterGroup& param )
: timeStepControl_()
, restart_factor_( param.getDefault("solver.restartfactor", double(0.1) ) )
, 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_( std::numeric_limits< double >::max() )
{
// valid are "pid" and "pid+iteration"
std::string control = param.getDefault("timestep.control", std::string("pid") );
const double tol = param.getDefault("timestep.control.tol", double(1e-3) );
if( control == "pid" )
timeStepControl_ = TimeStepControlType( new PIDTimeStepControl( tol ) );
else if ( control == "pid+iteration" )
{
const int iterations = param.getDefault("timestep.control.targetiteration", int(25) );
timeStepControl_ = TimeStepControlType( new PIDAndIterationCountTimeStepControl( iterations, tol ) );
}
else
OPM_THROW(std::runtime_error,"Unsupported time step control selected "<< control );
}
template <class Solver, class State, class WellState>
void AdaptiveTimeStepping::
step( Solver& solver, State& state, WellState& well_state,
const double time, const double timestep )
{
// create adaptive step timer with previously used sub step size
AdaptiveSimulatorTimer timer( time, time+timestep, last_timestep_ );
// copy states in case solver has to be restarted (to be revised)
State last_state( state );
WellState last_well_state( well_state );
// counter for solver restarts
int restarts = 0;
// sub step time loop
while( ! timer.done() )
{
// initialize time step control in case current state is needed later
timeStepControl_->initialize( state );
int linearIterations = -1;
try {
// (linearIterations < 0 means on convergence in solver)
linearIterations = solver.step(timer.currentStepLength(), state, well_state);
if( solver_verbose_ ) {
// report number of linear iterations
std::cout << "Overall linear iterations used: " << linearIterations << std::endl;
}
}
catch (Opm::NumericalProblem)
{
// since linearIterations is < 0 this will restart the solver
}
// (linearIterations < 0 means on convergence in solver)
if( linearIterations >= 0 )
{
// advance by current dt
++timer;
// compute new time step estimate
const double dtEstimate =
timeStepControl_->computeTimeStepSize( timer.currentStepLength(), linearIterations, state );
if( timestep_verbose_ )
std::cout << "Suggested time step size = " << dtEstimate/86400.0 << " (days)" << std::endl;
// set new time step length
timer.provideTimeStepEstimate( dtEstimate );
// update states
last_state = state ;
last_well_state = well_state;
}
else // in case of no convergence
{
// increase restart counter
if( restarts >= solver_restart_max_ ) {
OPM_THROW(Opm::NumericalProblem,"Solver failed to converge after " << restarts << " restarts.");
}
const double newTimeStep = restart_factor_ * timer.currentStepLength();
// we need to revise this
timer.provideTimeStepEstimate( newTimeStep );
if( solver_verbose_ )
std::cerr << "Solver convergence failed, restarting solver with new time step ("<< newTimeStep <<" days)." << std::endl;
// reset states
state = last_state;
well_state = last_well_state;
++restarts;
}
}
// store last small time step for next reportStep
last_timestep_ = timer.suggestedAverage();
if( timestep_verbose_ )
{
timer.report( std::cout );
std::cout << "Last suggested step size = " << last_timestep_/86400.0 << " (days)" << std::endl;
}
if( ! std::isfinite( last_timestep_ ) ) // check for NaN
last_timestep_ = timestep;
}
}
#endif

View File

@ -0,0 +1,51 @@
/*
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_TIMESTEPCONTROLINTERFACE_HEADER_INCLUDED
#define OPM_TIMESTEPCONTROLINTERFACE_HEADER_INCLUDED
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 () {}
};
}
#endif