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/ErrorMacros.hpp>
13 #include <opm/simulators/timestepping/SimulatorTimer.hpp>
14 #include <opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp>
15 #include <opm/simulators/timestepping/TimeStepControlInterface.hpp>
16 #include <opm/simulators/timestepping/TimeStepControl.hpp>
17 #include <opm/core/props/phaseUsageFromDeck.hpp>
18 #include <opm/input/eclipse/Schedule/ScheduleState.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(double max_next_tstep,
293  const Tuning& tuning,
294  const UnitSystem& unitSystem,
295  const bool terminalOutput = true)
296  : timeStepControl_()
297  , restartFactor_(tuning.TSFCNV)
298  , growthFactor_(tuning.TFDIFF)
299  , maxGrowth_(tuning.TSFMAX)
300  , maxTimeStep_(tuning.TSMAXZ) // 365.25
301  , minTimeStep_(tuning.TSFMIN) // 0.1;
303  , solverRestartMax_(EWOMS_GET_PARAM(TypeTag, int, SolverMaxRestarts)) // 10
304  , solverVerbose_(EWOMS_GET_PARAM(TypeTag, int, SolverVerbosity) > 0 && terminalOutput) // 2
305  , timestepVerbose_(EWOMS_GET_PARAM(TypeTag, int, TimeStepVerbosity) > 0 && terminalOutput) // 2
306  , suggestedNextTimestep_(max_next_tstep) // 1.0
307  , fullTimestepInitially_(EWOMS_GET_PARAM(TypeTag, bool, FullTimeStepInitially)) // false
308  , timestepAfterEvent_(tuning.TMAXWC) // 1e30
309  , useNewtonIteration_(false)
310  , minTimeStepBeforeShuttingProblematicWells_(EWOMS_GET_PARAM(TypeTag, double, MinTimeStepBeforeShuttingProblematicWellsInDays)*unit::day)
311  {
312  init_(unitSystem);
313  }
314 
315  static void registerParameters()
316  {
317  // TODO: make sure the help messages are correct (and useful)
318  EWOMS_REGISTER_PARAM(TypeTag, double, SolverRestartFactor,
319  "The factor time steps are elongated after restarts");
320  EWOMS_REGISTER_PARAM(TypeTag, double, SolverGrowthFactor,
321  "The factor time steps are elongated after a successful substep");
322  EWOMS_REGISTER_PARAM(TypeTag, double, SolverMaxGrowth,
323  "The maximum factor time steps are elongated after a report step");
324  EWOMS_REGISTER_PARAM(TypeTag, double, SolverMaxTimeStepInDays,
325  "The maximum size of a time step in days");
326  EWOMS_REGISTER_PARAM(TypeTag, double, SolverMinTimeStep,
327  "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");
328  EWOMS_REGISTER_PARAM(TypeTag, bool, SolverContinueOnConvergenceFailure,
329  "Continue instead of stop when minimum solver time step is reached");
330  EWOMS_REGISTER_PARAM(TypeTag, int, SolverMaxRestarts,
331  "The maximum number of breakdowns before a substep is given up and the simulator is terminated");
332  EWOMS_REGISTER_PARAM(TypeTag, int, SolverVerbosity,
333  "Specify the \"chattiness\" of the non-linear solver itself");
334  EWOMS_REGISTER_PARAM(TypeTag, int, TimeStepVerbosity,
335  "Specify the \"chattiness\" during the time integration");
336  EWOMS_REGISTER_PARAM(TypeTag, double, InitialTimeStepInDays,
337  "The size of the initial time step in days");
338  EWOMS_REGISTER_PARAM(TypeTag, bool, FullTimeStepInitially,
339  "Always attempt to finish a report step using a single substep");
340  EWOMS_REGISTER_PARAM(TypeTag, double, TimeStepAfterEventInDays,
341  "Time step size of the first time step after an event occurs during the simulation in days");
342  EWOMS_REGISTER_PARAM(TypeTag, std::string, TimeStepControl,
343  "The algorithm used to determine time-step sizes. valid options are: 'pid' (default), 'pid+iteration', 'pid+newtoniteration', 'iterationcount', 'newtoniterationcount' and 'hardcoded'");
344  EWOMS_REGISTER_PARAM(TypeTag, double, TimeStepControlTolerance,
345  "The tolerance used by the time step size control algorithm");
346  EWOMS_REGISTER_PARAM(TypeTag, int, TimeStepControlTargetIterations,
347  "The number of linear iterations which the time step control scheme should aim for (if applicable)");
348  EWOMS_REGISTER_PARAM(TypeTag, int, TimeStepControlTargetNewtonIterations,
349  "The number of Newton iterations which the time step control scheme should aim for (if applicable)");
350  EWOMS_REGISTER_PARAM(TypeTag, double, TimeStepControlDecayRate,
351  "The decay rate of the time step size of the number of target iterations is exceeded");
352  EWOMS_REGISTER_PARAM(TypeTag, double, TimeStepControlGrowthRate,
353  "The growth rate of the time step size of the number of target iterations is undercut");
354  EWOMS_REGISTER_PARAM(TypeTag, double, TimeStepControlDecayDampingFactor,
355  "The decay rate of the time step decrease when the target iterations is exceeded");
356  EWOMS_REGISTER_PARAM(TypeTag, double, TimeStepControlGrowthDampingFactor,
357  "The growth rate of the time step increase when the target iterations is undercut");
358  EWOMS_REGISTER_PARAM(TypeTag, std::string, TimeStepControlFileName,
359  "The name of the file which contains the hardcoded time steps sizes");
360  EWOMS_REGISTER_PARAM(TypeTag, double, MinTimeStepBeforeShuttingProblematicWellsInDays,
361  "The minimum time step size in days for which problematic wells are not shut");
362  EWOMS_REGISTER_PARAM(TypeTag, double, MinTimeStepBasedOnNewtonIterations,
363  "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");
364  }
365 
369  template <class Solver>
370  SimulatorReport step(const SimulatorTimer& simulatorTimer,
371  Solver& solver,
372  const bool isEvent,
373  const std::vector<int>* fipnum = nullptr)
374  {
375  SimulatorReport report;
376  const double timestep = simulatorTimer.currentStepLength();
377 
378  // init last time step as a fraction of the given time step
379  if (suggestedNextTimestep_ < 0) {
381  }
382 
384  suggestedNextTimestep_ = timestep;
385  }
386 
387  // use seperate time step after event
388  if (isEvent && timestepAfterEvent_ > 0) {
390  }
391 
392  auto& ebosSimulator = solver.model().ebosSimulator();
393  auto& ebosProblem = ebosSimulator.problem();
394 
395  // create adaptive step timer with previously used sub step size
396  AdaptiveSimulatorTimer substepTimer(simulatorTimer, suggestedNextTimestep_, maxTimeStep_);
397 
398  // counter for solver restarts
399  int restarts = 0;
400 
401  // sub step time loop
402  while (!substepTimer.done()) {
403  // get current delta t
404  const double dt = substepTimer.currentStepLength() ;
405  if (timestepVerbose_) {
406  logTimer(substepTimer);
407  }
408 
409  SimulatorReportSingle substepReport;
410  std::string causeOfFailure = "";
411  try {
412  substepReport = solver.step(substepTimer);
413  if (solverVerbose_) {
414  // report number of linear iterations
415  OpmLog::debug("Overall linear iterations used: " + std::to_string(substepReport.total_linear_iterations));
416  }
417  }
418  catch (const TooManyIterations& e) {
419  substepReport = solver.failureReport();
420  causeOfFailure = "Solver convergence failure - Iteration limit reached";
421 
422  logException_(e, solverVerbose_);
423  // since linearIterations is < 0 this will restart the solver
424  }
425  catch (const LinearSolverProblem& e) {
426  substepReport = solver.failureReport();
427  causeOfFailure = "Linear solver convergence failure";
428 
429  logException_(e, solverVerbose_);
430  // since linearIterations is < 0 this will restart the solver
431  }
432  catch (const NumericalIssue& e) {
433  substepReport = solver.failureReport();
434  causeOfFailure = "Solver convergence failure - Numerical problem encountered";
435 
436  logException_(e, solverVerbose_);
437  // since linearIterations is < 0 this will restart the solver
438  }
439  catch (const std::runtime_error& e) {
440  substepReport = solver.failureReport();
441 
442  logException_(e, solverVerbose_);
443  // also catch linear solver not converged
444  }
445  catch (const Dune::ISTLError& e) {
446  substepReport = solver.failureReport();
447 
448  logException_(e, solverVerbose_);
449  // also catch errors in ISTL AMG that occur when time step is too large
450  }
451  catch (const Dune::MatrixBlockError& e) {
452  substepReport = solver.failureReport();
453 
454  logException_(e, solverVerbose_);
455  // this can be thrown by ISTL's ILU0 in block mode, yet is not an ISTLError
456  }
457 
458  //Pass substep to eclwriter for summary output
459  ebosSimulator.problem().setSubStepReport(substepReport);
460 
461  report += substepReport;
462 
463  bool continue_on_uncoverged_solution = ignoreConvergenceFailure_ && !substepReport.converged && dt <= minTimeStep_;
464 
465  if (continue_on_uncoverged_solution) {
466  const auto msg = std::string("Solver failed to converge but timestep ")
467  + std::to_string(dt) + " is smaller or equal to "
468  + std::to_string(minTimeStep_) + "\n which is the minimum threshold given"
469  + "by option --solver-min-time-step= \n";
470  if (solverVerbose_) {
471  OpmLog::error(msg);
472  }
473  }
474 
475  if (substepReport.converged || continue_on_uncoverged_solution) {
476 
477  // advance by current dt
478  ++substepTimer;
479 
480  // create object to compute the time error, simply forwards the call to the model
481  SolutionTimeErrorSolverWrapperEbos<Solver> relativeChange(solver);
482 
483  // compute new time step estimate
484  const int iterations = useNewtonIteration_ ? substepReport.total_newton_iterations
485  : substepReport.total_linear_iterations;
486  double dtEstimate = timeStepControl_->computeTimeStepSize(dt, iterations, relativeChange,
487  substepTimer.simulationTimeElapsed());
488 
489  assert(dtEstimate > 0);
490  // limit the growth of the timestep size by the growth factor
491  dtEstimate = std::min(dtEstimate, double(maxGrowth_ * dt));
492  assert(dtEstimate > 0);
493  // further restrict time step size growth after convergence problems
494  if (restarts > 0) {
495  dtEstimate = std::min(growthFactor_ * dt, dtEstimate);
496  // solver converged, reset restarts counter
497  restarts = 0;
498  }
499 
500  // Further restrict time step size if we are in
501  // prediction mode with THP constraints.
502  if (solver.model().wellModel().hasTHPConstraints()) {
503  const double maxPredictionTHPTimestep = 16.0 * unit::day;
504  dtEstimate = std::min(dtEstimate, maxPredictionTHPTimestep);
505  }
506  assert(dtEstimate > 0);
507  if (timestepVerbose_) {
508  std::ostringstream ss;
509  substepReport.reportStep(ss);
510  OpmLog::info(ss.str());
511  }
512 
513  // write data if outputWriter was provided
514  // if the time step is done we do not need
515  // to write it as this will be done by the simulator
516  // anyway.
517  if (!substepTimer.done()) {
518  if (fipnum) {
519  solver.computeFluidInPlace(*fipnum);
520  }
521  time::StopWatch perfTimer;
522  perfTimer.start();
523 
524  ebosProblem.writeOutput();
525 
526  report.success.output_write_time += perfTimer.secsSinceStart();
527  }
528 
529  // set new time step length
530  substepTimer.provideTimeStepEstimate(dtEstimate);
531 
532  report.success.converged = substepTimer.done();
533  substepTimer.setLastStepFailed(false);
534 
535  }
536  else { // in case of no convergence
537  substepTimer.setLastStepFailed(true);
538 
539  // If we have restarted (i.e. cut the timestep) too
540  // many times, we have failed and throw an exception.
541  if (restarts >= solverRestartMax_) {
542  const auto msg = std::string("Solver failed to converge after cutting timestep ")
543  + std::to_string(restarts) + " times.";
544  if (solverVerbose_) {
545  OpmLog::error(msg);
546  }
547  OPM_THROW_NOLOG(NumericalIssue, msg);
548  }
549 
550  // The new, chopped timestep.
551  const double newTimeStep = restartFactor_ * dt;
552 
553 
554  // If we have restarted (i.e. cut the timestep) too
555  // much, we have failed and throw an exception.
556  if (newTimeStep < minTimeStep_) {
557  const auto msg = std::string("Solver failed to converge after cutting timestep to ")
558  + std::to_string(minTimeStep_) + "\n which is the minimum threshold given"
559  + "by option --solver-min-time-step= \n";
560  if (solverVerbose_) {
561  OpmLog::error(msg);
562  }
563  OPM_THROW_NOLOG(NumericalIssue, msg);
564  }
565 
566  // Define utility function for chopping timestep.
567  auto chopTimestep = [&]() {
568  substepTimer.provideTimeStepEstimate(newTimeStep);
569  if (solverVerbose_) {
570  std::string msg;
571  msg = causeOfFailure + "\nTimestep chopped to "
572  + std::to_string(unit::convert::to(substepTimer.currentStepLength(), unit::day)) + " days\n";
573  OpmLog::problem(msg);
574  }
575  ++restarts;
576  };
577 
578  const double minimumChoppedTimestep = minTimeStepBeforeShuttingProblematicWells_;
579  if (newTimeStep > minimumChoppedTimestep) {
580  chopTimestep();
581  } else {
582  // We are below the threshold, and will check if there are any
583  // wells we should close rather than chopping again.
584  std::set<std::string> failing_wells = consistentlyFailingWells(solver.model().stepReports());
585  if (failing_wells.empty()) {
586  // Found no wells to close, chop the timestep as above.
587  chopTimestep();
588  } else {
589  // Close all consistently failing wells.
590  int num_shut_wells = 0;
591  for (const auto& well : failing_wells) {
592  bool was_shut = solver.model().wellModel().forceShutWellByName(well, substepTimer.simulationTimeElapsed());
593  if (was_shut) {
594  ++num_shut_wells;
595  }
596  }
597  if (num_shut_wells == 0) {
598  // None of the problematic wells were shut.
599  // We must fall back to chopping again.
600  chopTimestep();
601  } else {
602  substepTimer.provideTimeStepEstimate(dt);
603  if (solverVerbose_) {
604  std::string msg;
605  msg = "\nProblematic well(s) were shut: ";
606  for (const auto& well : failing_wells) {
607  msg += well;
608  msg += " ";
609  }
610  msg += "(retrying timestep)\n";
611  OpmLog::problem(msg);
612  }
613  }
614  }
615  }
616  }
617  ebosProblem.setNextTimeStepSize(substepTimer.currentStepLength());
618  }
619 
620  // store estimated time step for next reportStep
621  suggestedNextTimestep_ = substepTimer.currentStepLength();
622  if (timestepVerbose_) {
623  std::ostringstream ss;
624  substepTimer.report(ss);
625  ss << "Suggested next step size = " << unit::convert::to(suggestedNextTimestep_, unit::day) << " (days)" << std::endl;
626  OpmLog::debug(ss.str());
627  }
628 
629  if (! std::isfinite(suggestedNextTimestep_)) { // check for NaN
630  suggestedNextTimestep_ = timestep;
631  }
632  return report;
633  }
634 
638  double suggestedNextStep() const
639  { return suggestedNextTimestep_; }
640 
641  void setSuggestedNextStep(const double x)
642  { suggestedNextTimestep_ = x; }
643 
644  void updateTUNING(double max_next_tstep, const Tuning& tuning)
645  {
646  restartFactor_ = tuning.TSFCNV;
647  growthFactor_ = tuning.TFDIFF;
648  maxGrowth_ = tuning.TSFMAX;
649  maxTimeStep_ = tuning.TSMAXZ;
650  suggestedNextTimestep_ = max_next_tstep;
651  timestepAfterEvent_ = tuning.TMAXWC;
652  }
653 
654 
655  protected:
656  void init_(const UnitSystem& unitSystem)
657  {
658  // valid are "pid" and "pid+iteration"
659  std::string control = EWOMS_GET_PARAM(TypeTag, std::string, TimeStepControl); // "pid"
660 
661  const double tol = EWOMS_GET_PARAM(TypeTag, double, TimeStepControlTolerance); // 1e-1
662  if (control == "pid") {
663  timeStepControl_ = TimeStepControlType(new PIDTimeStepControl(tol));
664  }
665  else if (control == "pid+iteration") {
666  const int iterations = EWOMS_GET_PARAM(TypeTag, int, TimeStepControlTargetIterations); // 30
667  const double decayDampingFactor = EWOMS_GET_PARAM(TypeTag, double, TimeStepControlDecayDampingFactor); // 1.0
668  const double growthDampingFactor = EWOMS_GET_PARAM(TypeTag, double, TimeStepControlGrowthDampingFactor); // 3.2
669  timeStepControl_ = TimeStepControlType(new PIDAndIterationCountTimeStepControl(iterations, decayDampingFactor, growthDampingFactor, tol));
670  }
671  else if (control == "pid+newtoniteration") {
672  const int iterations = EWOMS_GET_PARAM(TypeTag, int, TimeStepControlTargetNewtonIterations); // 8
673  const double decayDampingFactor = EWOMS_GET_PARAM(TypeTag, double, TimeStepControlDecayDampingFactor); // 1.0
674  const double growthDampingFactor = EWOMS_GET_PARAM(TypeTag, double, TimeStepControlGrowthDampingFactor); // 3.2
675  const double nonDimensionalMinTimeStepIterations = EWOMS_GET_PARAM(TypeTag, double, MinTimeStepBasedOnNewtonIterations); // 0.0 by default
676  // the min time step can be reduced by the newton iteration numbers
677  double minTimeStepReducedByIterations = unitSystem.to_si(UnitSystem::measure::time, nonDimensionalMinTimeStepIterations);
678  timeStepControl_ = TimeStepControlType(new PIDAndIterationCountTimeStepControl(iterations, decayDampingFactor,
679  growthDampingFactor, tol, minTimeStepReducedByIterations));
680  useNewtonIteration_ = true;
681  }
682  else if (control == "iterationcount") {
683  const int iterations = EWOMS_GET_PARAM(TypeTag, int, TimeStepControlTargetIterations); // 30
684  const double decayrate = EWOMS_GET_PARAM(TypeTag, double, TimeStepControlDecayRate); // 0.75
685  const double growthrate = EWOMS_GET_PARAM(TypeTag, double, TimeStepControlGrowthRate); // 1.25
686  timeStepControl_ = TimeStepControlType(new SimpleIterationCountTimeStepControl(iterations, decayrate, growthrate));
687  }
688  else if (control == "newtoniterationcount") {
689  const int iterations = EWOMS_GET_PARAM(TypeTag, int, TimeStepControlTargetNewtonIterations); // 8
690  const double decayrate = EWOMS_GET_PARAM(TypeTag, double, TimeStepControlDecayRate); // 0.75
691  const double growthrate = EWOMS_GET_PARAM(TypeTag, double, TimeStepControlGrowthRate); // 1.25
692  timeStepControl_ = TimeStepControlType(new SimpleIterationCountTimeStepControl(iterations, decayrate, growthrate));
693  useNewtonIteration_ = true;
694  }
695  else if (control == "hardcoded") {
696  const std::string filename = EWOMS_GET_PARAM(TypeTag, std::string, TimeStepControlFileName); // "timesteps"
697  timeStepControl_ = TimeStepControlType(new HardcodedTimeStepControl(filename));
698 
699  }
700  else
701  OPM_THROW(std::runtime_error,"Unsupported time step control selected "<< control);
702 
703  // make sure growth factor is something reasonable
704  assert(growthFactor_ >= 1.0);
705  }
706 
707  template <class ProblemType>
708  std::set<std::string> consistentlyFailingWells(const std::vector<ProblemType>& sr)
709  {
710  // If there are wells that cause repeated failures, we
711  // close them, and restart the un-chopped timestep.
712  std::ostringstream msg;
713  msg << " Excessive chopping detected in report step "
714  << sr.back().report_step << ", substep " << sr.back().current_step << "\n";
715 
716  std::set<std::string> failing_wells;
717 
718  // return empty set if no report exists
719  // well failures in assembly is not yet registred
720  if(sr.back().report.empty())
721  return failing_wells;
722 
723  const auto& wfs = sr.back().report.back().wellFailures();
724  for (const auto& wf : wfs) {
725  msg << " Well that failed: " << wf.wellName() << "\n";
726  }
727  msg.flush();
728  OpmLog::debug(msg.str());
729 
730  // Check the last few step reports.
731  const int num_steps = 3;
732  const int rep_step = sr.back().report_step;
733  const int sub_step = sr.back().current_step;
734  const int sr_size = sr.size();
735  if (sr_size >= num_steps) {
736  for (const auto& wf : wfs) {
737  failing_wells.insert(wf.wellName());
738  }
739  for (int s = 1; s < num_steps; ++s) {
740  const auto& srep = sr[sr_size - 1 - s];
741  // Report must be from same report step and substep, otherwise we have
742  // not chopped/retried enough times on this step.
743  if (srep.report_step != rep_step || srep.current_step != sub_step) {
744  break;
745  }
746  // Get the failing wells for this step, that also failed all other steps.
747  std::set<std::string> failing_wells_step;
748  for (const auto& wf : srep.report.back().wellFailures()) {
749  if (failing_wells.count(wf.wellName()) > 0) {
750  failing_wells_step.insert(wf.wellName());
751  }
752  }
753  failing_wells.swap(failing_wells_step);
754  }
755  }
756  return failing_wells;
757  }
758 
759  typedef std::unique_ptr<TimeStepControlInterface> TimeStepControlType;
760 
761  TimeStepControlType timeStepControl_;
762  double restartFactor_;
763  double growthFactor_;
764  double maxGrowth_;
765  double maxTimeStep_;
766  double minTimeStep_;
775  double minTimeStepBeforeShuttingProblematicWells_;
776  };
777 }
778 
779 #endif // OPM_ADAPTIVE_TIME_STEPPING_EBOS_HPP
Simulation timer for adaptive time stepping.
Definition: AdaptiveSimulatorTimer.hpp:41
bool done() const
Definition: AdaptiveSimulatorTimer.cpp:130
double currentStepLength() const
Definition: AdaptiveSimulatorTimer.cpp:112
void provideTimeStepEstimate(const double dt_estimate)
provide and estimate for new time step size
Definition: AdaptiveSimulatorTimer.cpp:76
void report(std::ostream &os) const
report start and end time as well as used steps so far
Definition: AdaptiveSimulatorTimer.cpp:157
double simulationTimeElapsed() const
Definition: AdaptiveSimulatorTimer.cpp:128
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:765
double growthFactor_
factor to multiply time step when solver recovered from failed convergence
Definition: AdaptiveTimeSteppingEbos.hpp:763
double minTimeStep_
minimal allowed time step size before throwing
Definition: AdaptiveTimeSteppingEbos.hpp:766
bool fullTimestepInitially_
beginning with the size of the time step from data file
Definition: AdaptiveTimeSteppingEbos.hpp:772
double suggestedNextStep() const
Returns the simulator report for the failed substeps of the last report step.
Definition: AdaptiveTimeSteppingEbos.hpp:638
int solverRestartMax_
how many restart of solver are allowed
Definition: AdaptiveTimeSteppingEbos.hpp:768
double timestepAfterEvent_
suggested size of timestep after an event
Definition: AdaptiveTimeSteppingEbos.hpp:773
TimeStepControlType timeStepControl_
time step control object
Definition: AdaptiveTimeSteppingEbos.hpp:761
AdaptiveTimeSteppingEbos(const UnitSystem &unitSystem, const bool terminalOutput=true)
contructor taking parameter object
Definition: AdaptiveTimeSteppingEbos.hpp:265
bool timestepVerbose_
timestep verbosity
Definition: AdaptiveTimeSteppingEbos.hpp:770
bool ignoreConvergenceFailure_
continue instead of stop when minimum time step is reached
Definition: AdaptiveTimeSteppingEbos.hpp:767
double restartFactor_
factor to multiply time step with when solver fails to converge
Definition: AdaptiveTimeSteppingEbos.hpp:762
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:370
bool solverVerbose_
solver verbosity
Definition: AdaptiveTimeSteppingEbos.hpp:769
AdaptiveTimeSteppingEbos(double max_next_tstep, const Tuning &tuning, const UnitSystem &unitSystem, const bool terminalOutput=true)
contructor taking parameter object
Definition: AdaptiveTimeSteppingEbos.hpp:292
bool useNewtonIteration_
use newton iteration count for adaptive time step control
Definition: AdaptiveTimeSteppingEbos.hpp:774
double suggestedNextTimestep_
suggested size of next timestep
Definition: AdaptiveTimeSteppingEbos.hpp:771
double maxGrowth_
factor that limits the maximum growth of a time step
Definition: AdaptiveTimeSteppingEbos.hpp:764
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: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:89
Definition: SimulatorReport.hpp:70