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