AdaptiveTimeStepping_impl.hpp
1 /*
2  Copyright 2014 IRIS AS
3 
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 */
19 #ifndef OPM_ADAPTIVETIMESTEPPING_IMPL_HEADER_INCLUDED
20 #define OPM_ADAPTIVETIMESTEPPING_IMPL_HEADER_INCLUDED
21 
22 #include <iostream>
23 #include <string>
24 #include <utility>
25 
26 #include <opm/simulators/timestepping/SimulatorTimer.hpp>
27 #include <opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp>
28 #include <opm/simulators/timestepping/TimeStepControl.hpp>
29 #include <opm/core/utility/StopWatch.hpp>
30 #include <opm/common/Exceptions.hpp>
31 #include <opm/common/OpmLog/OpmLog.hpp>
32 #include <dune/istl/istlexception.hh>
33 #include <dune/istl/ilu.hh> // For MatrixBlockException
34 #include <opm/parser/eclipse/EclipseState/Schedule/Tuning.hpp>
35 
36 namespace Opm {
37 
38  namespace detail
39  {
40  template <class Solver, class State>
42  {
43  const Solver& solver_;
44  const State& previous_;
45  const State& current_;
46  public:
47  SolutionTimeErrorSolverWrapper( const Solver& solver,
48  const State& previous,
49  const State& current )
50  : solver_( solver ),
51  previous_( previous ),
52  current_( current )
53  {}
54 
56  double relativeChange() const
57  {
58  return solver_.model().relativeChange( previous_, current_ );
59  }
60  };
61 
62  template<class E>
63  void logException(const E& exception, bool verbose)
64  {
65  if( verbose )
66  {
67  std::string message;
68  message = "Caught Exception: ";
69  message += exception.what();
70  OpmLog::debug(message);
71  }
72  }
73  }
74 
75  // AdaptiveTimeStepping
76  //---------------------
77  inline AdaptiveTimeStepping::AdaptiveTimeStepping( const Tuning& tuning,
78  size_t time_step,
79  const ParameterGroup& param,
80  const bool terminal_output )
81  : timeStepControl_()
82  , restart_factor_( tuning.getTSFCNV(time_step) )
83  , growth_factor_(tuning.getTFDIFF(time_step) )
84  , max_growth_( tuning.getTSFMAX(time_step) )
85  // default is 1 year, convert to seconds
86  , max_time_step_( tuning.getTSMAXZ(time_step) )
87  , solver_restart_max_( param.getDefault("solver.restart", int(10) ) )
88  , solver_verbose_( param.getDefault("solver.verbose", bool(true) ) && terminal_output )
89  , timestep_verbose_( param.getDefault("timestep.verbose", bool(true) ) && terminal_output )
90  , suggested_next_timestep_( tuning.getTSINIT(time_step) )
91  , full_timestep_initially_( param.getDefault("full_timestep_initially", bool(false) ) )
92  , timestep_after_event_( tuning.getTMAXWC(time_step))
93  , use_newton_iteration_(false)
94  {
95  init(param);
96 
97  }
98 
99  inline AdaptiveTimeStepping::AdaptiveTimeStepping( const ParameterGroup& param,
100  const bool terminal_output )
101  : timeStepControl_()
102  , restart_factor_( param.getDefault("solver.restartfactor", double(0.33) ) )
103  , growth_factor_( param.getDefault("solver.growthfactor", double(2) ) )
104  , max_growth_( param.getDefault("timestep.control.maxgrowth", double(3.0) ) )
105  // default is 1 year, convert to seconds
106  , max_time_step_( unit::convert::from(param.getDefault("timestep.max_timestep_in_days", 365.0 ), unit::day) )
107  , solver_restart_max_( param.getDefault("solver.restart", int(10) ) )
108  , solver_verbose_( param.getDefault("solver.verbose", bool(true) ) && terminal_output )
109  , timestep_verbose_( param.getDefault("timestep.verbose", bool(true) ) && terminal_output )
110  , suggested_next_timestep_( unit::convert::from(param.getDefault("timestep.initial_timestep_in_days", -1.0 ), unit::day) )
111  , full_timestep_initially_( param.getDefault("full_timestep_initially", bool(false) ) )
112  , timestep_after_event_( unit::convert::from(param.getDefault("timestep.timestep_in_days_after_event", -1.0 ), unit::day))
113  , use_newton_iteration_(false)
114  {
115  init(param);
116  }
117 
118  inline void AdaptiveTimeStepping::
119  init(const ParameterGroup& param)
120  {
121  // valid are "pid" and "pid+iteration"
122  std::string control = param.getDefault("timestep.control", std::string("pid") );
123  // iterations is the accumulation of all linear iterations over all newton steops per time step
124  const int defaultTargetIterations = 30;
125  const int defaultTargetNewtonIterations = 8;
126 
127  const double tol = param.getDefault("timestep.control.tol", double(1e-1) );
128  if( control == "pid" ) {
129  timeStepControl_ = TimeStepControlType( new PIDTimeStepControl( tol ) );
130  }
131  else if ( control == "pid+iteration" )
132  {
133  const int iterations = param.getDefault("timestep.control.targetiteration", defaultTargetIterations );
134  timeStepControl_ = TimeStepControlType( new PIDAndIterationCountTimeStepControl( iterations, tol ) );
135  }
136  else if ( control == "pid+newtoniteration" )
137  {
138  const int iterations = param.getDefault("timestep.control.targetiteration", defaultTargetNewtonIterations );
139  timeStepControl_ = TimeStepControlType( new PIDAndIterationCountTimeStepControl( iterations, tol ) );
140  use_newton_iteration_ = true;
141  }
142  else if ( control == "iterationcount" )
143  {
144  const int iterations = param.getDefault("timestep.control.targetiteration", defaultTargetIterations );
145  const double decayrate = param.getDefault("timestep.control.decayrate", double(0.75) );
146  const double growthrate = param.getDefault("timestep.control.growthrate", double(1.25) );
147  timeStepControl_ = TimeStepControlType( new SimpleIterationCountTimeStepControl( iterations, decayrate, growthrate ) );
148  } else if ( control == "hardcoded") {
149  const std::string filename = param.getDefault("timestep.control.filename", std::string("timesteps"));
150  timeStepControl_ = TimeStepControlType( new HardcodedTimeStepControl( filename ) );
151 
152  }
153  else
154  OPM_THROW(std::runtime_error,"Unsupported time step control selected "<< control );
155 
156  // make sure growth factor is something reasonable
157  assert( growth_factor_ >= 1.0 );
158  }
159 
160 
161 
162  template <class Solver, class State, class WellState>
163  SimulatorReport AdaptiveTimeStepping::
164  step( const SimulatorTimer& simulatorTimer, Solver& solver, State& state, WellState& well_state, const bool event )
165  {
166  return stepImpl( simulatorTimer, solver, state, well_state, event, nullptr, nullptr );
167  }
168 
169  template <class Solver, class State, class WellState, class Output>
170  SimulatorReport AdaptiveTimeStepping::
171  step( const SimulatorTimer& simulatorTimer,
172  Solver& solver, State& state, WellState& well_state,
173  const bool event,
174  Output& outputWriter,
175  const std::vector<int>* fipnum)
176  {
177  return stepImpl( simulatorTimer, solver, state, well_state, event, &outputWriter, fipnum );
178  }
179 
180 
181  // implementation of the step method
182  template <class Solver, class State, class WState, class Output >
183  SimulatorReport AdaptiveTimeStepping::
184  stepImpl( const SimulatorTimer& simulatorTimer,
185  Solver& solver, State& state, WState& well_state,
186  const bool event,
187  Output* outputWriter,
188  const std::vector<int>* fipnum)
189  {
190  SimulatorReport report;
191  const double timestep = simulatorTimer.currentStepLength();
192 
193  // init last time step as a fraction of the given time step
194  if( suggested_next_timestep_ < 0 ) {
196  }
197 
199  suggested_next_timestep_ = timestep;
200  }
201 
202  // use seperate time step after event
203  if (event && timestep_after_event_ > 0) {
205  }
206 
207  // create adaptive step timer with previously used sub step size
208  AdaptiveSimulatorTimer substepTimer( simulatorTimer, suggested_next_timestep_, max_time_step_ );
209 
210  // copy states in case solver has to be restarted (to be revised)
211  State last_state( state );
212  WState last_well_state( well_state );
213 
214  // reset the statistics for the failed substeps
215  failureReport_ = SimulatorReport();
216 
217  // counter for solver restarts
218  int restarts = 0;
219 
220  // sub step time loop
221  while( ! substepTimer.done() )
222  {
223  // get current delta t
224  const double dt = substepTimer.currentStepLength() ;
225  if( timestep_verbose_ )
226  {
227  std::ostringstream ss;
228  ss <<" Substep " << substepTimer.currentStepNum() << ", stepsize "
229  << unit::convert::to(substepTimer.currentStepLength(), unit::day) << " days.";
230  OpmLog::info(ss.str());
231  }
232 
233  SimulatorReport substepReport;
234  std::string cause_of_failure = "";
235  try {
236  substepReport = solver.step( substepTimer, state, well_state);
237  report += substepReport;
238 
239  if( solver_verbose_ ) {
240  // report number of linear iterations
241  OpmLog::note("Overall linear iterations used: " + std::to_string(substepReport.total_linear_iterations));
242  }
243  }
244  catch (const Opm::TooManyIterations& e) {
245  substepReport += solver.failureReport();
246  cause_of_failure = "Solver convergence failure - Iteration limit reached";
247 
248  detail::logException(e, solver_verbose_);
249  // since linearIterations is < 0 this will restart the solver
250  }
251  catch (const Opm::LinearSolverProblem& e) {
252  substepReport += solver.failureReport();
253  cause_of_failure = "Linear solver convergence failure";
254 
255  detail::logException(e, solver_verbose_);
256  // since linearIterations is < 0 this will restart the solver
257  }
258  catch (const Opm::NumericalProblem& e) {
259  substepReport += solver.failureReport();
260  cause_of_failure = "Solver convergence failure - Numerical problem encountered";
261 
262  detail::logException(e, solver_verbose_);
263  // since linearIterations is < 0 this will restart the solver
264  }
265  catch (const std::runtime_error& e) {
266  substepReport += solver.failureReport();
267 
268  detail::logException(e, solver_verbose_);
269  // also catch linear solver not converged
270  }
271  catch (const Dune::ISTLError& e) {
272  substepReport += solver.failureReport();
273 
274  detail::logException(e, solver_verbose_);
275  // also catch errors in ISTL AMG that occur when time step is too large
276  }
277  catch (const Dune::MatrixBlockError& e) {
278  substepReport += solver.failureReport();
279 
280  detail::logException(e, solver_verbose_);
281  // this can be thrown by ISTL's ILU0 in block mode, yet is not an ISTLError
282  }
283 
284  if( substepReport.converged )
285  {
286  // advance by current dt
287  ++substepTimer;
288 
289  // create object to compute the time error, simply forwards the call to the model
290  detail::SolutionTimeErrorSolverWrapper< Solver, State >
291  relativeChange( solver, last_state, state );
292 
293  // compute new time step estimate
294  const int iterations = use_newton_iteration_ ? substepReport.total_newton_iterations
295  : substepReport.total_linear_iterations;
296  double dtEstimate = timeStepControl_->computeTimeStepSize( dt, iterations, relativeChange,
297  substepTimer.simulationTimeElapsed());
298 
299  // limit the growth of the timestep size by the growth factor
300  dtEstimate = std::min( dtEstimate, double(max_growth_ * dt) );
301 
302  // further restrict time step size growth after convergence problems
303  if( restarts > 0 ) {
304  dtEstimate = std::min( growth_factor_ * dt, dtEstimate );
305  // solver converged, reset restarts counter
306  restarts = 0;
307  }
308 
309  if( timestep_verbose_ )
310  {
311  std::ostringstream ss;
312  ss << " Substep summary: ";
313  if (substepReport.total_well_iterations != 0) {
314  ss << "well its = " << std::setw(2) << substepReport.total_well_iterations << ", ";
315  }
316  ss << "newton its = " << std::setw(2) << substepReport.total_newton_iterations << ", "
317  << "linearizations = " << std::setw(2) << substepReport.total_linearizations
318  << " (" << std::fixed << std::setprecision(3) << std::setw(6) << substepReport.assemble_time << " sec), "
319  << "linear its = " << std::setw(3) << substepReport.total_linear_iterations
320  << " (" << std::fixed << std::setprecision(3) << std::setw(6) << substepReport.linear_solve_time << " sec)";
321  OpmLog::info(ss.str());
322  }
323 
324  // write data if outputWriter was provided
325  // if the time step is done we do not need
326  // to write it as this will be done by the simulator
327  // anyway.
328  if( outputWriter && !substepTimer.done() ) {
329  if (fipnum) {
330  solver.computeFluidInPlace(state, *fipnum);
331  }
332  Opm::time::StopWatch perfTimer;
333  perfTimer.start();
334  bool substep = true;
335  const auto& physicalModel = solver.model();
336  outputWriter->writeTimeStep( substepTimer, state, well_state, physicalModel, substep);
337  report.output_write_time += perfTimer.secsSinceStart();
338  }
339 
340  // set new time step length
341  substepTimer.provideTimeStepEstimate( dtEstimate );
342 
343  // update states
344  last_state = state ;
345  last_well_state = well_state;
346 
347  report.converged = substepTimer.done();
348  substepTimer.setLastStepFailed(false);
349 
350  }
351  else // in case of no convergence (linearIterations < 0)
352  {
353  substepTimer.setLastStepFailed(true);
354 
355  failureReport_ += substepReport;
356 
357  // increase restart counter
358  if( restarts >= solver_restart_max_ ) {
359  const auto msg = std::string("Solver failed to converge after cutting timestep ")
360  + std::to_string(restarts) + " times.";
361  if (solver_verbose_) {
362  OpmLog::error(msg);
363  }
364  OPM_THROW_NOLOG(Opm::NumericalProblem, msg);
365  }
366 
367  const double newTimeStep = restart_factor_ * dt;
368  // we need to revise this
369  substepTimer.provideTimeStepEstimate( newTimeStep );
370  if( solver_verbose_ ) {
371  std::string msg;
372  msg = cause_of_failure + "\nTimestep chopped to "
373  + std::to_string(unit::convert::to( substepTimer.currentStepLength(), unit::day )) + " days\n";
374  OpmLog::problem(msg);
375  }
376  // reset states
377  state = last_state;
378  well_state = last_well_state;
379 
380  ++restarts;
381  }
382  }
383 
384 
385  // store estimated time step for next reportStep
386  suggested_next_timestep_ = substepTimer.currentStepLength();
387  if( timestep_verbose_ )
388  {
389  std::ostringstream ss;
390  substepTimer.report(ss);
391  ss << "Suggested next step size = " << unit::convert::to( suggested_next_timestep_, unit::day ) << " (days)" << std::endl;
392  OpmLog::note(ss.str());
393  }
394 
395  if( ! std::isfinite( suggested_next_timestep_ ) ) { // check for NaN
396  suggested_next_timestep_ = timestep;
397  }
398  return report;
399  }
400 }
401 
402 #endif
TimeStepControlType timeStepControl_
time step control object
Definition: AdaptiveTimeStepping.hpp:113
double currentStepLength() const
Current step length.
Definition: SimulatorTimer.cpp:91
double suggested_next_timestep_
suggested size of next timestep
Definition: AdaptiveTimeStepping.hpp:121
const double max_time_step_
maximal allowed time step size
Definition: AdaptiveTimeStepping.hpp:117
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: AdditionalObjectDeleter.hpp:22
SimulatorReport step(const SimulatorTimer &timer, Solver &solver, State &state, WellState &well_state, const bool event)
step method that acts like the solver::step method in a sub cycle of time steps
Definition: AdaptiveTimeStepping_impl.hpp:164
SimulatorReport failureReport_
statistics for the failed substeps of the last timestep
Definition: AdaptiveTimeStepping.hpp:112
const double restart_factor_
factor to multiply time step with when solver fails to converge
Definition: AdaptiveTimeStepping.hpp:114
double relativeChange() const
return || u^n+1 - u^n || / || u^n+1 ||
Definition: AdaptiveTimeStepping_impl.hpp:56
const double max_growth_
factor that limits the maximum growth of a time step
Definition: AdaptiveTimeStepping.hpp:116
const double growth_factor_
factor to multiply time step when solver recovered from failed convergence
Definition: AdaptiveTimeStepping.hpp:115
AdaptiveTimeStepping(const ParameterGroup &param, const bool terminal_output=true)
contructor taking parameter object
Definition: AdaptiveTimeStepping_impl.hpp:99
const int solver_restart_max_
how many restart of solver are allowed
Definition: AdaptiveTimeStepping.hpp:118
const double timestep_after_event_
suggested size of timestep after an event
Definition: AdaptiveTimeStepping.hpp:123
const bool timestep_verbose_
timestep verbosity
Definition: AdaptiveTimeStepping.hpp:120
Definition: AdaptiveTimeStepping_impl.hpp:41
RelativeChangeInterface.
Definition: TimeStepControlInterface.hpp:31
const bool solver_verbose_
solver verbosity
Definition: AdaptiveTimeStepping.hpp:119
PID controller based adaptive time step control as suggested in: Turek and Kuzmin.
Definition: TimeStepControl.hpp:74
Definition: SimulatorTimer.hpp:34
bool use_newton_iteration_
use newton iteration count for adaptive time step control
Definition: AdaptiveTimeStepping.hpp:124
bool full_timestep_initially_
beginning with the size of the time step from data file
Definition: AdaptiveTimeStepping.hpp:122