Set timesteps after events

The time step after an event can either be set using
timestep_in_days_after_event or using the TUNING keyword in the deck.
This commit is contained in:
Tor Harald Sandve 2017-02-24 13:51:10 +01:00
parent aaa4af6c49
commit 19a16ceeca
6 changed files with 62 additions and 20 deletions

View File

@ -105,9 +105,9 @@ namespace Opm
std::ofstream tstep_os(tstep_filename.c_str());
const auto& schedule = eclipse_state_->getSchedule();
const auto& events = schedule.getEvents();
// adaptive time stepping
const auto& events = schedule.getEvents();
std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping;
if( param_.getDefault("timestep.adaptive", true ) )
{
@ -234,8 +234,12 @@ namespace Opm
//
// \Note: The report steps are met in any case
// \Note: The sub stepping will require a copy of the state variables
if( adaptiveTimeStepping ) {
report += adaptiveTimeStepping->step( timer, *solver, state, well_state, output_writer_,
if( adaptiveTimeStepping ) {
bool event = events.hasEvent(ScheduleEvents::NEW_WELL, timer.currentStepNum()) ||
events.hasEvent(ScheduleEvents::PRODUCTION_UPDATE, timer.currentStepNum()) ||
events.hasEvent(ScheduleEvents::INJECTION_UPDATE, timer.currentStepNum()) ||
events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE, timer.currentStepNum());
report += adaptiveTimeStepping->step( timer, *solver, state, well_state, event, output_writer_,
output_writer_.requireFIPNUM() ? &fipnum : nullptr );
}
else {

View File

@ -162,6 +162,7 @@ public:
const auto& schedule = eclState().getSchedule();
// adaptive time stepping
const auto& events = schedule.getEvents();
std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping;
if( param_.getDefault("timestep.adaptive", true ) )
{
@ -293,7 +294,11 @@ public:
// \Note: The report steps are met in any case
// \Note: The sub stepping will require a copy of the state variables
if( adaptiveTimeStepping ) {
report += adaptiveTimeStepping->step( timer, *solver, state, well_state, output_writer_,
bool event = events.hasEvent(ScheduleEvents::NEW_WELL, timer.currentStepNum()) ||
events.hasEvent(ScheduleEvents::PRODUCTION_UPDATE, timer.currentStepNum()) ||
events.hasEvent(ScheduleEvents::INJECTION_UPDATE, timer.currentStepNum()) ||
events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE, timer.currentStepNum());
report += adaptiveTimeStepping->step( timer, *solver, state, well_state, event, output_writer_,
output_writer_.requireFIPNUM() ? &fipnum : nullptr );
}
else {

View File

@ -66,6 +66,7 @@ namespace Opm
std::ofstream tstep_os(tstep_filename.c_str());
// adaptive time stepping
const auto& events = eclipse_state_->getSchedule().getEvents();
std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping;
if( param_.getDefault("timestep.adaptive", true ) )
{
@ -173,7 +174,11 @@ namespace Opm
// \Note: The report steps are met in any case
// \Note: The sub stepping will require a copy of the state variables
if( adaptiveTimeStepping ) {
adaptiveTimeStepping->step( timer, *solver, state, well_state, output_writer_ );
bool event = events.hasEvent(ScheduleEvents::NEW_WELL, timer.currentStepNum()) ||
events.hasEvent(ScheduleEvents::PRODUCTION_UPDATE, timer.currentStepNum()) ||
events.hasEvent(ScheduleEvents::INJECTION_UPDATE, timer.currentStepNum()) ||
events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE, timer.currentStepNum());
adaptiveTimeStepping->step( timer, *solver, state, well_state, event, output_writer_);
}
else {
// solve for complete report step

View File

@ -63,10 +63,12 @@ namespace Opm {
\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 event event status for possible tuning
*/
template <class Solver, class State, class WellState>
SimulatorReport step( const SimulatorTimer& timer,
Solver& solver, State& state, WellState& well_state );
Solver& solver, State& state, WellState& well_state,
const bool event);
/** \brief step method that acts like the solver::step method
in a sub cycle of time steps
@ -76,11 +78,13 @@ namespace Opm {
\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 event event status for possible tuning
\param outputWriter writer object to write sub steps
*/
template <class Solver, class State, class WellState, class Output>
SimulatorReport step( const SimulatorTimer& timer,
Solver& solver, State& state, WellState& well_state,
const bool event,
Output& outputWriter,
const std::vector<int>* fipnum = nullptr);
@ -92,6 +96,7 @@ namespace Opm {
template <class Solver, class State, class WellState, class Output>
SimulatorReport stepImpl( const SimulatorTimer& timer,
Solver& solver, State& state, WellState& well_state,
const bool event,
Output* outputWriter,
const std::vector<int>* fipnum);
@ -109,6 +114,8 @@ namespace Opm {
const bool timestep_verbose_; //!< timestep verbosity
double suggested_next_timestep_; //!< suggested size of next timestep
bool full_timestep_initially_; //!< beginning with the size of the time step from data file
const double timestep_after_event_; //!< suggested size of timestep after an event
bool use_newton_iteration_; //!< use newton iteration count for adaptive time step control
};
}

View File

@ -89,6 +89,8 @@ namespace Opm {
, timestep_verbose_( param.getDefault("timestep.verbose", bool(true) ) && terminal_output )
, suggested_next_timestep_( tuning.getTSINIT(time_step) )
, full_timestep_initially_( param.getDefault("full_timestep_initially", bool(false) ) )
, timestep_after_event_( tuning.getTMAXWC(time_step))
, use_newton_iteration_(false)
{
init(param);
@ -107,6 +109,8 @@ namespace Opm {
, timestep_verbose_( param.getDefault("timestep.verbose", bool(true) ) && terminal_output )
, suggested_next_timestep_( unit::convert::from(param.getDefault("timestep.initial_timestep_in_days", -1.0 ), unit::day) )
, full_timestep_initially_( param.getDefault("full_timestep_initially", bool(false) ) )
, timestep_after_event_( unit::convert::from(param.getDefault("timestep.timestep_in_days_after_event", -1.0 ), unit::day))
, use_newton_iteration_(false)
{
init(param);
}
@ -118,6 +122,7 @@ namespace Opm {
std::string control = param.getDefault("timestep.control", std::string("pid") );
// iterations is the accumulation of all linear iterations over all newton steops per time step
const int defaultTargetIterations = 30;
const int defaultTargetNewtonIterations = 8;
const double tol = param.getDefault("timestep.control.tol", double(1e-1) );
if( control == "pid" ) {
@ -128,6 +133,12 @@ namespace Opm {
const int iterations = param.getDefault("timestep.control.targetiteration", defaultTargetIterations );
timeStepControl_ = TimeStepControlType( new PIDAndIterationCountTimeStepControl( iterations, tol ) );
}
else if ( control == "pid+newtoniteration" )
{
const int iterations = param.getDefault("timestep.control.targetiteration", defaultTargetNewtonIterations );
timeStepControl_ = TimeStepControlType( new PIDAndIterationCountTimeStepControl( iterations, tol ) );
use_newton_iteration_ = true;
}
else if ( control == "iterationcount" )
{
const int iterations = param.getDefault("timestep.control.targetiteration", defaultTargetIterations );
@ -150,19 +161,20 @@ namespace Opm {
template <class Solver, class State, class WellState>
SimulatorReport AdaptiveTimeStepping::
step( const SimulatorTimer& simulatorTimer, Solver& solver, State& state, WellState& well_state )
step( const SimulatorTimer& simulatorTimer, Solver& solver, State& state, WellState& well_state, const bool event )
{
return stepImpl( simulatorTimer, solver, state, well_state, nullptr, nullptr );
return stepImpl( simulatorTimer, solver, state, well_state, event, nullptr, nullptr );
}
template <class Solver, class State, class WellState, class Output>
SimulatorReport AdaptiveTimeStepping::
step( const SimulatorTimer& simulatorTimer,
Solver& solver, State& state, WellState& well_state,
const bool event,
Output& outputWriter,
const std::vector<int>* fipnum)
{
return stepImpl( simulatorTimer, solver, state, well_state, &outputWriter, fipnum );
return stepImpl( simulatorTimer, solver, state, well_state, event, &outputWriter, fipnum );
}
@ -171,6 +183,7 @@ namespace Opm {
SimulatorReport AdaptiveTimeStepping::
stepImpl( const SimulatorTimer& simulatorTimer,
Solver& solver, State& state, WState& well_state,
const bool event,
Output* outputWriter,
const std::vector<int>* fipnum)
{
@ -186,8 +199,10 @@ namespace Opm {
suggested_next_timestep_ = timestep;
}
// TODO
// take change in well state into account
// use seperate time step after event
if (event && timestep_after_event_ > 0) {
suggested_next_timestep_ = timestep_after_event_;
}
// create adaptive step timer with previously used sub step size
AdaptiveSimulatorTimer substepTimer( simulatorTimer, suggested_next_timestep_, max_time_step_ );
@ -249,8 +264,10 @@ namespace Opm {
relativeChange( solver, last_state, state );
// compute new time step estimate
double dtEstimate =
timeStepControl_->computeTimeStepSize( dt, substepReport.total_linear_iterations, relativeChange, substepTimer.simulationTimeElapsed());
const int iterations = use_newton_iteration_ ? substepReport.total_newton_iterations
: substepReport.total_linear_iterations;
double dtEstimate = timeStepControl_->computeTimeStepSize( dt, iterations, relativeChange,
substepTimer.simulationTimeElapsed());
// limit the growth of the timestep size by the growth factor
dtEstimate = std::min( dtEstimate, double(max_growth_ * dt) );

View File

@ -177,16 +177,20 @@ namespace Opm
double PIDAndIterationCountTimeStepControl::
computeTimeStepSize( const double dt, const int iterations, const RelativeChangeInterface& relChange, const double simulationTimeElapsed ) const
{
double dtEstimate = BaseType :: computeTimeStepSize( dt, iterations, relChange, simulationTimeElapsed);
double dtEstimatePID = BaseType :: computeTimeStepSize( dt, iterations, relChange, simulationTimeElapsed);
// further reduce step size if to many iterations were used
if( iterations > target_iterations_ )
{
// if iterations was the same or dts were the same, do some magic
dtEstimate *= double( target_iterations_ ) / double(iterations);
// adjust timesteps based on target iteration
double dtEstimateIter;
if (iterations > target_iterations_) {
double off_target_fraction = double(iterations - target_iterations_) / target_iterations_;
dtEstimateIter = dt / (1.0 + off_target_fraction);
} else {
double off_target_fraction = double(target_iterations_ - iterations) / target_iterations_;
// Be a bit more careful when increasing. The 1.2 factor is from ebos.
dtEstimateIter = dt * (1.0 + off_target_fraction / 1.2);
}
return dtEstimate;
return std::min(dtEstimatePID, dtEstimateIter);
}
} // end namespace Opm