92 template <
class Solver>
95 const Solver& solver_;
97 SolutionTimeErrorSolverWrapper(
const Solver& solver)
102 double relativeChange()
const
103 {
return solver_.model().relativeChange(); }
107 void logException_(
const E& exception,
bool verbose)
111 message =
"Caught Exception: ";
112 message += exception.what();
113 OpmLog::debug(message);
122 const double max_next_tstep = -1.0,
123 const bool terminalOutput =
true)
125 ,
restartFactor_(Parameters::Get<Parameters::SolverRestartFactor<Scalar>>())
126 ,
growthFactor_(Parameters::Get<Parameters::SolverGrowthFactor<Scalar>>())
127 ,
maxGrowth_(Parameters::Get<Parameters::SolverMaxGrowth<Scalar>>())
128 ,
maxTimeStep_(Parameters::Get<Parameters::SolverMaxTimeStepInDays<Scalar>>() * 24 * 60 * 60)
129 ,
minTimeStep_(unitSystem.to_si(UnitSystem::measure::time, Parameters::Get<Parameters::SolverMinTimeStep<Scalar>>()))
132 ,
solverVerbose_(Parameters::Get<Parameters::SolverVerbosity>() > 0 && terminalOutput)
133 ,
timestepVerbose_(Parameters::Get<Parameters::TimeStepVerbosity>() > 0 && terminalOutput)
134 ,
suggestedNextTimestep_((max_next_tstep <= 0 ? Parameters::Get<Parameters::InitialTimeStepInDays>() : max_next_tstep) * 24 * 60 * 60)
136 ,
timestepAfterEvent_(Parameters::Get<Parameters::TimeStepAfterEventInDays<Scalar>>() * 24 * 60 * 60)
138 , minTimeStepBeforeShuttingProblematicWells_(Parameters::Get<Parameters::MinTimeStepBeforeShuttingProblematicWellsInDays>() * unit::day)
150 const Tuning& tuning,
151 const UnitSystem& unitSystem,
152 const bool terminalOutput =
true)
161 ,
solverVerbose_(Parameters::Get<Parameters::SolverVerbosity>() > 0 && terminalOutput)
162 ,
timestepVerbose_(Parameters::Get<Parameters::TimeStepVerbosity>() > 0 && terminalOutput)
163 ,
suggestedNextTimestep_(max_next_tstep <= 0 ? Parameters::Get<Parameters::InitialTimeStepInDays>() * 24 * 60 * 60 : max_next_tstep)
167 , minTimeStepBeforeShuttingProblematicWells_(Parameters::Get<Parameters::MinTimeStepBeforeShuttingProblematicWellsInDays>() * unit::day)
172 static void registerParameters()
174 registerEclTimeSteppingParameters<Scalar>();
175 detail::registerAdaptiveParameters();
183 template <
class Solver>
187 const std::vector<int>* fipnum =
nullptr,
188 const std::function<
bool()> tuningUpdater = [](){
return false;})
209 auto& simulator = solver.model().simulator();
210 auto& problem = simulator.problem();
219 while (!substepTimer.done()) {
223 if (tuningUpdater()) {
228 const double dt = substepTimer.currentStepLength();
230 detail::logTimer(substepTimer);
233 SimulatorReportSingle substepReport;
234 std::string causeOfFailure;
236 substepReport = solver.step(substepTimer);
240 OpmLog::debug(
"Overall linear iterations used: " + std::to_string(substepReport.total_linear_iterations));
243 catch (
const TooManyIterations& e) {
244 substepReport = solver.failureReport();
245 causeOfFailure =
"Solver convergence failure - Iteration limit reached";
250 catch (
const ConvergenceMonitorFailure& e) {
251 substepReport = solver.failureReport();
252 causeOfFailure =
"Convergence monitor failure";
254 catch (
const LinearSolverProblem& e) {
255 substepReport = solver.failureReport();
256 causeOfFailure =
"Linear solver convergence failure";
261 catch (
const NumericalProblem& e) {
262 substepReport = solver.failureReport();
263 causeOfFailure =
"Solver convergence failure - Numerical problem encountered";
268 catch (
const std::runtime_error& e) {
269 substepReport = solver.failureReport();
274 catch (
const Dune::ISTLError& e) {
275 substepReport = solver.failureReport();
280 catch (
const Dune::MatrixBlockError& e) {
281 substepReport = solver.failureReport();
288 simulator.problem().setSubStepReport(substepReport);
290 report += substepReport;
293 !substepReport.converged &&
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",
304 OpmLog::problem(msg);
307 if (substepReport.converged || continue_on_uncoverged_solution) {
313 SolutionTimeErrorSolverWrapper<Solver> relativeChange(solver);
317 : substepReport.total_linear_iterations;
318 double dtEstimate =
timeStepControl_->computeTimeStepSize(dt, iterations, relativeChange,
319 substepTimer.simulationTimeElapsed());
321 assert(dtEstimate > 0);
323 dtEstimate = std::min(dtEstimate,
double(
maxGrowth_ * dt));
324 assert(dtEstimate > 0);
334 if (solver.model().wellModel().hasTHPConstraints()) {
335 const double maxPredictionTHPTimestep = 16.0 * unit::day;
336 dtEstimate = std::min(dtEstimate, maxPredictionTHPTimestep);
338 assert(dtEstimate > 0);
340 std::ostringstream ss;
341 substepReport.reportStep(ss);
342 OpmLog::info(ss.str());
349 if (!substepTimer.done()) {
351 solver.computeFluidInPlace(*fipnum);
353 time::StopWatch perfTimer;
356 problem.writeOutput();
358 report.success.output_write_time += perfTimer.secsSinceStart();
362 substepTimer.provideTimeStepEstimate(dtEstimate);
364 report.success.converged = substepTimer.done();
365 substepTimer.setLastStepFailed(
false);
369 substepTimer.setLastStepFailed(
true);
374 const auto msg = fmt::format(
375 "Solver failed to converge after cutting timestep {} times.",
382 throw TimeSteppingBreakdown{msg};
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",
401 throw TimeSteppingBreakdown{msg};
405 auto chopTimestep = [&]() {
406 substepTimer.provideTimeStepEstimate(newTimeStep);
408 const auto msg = fmt::format(
409 "{}\nTimestep chopped to {} days\n",
411 std::to_string(unit::convert::to(substepTimer.currentStepLength(), unit::day))
413 OpmLog::problem(msg);
418 const double minimumChoppedTimestep = minTimeStepBeforeShuttingProblematicWells_;
419 if (newTimeStep > minimumChoppedTimestep) {
424 std::set<std::string> failing_wells = detail::consistentlyFailingWells(solver.model().stepReports());
425 if (failing_wells.empty()) {
430 int num_shut_wells = 0;
431 for (
const auto& well : failing_wells) {
432 bool was_shut = solver.model().wellModel().forceShutWellByName(well, substepTimer.simulationTimeElapsed());
437 if (num_shut_wells == 0) {
442 substepTimer.provideTimeStepEstimate(dt);
445 msg =
"\nProblematic well(s) were shut: ";
446 for (
const auto& well : failing_wells) {
450 msg +=
"(retrying timestep)\n";
451 OpmLog::problem(msg);
457 problem.setNextTimeStepSize(substepTimer.currentStepLength());
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());
481 void setSuggestedNextStep(
const double x)
484 void updateTUNING(
double max_next_tstep,
const Tuning& tuning)
490 updateNEXTSTEP(max_next_tstep);
494 void updateNEXTSTEP(
double max_next_tstep)
497 if (max_next_tstep > 0) {
502 template<
class Serializer>
503 void serializeOp(Serializer& serializer)
507 case TimeStepControlType::HardCodedTimeStep:
508 allocAndSerialize<HardcodedTimeStepControl>(serializer);
510 case TimeStepControlType::PIDAndIterationCount:
511 allocAndSerialize<PIDAndIterationCountTimeStepControl>(serializer);
513 case TimeStepControlType::SimpleIterationCount:
514 allocAndSerialize<SimpleIterationCountTimeStepControl>(serializer);
516 case TimeStepControlType::PID:
517 allocAndSerialize<PIDTimeStepControl>(serializer);
533 serializer(minTimeStepBeforeShuttingProblematicWells_);
536 static AdaptiveTimeStepping<TypeTag> serializationTestObjectHardcoded()
538 return serializationTestObject_<HardcodedTimeStepControl>();
541 static AdaptiveTimeStepping<TypeTag> serializationTestObjectPID()
543 return serializationTestObject_<PIDTimeStepControl>();
546 static AdaptiveTimeStepping<TypeTag> serializationTestObjectPIDIt()
548 return serializationTestObject_<PIDAndIterationCountTimeStepControl>();
551 static AdaptiveTimeStepping<TypeTag> serializationTestObjectSimple()
553 return serializationTestObject_<SimpleIterationCountTimeStepControl>();
556 bool operator==(
const AdaptiveTimeStepping<TypeTag>& rhs)
566 case TimeStepControlType::HardCodedTimeStep:
567 result = castAndComp<HardcodedTimeStepControl>(rhs);
569 case TimeStepControlType::PIDAndIterationCount:
570 result = castAndComp<PIDAndIterationCountTimeStepControl>(rhs);
572 case TimeStepControlType::SimpleIterationCount:
573 result = castAndComp<SimpleIterationCountTimeStepControl>(rhs);
575 case TimeStepControlType::PID:
576 result = castAndComp<PIDTimeStepControl>(rhs);
592 this->minTimeStepBeforeShuttingProblematicWells_ ==
593 rhs.minTimeStepBeforeShuttingProblematicWells_;
597 template<
class Controller>
598 static AdaptiveTimeStepping<TypeTag> serializationTestObject_()
600 AdaptiveTimeStepping<TypeTag> result;
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());
620 template<
class T,
class Serializer>
621 void allocAndSerialize(Serializer& serializer)
623 if (!serializer.isSerializing()) {
630 bool castAndComp(
const AdaptiveTimeStepping<TypeTag>& Rhs)
const
633 const T* rhs =
static_cast<const T*
>(Rhs.timeStepControl_.get());
638 void init_(
const UnitSystem& unitSystem)
641 std::string control = Parameters::Get<Parameters::TimeStepControl>();
643 const double tol = Parameters::Get<Parameters::TimeStepControlTolerance>();
644 if (control ==
"pid") {
648 else if (control ==
"pid+iteration") {
649 const int iterations = Parameters::Get<Parameters::TimeStepControlTargetIterations>();
650 const double decayDampingFactor = Parameters::Get<Parameters::TimeStepControlDecayDampingFactor>();
651 const double growthDampingFactor = Parameters::Get<Parameters::TimeStepControlGrowthDampingFactor>();
652 timeStepControl_ = std::make_unique<PIDAndIterationCountTimeStepControl>(iterations, decayDampingFactor, growthDampingFactor, tol);
655 else if (control ==
"pid+newtoniteration") {
656 const int iterations = Parameters::Get<Parameters::TimeStepControlTargetNewtonIterations>();
657 const double decayDampingFactor = Parameters::Get<Parameters::TimeStepControlDecayDampingFactor>();
658 const double growthDampingFactor = Parameters::Get<Parameters::TimeStepControlGrowthDampingFactor>();
659 const double nonDimensionalMinTimeStepIterations = Parameters::Get<Parameters::MinTimeStepBasedOnNewtonIterations>();
661 double minTimeStepReducedByIterations = unitSystem.to_si(UnitSystem::measure::time, nonDimensionalMinTimeStepIterations);
662 timeStepControl_ = std::make_unique<PIDAndIterationCountTimeStepControl>(iterations, decayDampingFactor,
663 growthDampingFactor, tol, minTimeStepReducedByIterations);
667 else if (control ==
"iterationcount") {
668 const int iterations = Parameters::Get<Parameters::TimeStepControlTargetIterations>();
669 const double decayrate = Parameters::Get<Parameters::TimeStepControlDecayRate>();
670 const double growthrate = Parameters::Get<Parameters::TimeStepControlGrowthRate>();
671 timeStepControl_ = std::make_unique<SimpleIterationCountTimeStepControl>(iterations, decayrate, growthrate);
674 else if (control ==
"newtoniterationcount") {
675 const int iterations = Parameters::Get<Parameters::TimeStepControlTargetNewtonIterations>();
676 const double decayrate = Parameters::Get<Parameters::TimeStepControlDecayRate>();
677 const double growthrate = Parameters::Get<Parameters::TimeStepControlGrowthRate>();
678 timeStepControl_ = std::make_unique<SimpleIterationCountTimeStepControl>(iterations, decayrate, growthrate);
682 else if (control ==
"hardcoded") {
683 const std::string filename = Parameters::Get<Parameters::TimeStepControlFileName>();
688 OPM_THROW(std::runtime_error,
689 "Unsupported time step control selected " + control);
695 using TimeStepController = std::unique_ptr<TimeStepControlInterface>;
712 double minTimeStepBeforeShuttingProblematicWells_;