2014-10-06 07:06:07 -05:00
|
|
|
#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_()
|
2014-10-13 04:17:37 -05:00
|
|
|
, initial_fraction_( param.getDefault("solver.initialfraction", double(0.25) ) )
|
2014-10-06 07:06:07 -05:00
|
|
|
, 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) ) )
|
2014-10-13 04:17:37 -05:00
|
|
|
, last_timestep_( -1.0 )
|
2014-10-06 07:06:07 -05:00
|
|
|
{
|
|
|
|
// 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 )
|
|
|
|
{
|
2014-10-13 04:17:37 -05:00
|
|
|
// init last time step as a fraction of the given time step
|
|
|
|
if( last_timestep_ < 0 )
|
|
|
|
last_timestep_ = initial_fraction_ * timestep ;
|
|
|
|
|
2014-10-06 07:06:07 -05:00
|
|
|
// 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() )
|
|
|
|
{
|
2014-10-10 06:55:28 -05:00
|
|
|
// get current delta t
|
|
|
|
const double dt = timer.currentStepLength() ;
|
|
|
|
|
2014-10-06 07:06:07 -05:00
|
|
|
// 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)
|
2014-10-10 06:55:28 -05:00
|
|
|
linearIterations = solver.step( dt, state, well_state);
|
2014-10-06 07:06:07 -05:00
|
|
|
|
|
|
|
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 =
|
2014-10-10 06:55:28 -05:00
|
|
|
timeStepControl_->computeTimeStepSize( dt, linearIterations, state );
|
2014-10-06 07:06:07 -05:00
|
|
|
if( timestep_verbose_ )
|
2014-10-07 02:48:57 -05:00
|
|
|
std::cout << "Suggested time step size = " << unit::convert::to(dtEstimate, unit::day) << " (days)" << std::endl;
|
2014-10-06 07:06:07 -05:00
|
|
|
|
|
|
|
// 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.");
|
|
|
|
}
|
|
|
|
|
2014-10-10 06:55:28 -05:00
|
|
|
const double newTimeStep = restart_factor_ * dt;
|
2014-10-06 07:06:07 -05:00
|
|
|
// we need to revise this
|
|
|
|
timer.provideTimeStepEstimate( newTimeStep );
|
|
|
|
if( solver_verbose_ )
|
2014-10-13 04:17:37 -05:00
|
|
|
std::cerr << "Solver convergence failed, restarting solver with new time step ("
|
|
|
|
<< unit::convert::to( newTimeStep, unit::day ) <<" days)." << std::endl;
|
2014-10-06 07:06:07 -05:00
|
|
|
|
|
|
|
// 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 );
|
2014-10-07 02:48:57 -05:00
|
|
|
std::cout << "Last suggested step size = " << unit::convert::to( last_timestep_, unit::day ) << " (days)" << std::endl;
|
2014-10-06 07:06:07 -05:00
|
|
|
}
|
|
|
|
|
|
|
|
if( ! std::isfinite( last_timestep_ ) ) // check for NaN
|
|
|
|
last_timestep_ = timestep;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
#endif
|