My Project
AdaptiveTimeSteppingEbos.hpp
1 /*
2 */
3 #ifndef OPM_ADAPTIVE_TIME_STEPPING_EBOS_HPP
4 #define OPM_ADAPTIVE_TIME_STEPPING_EBOS_HPP
5 
6 #include <iostream>
7 #include <utility>
8 
9 #include <opm/simulators/timestepping/SimulatorReport.hpp>
10 #include <opm/grid/utility/StopWatch.hpp>
11 #include <opm/common/OpmLog/OpmLog.hpp>
12 #include <opm/common/utility/parameters/ParameterGroup.hpp>
13 #include <opm/common/ErrorMacros.hpp>
14 #include <opm/simulators/timestepping/SimulatorTimer.hpp>
15 #include <opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp>
16 #include <opm/simulators/timestepping/TimeStepControlInterface.hpp>
17 #include <opm/simulators/timestepping/TimeStepControl.hpp>
18 #include <opm/core/props/phaseUsageFromDeck.hpp>
19 #include <opm/input/eclipse/Schedule/ScheduleState.hpp>
20 #include <opm/common/Exceptions.hpp>
21 
22 namespace Opm::Properties {
23 
24 namespace TTag {
26 }
27 
28 template<class TypeTag, class MyTypeTag>
30  using type = UndefinedProperty;
31 };
32 template<class TypeTag, class MyTypeTag>
34  using type = UndefinedProperty;
35 };
36 template<class TypeTag, class MyTypeTag>
38  using type = UndefinedProperty;
39 };
40 template<class TypeTag, class MyTypeTag>
42  using type = UndefinedProperty;
43 };
44 template<class TypeTag, class MyTypeTag>
46  using type = UndefinedProperty;
47 };
48 template<class TypeTag, class MyTypeTag>
50  using type = UndefinedProperty;
51 };
52 template<class TypeTag, class MyTypeTag>
54  using type = UndefinedProperty;
55 };
56 template<class TypeTag, class MyTypeTag>
58  using type = UndefinedProperty;
59 };
60 template<class TypeTag, class MyTypeTag>
62  using type = UndefinedProperty;
63 };
64 template<class TypeTag, class MyTypeTag>
66  using type = UndefinedProperty;
67 };
68 template<class TypeTag, class MyTypeTag>
70  using type = UndefinedProperty;
71 };
72 template<class TypeTag, class MyTypeTag>
74  using type = UndefinedProperty;
75 };
76 template<class TypeTag, class MyTypeTag>
78  using type = UndefinedProperty;
79 };
80 template<class TypeTag, class MyTypeTag>
82  using type = UndefinedProperty;
83 };
84 template<class TypeTag, class MyTypeTag>
86  using type = UndefinedProperty;
87 };
88 template<class TypeTag, class MyTypeTag>
90  using type = UndefinedProperty;
91 };
92 template<class TypeTag, class MyTypeTag>
94  using type = UndefinedProperty;
95 };
96 template<class TypeTag, class MyTypeTag>
98  using type = UndefinedProperty;
99 };
100 template<class TypeTag, class MyTypeTag>
102  using type = UndefinedProperty;
103 };
104 template<class TypeTag, class MyTypeTag>
106  using type = UndefinedProperty;
107 };
108 template<class TypeTag, class MyTypeTag>
110  using type = UndefinedProperty;
111 };
112 template<class TypeTag, class MyTypeTag>
114  using type = UndefinedProperty;
115 };
116 template<class TypeTag, class MyTypeTag>
118  using type = UndefinedProperty;
119 };
120 
121 template<class TypeTag>
122 struct SolverRestartFactor<TypeTag, TTag::FlowTimeSteppingParameters> {
123  using type = GetPropType<TypeTag, Scalar>;
124  static constexpr type value = 0.33;
125 };
126 template<class TypeTag>
127 struct SolverGrowthFactor<TypeTag, TTag::FlowTimeSteppingParameters> {
128  using type = GetPropType<TypeTag, Scalar>;
129  static constexpr type value = 2.0;
130 };
131 template<class TypeTag>
132 struct SolverMaxGrowth<TypeTag, TTag::FlowTimeSteppingParameters> {
133  using type = GetPropType<TypeTag, Scalar>;
134  static constexpr type value = 3.0;
135 };
136 template<class TypeTag>
137 struct SolverMaxTimeStepInDays<TypeTag, TTag::FlowTimeSteppingParameters> {
138  using type = GetPropType<TypeTag, Scalar>;
139  static constexpr type value = 365.0;
140 };
141 template<class TypeTag>
142 struct SolverMinTimeStep<TypeTag, TTag::FlowTimeSteppingParameters> {
143  using type = GetPropType<TypeTag, Scalar>;
144  static constexpr type value = 1.0e-12;
145 };
146 template<class TypeTag>
147 struct SolverContinueOnConvergenceFailure<TypeTag, TTag::FlowTimeSteppingParameters> {
148  static constexpr bool value = false;
149 };
150 template<class TypeTag>
151 struct SolverMaxRestarts<TypeTag, TTag::FlowTimeSteppingParameters> {
152  static constexpr int value = 10;
153 };
154 template<class TypeTag>
155 struct SolverVerbosity<TypeTag, TTag::FlowTimeSteppingParameters> {
156  static constexpr int value = 1;
157 };
158 template<class TypeTag>
159 struct TimeStepVerbosity<TypeTag, TTag::FlowTimeSteppingParameters> {
160  static constexpr int value = 1;
161 };
162 template<class TypeTag>
163 struct InitialTimeStepInDays<TypeTag, TTag::FlowTimeSteppingParameters> {
164  using type = GetPropType<TypeTag, Scalar>;
165  static constexpr type value = 1.0;
166 };
167 template<class TypeTag>
168 struct FullTimeStepInitially<TypeTag, TTag::FlowTimeSteppingParameters> {
169  static constexpr bool value = false;
170 };
171 template<class TypeTag>
172 struct TimeStepAfterEventInDays<TypeTag, TTag::FlowTimeSteppingParameters> {
173  using type = GetPropType<TypeTag, Scalar>;
174  static constexpr type value = -1.0;
175 };
176 template<class TypeTag>
177 struct TimeStepControl<TypeTag, TTag::FlowTimeSteppingParameters> {
178  static constexpr auto value = "pid+newtoniteration";
179 };
180 template<class TypeTag>
181 struct TimeStepControlTolerance<TypeTag, TTag::FlowTimeSteppingParameters> {
182  using type = GetPropType<TypeTag, Scalar>;
183  static constexpr type value = 1e-1;
184 };
185 template<class TypeTag>
186 struct TimeStepControlTargetIterations<TypeTag, TTag::FlowTimeSteppingParameters> {
187  static constexpr int value = 30;
188 };
189 template<class TypeTag>
190 struct TimeStepControlTargetNewtonIterations<TypeTag, TTag::FlowTimeSteppingParameters> {
191  static constexpr int value = 8;
192 };
193 template<class TypeTag>
194 struct TimeStepControlDecayRate<TypeTag, TTag::FlowTimeSteppingParameters> {
195  using type = GetPropType<TypeTag, Scalar>;
196  static constexpr type value = 0.75;
197 };
198 template<class TypeTag>
199 struct TimeStepControlGrowthRate<TypeTag, TTag::FlowTimeSteppingParameters> {
200  using type = GetPropType<TypeTag, Scalar>;
201  static constexpr type value = 1.25;
202 };
203 template<class TypeTag>
204 struct TimeStepControlDecayDampingFactor<TypeTag, TTag::FlowTimeSteppingParameters> {
205  using type = GetPropType<TypeTag, Scalar>;
206  static constexpr type value = 1.0;
207 };
208 template<class TypeTag>
209 struct TimeStepControlGrowthDampingFactor<TypeTag, TTag::FlowTimeSteppingParameters> {
210  using type = GetPropType<TypeTag, Scalar>;
211  static constexpr type value = 3.2;
212 };
213 template<class TypeTag>
214 struct TimeStepControlFileName<TypeTag, TTag::FlowTimeSteppingParameters> {
215  static constexpr auto value = "timesteps";
216 };
217 template<class TypeTag>
218 struct MinTimeStepBeforeShuttingProblematicWellsInDays<TypeTag, TTag::FlowTimeSteppingParameters> {
219  using type = GetPropType<TypeTag, Scalar>;
220  static constexpr type value = 0.01;
221 };
222 
223 template<class TypeTag>
224 struct MinTimeStepBasedOnNewtonIterations<TypeTag, TTag::FlowTimeSteppingParameters> {
225  using type = GetPropType<TypeTag, Scalar>;
226  static constexpr type value = 0.0;
227 };
228 
229 } // namespace Opm::Properties
230 
231 namespace Opm {
232  // AdaptiveTimeStepping
233  //---------------------
234  void logTimer(const AdaptiveSimulatorTimer& substepTimer);
235 
236  template<class TypeTag>
238  {
239  template <class Solver>
240  class SolutionTimeErrorSolverWrapperEbos : public RelativeChangeInterface
241  {
242  const Solver& solver_;
243  public:
244  SolutionTimeErrorSolverWrapperEbos(const Solver& solver)
245  : solver_(solver)
246  {}
247 
249  double relativeChange() const
250  { return solver_.model().relativeChange(); }
251  };
252 
253  template<class E>
254  void logException_(const E& exception, bool verbose)
255  {
256  if (verbose) {
257  std::string message;
258  message = "Caught Exception: ";
259  message += exception.what();
260  OpmLog::debug(message);
261  }
262  }
263 
264  public:
266  AdaptiveTimeSteppingEbos(const UnitSystem& unitSystem,
267  const bool terminalOutput = true)
268  : timeStepControl_()
269  , restartFactor_(EWOMS_GET_PARAM(TypeTag, double, SolverRestartFactor)) // 0.33
270  , growthFactor_(EWOMS_GET_PARAM(TypeTag, double, SolverGrowthFactor)) // 2.0
271  , maxGrowth_(EWOMS_GET_PARAM(TypeTag, double, SolverMaxGrowth)) // 3.0
272  , maxTimeStep_(EWOMS_GET_PARAM(TypeTag, double, SolverMaxTimeStepInDays)*24*60*60) // 365.25
273  , minTimeStep_(unitSystem.to_si(UnitSystem::measure::time, EWOMS_GET_PARAM(TypeTag, double, SolverMinTimeStep))) // 1e-12;
274  , ignoreConvergenceFailure_(EWOMS_GET_PARAM(TypeTag, bool, SolverContinueOnConvergenceFailure)) // false;
275  , solverRestartMax_(EWOMS_GET_PARAM(TypeTag, int, SolverMaxRestarts)) // 10
276  , solverVerbose_(EWOMS_GET_PARAM(TypeTag, int, SolverVerbosity) > 0 && terminalOutput) // 2
277  , timestepVerbose_(EWOMS_GET_PARAM(TypeTag, int, TimeStepVerbosity) > 0 && terminalOutput) // 2
278  , suggestedNextTimestep_(EWOMS_GET_PARAM(TypeTag, double, InitialTimeStepInDays)*24*60*60) // 1.0
279  , fullTimestepInitially_(EWOMS_GET_PARAM(TypeTag, bool, FullTimeStepInitially)) // false
280  , timestepAfterEvent_(EWOMS_GET_PARAM(TypeTag, double, TimeStepAfterEventInDays)*24*60*60) // 1e30
281  , useNewtonIteration_(false)
282  , minTimeStepBeforeShuttingProblematicWells_(EWOMS_GET_PARAM(TypeTag, double, MinTimeStepBeforeShuttingProblematicWellsInDays)*unit::day)
283 
284  {
285  init_(unitSystem);
286  }
287 
288 
289 
293  AdaptiveTimeSteppingEbos(double max_next_tstep,
294  const Tuning& tuning,
295  const UnitSystem& unitSystem,
296  const bool terminalOutput = true)
297  : timeStepControl_()
298  , restartFactor_(tuning.TSFCNV)
299  , growthFactor_(tuning.TFDIFF)
300  , maxGrowth_(tuning.TSFMAX)
301  , maxTimeStep_(tuning.TSMAXZ) // 365.25
302  , minTimeStep_(tuning.TSFMIN) // 0.1;
304  , solverRestartMax_(EWOMS_GET_PARAM(TypeTag, int, SolverMaxRestarts)) // 10
305  , solverVerbose_(EWOMS_GET_PARAM(TypeTag, int, SolverVerbosity) > 0 && terminalOutput) // 2
306  , timestepVerbose_(EWOMS_GET_PARAM(TypeTag, int, TimeStepVerbosity) > 0 && terminalOutput) // 2
307  , suggestedNextTimestep_(max_next_tstep) // 1.0
308  , fullTimestepInitially_(EWOMS_GET_PARAM(TypeTag, bool, FullTimeStepInitially)) // false
309  , timestepAfterEvent_(tuning.TMAXWC) // 1e30
310  , useNewtonIteration_(false)
311  , minTimeStepBeforeShuttingProblematicWells_(EWOMS_GET_PARAM(TypeTag, double, MinTimeStepBeforeShuttingProblematicWellsInDays)*unit::day)
312  {
313  init_(unitSystem);
314  }
315 
316  static void registerParameters()
317  {
318  // TODO: make sure the help messages are correct (and useful)
319  EWOMS_REGISTER_PARAM(TypeTag, double, SolverRestartFactor,
320  "The factor time steps are elongated after restarts");
321  EWOMS_REGISTER_PARAM(TypeTag, double, SolverGrowthFactor,
322  "The factor time steps are elongated after a successful substep");
323  EWOMS_REGISTER_PARAM(TypeTag, double, SolverMaxGrowth,
324  "The maximum factor time steps are elongated after a report step");
325  EWOMS_REGISTER_PARAM(TypeTag, double, SolverMaxTimeStepInDays,
326  "The maximum size of a time step in days");
327  EWOMS_REGISTER_PARAM(TypeTag, double, SolverMinTimeStep,
328  "The minimum size of a time step in days for field and metric and hours for lab. If a step cannot converge without getting cut below this step size the simulator will stop");
329  EWOMS_REGISTER_PARAM(TypeTag, bool, SolverContinueOnConvergenceFailure,
330  "Continue instead of stop when minimum solver time step is reached");
331  EWOMS_REGISTER_PARAM(TypeTag, int, SolverMaxRestarts,
332  "The maximum number of breakdowns before a substep is given up and the simulator is terminated");
333  EWOMS_REGISTER_PARAM(TypeTag, int, SolverVerbosity,
334  "Specify the \"chattiness\" of the non-linear solver itself");
335  EWOMS_REGISTER_PARAM(TypeTag, int, TimeStepVerbosity,
336  "Specify the \"chattiness\" during the time integration");
337  EWOMS_REGISTER_PARAM(TypeTag, double, InitialTimeStepInDays,
338  "The size of the initial time step in days");
339  EWOMS_REGISTER_PARAM(TypeTag, bool, FullTimeStepInitially,
340  "Always attempt to finish a report step using a single substep");
341  EWOMS_REGISTER_PARAM(TypeTag, double, TimeStepAfterEventInDays,
342  "Time step size of the first time step after an event occurs during the simulation in days");
343  EWOMS_REGISTER_PARAM(TypeTag, std::string, TimeStepControl,
344  "The algorithm used to determine time-step sizes. valid options are: 'pid' (default), 'pid+iteration', 'pid+newtoniteration', 'iterationcount', 'newtoniterationcount' and 'hardcoded'");
345  EWOMS_REGISTER_PARAM(TypeTag, double, TimeStepControlTolerance,
346  "The tolerance used by the time step size control algorithm");
347  EWOMS_REGISTER_PARAM(TypeTag, int, TimeStepControlTargetIterations,
348  "The number of linear iterations which the time step control scheme should aim for (if applicable)");
349  EWOMS_REGISTER_PARAM(TypeTag, int, TimeStepControlTargetNewtonIterations,
350  "The number of Newton iterations which the time step control scheme should aim for (if applicable)");
351  EWOMS_REGISTER_PARAM(TypeTag, double, TimeStepControlDecayRate,
352  "The decay rate of the time step size of the number of target iterations is exceeded");
353  EWOMS_REGISTER_PARAM(TypeTag, double, TimeStepControlGrowthRate,
354  "The growth rate of the time step size of the number of target iterations is undercut");
355  EWOMS_REGISTER_PARAM(TypeTag, double, TimeStepControlDecayDampingFactor,
356  "The decay rate of the time step decrease when the target iterations is exceeded");
357  EWOMS_REGISTER_PARAM(TypeTag, double, TimeStepControlGrowthDampingFactor,
358  "The growth rate of the time step increase when the target iterations is undercut");
359  EWOMS_REGISTER_PARAM(TypeTag, std::string, TimeStepControlFileName,
360  "The name of the file which contains the hardcoded time steps sizes");
361  EWOMS_REGISTER_PARAM(TypeTag, double, MinTimeStepBeforeShuttingProblematicWellsInDays,
362  "The minimum time step size in days for which problematic wells are not shut");
363  EWOMS_REGISTER_PARAM(TypeTag, double, MinTimeStepBasedOnNewtonIterations,
364  "The minimum time step size (in days for field and metric unit and hours for lab unit) can be reduced to based on newton iteration counts");
365  }
366 
370  template <class Solver>
371  SimulatorReport step(const SimulatorTimer& simulatorTimer,
372  Solver& solver,
373  const bool isEvent,
374  const std::vector<int>* fipnum = nullptr)
375  {
376  SimulatorReport report;
377  const double timestep = simulatorTimer.currentStepLength();
378 
379  // init last time step as a fraction of the given time step
380  if (suggestedNextTimestep_ < 0) {
382  }
383 
385  suggestedNextTimestep_ = timestep;
386  }
387 
388  // use seperate time step after event
389  if (isEvent && timestepAfterEvent_ > 0) {
391  }
392 
393  auto& ebosSimulator = solver.model().ebosSimulator();
394  auto& ebosProblem = ebosSimulator.problem();
395 
396  // create adaptive step timer with previously used sub step size
397  AdaptiveSimulatorTimer substepTimer(simulatorTimer, suggestedNextTimestep_, maxTimeStep_);
398 
399  // counter for solver restarts
400  int restarts = 0;
401 
402  // sub step time loop
403  while (!substepTimer.done()) {
404  // get current delta t
405  const double dt = substepTimer.currentStepLength() ;
406  if (timestepVerbose_) {
407  logTimer(substepTimer);
408  }
409 
410  SimulatorReportSingle substepReport;
411  std::string causeOfFailure = "";
412  try {
413  substepReport = solver.step(substepTimer);
414  if (solverVerbose_) {
415  // report number of linear iterations
416  OpmLog::debug("Overall linear iterations used: " + std::to_string(substepReport.total_linear_iterations));
417  }
418  }
419  catch (const TooManyIterations& e) {
420  substepReport = solver.failureReport();
421  causeOfFailure = "Solver convergence failure - Iteration limit reached";
422 
423  logException_(e, solverVerbose_);
424  // since linearIterations is < 0 this will restart the solver
425  }
426  catch (const LinearSolverProblem& e) {
427  substepReport = solver.failureReport();
428  causeOfFailure = "Linear solver convergence failure";
429 
430  logException_(e, solverVerbose_);
431  // since linearIterations is < 0 this will restart the solver
432  }
433  catch (const NumericalIssue& e) {
434  substepReport = solver.failureReport();
435  causeOfFailure = "Solver convergence failure - Numerical problem encountered";
436 
437  logException_(e, solverVerbose_);
438  // since linearIterations is < 0 this will restart the solver
439  }
440  catch (const std::runtime_error& e) {
441  substepReport = solver.failureReport();
442 
443  logException_(e, solverVerbose_);
444  // also catch linear solver not converged
445  }
446  catch (const Dune::ISTLError& e) {
447  substepReport = solver.failureReport();
448 
449  logException_(e, solverVerbose_);
450  // also catch errors in ISTL AMG that occur when time step is too large
451  }
452  catch (const Dune::MatrixBlockError& e) {
453  substepReport = solver.failureReport();
454 
455  logException_(e, solverVerbose_);
456  // this can be thrown by ISTL's ILU0 in block mode, yet is not an ISTLError
457  }
458 
459  //Pass substep to eclwriter for summary output
460  ebosSimulator.problem().setSubStepReport(substepReport);
461 
462  report += substepReport;
463 
464  bool continue_on_uncoverged_solution = ignoreConvergenceFailure_ && !substepReport.converged && dt <= minTimeStep_;
465 
466  if (continue_on_uncoverged_solution) {
467  const auto msg = std::string("Solver failed to converge but timestep ")
468  + std::to_string(dt) + " is smaller or equal to "
469  + std::to_string(minTimeStep_) + "\n which is the minimum threshold given"
470  + "by option --solver-min-time-step= \n";
471  if (solverVerbose_) {
472  OpmLog::error(msg);
473  }
474  }
475 
476  if (substepReport.converged || continue_on_uncoverged_solution) {
477 
478  // advance by current dt
479  ++substepTimer;
480 
481  // create object to compute the time error, simply forwards the call to the model
482  SolutionTimeErrorSolverWrapperEbos<Solver> relativeChange(solver);
483 
484  // compute new time step estimate
485  const int iterations = useNewtonIteration_ ? substepReport.total_newton_iterations
486  : substepReport.total_linear_iterations;
487  double dtEstimate = timeStepControl_->computeTimeStepSize(dt, iterations, relativeChange,
488  substepTimer.simulationTimeElapsed());
489 
490  assert(dtEstimate > 0);
491  // limit the growth of the timestep size by the growth factor
492  dtEstimate = std::min(dtEstimate, double(maxGrowth_ * dt));
493  assert(dtEstimate > 0);
494  // further restrict time step size growth after convergence problems
495  if (restarts > 0) {
496  dtEstimate = std::min(growthFactor_ * dt, dtEstimate);
497  // solver converged, reset restarts counter
498  restarts = 0;
499  }
500 
501  // Further restrict time step size if we are in
502  // prediction mode with THP constraints.
503  if (solver.model().wellModel().hasTHPConstraints()) {
504  const double maxPredictionTHPTimestep = 16.0 * unit::day;
505  dtEstimate = std::min(dtEstimate, maxPredictionTHPTimestep);
506  }
507  assert(dtEstimate > 0);
508  if (timestepVerbose_) {
509  std::ostringstream ss;
510  substepReport.reportStep(ss);
511  OpmLog::info(ss.str());
512  }
513 
514  // write data if outputWriter was provided
515  // if the time step is done we do not need
516  // to write it as this will be done by the simulator
517  // anyway.
518  if (!substepTimer.done()) {
519  if (fipnum) {
520  solver.computeFluidInPlace(*fipnum);
521  }
522  time::StopWatch perfTimer;
523  perfTimer.start();
524 
525  ebosProblem.writeOutput();
526 
527  report.success.output_write_time += perfTimer.secsSinceStart();
528  }
529 
530  // set new time step length
531  substepTimer.provideTimeStepEstimate(dtEstimate);
532 
533  report.success.converged = substepTimer.done();
534  substepTimer.setLastStepFailed(false);
535 
536  }
537  else { // in case of no convergence
538  substepTimer.setLastStepFailed(true);
539 
540  // If we have restarted (i.e. cut the timestep) too
541  // many times, we have failed and throw an exception.
542  if (restarts >= solverRestartMax_) {
543  const auto msg = std::string("Solver failed to converge after cutting timestep ")
544  + std::to_string(restarts) + " times.";
545  if (solverVerbose_) {
546  OpmLog::error(msg);
547  }
548  OPM_THROW_NOLOG(NumericalIssue, msg);
549  }
550 
551  // The new, chopped timestep.
552  const double newTimeStep = restartFactor_ * dt;
553 
554 
555  // If we have restarted (i.e. cut the timestep) too
556  // much, we have failed and throw an exception.
557  if (newTimeStep < minTimeStep_) {
558  const auto msg = std::string("Solver failed to converge after cutting timestep to ")
559  + std::to_string(minTimeStep_) + "\n which is the minimum threshold given"
560  + "by option --solver-min-time-step= \n";
561  if (solverVerbose_) {
562  OpmLog::error(msg);
563  }
564  OPM_THROW_NOLOG(NumericalIssue, msg);
565  }
566 
567  // Define utility function for chopping timestep.
568  auto chopTimestep = [&]() {
569  substepTimer.provideTimeStepEstimate(newTimeStep);
570  if (solverVerbose_) {
571  std::string msg;
572  msg = causeOfFailure + "\nTimestep chopped to "
573  + std::to_string(unit::convert::to(substepTimer.currentStepLength(), unit::day)) + " days\n";
574  OpmLog::problem(msg);
575  }
576  ++restarts;
577  };
578 
579  const double minimumChoppedTimestep = minTimeStepBeforeShuttingProblematicWells_;
580  if (newTimeStep > minimumChoppedTimestep) {
581  chopTimestep();
582  } else {
583  // We are below the threshold, and will check if there are any
584  // wells we should close rather than chopping again.
585  std::set<std::string> failing_wells = consistentlyFailingWells(solver.model().stepReports());
586  if (failing_wells.empty()) {
587  // Found no wells to close, chop the timestep as above.
588  chopTimestep();
589  } else {
590  // Close all consistently failing wells.
591  int num_shut_wells = 0;
592  for (const auto& well : failing_wells) {
593  bool was_shut = solver.model().wellModel().forceShutWellByName(well, substepTimer.simulationTimeElapsed());
594  if (was_shut) {
595  ++num_shut_wells;
596  }
597  }
598  if (num_shut_wells == 0) {
599  // None of the problematic wells were shut.
600  // We must fall back to chopping again.
601  chopTimestep();
602  } else {
603  substepTimer.provideTimeStepEstimate(dt);
604  if (solverVerbose_) {
605  std::string msg;
606  msg = "\nProblematic well(s) were shut: ";
607  for (const auto& well : failing_wells) {
608  msg += well;
609  msg += " ";
610  }
611  msg += "(retrying timestep)\n";
612  OpmLog::problem(msg);
613  }
614  }
615  }
616  }
617  }
618  ebosProblem.setNextTimeStepSize(substepTimer.currentStepLength());
619  }
620 
621  // store estimated time step for next reportStep
622  suggestedNextTimestep_ = substepTimer.currentStepLength();
623  if (timestepVerbose_) {
624  std::ostringstream ss;
625  substepTimer.report(ss);
626  ss << "Suggested next step size = " << unit::convert::to(suggestedNextTimestep_, unit::day) << " (days)" << std::endl;
627  OpmLog::debug(ss.str());
628  }
629 
630  if (! std::isfinite(suggestedNextTimestep_)) { // check for NaN
631  suggestedNextTimestep_ = timestep;
632  }
633  return report;
634  }
635 
639  double suggestedNextStep() const
640  { return suggestedNextTimestep_; }
641 
642  void setSuggestedNextStep(const double x)
643  { suggestedNextTimestep_ = x; }
644 
645  void updateTUNING(double max_next_tstep, const Tuning& tuning)
646  {
647  restartFactor_ = tuning.TSFCNV;
648  growthFactor_ = tuning.TFDIFF;
649  maxGrowth_ = tuning.TSFMAX;
650  maxTimeStep_ = tuning.TSMAXZ;
651  suggestedNextTimestep_ = max_next_tstep;
652  timestepAfterEvent_ = tuning.TMAXWC;
653  }
654 
655 
656  protected:
657  void init_(const UnitSystem& unitSystem)
658  {
659  // valid are "pid" and "pid+iteration"
660  std::string control = EWOMS_GET_PARAM(TypeTag, std::string, TimeStepControl); // "pid"
661 
662  const double tol = EWOMS_GET_PARAM(TypeTag, double, TimeStepControlTolerance); // 1e-1
663  if (control == "pid") {
664  timeStepControl_ = TimeStepControlType(new PIDTimeStepControl(tol));
665  }
666  else if (control == "pid+iteration") {
667  const int iterations = EWOMS_GET_PARAM(TypeTag, int, TimeStepControlTargetIterations); // 30
668  const double decayDampingFactor = EWOMS_GET_PARAM(TypeTag, double, TimeStepControlDecayDampingFactor); // 1.0
669  const double growthDampingFactor = EWOMS_GET_PARAM(TypeTag, double, TimeStepControlGrowthDampingFactor); // 3.2
670  timeStepControl_ = TimeStepControlType(new PIDAndIterationCountTimeStepControl(iterations, decayDampingFactor, growthDampingFactor, tol));
671  }
672  else if (control == "pid+newtoniteration") {
673  const int iterations = EWOMS_GET_PARAM(TypeTag, int, TimeStepControlTargetNewtonIterations); // 8
674  const double decayDampingFactor = EWOMS_GET_PARAM(TypeTag, double, TimeStepControlDecayDampingFactor); // 1.0
675  const double growthDampingFactor = EWOMS_GET_PARAM(TypeTag, double, TimeStepControlGrowthDampingFactor); // 3.2
676  const double nonDimensionalMinTimeStepIterations = EWOMS_GET_PARAM(TypeTag, double, MinTimeStepBasedOnNewtonIterations); // 0.0 by default
677  // the min time step can be reduced by the newton iteration numbers
678  double minTimeStepReducedByIterations = unitSystem.to_si(UnitSystem::measure::time, nonDimensionalMinTimeStepIterations);
679  timeStepControl_ = TimeStepControlType(new PIDAndIterationCountTimeStepControl(iterations, decayDampingFactor,
680  growthDampingFactor, tol, minTimeStepReducedByIterations));
681  useNewtonIteration_ = true;
682  }
683  else if (control == "iterationcount") {
684  const int iterations = EWOMS_GET_PARAM(TypeTag, int, TimeStepControlTargetIterations); // 30
685  const double decayrate = EWOMS_GET_PARAM(TypeTag, double, TimeStepControlDecayRate); // 0.75
686  const double growthrate = EWOMS_GET_PARAM(TypeTag, double, TimeStepControlGrowthRate); // 1.25
687  timeStepControl_ = TimeStepControlType(new SimpleIterationCountTimeStepControl(iterations, decayrate, growthrate));
688  }
689  else if (control == "newtoniterationcount") {
690  const int iterations = EWOMS_GET_PARAM(TypeTag, int, TimeStepControlTargetNewtonIterations); // 8
691  const double decayrate = EWOMS_GET_PARAM(TypeTag, double, TimeStepControlDecayRate); // 0.75
692  const double growthrate = EWOMS_GET_PARAM(TypeTag, double, TimeStepControlGrowthRate); // 1.25
693  timeStepControl_ = TimeStepControlType(new SimpleIterationCountTimeStepControl(iterations, decayrate, growthrate));
694  useNewtonIteration_ = true;
695  }
696  else if (control == "hardcoded") {
697  const std::string filename = EWOMS_GET_PARAM(TypeTag, std::string, TimeStepControlFileName); // "timesteps"
698  timeStepControl_ = TimeStepControlType(new HardcodedTimeStepControl(filename));
699 
700  }
701  else
702  OPM_THROW(std::runtime_error,"Unsupported time step control selected "<< control);
703 
704  // make sure growth factor is something reasonable
705  assert(growthFactor_ >= 1.0);
706  }
707 
708  template <class ProblemType>
709  std::set<std::string> consistentlyFailingWells(const std::vector<ProblemType>& sr)
710  {
711  // If there are wells that cause repeated failures, we
712  // close them, and restart the un-chopped timestep.
713  std::ostringstream msg;
714  msg << " Excessive chopping detected in report step "
715  << sr.back().report_step << ", substep " << sr.back().current_step << "\n";
716 
717  std::set<std::string> failing_wells;
718 
719  // return empty set if no report exists
720  // well failures in assembly is not yet registred
721  if(sr.back().report.empty())
722  return failing_wells;
723 
724  const auto& wfs = sr.back().report.back().wellFailures();
725  for (const auto& wf : wfs) {
726  msg << " Well that failed: " << wf.wellName() << "\n";
727  }
728  msg.flush();
729  OpmLog::debug(msg.str());
730 
731  // Check the last few step reports.
732  const int num_steps = 3;
733  const int rep_step = sr.back().report_step;
734  const int sub_step = sr.back().current_step;
735  const int sr_size = sr.size();
736  if (sr_size >= num_steps) {
737  for (const auto& wf : wfs) {
738  failing_wells.insert(wf.wellName());
739  }
740  for (int s = 1; s < num_steps; ++s) {
741  const auto& srep = sr[sr_size - 1 - s];
742  // Report must be from same report step and substep, otherwise we have
743  // not chopped/retried enough times on this step.
744  if (srep.report_step != rep_step || srep.current_step != sub_step) {
745  break;
746  }
747  // Get the failing wells for this step, that also failed all other steps.
748  std::set<std::string> failing_wells_step;
749  for (const auto& wf : srep.report.back().wellFailures()) {
750  if (failing_wells.count(wf.wellName()) > 0) {
751  failing_wells_step.insert(wf.wellName());
752  }
753  }
754  failing_wells.swap(failing_wells_step);
755  }
756  }
757  return failing_wells;
758  }
759 
760  typedef std::unique_ptr<TimeStepControlInterface> TimeStepControlType;
761 
762  TimeStepControlType timeStepControl_;
763  double restartFactor_;
764  double growthFactor_;
765  double maxGrowth_;
766  double maxTimeStep_;
767  double minTimeStep_;
776  double minTimeStepBeforeShuttingProblematicWells_;
777  };
778 }
779 
780 #endif // OPM_ADAPTIVE_TIME_STEPPING_EBOS_HPP
Simulation timer for adaptive time stepping.
Definition: AdaptiveSimulatorTimer.hpp:41
bool done() const
Definition: AdaptiveSimulatorTimer.cpp:126
double currentStepLength() const
Definition: AdaptiveSimulatorTimer.cpp:108
void provideTimeStepEstimate(const double dt_estimate)
provide and estimate for new time step size
Definition: AdaptiveSimulatorTimer.cpp:72
void report(std::ostream &os) const
report start and end time as well as used steps so far
Definition: AdaptiveSimulatorTimer.cpp:153
double simulationTimeElapsed() const
Definition: AdaptiveSimulatorTimer.cpp:124
void setLastStepFailed(bool lastStepFailed)
tell the timestepper whether timestep failed or not
Definition: AdaptiveSimulatorTimer.hpp:104
Definition: AdaptiveTimeSteppingEbos.hpp:238
double maxTimeStep_
maximal allowed time step size in days
Definition: AdaptiveTimeSteppingEbos.hpp:766
double growthFactor_
factor to multiply time step when solver recovered from failed convergence
Definition: AdaptiveTimeSteppingEbos.hpp:764
double minTimeStep_
minimal allowed time step size before throwing
Definition: AdaptiveTimeSteppingEbos.hpp:767
bool fullTimestepInitially_
beginning with the size of the time step from data file
Definition: AdaptiveTimeSteppingEbos.hpp:773
double suggestedNextStep() const
Returns the simulator report for the failed substeps of the last report step.
Definition: AdaptiveTimeSteppingEbos.hpp:639
int solverRestartMax_
how many restart of solver are allowed
Definition: AdaptiveTimeSteppingEbos.hpp:769
double timestepAfterEvent_
suggested size of timestep after an event
Definition: AdaptiveTimeSteppingEbos.hpp:774
TimeStepControlType timeStepControl_
time step control object
Definition: AdaptiveTimeSteppingEbos.hpp:762
AdaptiveTimeSteppingEbos(const UnitSystem &unitSystem, const bool terminalOutput=true)
contructor taking parameter object
Definition: AdaptiveTimeSteppingEbos.hpp:266
bool timestepVerbose_
timestep verbosity
Definition: AdaptiveTimeSteppingEbos.hpp:771
bool ignoreConvergenceFailure_
continue instead of stop when minimum time step is reached
Definition: AdaptiveTimeSteppingEbos.hpp:768
double restartFactor_
factor to multiply time step with when solver fails to converge
Definition: AdaptiveTimeSteppingEbos.hpp:763
SimulatorReport step(const SimulatorTimer &simulatorTimer, Solver &solver, const bool isEvent, const std::vector< int > *fipnum=nullptr)
step method that acts like the solver::step method in a sub cycle of time steps
Definition: AdaptiveTimeSteppingEbos.hpp:371
bool solverVerbose_
solver verbosity
Definition: AdaptiveTimeSteppingEbos.hpp:770
AdaptiveTimeSteppingEbos(double max_next_tstep, const Tuning &tuning, const UnitSystem &unitSystem, const bool terminalOutput=true)
contructor taking parameter object
Definition: AdaptiveTimeSteppingEbos.hpp:293
bool useNewtonIteration_
use newton iteration count for adaptive time step control
Definition: AdaptiveTimeSteppingEbos.hpp:775
double suggestedNextTimestep_
suggested size of next timestep
Definition: AdaptiveTimeSteppingEbos.hpp:772
double maxGrowth_
factor that limits the maximum growth of a time step
Definition: AdaptiveTimeSteppingEbos.hpp:765
RelativeChangeInterface.
Definition: TimeStepControlInterface.hpp:32
Definition: SimulatorTimer.hpp:37
double currentStepLength() const override
Current step length.
Definition: SimulatorTimer.cpp:94
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27
Definition: AdaptiveTimeSteppingEbos.hpp:69
Definition: AdaptiveTimeSteppingEbos.hpp:65
Definition: AdaptiveTimeSteppingEbos.hpp:117
Definition: AdaptiveTimeSteppingEbos.hpp:113
Definition: AdaptiveTimeSteppingEbos.hpp:49
Definition: AdaptiveTimeSteppingEbos.hpp:33
Definition: AdaptiveTimeSteppingEbos.hpp:37
Definition: AdaptiveTimeSteppingEbos.hpp:53
Definition: AdaptiveTimeSteppingEbos.hpp:41
Definition: AdaptiveTimeSteppingEbos.hpp:45
Definition: AdaptiveTimeSteppingEbos.hpp:29
Definition: AdaptiveTimeSteppingEbos.hpp:57
Definition: AdaptiveTimeSteppingEbos.hpp:25
Definition: AdaptiveTimeSteppingEbos.hpp:73
Definition: AdaptiveTimeSteppingEbos.hpp:101
Definition: AdaptiveTimeSteppingEbos.hpp:93
Definition: AdaptiveTimeSteppingEbos.hpp:109
Definition: AdaptiveTimeSteppingEbos.hpp:105
Definition: AdaptiveTimeSteppingEbos.hpp:97
Definition: AdaptiveTimeSteppingEbos.hpp:85
Definition: AdaptiveTimeSteppingEbos.hpp:89
Definition: AdaptiveTimeSteppingEbos.hpp:81
Definition: AdaptiveTimeSteppingEbos.hpp:77
Definition: AdaptiveTimeSteppingEbos.hpp:61
A struct for returning timing data from a simulator to its caller.
Definition: SimulatorReport.hpp:31
void reportStep(std::ostringstream &os) const
Print a report suitable for a single simulation step.
Definition: SimulatorReport.cpp:88
Definition: SimulatorReport.hpp:69