My Project
Loading...
Searching...
No Matches
AdaptiveTimeStepping.hpp
1/*
2*/
3#ifndef OPM_ADAPTIVE_TIME_STEPPING_HPP
4#define OPM_ADAPTIVE_TIME_STEPPING_HPP
5
6#include <dune/common/version.hh>
7#if DUNE_VERSION_GTE(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/grid/utility/StopWatch.hpp>
18
19#include <opm/input/eclipse/Units/Units.hpp>
20#include <opm/input/eclipse/Units/UnitSystem.hpp>
21
22#include <opm/input/eclipse/Schedule/Tuning.hpp>
23
27
28#include <opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp>
29#include <opm/simulators/timestepping/EclTimeSteppingParams.hpp>
30#include <opm/simulators/timestepping/SimulatorReport.hpp>
31#include <opm/simulators/timestepping/SimulatorTimer.hpp>
32#include <opm/simulators/timestepping/TimeStepControl.hpp>
33#include <opm/simulators/timestepping/TimeStepControlInterface.hpp>
34
35#include <opm/simulators/utils/phaseUsageFromDeck.hpp>
36
37#include <fmt/format.h>
38
39#include <algorithm>
40#include <cassert>
41#include <cmath>
42#include <functional>
43#include <memory>
44#include <set>
45#include <sstream>
46#include <stdexcept>
47#include <string>
48#include <vector>
49
50namespace Opm::Parameters {
51
52struct SolverContinueOnConvergenceFailure { static constexpr bool value = false; };
53struct SolverMaxRestarts { static constexpr int value = 10; };
54struct SolverVerbosity { static constexpr int value = 1; };
55struct TimeStepVerbosity { static constexpr int value = 1; };
56struct InitialTimeStepInDays { static constexpr double value = 1.0; };
57struct FullTimeStepInitially { static constexpr bool value = false; };
58struct TimeStepControl { static constexpr auto value = "pid+newtoniteration"; };
59struct TimeStepControlTolerance { static constexpr double value = 1e-1; };
60struct TimeStepControlTargetIterations { static constexpr int value = 30; };
61struct TimeStepControlTargetNewtonIterations { static constexpr int value = 8; };
62struct TimeStepControlDecayRate { static constexpr double value = 0.75; };
63struct TimeStepControlGrowthRate { static constexpr double value = 1.25; };
64struct TimeStepControlDecayDampingFactor { static constexpr double value = 1.0; };
65struct TimeStepControlGrowthDampingFactor { static constexpr double value = 3.2; };
66struct TimeStepControlFileName { static constexpr auto value = "timesteps"; };
67struct MinTimeStepBeforeShuttingProblematicWellsInDays { static constexpr double value = 0.01; };
68struct MinTimeStepBasedOnNewtonIterations { static constexpr double value = 0.0; };
69
70} // namespace Opm::Parameters
71
72namespace Opm {
73
74struct StepReport;
75
76namespace detail {
77
78void logTimer(const AdaptiveSimulatorTimer& substepTimer);
79
80std::set<std::string> consistentlyFailingWells(const std::vector<StepReport>& sr);
81
82void registerAdaptiveParameters();
83
84}
85
86 // AdaptiveTimeStepping
87 //---------------------
88 template<class TypeTag>
90 {
92 template <class Solver>
93 class SolutionTimeErrorSolverWrapper : public RelativeChangeInterface
94 {
95 const Solver& solver_;
96 public:
97 SolutionTimeErrorSolverWrapper(const Solver& solver)
98 : solver_(solver)
99 {}
100
102 double relativeChange() const
103 { return solver_.model().relativeChange(); }
104 };
105
106 template<class E>
107 void logException_(const E& exception, bool verbose)
108 {
109 if (verbose) {
110 std::string message;
111 message = "Caught Exception: ";
112 message += exception.what();
113 OpmLog::debug(message);
114 }
115 }
116
117 public:
118 AdaptiveTimeStepping() = default;
119
121 AdaptiveTimeStepping(const UnitSystem& unitSystem,
122 const double max_next_tstep = -1.0,
123 const bool terminalOutput = true)
125 , restartFactor_(Parameters::Get<Parameters::SolverRestartFactor<Scalar>>()) // 0.33
126 , growthFactor_(Parameters::Get<Parameters::SolverGrowthFactor<Scalar>>()) // 2.0
127 , maxGrowth_(Parameters::Get<Parameters::SolverMaxGrowth<Scalar>>()) // 3.0
128 , maxTimeStep_(Parameters::Get<Parameters::SolverMaxTimeStepInDays<Scalar>>() * 24 * 60 * 60) // 365.25
129 , minTimeStep_(unitSystem.to_si(UnitSystem::measure::time, Parameters::Get<Parameters::SolverMinTimeStep<Scalar>>())) // 1e-12;
130 , ignoreConvergenceFailure_(Parameters::Get<Parameters::SolverContinueOnConvergenceFailure>()) // false;
131 , solverRestartMax_(Parameters::Get<Parameters::SolverMaxRestarts>()) // 10
132 , solverVerbose_(Parameters::Get<Parameters::SolverVerbosity>() > 0 && terminalOutput) // 2
133 , timestepVerbose_(Parameters::Get<Parameters::TimeStepVerbosity>() > 0 && terminalOutput) // 2
134 , suggestedNextTimestep_((max_next_tstep <= 0 ? Parameters::Get<Parameters::InitialTimeStepInDays>() : max_next_tstep) * 24 * 60 * 60) // 1.0
135 , fullTimestepInitially_(Parameters::Get<Parameters::FullTimeStepInitially>()) // false
136 , timestepAfterEvent_(Parameters::Get<Parameters::TimeStepAfterEventInDays<Scalar>>() * 24 * 60 * 60) // 1e30
137 , useNewtonIteration_(false)
138 , minTimeStepBeforeShuttingProblematicWells_(Parameters::Get<Parameters::MinTimeStepBeforeShuttingProblematicWellsInDays>() * unit::day)
139
140 {
141 init_(unitSystem);
142 }
143
144
145
149 AdaptiveTimeStepping(double max_next_tstep,
150 const Tuning& tuning,
151 const UnitSystem& unitSystem,
152 const bool terminalOutput = true)
154 , restartFactor_(tuning.TSFCNV)
155 , growthFactor_(tuning.TFDIFF)
156 , maxGrowth_(tuning.TSFMAX)
157 , maxTimeStep_(tuning.TSMAXZ) // 365.25
158 , minTimeStep_(tuning.TSFMIN) // 0.1;
160 , solverRestartMax_(Parameters::Get<Parameters::SolverMaxRestarts>()) // 10
161 , solverVerbose_(Parameters::Get<Parameters::SolverVerbosity>() > 0 && terminalOutput) // 2
162 , timestepVerbose_(Parameters::Get<Parameters::TimeStepVerbosity>() > 0 && terminalOutput) // 2
163 , suggestedNextTimestep_(max_next_tstep <= 0 ? Parameters::Get<Parameters::InitialTimeStepInDays>() * 24 * 60 * 60 : max_next_tstep) // 1.0
164 , fullTimestepInitially_(Parameters::Get<Parameters::FullTimeStepInitially>()) // false
165 , timestepAfterEvent_(tuning.TMAXWC) // 1e30
166 , useNewtonIteration_(false)
167 , minTimeStepBeforeShuttingProblematicWells_(Parameters::Get<Parameters::MinTimeStepBeforeShuttingProblematicWellsInDays>() * unit::day)
168 {
169 init_(unitSystem);
170 }
171
172 static void registerParameters()
173 {
174 registerEclTimeSteppingParameters<Scalar>();
175 detail::registerAdaptiveParameters();
176 }
177
183 template <class Solver>
184 SimulatorReport step(const SimulatorTimer& simulatorTimer,
185 Solver& solver,
186 const bool isEvent,
187 const std::vector<int>* fipnum = nullptr,
188 const std::function<bool()> tuningUpdater = [](){return false;})
189 {
190 // Maybe update tuning
191 tuningUpdater();
192 SimulatorReport report;
193 const double timestep = simulatorTimer.currentStepLength();
194
195 // init last time step as a fraction of the given time step
196 if (suggestedNextTimestep_ < 0) {
198 }
199
201 suggestedNextTimestep_ = timestep;
202 }
203
204 // use seperate time step after event
205 if (isEvent && timestepAfterEvent_ > 0) {
207 }
208
209 auto& simulator = solver.model().simulator();
210 auto& problem = simulator.problem();
211
212 // create adaptive step timer with previously used sub step size
213 AdaptiveSimulatorTimer substepTimer(simulatorTimer, suggestedNextTimestep_, maxTimeStep_);
214
215 // counter for solver restarts
216 int restarts = 0;
217
218 // sub step time loop
219 while (!substepTimer.done()) {
220 // Maybe update tuning
221 // get current delta t
222 auto oldValue = suggestedNextTimestep_;
223 if (tuningUpdater()) {
224 // Use provideTimeStepEstimate to make we sure don't simulate longer than the report step is.
225 substepTimer.provideTimeStepEstimate(suggestedNextTimestep_);
226 suggestedNextTimestep_ = oldValue;
227 }
228 const double dt = substepTimer.currentStepLength();
229 if (timestepVerbose_) {
230 detail::logTimer(substepTimer);
231 }
232
233 SimulatorReportSingle substepReport;
234 std::string causeOfFailure;
235 try {
236 substepReport = solver.step(substepTimer);
237
238 if (solverVerbose_) {
239 // report number of linear iterations
240 OpmLog::debug("Overall linear iterations used: " + std::to_string(substepReport.total_linear_iterations));
241 }
242 }
243 catch (const TooManyIterations& e) {
244 substepReport = solver.failureReport();
245 causeOfFailure = "Solver convergence failure - Iteration limit reached";
246
247 logException_(e, solverVerbose_);
248 // since linearIterations is < 0 this will restart the solver
249 }
250 catch (const ConvergenceMonitorFailure& e) {
251 substepReport = solver.failureReport();
252 causeOfFailure = "Convergence monitor failure";
253 }
254 catch (const LinearSolverProblem& e) {
255 substepReport = solver.failureReport();
256 causeOfFailure = "Linear solver convergence failure";
257
258 logException_(e, solverVerbose_);
259 // since linearIterations is < 0 this will restart the solver
260 }
261 catch (const NumericalProblem& e) {
262 substepReport = solver.failureReport();
263 causeOfFailure = "Solver convergence failure - Numerical problem encountered";
264
265 logException_(e, solverVerbose_);
266 // since linearIterations is < 0 this will restart the solver
267 }
268 catch (const std::runtime_error& e) {
269 substepReport = solver.failureReport();
270
271 logException_(e, solverVerbose_);
272 // also catch linear solver not converged
273 }
274 catch (const Dune::ISTLError& e) {
275 substepReport = solver.failureReport();
276
277 logException_(e, solverVerbose_);
278 // also catch errors in ISTL AMG that occur when time step is too large
279 }
280 catch (const Dune::MatrixBlockError& e) {
281 substepReport = solver.failureReport();
282
283 logException_(e, solverVerbose_);
284 // this can be thrown by ISTL's ILU0 in block mode, yet is not an ISTLError
285 }
286
287 //Pass substep to eclwriter for summary output
288 simulator.problem().setSubStepReport(substepReport);
289
290 report += substepReport;
291
292 bool continue_on_uncoverged_solution = ignoreConvergenceFailure_ &&
293 !substepReport.converged &&
294 dt <= minTimeStep_;
295
296 if (continue_on_uncoverged_solution && solverVerbose_) {
297 const auto msg = fmt::format(
298 "Solver failed to converge but timestep "
299 "{} is smaller or equal to {}\n"
300 "which is the minimum threshold given "
301 "by option --solver-min-time-step\n",
302 dt, minTimeStep_
303 );
304 OpmLog::problem(msg);
305 }
306
307 if (substepReport.converged || continue_on_uncoverged_solution) {
308
309 // advance by current dt
310 ++substepTimer;
311
312 // create object to compute the time error, simply forwards the call to the model
313 SolutionTimeErrorSolverWrapper<Solver> relativeChange(solver);
314
315 // compute new time step estimate
316 const int iterations = useNewtonIteration_ ? substepReport.total_newton_iterations
317 : substepReport.total_linear_iterations;
318 double dtEstimate = timeStepControl_->computeTimeStepSize(dt, iterations, relativeChange,
319 substepTimer.simulationTimeElapsed());
320
321 assert(dtEstimate > 0);
322 // limit the growth of the timestep size by the growth factor
323 dtEstimate = std::min(dtEstimate, double(maxGrowth_ * dt));
324 assert(dtEstimate > 0);
325 // further restrict time step size growth after convergence problems
326 if (restarts > 0) {
327 dtEstimate = std::min(growthFactor_ * dt, dtEstimate);
328 // solver converged, reset restarts counter
329 restarts = 0;
330 }
331
332 // Further restrict time step size if we are in
333 // prediction mode with THP constraints.
334 if (solver.model().wellModel().hasTHPConstraints()) {
335 const double maxPredictionTHPTimestep = 16.0 * unit::day;
336 dtEstimate = std::min(dtEstimate, maxPredictionTHPTimestep);
337 }
338 assert(dtEstimate > 0);
339 if (timestepVerbose_) {
340 std::ostringstream ss;
341 substepReport.reportStep(ss);
342 OpmLog::info(ss.str());
343 }
344
345 // write data if outputWriter was provided
346 // if the time step is done we do not need
347 // to write it as this will be done by the simulator
348 // anyway.
349 if (!substepTimer.done()) {
350 if (fipnum) {
351 solver.computeFluidInPlace(*fipnum);
352 }
353 time::StopWatch perfTimer;
354 perfTimer.start();
355
356 problem.writeOutput();
357
358 report.success.output_write_time += perfTimer.secsSinceStart();
359 }
360
361 // set new time step length
362 substepTimer.provideTimeStepEstimate(dtEstimate);
363
364 report.success.converged = substepTimer.done();
365 substepTimer.setLastStepFailed(false);
366
367 }
368 else { // in case of no convergence
369 substepTimer.setLastStepFailed(true);
370
371 // If we have restarted (i.e. cut the timestep) too
372 // many times, we have failed and throw an exception.
373 if (restarts >= solverRestartMax_) {
374 const auto msg = fmt::format(
375 "Solver failed to converge after cutting timestep {} times.",
376 restarts
377 );
378 if (solverVerbose_) {
379 OpmLog::error(msg);
380 }
381 // Use throw directly to prevent file and line
382 throw TimeSteppingBreakdown{msg};
383 }
384
385 // The new, chopped timestep.
386 const double newTimeStep = restartFactor_ * dt;
387
388
389 // If we have restarted (i.e. cut the timestep) too
390 // much, we have failed and throw an exception.
391 if (newTimeStep < minTimeStep_) {
392 const auto msg = fmt::format(
393 "Solver failed to converge after cutting timestep to {}\n"
394 "which is the minimum threshold given by option --solver-min-time-step\n",
396 );
397 if (solverVerbose_) {
398 OpmLog::error(msg);
399 }
400 // Use throw directly to prevent file and line
401 throw TimeSteppingBreakdown{msg};
402 }
403
404 // Define utility function for chopping timestep.
405 auto chopTimestep = [&]() {
406 substepTimer.provideTimeStepEstimate(newTimeStep);
407 if (solverVerbose_) {
408 const auto msg = fmt::format(
409 "{}\nTimestep chopped to {} days\n",
410 causeOfFailure,
411 std::to_string(unit::convert::to(substepTimer.currentStepLength(), unit::day))
412 );
413 OpmLog::problem(msg);
414 }
415 ++restarts;
416 };
417
418 const double minimumChoppedTimestep = minTimeStepBeforeShuttingProblematicWells_;
419 if (newTimeStep > minimumChoppedTimestep) {
420 chopTimestep();
421 } else {
422 // We are below the threshold, and will check if there are any
423 // wells we should close rather than chopping again.
424 std::set<std::string> failing_wells = detail::consistentlyFailingWells(solver.model().stepReports());
425 if (failing_wells.empty()) {
426 // Found no wells to close, chop the timestep as above.
427 chopTimestep();
428 } else {
429 // Close all consistently failing wells.
430 int num_shut_wells = 0;
431 for (const auto& well : failing_wells) {
432 bool was_shut = solver.model().wellModel().forceShutWellByName(well, substepTimer.simulationTimeElapsed());
433 if (was_shut) {
434 ++num_shut_wells;
435 }
436 }
437 if (num_shut_wells == 0) {
438 // None of the problematic wells were shut.
439 // We must fall back to chopping again.
440 chopTimestep();
441 } else {
442 substepTimer.provideTimeStepEstimate(dt);
443 if (solverVerbose_) {
444 std::string msg;
445 msg = "\nProblematic well(s) were shut: ";
446 for (const auto& well : failing_wells) {
447 msg += well;
448 msg += " ";
449 }
450 msg += "(retrying timestep)\n";
451 OpmLog::problem(msg);
452 }
453 }
454 }
455 }
456 }
457 problem.setNextTimeStepSize(substepTimer.currentStepLength());
458 }
459
460 // store estimated time step for next reportStep
461 suggestedNextTimestep_ = substepTimer.currentStepLength();
462 if (timestepVerbose_) {
463 std::ostringstream ss;
464 substepTimer.report(ss);
465 ss << "Suggested next step size = " << unit::convert::to(suggestedNextTimestep_, unit::day) << " (days)" << std::endl;
466 OpmLog::debug(ss.str());
467 }
468
469 if (! std::isfinite(suggestedNextTimestep_)) { // check for NaN
470 suggestedNextTimestep_ = timestep;
471 }
472 return report;
473 }
474
478 double suggestedNextStep() const
479 { return suggestedNextTimestep_; }
480
481 void setSuggestedNextStep(const double x)
483
484 void updateTUNING(double max_next_tstep, const Tuning& tuning)
485 {
486 restartFactor_ = tuning.TSFCNV;
487 growthFactor_ = tuning.TFDIFF;
488 maxGrowth_ = tuning.TSFMAX;
489 maxTimeStep_ = tuning.TSMAXZ;
490 updateNEXTSTEP(max_next_tstep);
491 timestepAfterEvent_ = tuning.TMAXWC;
492 }
493
494 void updateNEXTSTEP(double max_next_tstep)
495 {
496 // \Note Only update next suggested step if TSINIT was explicitly set in TUNING or NEXTSTEP is active.
497 if (max_next_tstep > 0) {
498 suggestedNextTimestep_ = max_next_tstep;
499 }
500 }
501
502 template<class Serializer>
503 void serializeOp(Serializer& serializer)
504 {
505 serializer(timeStepControlType_);
506 switch (timeStepControlType_) {
507 case TimeStepControlType::HardCodedTimeStep:
508 allocAndSerialize<HardcodedTimeStepControl>(serializer);
509 break;
510 case TimeStepControlType::PIDAndIterationCount:
511 allocAndSerialize<PIDAndIterationCountTimeStepControl>(serializer);
512 break;
513 case TimeStepControlType::SimpleIterationCount:
514 allocAndSerialize<SimpleIterationCountTimeStepControl>(serializer);
515 break;
516 case TimeStepControlType::PID:
517 allocAndSerialize<PIDTimeStepControl>(serializer);
518 break;
519 }
520 serializer(restartFactor_);
521 serializer(growthFactor_);
522 serializer(maxGrowth_);
523 serializer(maxTimeStep_);
524 serializer(minTimeStep_);
525 serializer(ignoreConvergenceFailure_);
526 serializer(solverRestartMax_);
527 serializer(solverVerbose_);
528 serializer(timestepVerbose_);
529 serializer(suggestedNextTimestep_);
530 serializer(fullTimestepInitially_);
531 serializer(timestepAfterEvent_);
532 serializer(useNewtonIteration_);
533 serializer(minTimeStepBeforeShuttingProblematicWells_);
534 }
535
536 static AdaptiveTimeStepping<TypeTag> serializationTestObjectHardcoded()
537 {
538 return serializationTestObject_<HardcodedTimeStepControl>();
539 }
540
541 static AdaptiveTimeStepping<TypeTag> serializationTestObjectPID()
542 {
543 return serializationTestObject_<PIDTimeStepControl>();
544 }
545
546 static AdaptiveTimeStepping<TypeTag> serializationTestObjectPIDIt()
547 {
548 return serializationTestObject_<PIDAndIterationCountTimeStepControl>();
549 }
550
551 static AdaptiveTimeStepping<TypeTag> serializationTestObjectSimple()
552 {
553 return serializationTestObject_<SimpleIterationCountTimeStepControl>();
554 }
555
556 bool operator==(const AdaptiveTimeStepping<TypeTag>& rhs)
557 {
558 if (timeStepControlType_ != rhs.timeStepControlType_ ||
559 (timeStepControl_ && !rhs.timeStepControl_) ||
560 (!timeStepControl_ && rhs.timeStepControl_)) {
561 return false;
562 }
563
564 bool result = false;
565 switch (timeStepControlType_) {
566 case TimeStepControlType::HardCodedTimeStep:
567 result = castAndComp<HardcodedTimeStepControl>(rhs);
568 break;
569 case TimeStepControlType::PIDAndIterationCount:
570 result = castAndComp<PIDAndIterationCountTimeStepControl>(rhs);
571 break;
572 case TimeStepControlType::SimpleIterationCount:
573 result = castAndComp<SimpleIterationCountTimeStepControl>(rhs);
574 break;
575 case TimeStepControlType::PID:
576 result = castAndComp<PIDTimeStepControl>(rhs);
577 break;
578 }
579
580 return result &&
581 this->restartFactor_ == rhs.restartFactor_ &&
582 this->growthFactor_ == rhs.growthFactor_ &&
583 this->maxGrowth_ == rhs.maxGrowth_ &&
584 this->maxTimeStep_ == rhs.maxTimeStep_ &&
585 this->minTimeStep_ == rhs.minTimeStep_ &&
586 this->ignoreConvergenceFailure_ == rhs.ignoreConvergenceFailure_ &&
587 this->solverRestartMax_== rhs.solverRestartMax_ &&
588 this->solverVerbose_ == rhs.solverVerbose_ &&
589 this->fullTimestepInitially_ == rhs.fullTimestepInitially_ &&
590 this->timestepAfterEvent_ == rhs.timestepAfterEvent_ &&
591 this->useNewtonIteration_ == rhs.useNewtonIteration_ &&
592 this->minTimeStepBeforeShuttingProblematicWells_ ==
593 rhs.minTimeStepBeforeShuttingProblematicWells_;
594 }
595
596 private:
597 template<class Controller>
598 static AdaptiveTimeStepping<TypeTag> serializationTestObject_()
599 {
600 AdaptiveTimeStepping<TypeTag> result;
601
602 result.restartFactor_ = 1.0;
603 result.growthFactor_ = 2.0;
604 result.maxGrowth_ = 3.0;
605 result.maxTimeStep_ = 4.0;
606 result.minTimeStep_ = 5.0;
607 result.ignoreConvergenceFailure_ = true;
608 result.solverRestartMax_ = 6;
609 result.solverVerbose_ = true;
610 result.timestepVerbose_ = true;
611 result.suggestedNextTimestep_ = 7.0;
612 result.fullTimestepInitially_ = true;
613 result.useNewtonIteration_ = true;
614 result.minTimeStepBeforeShuttingProblematicWells_ = 9.0;
615 result.timeStepControlType_ = Controller::Type;
616 result.timeStepControl_ = std::make_unique<Controller>(Controller::serializationTestObject());
617
618 return result;
619 }
620 template<class T, class Serializer>
621 void allocAndSerialize(Serializer& serializer)
622 {
623 if (!serializer.isSerializing()) {
624 timeStepControl_ = std::make_unique<T>();
625 }
626 serializer(*static_cast<T*>(timeStepControl_.get()));
627 }
628
629 template<class T>
630 bool castAndComp(const AdaptiveTimeStepping<TypeTag>& Rhs) const
631 {
632 const T* lhs = static_cast<const T*>(timeStepControl_.get());
633 const T* rhs = static_cast<const T*>(Rhs.timeStepControl_.get());
634 return *lhs == *rhs;
635 }
636
637 protected:
638 void init_(const UnitSystem& unitSystem)
639 {
640 // valid are "pid" and "pid+iteration"
641 std::string control = Parameters::Get<Parameters::TimeStepControl>(); // "pid"
642
643 const double tol = Parameters::Get<Parameters::TimeStepControlTolerance>(); // 1e-1
644 if (control == "pid") {
645 timeStepControl_ = std::make_unique<PIDTimeStepControl>(tol);
646 timeStepControlType_ = TimeStepControlType::PID;
647 }
648 else if (control == "pid+iteration") {
649 const int iterations = Parameters::Get<Parameters::TimeStepControlTargetIterations>(); // 30
650 const double decayDampingFactor = Parameters::Get<Parameters::TimeStepControlDecayDampingFactor>(); // 1.0
651 const double growthDampingFactor = Parameters::Get<Parameters::TimeStepControlGrowthDampingFactor>(); // 3.2
652 timeStepControl_ = std::make_unique<PIDAndIterationCountTimeStepControl>(iterations, decayDampingFactor, growthDampingFactor, tol);
653 timeStepControlType_ = TimeStepControlType::PIDAndIterationCount;
654 }
655 else if (control == "pid+newtoniteration") {
656 const int iterations = Parameters::Get<Parameters::TimeStepControlTargetNewtonIterations>(); // 8
657 const double decayDampingFactor = Parameters::Get<Parameters::TimeStepControlDecayDampingFactor>(); // 1.0
658 const double growthDampingFactor = Parameters::Get<Parameters::TimeStepControlGrowthDampingFactor>(); // 3.2
659 const double nonDimensionalMinTimeStepIterations = Parameters::Get<Parameters::MinTimeStepBasedOnNewtonIterations>(); // 0.0 by default
660 // the min time step can be reduced by the newton iteration numbers
661 double minTimeStepReducedByIterations = unitSystem.to_si(UnitSystem::measure::time, nonDimensionalMinTimeStepIterations);
662 timeStepControl_ = std::make_unique<PIDAndIterationCountTimeStepControl>(iterations, decayDampingFactor,
663 growthDampingFactor, tol, minTimeStepReducedByIterations);
664 timeStepControlType_ = TimeStepControlType::PIDAndIterationCount;
665 useNewtonIteration_ = true;
666 }
667 else if (control == "iterationcount") {
668 const int iterations = Parameters::Get<Parameters::TimeStepControlTargetIterations>(); // 30
669 const double decayrate = Parameters::Get<Parameters::TimeStepControlDecayRate>(); // 0.75
670 const double growthrate = Parameters::Get<Parameters::TimeStepControlGrowthRate>(); // 1.25
671 timeStepControl_ = std::make_unique<SimpleIterationCountTimeStepControl>(iterations, decayrate, growthrate);
672 timeStepControlType_ = TimeStepControlType::SimpleIterationCount;
673 }
674 else if (control == "newtoniterationcount") {
675 const int iterations = Parameters::Get<Parameters::TimeStepControlTargetNewtonIterations>(); // 8
676 const double decayrate = Parameters::Get<Parameters::TimeStepControlDecayRate>(); // 0.75
677 const double growthrate = Parameters::Get<Parameters::TimeStepControlGrowthRate>(); // 1.25
678 timeStepControl_ = std::make_unique<SimpleIterationCountTimeStepControl>(iterations, decayrate, growthrate);
679 useNewtonIteration_ = true;
680 timeStepControlType_ = TimeStepControlType::SimpleIterationCount;
681 }
682 else if (control == "hardcoded") {
683 const std::string filename = Parameters::Get<Parameters::TimeStepControlFileName>(); // "timesteps"
684 timeStepControl_ = std::make_unique<HardcodedTimeStepControl>(filename);
685 timeStepControlType_ = TimeStepControlType::HardCodedTimeStep;
686 }
687 else
688 OPM_THROW(std::runtime_error,
689 "Unsupported time step control selected " + control);
690
691 // make sure growth factor is something reasonable
692 assert(growthFactor_ >= 1.0);
693 }
694
695 using TimeStepController = std::unique_ptr<TimeStepControlInterface>;
696
697 TimeStepControlType timeStepControlType_;
698 TimeStepController timeStepControl_;
701 double maxGrowth_;
712 double minTimeStepBeforeShuttingProblematicWells_;
713 };
714}
715
716#endif // OPM_ADAPTIVE_TIME_STEPPING_HPP
Defines a type tags and some fundamental properties all models.
Simulation timer for adaptive time stepping.
Definition AdaptiveSimulatorTimer.hpp:41
Definition AdaptiveTimeStepping.hpp:90
double minTimeStep_
minimal allowed time step size before throwing
Definition AdaptiveTimeStepping.hpp:703
double suggestedNextTimestep_
suggested size of next timestep
Definition AdaptiveTimeStepping.hpp:708
bool solverVerbose_
solver verbosity
Definition AdaptiveTimeStepping.hpp:706
double suggestedNextStep() const
Returns the simulator report for the failed substeps of the last report step.
Definition AdaptiveTimeStepping.hpp:478
double growthFactor_
factor to multiply time step when solver recovered from failed convergence
Definition AdaptiveTimeStepping.hpp:700
AdaptiveTimeStepping(const UnitSystem &unitSystem, const double max_next_tstep=-1.0, const bool terminalOutput=true)
contructor taking parameter object
Definition AdaptiveTimeStepping.hpp:121
TimeStepController timeStepControl_
time step control object
Definition AdaptiveTimeStepping.hpp:698
bool fullTimestepInitially_
beginning with the size of the time step from data file
Definition AdaptiveTimeStepping.hpp:709
AdaptiveTimeStepping(double max_next_tstep, const Tuning &tuning, const UnitSystem &unitSystem, const bool terminalOutput=true)
contructor taking parameter object
Definition AdaptiveTimeStepping.hpp:149
double timestepAfterEvent_
suggested size of timestep after an event
Definition AdaptiveTimeStepping.hpp:710
bool useNewtonIteration_
use newton iteration count for adaptive time step control
Definition AdaptiveTimeStepping.hpp:711
bool timestepVerbose_
timestep verbosity
Definition AdaptiveTimeStepping.hpp:707
int solverRestartMax_
how many restart of solver are allowed
Definition AdaptiveTimeStepping.hpp:705
TimeStepControlType timeStepControlType_
type of time step control object
Definition AdaptiveTimeStepping.hpp:697
SimulatorReport step(const SimulatorTimer &simulatorTimer, Solver &solver, const bool isEvent, const std::vector< int > *fipnum=nullptr, const std::function< bool()> tuningUpdater=[](){return false;})
step method that acts like the solver::step method in a sub cycle of time steps
Definition AdaptiveTimeStepping.hpp:184
bool ignoreConvergenceFailure_
continue instead of stop when minimum time step is reached
Definition AdaptiveTimeStepping.hpp:704
double maxTimeStep_
maximal allowed time step size in days
Definition AdaptiveTimeStepping.hpp:702
double restartFactor_
factor to multiply time step with when solver fails to converge
Definition AdaptiveTimeStepping.hpp:699
double maxGrowth_
factor that limits the maximum growth of a time step
Definition AdaptiveTimeStepping.hpp:701
RelativeChangeInterface.
Definition TimeStepControlInterface.hpp:32
Definition SimulatorTimer.hpp:39
double currentStepLength() const override
Current step length.
Definition SimulatorTimer.cpp:109
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition propertysystem.hh:235
This file provides the infrastructure to retrieve run-time parameters.
The Opm property system, traits with inheritance.
Definition AdaptiveTimeStepping.hpp:57
Definition AdaptiveTimeStepping.hpp:56
Definition AdaptiveTimeStepping.hpp:68
Definition AdaptiveTimeStepping.hpp:52
Definition AdaptiveTimeStepping.hpp:53
Definition AdaptiveTimeStepping.hpp:54
Definition AdaptiveTimeStepping.hpp:64
Definition AdaptiveTimeStepping.hpp:62
Definition AdaptiveTimeStepping.hpp:66
Definition AdaptiveTimeStepping.hpp:65
Definition AdaptiveTimeStepping.hpp:63
Definition AdaptiveTimeStepping.hpp:60
Definition AdaptiveTimeStepping.hpp:61
Definition AdaptiveTimeStepping.hpp:59
Definition AdaptiveTimeStepping.hpp:58
Definition AdaptiveTimeStepping.hpp:55
Definition SimulatorReport.hpp:100
Definition ConvergenceReport.hpp:447