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