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