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