224 template <
class Solver>
227 const Solver& solver_;
229 SolutionTimeErrorSolverWrapper(
const Solver& solver)
234 double relativeChange()
const
235 {
return solver_.model().relativeChange(); }
239 void logException_(
const E& exception,
bool verbose)
243 message =
"Caught Exception: ";
244 message += exception.what();
245 OpmLog::debug(message);
254 const double max_next_tstep = -1.0,
255 const bool terminalOutput =
true)
257 ,
restartFactor_(Parameters::get<TypeTag, Properties::SolverRestartFactor>())
258 ,
growthFactor_(Parameters::get<TypeTag, Properties::SolverGrowthFactor>())
259 ,
maxGrowth_(Parameters::get<TypeTag, Properties::SolverMaxGrowth>())
260 ,
maxTimeStep_(Parameters::get<TypeTag, Properties::SolverMaxTimeStepInDays>() * 24 * 60 * 60)
261 ,
minTimeStep_(unitSystem.to_si(UnitSystem::measure::time, Parameters::get<TypeTag, Properties::SolverMinTimeStep>()))
264 ,
solverVerbose_(Parameters::get<TypeTag, Properties::SolverVerbosity>() > 0 && terminalOutput)
265 ,
timestepVerbose_(Parameters::get<TypeTag, Properties::TimeStepVerbosity>() > 0 && terminalOutput)
266 ,
suggestedNextTimestep_((max_next_tstep <= 0 ? Parameters::get<TypeTag, Properties::InitialTimeStepInDays>() : max_next_tstep) * 24 * 60 * 60)
268 ,
timestepAfterEvent_(Parameters::get<TypeTag, Properties::TimeStepAfterEventInDays>() * 24 * 60 * 60)
270 , minTimeStepBeforeShuttingProblematicWells_(Parameters::get<TypeTag, Properties::MinTimeStepBeforeShuttingProblematicWellsInDays>() * unit::day)
282 const Tuning& tuning,
283 const UnitSystem& unitSystem,
284 const bool terminalOutput =
true)
293 ,
solverVerbose_(Parameters::get<TypeTag, Properties::SolverVerbosity>() > 0 && terminalOutput)
294 ,
timestepVerbose_(Parameters::get<TypeTag, Properties::TimeStepVerbosity>() > 0 && terminalOutput)
295 ,
suggestedNextTimestep_(max_next_tstep <= 0 ? Parameters::get<TypeTag, Properties::InitialTimeStepInDays>() * 24 * 60 * 60 : max_next_tstep)
299 , minTimeStepBeforeShuttingProblematicWells_(Parameters::get<TypeTag, Properties::MinTimeStepBeforeShuttingProblematicWellsInDays>() * unit::day)
304 static void registerParameters()
306 registerEclTimeSteppingParameters<TypeTag>();
308 Parameters::registerParam<TypeTag, Properties::SolverContinueOnConvergenceFailure>
309 (
"Continue instead of stop when minimum solver time step is reached");
310 Parameters::registerParam<TypeTag, Properties::SolverMaxRestarts>
311 (
"The maximum number of breakdowns before a substep is given up and "
312 "the simulator is terminated");
313 Parameters::registerParam<TypeTag, Properties::SolverVerbosity>
314 (
"Specify the \"chattiness\" of the non-linear solver itself");
315 Parameters::registerParam<TypeTag, Properties::TimeStepVerbosity>
316 (
"Specify the \"chattiness\" during the time integration");
317 Parameters::registerParam<TypeTag, Properties::InitialTimeStepInDays>
318 (
"The size of the initial time step in days");
319 Parameters::registerParam<TypeTag, Properties::FullTimeStepInitially>
320 (
"Always attempt to finish a report step using a single substep");
321 Parameters::registerParam<TypeTag, Properties::TimeStepControl>
322 (
"The algorithm used to determine time-step sizes. "
323 "Valid options are: "
326 "'pid+newtoniteration', "
328 "'newtoniterationcount' "
330 Parameters::registerParam<TypeTag, Properties::TimeStepControlTolerance>
331 (
"The tolerance used by the time step size control algorithm");
332 Parameters::registerParam<TypeTag, Properties::TimeStepControlTargetIterations>
333 (
"The number of linear iterations which the time step control scheme "
334 "should aim for (if applicable)");
335 Parameters::registerParam<TypeTag, Properties::TimeStepControlTargetNewtonIterations>
336 (
"The number of Newton iterations which the time step control scheme "
337 "should aim for (if applicable)");
338 Parameters::registerParam<TypeTag, Properties::TimeStepControlDecayRate>
339 (
"The decay rate of the time step size of the number of "
340 "target iterations is exceeded");
341 Parameters::registerParam<TypeTag, Properties::TimeStepControlGrowthRate>
342 (
"The growth rate of the time step size of the number of "
343 "target iterations is undercut");
344 Parameters::registerParam<TypeTag, Properties::TimeStepControlDecayDampingFactor>
345 (
"The decay rate of the time step decrease when the "
346 "target iterations is exceeded");
347 Parameters::registerParam<TypeTag, Properties::TimeStepControlGrowthDampingFactor>
348 (
"The growth rate of the time step increase when the "
349 "target iterations is undercut");
350 Parameters::registerParam<TypeTag, Properties::TimeStepControlFileName>
351 (
"The name of the file which contains the hardcoded time steps sizes");
352 Parameters::registerParam<TypeTag, Properties::MinTimeStepBeforeShuttingProblematicWellsInDays>
353 (
"The minimum time step size in days for which problematic wells are not shut");
354 Parameters::registerParam<TypeTag, Properties::MinTimeStepBasedOnNewtonIterations>
355 (
"The minimum time step size (in days for field and metric unit and hours for lab unit) "
356 "can be reduced to based on newton iteration counts");
364 template <
class Solver>
368 const std::vector<int>* fipnum =
nullptr,
369 const std::function<
bool()> tuningUpdater = [](){
return false;})
390 auto& simulator = solver.model().simulator();
391 auto& problem = simulator.problem();
400 while (!substepTimer.done()) {
404 if (tuningUpdater()) {
409 const double dt = substepTimer.currentStepLength();
411 detail::logTimer(substepTimer);
414 SimulatorReportSingle substepReport;
415 std::string causeOfFailure;
417 substepReport = solver.step(substepTimer);
421 OpmLog::debug(
"Overall linear iterations used: " + std::to_string(substepReport.total_linear_iterations));
424 catch (
const TooManyIterations& e) {
425 substepReport = solver.failureReport();
426 causeOfFailure =
"Solver convergence failure - Iteration limit reached";
431 catch (
const LinearSolverProblem& e) {
432 substepReport = solver.failureReport();
433 causeOfFailure =
"Linear solver convergence failure";
438 catch (
const NumericalProblem& e) {
439 substepReport = solver.failureReport();
440 causeOfFailure =
"Solver convergence failure - Numerical problem encountered";
445 catch (
const std::runtime_error& e) {
446 substepReport = solver.failureReport();
451 catch (
const Dune::ISTLError& e) {
452 substepReport = solver.failureReport();
457 catch (
const Dune::MatrixBlockError& e) {
458 substepReport = solver.failureReport();
465 simulator.problem().setSubStepReport(substepReport);
467 report += substepReport;
470 !substepReport.converged &&
474 const auto msg = fmt::format(
475 "Solver failed to converge but timestep "
476 "{} is smaller or equal to {}\n"
477 "which is the minimum threshold given "
478 "by option --solver-min-time-step\n",
481 OpmLog::problem(msg);
484 if (substepReport.converged || continue_on_uncoverged_solution) {
490 SolutionTimeErrorSolverWrapper<Solver> relativeChange(solver);
494 : substepReport.total_linear_iterations;
495 double dtEstimate =
timeStepControl_->computeTimeStepSize(dt, iterations, relativeChange,
496 substepTimer.simulationTimeElapsed());
498 assert(dtEstimate > 0);
500 dtEstimate = std::min(dtEstimate,
double(
maxGrowth_ * dt));
501 assert(dtEstimate > 0);
511 if (solver.model().wellModel().hasTHPConstraints()) {
512 const double maxPredictionTHPTimestep = 16.0 * unit::day;
513 dtEstimate = std::min(dtEstimate, maxPredictionTHPTimestep);
515 assert(dtEstimate > 0);
517 std::ostringstream ss;
518 substepReport.reportStep(ss);
519 OpmLog::info(ss.str());
526 if (!substepTimer.done()) {
528 solver.computeFluidInPlace(*fipnum);
530 time::StopWatch perfTimer;
533 problem.writeOutput(simulatorTimer);
535 report.success.output_write_time += perfTimer.secsSinceStart();
539 substepTimer.provideTimeStepEstimate(dtEstimate);
541 report.success.converged = substepTimer.done();
542 substepTimer.setLastStepFailed(
false);
546 substepTimer.setLastStepFailed(
true);
551 const auto msg = fmt::format(
552 "Solver failed to converge after cutting timestep {} times.",
559 throw TimeSteppingBreakdown{msg};
569 const auto msg = fmt::format(
570 "Solver failed to converge after cutting timestep to {}\n"
571 "which is the minimum threshold given by option --solver-min-time-step\n",
578 throw TimeSteppingBreakdown{msg};
582 auto chopTimestep = [&]() {
583 substepTimer.provideTimeStepEstimate(newTimeStep);
585 const auto msg = fmt::format(
586 "{}\nTimestep chopped to {} days\n",
588 std::to_string(unit::convert::to(substepTimer.currentStepLength(), unit::day))
590 OpmLog::problem(msg);
595 const double minimumChoppedTimestep = minTimeStepBeforeShuttingProblematicWells_;
596 if (newTimeStep > minimumChoppedTimestep) {
601 std::set<std::string> failing_wells = detail::consistentlyFailingWells(solver.model().stepReports());
602 if (failing_wells.empty()) {
607 int num_shut_wells = 0;
608 for (
const auto& well : failing_wells) {
609 bool was_shut = solver.model().wellModel().forceShutWellByName(well, substepTimer.simulationTimeElapsed());
614 if (num_shut_wells == 0) {
619 substepTimer.provideTimeStepEstimate(dt);
622 msg =
"\nProblematic well(s) were shut: ";
623 for (
const auto& well : failing_wells) {
627 msg +=
"(retrying timestep)\n";
628 OpmLog::problem(msg);
634 problem.setNextTimeStepSize(substepTimer.currentStepLength());
640 std::ostringstream ss;
641 substepTimer.report(ss);
642 ss <<
"Suggested next step size = " << unit::convert::to(
suggestedNextTimestep_, unit::day) <<
" (days)" << std::endl;
643 OpmLog::debug(ss.str());
658 void setSuggestedNextStep(
const double x)
661 void updateTUNING(
double max_next_tstep,
const Tuning& tuning)
667 updateNEXTSTEP(max_next_tstep);
671 void updateNEXTSTEP(
double max_next_tstep)
674 if (max_next_tstep > 0) {
679 template<
class Serializer>
680 void serializeOp(Serializer& serializer)
684 case TimeStepControlType::HardCodedTimeStep:
685 allocAndSerialize<HardcodedTimeStepControl>(serializer);
687 case TimeStepControlType::PIDAndIterationCount:
688 allocAndSerialize<PIDAndIterationCountTimeStepControl>(serializer);
690 case TimeStepControlType::SimpleIterationCount:
691 allocAndSerialize<SimpleIterationCountTimeStepControl>(serializer);
693 case TimeStepControlType::PID:
694 allocAndSerialize<PIDTimeStepControl>(serializer);
710 serializer(minTimeStepBeforeShuttingProblematicWells_);
713 static AdaptiveTimeStepping<TypeTag> serializationTestObjectHardcoded()
715 return serializationTestObject_<HardcodedTimeStepControl>();
718 static AdaptiveTimeStepping<TypeTag> serializationTestObjectPID()
720 return serializationTestObject_<PIDTimeStepControl>();
723 static AdaptiveTimeStepping<TypeTag> serializationTestObjectPIDIt()
725 return serializationTestObject_<PIDAndIterationCountTimeStepControl>();
728 static AdaptiveTimeStepping<TypeTag> serializationTestObjectSimple()
730 return serializationTestObject_<SimpleIterationCountTimeStepControl>();
733 bool operator==(
const AdaptiveTimeStepping<TypeTag>& rhs)
743 case TimeStepControlType::HardCodedTimeStep:
744 result = castAndComp<HardcodedTimeStepControl>(rhs);
746 case TimeStepControlType::PIDAndIterationCount:
747 result = castAndComp<PIDAndIterationCountTimeStepControl>(rhs);
749 case TimeStepControlType::SimpleIterationCount:
750 result = castAndComp<SimpleIterationCountTimeStepControl>(rhs);
752 case TimeStepControlType::PID:
753 result = castAndComp<PIDTimeStepControl>(rhs);
769 this->minTimeStepBeforeShuttingProblematicWells_ ==
770 rhs.minTimeStepBeforeShuttingProblematicWells_;
774 template<
class Controller>
775 static AdaptiveTimeStepping<TypeTag> serializationTestObject_()
777 AdaptiveTimeStepping<TypeTag> result;
779 result.restartFactor_ = 1.0;
780 result.growthFactor_ = 2.0;
781 result.maxGrowth_ = 3.0;
782 result.maxTimeStep_ = 4.0;
783 result.minTimeStep_ = 5.0;
784 result.ignoreConvergenceFailure_ =
true;
785 result.solverRestartMax_ = 6;
786 result.solverVerbose_ =
true;
787 result.timestepVerbose_ =
true;
788 result.suggestedNextTimestep_ = 7.0;
789 result.fullTimestepInitially_ =
true;
790 result.useNewtonIteration_ =
true;
791 result.minTimeStepBeforeShuttingProblematicWells_ = 9.0;
792 result.timeStepControlType_ = Controller::Type;
793 result.timeStepControl_ = std::make_unique<Controller>(Controller::serializationTestObject());
797 template<
class T,
class Serializer>
798 void allocAndSerialize(Serializer& serializer)
800 if (!serializer.isSerializing()) {
807 bool castAndComp(
const AdaptiveTimeStepping<TypeTag>& Rhs)
const
810 const T* rhs =
static_cast<const T*
>(Rhs.timeStepControl_.get());
815 void init_(
const UnitSystem& unitSystem)
818 std::string control = Parameters::get<TypeTag, Properties::TimeStepControl>();
820 const double tol = Parameters::get<TypeTag, Properties::TimeStepControlTolerance>();
821 if (control ==
"pid") {
825 else if (control ==
"pid+iteration") {
826 const int iterations = Parameters::get<TypeTag, Properties::TimeStepControlTargetIterations>();
827 const double decayDampingFactor = Parameters::get<TypeTag, Properties::TimeStepControlDecayDampingFactor>();
828 const double growthDampingFactor = Parameters::get<TypeTag, Properties::TimeStepControlGrowthDampingFactor>();
829 timeStepControl_ = std::make_unique<PIDAndIterationCountTimeStepControl>(iterations, decayDampingFactor, growthDampingFactor, tol);
832 else if (control ==
"pid+newtoniteration") {
833 const int iterations = Parameters::get<TypeTag, Properties::TimeStepControlTargetNewtonIterations>();
834 const double decayDampingFactor = Parameters::get<TypeTag, Properties::TimeStepControlDecayDampingFactor>();
835 const double growthDampingFactor = Parameters::get<TypeTag, Properties::TimeStepControlGrowthDampingFactor>();
836 const double nonDimensionalMinTimeStepIterations = Parameters::get<TypeTag, Properties::MinTimeStepBasedOnNewtonIterations>();
838 double minTimeStepReducedByIterations = unitSystem.to_si(UnitSystem::measure::time, nonDimensionalMinTimeStepIterations);
839 timeStepControl_ = std::make_unique<PIDAndIterationCountTimeStepControl>(iterations, decayDampingFactor,
840 growthDampingFactor, tol, minTimeStepReducedByIterations);
844 else if (control ==
"iterationcount") {
845 const int iterations = Parameters::get<TypeTag, Properties::TimeStepControlTargetIterations>();
846 const double decayrate = Parameters::get<TypeTag, Properties::TimeStepControlDecayRate>();
847 const double growthrate = Parameters::get<TypeTag, Properties::TimeStepControlGrowthRate>();
848 timeStepControl_ = std::make_unique<SimpleIterationCountTimeStepControl>(iterations, decayrate, growthrate);
851 else if (control ==
"newtoniterationcount") {
852 const int iterations = Parameters::get<TypeTag, Properties::TimeStepControlTargetNewtonIterations>();
853 const double decayrate = Parameters::get<TypeTag, Properties::TimeStepControlDecayRate>();
854 const double growthrate = Parameters::get<TypeTag, Properties::TimeStepControlGrowthRate>();
855 timeStepControl_ = std::make_unique<SimpleIterationCountTimeStepControl>(iterations, decayrate, growthrate);
859 else if (control ==
"hardcoded") {
860 const std::string filename = Parameters::get<TypeTag, Properties::TimeStepControlFileName>();
865 OPM_THROW(std::runtime_error,
866 "Unsupported time step control selected " + control);
872 using TimeStepController = std::unique_ptr<TimeStepControlInterface>;
889 double minTimeStepBeforeShuttingProblematicWells_;